Solve CPD

        Args:
            X (numpy.ndarray): Target tensor of CPD.
            rank (int): Rank of the approximate tensor.
            max_iter (int): Max iteration of the ALS.
            stopping_criterion (float): Threshold for stopping the ALS.
                If the value is negative, the convergence check is ignored;
                in other words, it may reduce the computation time.
            lambda_reg (float): regularization parameter. Larger lambda_reg
                means larger regularization.
            dtype (numpy.dtype): Data type

        Returns:
            list of numpy.ndarray: Decomposed matrices.
            numpy.ndarray: Lambda of the CPD.

def solve(self, X, rank, max_iter=500, stopping_criterion=1e-5,
              lambda_reg=0.0, dtype=np.float32, rng=None):
        """Solve CPD

        Args:
            X (numpy.ndarray): Target tensor of CPD.
            rank (int): Rank of the approximate tensor.
            max_iter (int): Max iteration of the ALS.
            stopping_criterion (float): Threshold for stopping the ALS.
                If the value is negative, the convergence check is ignored;
                in other words, it may reduce the computation time.
            lambda_reg (float): regularization parameter. Larger lambda_reg
                means larger regularization.
            dtype (numpy.dtype): Data type

        Returns:
            list of numpy.ndarray: Decomposed matrices.
            numpy.ndarray: Lambda of the CPD.

        """

        N = X.ndim  # Tensor dimensions
        squared_norm_X = np.sum(X ** 2)  # Frobenious norm square

        # Initialize
        if rng is None:
            rng = np.random.RandomState(313)
        A = [None for _ in range(N)]
        for n in range(1, N):
            A[n] = np.array(rng.rand(X.shape[n], rank), dtype=dtype)

        # Solve ALS problem
        criterion = 0
        for itr in range(max_iter):
            criterion_prev = criterion

            # Fix one dimension
            for n in range(N):
                # Solve sub problem
                V = self.hadamard_products_of_gramians(A, n)
                P = self.khatrirao_products(A, n)
                A_n = np.tensordot(
                    X,
                    P,
                    axes=(
                        [o for o in range(N) if o != n],
                        np.arange(N-1)[::-1]))
                # L2 regularization
                V = V+np.eye(rank)*lambda_reg

                A_n = A_n.dot(pinv(V))

                # Normalize
                if itr == 0:
                    lmbda = np.sqrt((A_n ** 2).sum(axis=0))
                else:
                    lmbda = A_n.max(axis=0)
                    lmbda[lmbda < 1] = 1
                A[n] = A_n / lmbda

            # Check convergence
            if stopping_criterion < 0:
                continue
            X_approx = self.approximate_tensor(A, lmbda)
            squared_norm_residual = squared_norm_X + \
                np.sum(X_approx**2) - 2 * np.sum(X*X_approx)
            criterion = 1.0 - (squared_norm_residual / squared_norm_X)
            criterion_change = abs(criterion_prev - criterion)

            if itr > 0 and criterion_change < stopping_criterion:
                break

        return A, lmbda