Gradient Descent optimization process
        
        Tunes the embeddings (Y) so that their pairwise distance distribution
        matches the input high-dimensional data (X) pairwise distance distribution.
        In other words, minimizes the KL divergence cost.

def fit(self, X):
        """Gradient Descent optimization process
        
        Tunes the embeddings (Y) so that their pairwise distance distribution
        matches the input high-dimensional data (X) pairwise distance distribution.
        In other words, minimizes the KL divergence cost.

        """

        # compute the std. deviation with mean vector x in X
        self._compute_std_dev(X)

        # Kullback–Leibler divergence
        kl_cost = KL_Divergence()

        # compute high-dimensional affinities (Gaussian Distribution)
        high_dim_dist = self._get_high_dim_dist(X)
        # sample initial solutions
        Y = np.random.randn(X.shape[0], self.n_components)

        prev_Ys = [Y, Y]

        for iteration in range(1, self.n_iter+1):
            # compute low-dimensional affinities (Student t-Distribution)
            low_dim_dist = self._get_low_dim_dist(Y)
    
            for i in range(Y.shape[0]):
                # compute gradient
                grad = kl_cost.gradient(high_dim_dist, low_dim_dist, Y, i)
                # set new Y[i]
                Y[i] = prev_Ys[1][i] + self.learning_rate * grad + self.momentum * (prev_Ys[1][i] - prev_Ys[0][i])

            prev_Ys = [prev_Ys[1], Y]

            if iteration % 100 == 0 and self.verbose:
                low_dim_dist = self._get_low_dim_dist(Y)
                print(f"ITERATION: {iteration}{3*' '}|||{3*' '}KL divergence: {kl_cost(high_dim_dist, low_dim_dist)}")

        self.embeddings = Y
        return self