org.apache.commons.math3.linear.SingularValueDecomposition Java Examples

/**
 * Find a solution of the linear equation A x = b where
 * <ul>
 * <li>A is an n x m - matrix given as double[n][m]</li>
 * <li>b is an m - vector given as double[m],</li>
 * <li>x is an n - vector given as double[n],</li>
 * </ul>
 *
 * @param matrixA The matrix A (left hand side of the linear equation).
 * @param b The vector (right hand of the linear equation).
 * @return A solution x to A x = b.
 */
public static double[] solveLinearEquationSVD(final double[][] matrixA, final double[] b) {

	if(isSolverUseApacheCommonsMath) {
		final Array2DRowRealMatrix matrix = new Array2DRowRealMatrix(matrixA);

		// Using SVD - very slow
		final DecompositionSolver solver = new SingularValueDecomposition(matrix).getSolver();

		return solver.solve(new Array2DRowRealMatrix(b)).getColumn(0);
	}
	else {
		return org.jblas.Solve.solve(new org.jblas.DoubleMatrix(matrixA), new org.jblas.DoubleMatrix(b)).data;

		// For use of colt:
		// cern.colt.matrix.linalg.Algebra linearAlgebra = new cern.colt.matrix.linalg.Algebra();
		// return linearAlgebra.solve(new DenseDoubleMatrix2D(A), linearAlgebra.transpose(new DenseDoubleMatrix2D(new double[][] { b }))).viewColumn(0).toArray();

		// For use of parallel colt:
		// return new cern.colt.matrix.tdouble.algo.decomposition.DenseDoubleLUDecomposition(new cern.colt.matrix.tdouble.impl.DenseDoubleMatrix2D(A)).solve(new cern.colt.matrix.tdouble.impl.DenseDoubleMatrix1D(b)).toArray();
	}
}
 

/**
 * Find a solution of the linear equation A x = b where
 * <ul>
 * <li>A is an n x m - matrix given as double[n][m]</li>
 * <li>b is an m - vector given as double[m],</li>
 * <li>x is an n - vector given as double[n],</li>
 * </ul>
 *
 * @param matrixA The matrix A (left hand side of the linear equation).
 * @param b The vector (right hand of the linear equation).
 * @return A solution x to A x = b.
 */
public static double[] solveLinearEquationSVD(final double[][] matrixA, final double[] b) {

	if(isSolverUseApacheCommonsMath) {
		final Array2DRowRealMatrix matrix = new Array2DRowRealMatrix(matrixA);

		// Using SVD - very slow
		final DecompositionSolver solver = new SingularValueDecomposition(matrix).getSolver();

		return solver.solve(new Array2DRowRealMatrix(b)).getColumn(0);
	}
	else {
		return org.jblas.Solve.solve(new org.jblas.DoubleMatrix(matrixA), new org.jblas.DoubleMatrix(b)).data;

		// For use of colt:
		// cern.colt.matrix.linalg.Algebra linearAlgebra = new cern.colt.matrix.linalg.Algebra();
		// return linearAlgebra.solve(new DenseDoubleMatrix2D(A), linearAlgebra.transpose(new DenseDoubleMatrix2D(new double[][] { b }))).viewColumn(0).toArray();

		// For use of parallel colt:
		// return new cern.colt.matrix.tdouble.algo.decomposition.DenseDoubleLUDecomposition(new cern.colt.matrix.tdouble.impl.DenseDoubleMatrix2D(A)).solve(new cern.colt.matrix.tdouble.impl.DenseDoubleMatrix1D(b)).toArray();
	}
}
 

@Test
public void testPower1() {
    RealMatrix A = new Array2DRowRealMatrix(
        new double[][] {new double[] {1, 2, 3}, new double[] {4, 5, 6}});

    double[] x = new double[3];
    x[0] = Math.random();
    x[1] = Math.random();
    x[2] = Math.random();

    double[] u = new double[2];
    double[] v = new double[3];

    double s = MatrixUtils.power1(A, x, 2, u, v);

    SingularValueDecomposition svdA = new SingularValueDecomposition(A);

    Assert.assertArrayEquals(svdA.getU().getColumn(0), u, 0.001d);
    Assert.assertArrayEquals(svdA.getV().getColumn(0), v, 0.001d);
    Assert.assertEquals(svdA.getSingularValues()[0], s, 0.001d);
}
 

/**
 * Singular Value Decomposition (SVD) based naive scoring.
 */
private double computeScoreSVD(@Nonnull final RealMatrix H, @Nonnull final RealMatrix G) {
    SingularValueDecomposition svdH = new SingularValueDecomposition(H);
    RealMatrix UT = svdH.getUT();

    SingularValueDecomposition svdG = new SingularValueDecomposition(G);
    RealMatrix Q = svdG.getU();

    // find the largest singular value for the r principal components
    RealMatrix UTQ = UT.getSubMatrix(0, r - 1, 0, window - 1)
                       .multiply(Q.getSubMatrix(0, window - 1, 0, r - 1));
    SingularValueDecomposition svdUTQ = new SingularValueDecomposition(UTQ);
    double[] s = svdUTQ.getSingularValues();

    return 1.d - s[0];
}
 

/**
 * Return the solution x of XTX x = XT y for a given y.
 * @TODO Performance upon repeated call can be optimized by caching XTX.
 *
 * @param dependents The sample vector of the random variable y.
 * @return The solution x of XTX x = XT y.
 */
public double[] getLinearRegressionParameters(final RandomVariable dependents) {

	final RandomVariable[] basisFunctions = basisFunctionsEstimator.getBasisFunctions();

	synchronized (solverLock) {
		if(solver == null) {
			// Build XTX - the symmetric matrix consisting of the scalar products of the basis functions.
			final double[][] XTX = new double[basisFunctions.length][basisFunctions.length];
			for(int i=0; i<basisFunctions.length; i++) {
				for(int j=i; j<basisFunctions.length; j++) {
					XTX[i][j] = basisFunctions[i].mult(basisFunctions[j]).getAverage();	// Scalar product
					XTX[j][i] = XTX[i][j];												// Symmetric matrix
				}
			}

			solver = new SingularValueDecomposition(new Array2DRowRealMatrix(XTX, false)).getSolver();
		}
	}

	// Build XTy - the projection of the dependents random variable on the basis functions.
	final double[] XTy = new double[basisFunctions.length];
	for(int i=0; i<basisFunctions.length; i++) {
		XTy[i] = dependents.mult(basisFunctions[i]).getAverage();				// Scalar product
	}

	// Solve X^T X x = X^T y - which gives us the regression coefficients x = linearRegressionParameters
	final double[] linearRegressionParameters = solver.solve(new ArrayRealVector(XTy)).toArray();

	return linearRegressionParameters;
}
 

/**
 * Performs Singular Value Decomposition. Calls Apache Commons Math SVD.
 * X = U * Sigma * Vt, where X is the input matrix,
 * U is the left singular matrix, Sigma is the singular values matrix returned as a
 * column matrix and Vt is the transpose of the right singular matrix V.
 * However, the returned array has  { U, Sigma, V}
 * 
 * @param in Input matrix
 * @return An array containing U, Sigma & V
 */
private static MatrixBlock[] computeSvd(MatrixBlock in) {
	Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);

	SingularValueDecomposition svd = new SingularValueDecomposition(matrixInput);
	double[] sigma = svd.getSingularValues();
	RealMatrix u = svd.getU();
	RealMatrix v = svd.getV();
	MatrixBlock U = DataConverter.convertToMatrixBlock(u.getData());
	MatrixBlock Sigma = DataConverter.convertToMatrixBlock(sigma, true);
	Sigma = LibMatrixReorg.diag(Sigma, new MatrixBlock(Sigma.rlen, Sigma.rlen, true));
	MatrixBlock V = DataConverter.convertToMatrixBlock(v.getData());

	return new MatrixBlock[] { U, Sigma, V };
}
 

/** Create a SVD instance using Apache Commons Math.
 *
 * @param m matrix that is not {@code null}
 * @return SVD instance that is never {@code null}
 */
@Override
public SVD createSVD(final RealMatrix m) {

    Utils.nonNull(m, "Cannot create SVD on a null matrix.");

    final SingularValueDecomposition svd = new SingularValueDecomposition(m);
    final RealMatrix pinv = svd.getSolver().getInverse();
    return new SimpleSVD(svd.getU(), svd.getSingularValues(), svd.getV(), pinv);
}
 

@Override
public double getCondition(Matrix m) {
  ArgChecker.notNull(m, "m");
  if (m instanceof DoubleMatrix) {
    RealMatrix temp = CommonsMathWrapper.wrap((DoubleMatrix) m);
    SingularValueDecomposition svd = new SingularValueDecomposition(temp);
    return svd.getConditionNumber();
  }
  throw new IllegalArgumentException("Can only find condition number of DoubleMatrix; have " + m.getClass());
}
 

@Override
public DoubleMatrix getInverse(Matrix m) {
  ArgChecker.notNull(m, "matrix was null");
  if (m instanceof DoubleMatrix) {
    RealMatrix temp = CommonsMathWrapper.wrap((DoubleMatrix) m);
    SingularValueDecomposition sv = new SingularValueDecomposition(temp);
    RealMatrix inv = sv.getSolver().getInverse();
    return CommonsMathWrapper.unwrap(inv);
  }
  throw new IllegalArgumentException("Can only find inverse of DoubleMatrix; have " + m.getClass());
}
 

@Override
public SVDecompositionResult apply(DoubleMatrix x) {
  ArgChecker.notNull(x, "x");
  MatrixValidate.notNaNOrInfinite(x);
  RealMatrix commonsMatrix = CommonsMathWrapper.wrap(x);
  SingularValueDecomposition svd = new SingularValueDecomposition(commonsMatrix);
  return new SVDecompositionCommonsResult(svd);
}
 

/**
 * Creates an instance.
 * 
 * @param svd The result of the SV decomposition, not null
 */
public SVDecompositionCommonsResult(SingularValueDecomposition svd) {
  ArgChecker.notNull(svd, "svd");
  _condition = svd.getConditionNumber();
  _norm = svd.getNorm();
  _rank = svd.getRank();
  _s = CommonsMathWrapper.unwrap(svd.getS());
  _singularValues = svd.getSingularValues();
  _u = CommonsMathWrapper.unwrap(svd.getU());
  _uTranspose = CommonsMathWrapper.unwrap(svd.getUT());
  _v = CommonsMathWrapper.unwrap(svd.getV());
  _vTranspose = CommonsMathWrapper.unwrap(svd.getVT());
  _solver = svd.getSolver();
}
 

/**
 * Performs Singular Value Decomposition. Calls Apache Commons Math SVD.
 * X = U * Sigma * Vt, where X is the input matrix,
 * U is the left singular matrix, Sigma is the singular values matrix returned as a
 * column matrix and Vt is the transpose of the right singular matrix V.
 * However, the returned array has  { U, Sigma, V}
 * 
 * @param in Input matrix
 * @return An array containing U, Sigma & V
 */
private static MatrixBlock[] computeSvd(MatrixBlock in) {
	Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in);

	SingularValueDecomposition svd = new SingularValueDecomposition(matrixInput);
	double[] sigma = svd.getSingularValues();
	RealMatrix u = svd.getU();
	RealMatrix v = svd.getV();
	MatrixBlock U = DataConverter.convertToMatrixBlock(u.getData());
	MatrixBlock Sigma = DataConverter.convertToMatrixBlock(sigma, true);
	Sigma = LibMatrixReorg.diag(Sigma, new MatrixBlock(Sigma.rlen, Sigma.rlen, true));
	MatrixBlock V = DataConverter.convertToMatrixBlock(v.getData());

	return new MatrixBlock[] { U, Sigma, V };
}
 

/**
 * Return the solution x of XTX x = XT y for a given y.
 * @TODO Performance upon repeated call can be optimized by caching XTX.
 *
 * @param dependents The sample vector of the random variable y.
 * @return The solution x of XTX x = XT y.
 */
@Override
public double[] getLinearRegressionParameters(RandomVariable dependents) {

	final RandomVariable localizerWeights = dependents.squared().sub(Math.pow(dependents.getStandardDeviation()*standardDeviations,2.0)).choose(new Scalar(0.0), new Scalar(1.0));

	// Localize basis functions
	final RandomVariable[] basisFunctionsNonLocalized = getBasisFunctionsEstimator().getBasisFunctions();
	final RandomVariable[] basisFunctions = new RandomVariable[basisFunctionsNonLocalized.length];
	for(int i=0; i<basisFunctions.length; i++) {
		basisFunctions[i] = basisFunctionsNonLocalized[i].mult(localizerWeights);
	}

	// Localize dependents
	dependents = dependents.mult(localizerWeights);

	// Build XTX - the symmetric matrix consisting of the scalar products of the basis functions.
	final double[][] XTX = new double[basisFunctions.length][basisFunctions.length];
	for(int i=0; i<basisFunctions.length; i++) {
		for(int j=i; j<basisFunctions.length; j++) {
			XTX[i][j] = basisFunctions[i].mult(basisFunctions[j]).getAverage();	// Scalar product
			XTX[j][i] = XTX[i][j];												// Symmetric matrix
		}
	}

	final DecompositionSolver solver = new SingularValueDecomposition(new Array2DRowRealMatrix(XTX, false)).getSolver();

	// Build XTy - the projection of the dependents random variable on the basis functions.
	final double[] XTy = new double[basisFunctions.length];
	for(int i=0; i<basisFunctions.length; i++) {
		XTy[i] = dependents.mult(basisFunctions[i]).getAverage();				// Scalar product
	}

	// Solve X^T X x = X^T y - which gives us the regression coefficients x = linearRegressionParameters
	final double[] linearRegressionParameters = solver.solve(new ArrayRealVector(XTy)).toArray();

	return linearRegressionParameters;
}
 

/**
 * Pseudo-Inverse of a matrix calculated in the least square sense.
 *
 * @param matrix The given matrix A.
 * @return pseudoInverse The pseudo-inverse matrix P, such that A*P*A = A and P*A*P = P
 */
public static double[][] pseudoInverse(final double[][] matrix){
	if(isSolverUseApacheCommonsMath) {
		// Use LU from common math
		final SingularValueDecomposition svd = new SingularValueDecomposition(new Array2DRowRealMatrix(matrix));
		final double[][] matrixInverse = svd.getSolver().getInverse().getData();

		return matrixInverse;
	}
	else {
		return org.jblas.Solve.pinv(new org.jblas.DoubleMatrix(matrix)).toArray2();
	}
}
 

/**
 * Return the solution x of XTX x = XT y for a given y.
 * @TODO Performance upon repeated call can be optimized by caching XTX.
 *
 * @param dependents The sample vector of the random variable y.
 * @return The solution x of XTX x = XT y.
 */
public double[] getLinearRegressionParameters(final RandomVariable dependents) {

	final RandomVariable[] basisFunctions = basisFunctionsEstimator.getBasisFunctions();

	synchronized (solverLock) {
		if(solver == null) {
			// Build XTX - the symmetric matrix consisting of the scalar products of the basis functions.
			final double[][] XTX = new double[basisFunctions.length][basisFunctions.length];
			for(int i=0; i<basisFunctions.length; i++) {
				for(int j=i; j<basisFunctions.length; j++) {
					XTX[i][j] = basisFunctions[i].mult(basisFunctions[j]).getAverage();	// Scalar product
					XTX[j][i] = XTX[i][j];												// Symmetric matrix
				}
			}

			solver = new SingularValueDecomposition(new Array2DRowRealMatrix(XTX, false)).getSolver();
		}
	}

	// Build XTy - the projection of the dependents random variable on the basis functions.
	final double[] XTy = new double[basisFunctions.length];
	for(int i=0; i<basisFunctions.length; i++) {
		XTy[i] = dependents.mult(basisFunctions[i]).getAverage();				// Scalar product
	}

	// Solve X^T X x = X^T y - which gives us the regression coefficients x = linearRegressionParameters
	final double[] linearRegressionParameters = solver.solve(new ArrayRealVector(XTy)).toArray();

	return linearRegressionParameters;
}
 

/**
 * Return the solution x of XTX x = XT y for a given y.
 * @TODO Performance upon repeated call can be optimized by caching XTX.
 *
 * @param dependents The sample vector of the random variable y.
 * @return The solution x of XTX x = XT y.
 */
@Override
public double[] getLinearRegressionParameters(RandomVariable dependents) {

	final RandomVariable localizerWeights = dependents.squared().sub(Math.pow(dependents.getStandardDeviation()*standardDeviations,2.0)).choose(new Scalar(0.0), new Scalar(1.0));

	// Localize basis functions
	final RandomVariable[] basisFunctionsNonLocalized = getBasisFunctionsEstimator().getBasisFunctions();
	final RandomVariable[] basisFunctions = new RandomVariable[basisFunctionsNonLocalized.length];
	for(int i=0; i<basisFunctions.length; i++) {
		basisFunctions[i] = basisFunctionsNonLocalized[i].mult(localizerWeights);
	}

	// Localize dependents
	dependents = dependents.mult(localizerWeights);

	// Build XTX - the symmetric matrix consisting of the scalar products of the basis functions.
	final double[][] XTX = new double[basisFunctions.length][basisFunctions.length];
	for(int i=0; i<basisFunctions.length; i++) {
		for(int j=i; j<basisFunctions.length; j++) {
			XTX[i][j] = basisFunctions[i].mult(basisFunctions[j]).getAverage();	// Scalar product
			XTX[j][i] = XTX[i][j];												// Symmetric matrix
		}
	}

	final DecompositionSolver solver = new SingularValueDecomposition(new Array2DRowRealMatrix(XTX, false)).getSolver();

	// Build XTy - the projection of the dependents random variable on the basis functions.
	final double[] XTy = new double[basisFunctions.length];
	for(int i=0; i<basisFunctions.length; i++) {
		XTy[i] = dependents.mult(basisFunctions[i]).getAverage();				// Scalar product
	}

	// Solve X^T X x = X^T y - which gives us the regression coefficients x = linearRegressionParameters
	final double[] linearRegressionParameters = solver.solve(new ArrayRealVector(XTy)).toArray();

	return linearRegressionParameters;
}
 

/**
 * Pseudo-Inverse of a matrix calculated in the least square sense.
 *
 * @param matrix The given matrix A.
 * @return pseudoInverse The pseudo-inverse matrix P, such that A*P*A = A and P*A*P = P
 */
public static double[][] pseudoInverse(final double[][] matrix){
	if(isSolverUseApacheCommonsMath) {
		// Use LU from common math
		final SingularValueDecomposition svd = new SingularValueDecomposition(new Array2DRowRealMatrix(matrix));
		final double[][] matrixInverse = svd.getSolver().getInverse().getData();

		return matrixInverse;
	}
	else {
		return org.jblas.Solve.pinv(new org.jblas.DoubleMatrix(matrix)).toArray2();
	}
}
 

/**
 * Calculates the compact Singular Value Decomposition of a matrix. The
 * Singular Value Decomposition of matrix A is a set of three matrices: U, Σ
 * and V such that A = U × Σ × VT. Let A be a m × n matrix, then U is a m ×
 * p orthogonal matrix, Σ is a p × p diagonal matrix with positive or null
 * elements, V is a p × n orthogonal matrix (hence VT is also orthogonal)
 * where p=min(m,n).
 *
 * @param a Given matrix.
 * @return Result U/S/V arrays.
 */
public static Array[] svd(Array a) {
    int m = a.getShape()[0];
    int n = a.getShape()[1];
    int k = Math.min(m, n);
    Array Ua = Array.factory(DataType.DOUBLE, new int[]{m, k});
    Array Va = Array.factory(DataType.DOUBLE, new int[]{k, n});
    Array Sa = Array.factory(DataType.DOUBLE, new int[]{k});
    double[][] aa = (double[][]) ArrayUtil.copyToNDJavaArray_Double(a);
    RealMatrix matrix = new Array2DRowRealMatrix(aa, false);
    SingularValueDecomposition decomposition = new SingularValueDecomposition(matrix);
    RealMatrix U = decomposition.getU();
    RealMatrix V = decomposition.getVT();
    double[] sv = decomposition.getSingularValues();
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < k; j++) {
            Ua.setDouble(i * k + j, U.getEntry(i, j));
        }
    }
    for (int i = 0; i < k; i++) {
        for (int j = 0; j < n; j++) {
            Va.setDouble(i * n + j, V.getEntry(i, j));
        }
    }
    for (int i = 0; i < k; i++) {
        Sa.setDouble(i, sv[i]);
    }

    return new Array[]{Ua, Sa, Va};
}
 

private void updateInverseCovariance() {
    try {
        inverseCov = new LUDecomposition(cov).getSolver().getInverse();
    } catch (SingularMatrixException e) {
        singularCovariances.inc();
        inverseCov = new SingularValueDecomposition(cov).getSolver().getInverse();
    }
}
 

private double computeL(double mu) {
	double LPenalty = lpenalty * mu;
	SingularValueDecomposition svd = new SingularValueDecomposition(X.subtract(S));
	double[] penalizedD = softThreshold(svd.getSingularValues(), LPenalty);
	RealMatrix D_matrix = MatrixUtils.createRealDiagonalMatrix(penalizedD);
	L = svd.getU().multiply(D_matrix).multiply(svd.getVT());
	return sum(penalizedD) * LPenalty;
}
 

public void updateCoefficients(double l2penalty) {
       if (this.X_svd == null) {
       	this.X_svd = new SingularValueDecomposition(X);
       }
    RealMatrix V = this.X_svd.getV();
    double[] s = this.X_svd.getSingularValues();
    RealMatrix U = this.X_svd.getU();
    
    for (int i = 0; i < s.length; i++) {
    	s[i] = s[i] / (s[i]*s[i] + l2penalty);
    }
    RealMatrix S = MatrixUtils.createRealDiagonalMatrix(s);
    
    RealMatrix Z = V.multiply(S).multiply(U.transpose());
    
    this.coefficients = Z.operate(this.Y);
    
    this.fitted = this.X.operate(this.coefficients);
    double errorVariance = 0;
    for (int i = 0; i < residuals.length; i++) {
    	this.residuals[i] = this.Y[i] - this.fitted[i];
    	errorVariance += this.residuals[i] * this.residuals[i];
    }
    errorVariance = errorVariance / (X.getRowDimension() - X.getColumnDimension());
    
    RealMatrix errorVarianceMatrix = MatrixUtils.createRealIdentityMatrix(this.Y.length).scalarMultiply(errorVariance);
    RealMatrix coefficientsCovarianceMatrix = Z.multiply(errorVarianceMatrix).multiply(Z.transpose());
    this.standarderrors = getDiagonal(coefficientsCovarianceMatrix);
}
 

/**
 * L = A x R
 *
 * @return a matrix A that minimizes A x R - L
 */
@Nonnull
public static RealMatrix solve(@Nonnull final RealMatrix L, @Nonnull final RealMatrix R,
        final boolean exact) throws SingularMatrixException {
    LUDecomposition LU = new LUDecomposition(R);
    DecompositionSolver solver = LU.getSolver();
    final RealMatrix A;
    if (exact || solver.isNonSingular()) {
        A = LU.getSolver().solve(L);
    } else {
        SingularValueDecomposition SVD = new SingularValueDecomposition(R);
        A = SVD.getSolver().solve(L);
    }
    return A;
}
 

@Nonnull
public static RealMatrix inverse(@Nonnull final RealMatrix m, final boolean exact)
        throws SingularMatrixException {
    LUDecomposition LU = new LUDecomposition(m);
    DecompositionSolver solver = LU.getSolver();
    final RealMatrix inv;
    if (exact || solver.isNonSingular()) {
        inv = solver.getInverse();
    } else {
        SingularValueDecomposition SVD = new SingularValueDecomposition(m);
        inv = SVD.getSolver().getInverse();
    }
    return inv;
}
 

/**
 * pdf(x, x_hat) = exp(-0.5 * (x-x_hat) * inv(Σ) * (x-x_hat)T) / ( 2π^0.5d * det(Σ)^0.5)
 * 
 * @return value of probabilistic density function
 * @link https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Density_function
 */
public static double pdf(@Nonnull final RealVector x, @Nonnull final RealVector x_hat,
        @Nonnull final RealMatrix sigma) {
    final int dim = x.getDimension();
    Preconditions.checkArgument(x_hat.getDimension() == dim,
        "|x| != |x_hat|, |x|=" + dim + ", |x_hat|=" + x_hat.getDimension());
    Preconditions.checkArgument(sigma.getRowDimension() == dim,
        "|x| != |sigma|, |x|=" + dim + ", |sigma|=" + sigma.getRowDimension());
    Preconditions.checkArgument(sigma.isSquare(), "Sigma is not square matrix");

    LUDecomposition LU = new LUDecomposition(sigma);
    final double detSigma = LU.getDeterminant();
    double denominator = Math.pow(2.d * Math.PI, 0.5d * dim) * Math.pow(detSigma, 0.5d);
    if (denominator == 0.d) { // avoid divide by zero
        return 0.d;
    }

    final RealMatrix invSigma;
    DecompositionSolver solver = LU.getSolver();
    if (solver.isNonSingular() == false) {
        SingularValueDecomposition svd = new SingularValueDecomposition(sigma);
        invSigma = svd.getSolver().getInverse(); // least square solution
    } else {
        invSigma = solver.getInverse();
    }
    //EigenDecomposition eigen = new EigenDecomposition(sigma);
    //double detSigma = eigen.getDeterminant();
    //RealMatrix invSigma = eigen.getSolver().getInverse();

    RealVector diff = x.subtract(x_hat);
    RealVector premultiplied = invSigma.preMultiply(diff);
    double sum = premultiplied.dotProduct(diff);
    double numerator = Math.exp(-0.5d * sum);

    return numerator / denominator;
}
 

/**
 * Constructor
 * @param x     the input matrix
 */
XSVD(RealMatrix x) {
    final SingularValueDecomposition svd = new SingularValueDecomposition(x);
    this.rank = svd.getRank();
    this.u = toDataFrame(svd.getU());
    this.v = toDataFrame(svd.getV());
    this.s = toDataFrame(svd.getS());
    this.singularValues = Array.of(svd.getSingularValues());
}
 

@Override
public PCA fit(RealMatrix X) {
	synchronized(fitLock) {
		this.centerer = new MeanCenterer().fit(X);
		this.m = X.getRowDimension();
		this.n = X.getColumnDimension();
		
		// ensure n_components not too large
		if(this.n_components > n)
			this.n_components = n;
		
		final RealMatrix data = this.centerer.transform(X);
		SingularValueDecomposition svd = new SingularValueDecomposition(data);
		RealMatrix U = svd.getU(), S = svd.getS(), V = svd.getV().transpose();
		
		// flip Eigenvectors' sign to enforce deterministic output
		EntryPair<RealMatrix, RealMatrix> uv_sign_swap = eigenSignFlip(U, V);
		
		U = uv_sign_swap.getKey();
		V = uv_sign_swap.getValue();
		RealMatrix components_ = V;
		
		
		// get variance explained by singular value
		final double[] s = MatUtils.diagFromSquare(S.getData());
		this.variabilities = new double[s.length];
		for(int i= 0; i < s.length; i++) {
			variabilities[i] = (s[i]*s[i]) / (double)m;
			total_var += variabilities[i];
		}
		
		
		// get variability ratio
		this.variability_ratio = new double[s.length];
		for(int i = 0; i < s.length; i++) {
			variability_ratio[i] = variabilities[i] / total_var;
		}
		
		
		// post-process number of components if in var_mode
		double[] ratio_cumsum = VecUtils.cumsum(variability_ratio);
		if(this.var_mode) {
			for(int i = 0; i < ratio_cumsum.length; i++) {
				if(ratio_cumsum[i] >= this.variability) {
					this.n_components = i + 1;
					break;
				}
				
				// if it never hits the if block, the n_components is
				// equal to the number of columns in its entirety
			}
		}
		
		
		// get noise variance
		if(n_components < FastMath.min(n, m)) {
			this.noise_variance = VecUtils.mean(VecUtils.slice(variabilities, n_components, s.length));
		} else {
			this.noise_variance = 0.0;
		}
		
		
		// Set the components and other sliced variables
		this.components = new Array2DRowRealMatrix(MatUtils.slice(components_.getData(), 0, n_components), false);
		this.variabilities = VecUtils.slice(variabilities, 0, n_components);
		this.variability_ratio = VecUtils.slice(variability_ratio, 0, n_components);
		
		if(retain) {
			this.U = new Array2DRowRealMatrix(MatUtils.slice(U.getData(), 0, n_components), false);;
			this.S = new Array2DRowRealMatrix(MatUtils.slice(S.getData(), 0, n_components), false);;
		}
		
		
		return this;
	}
}
 

/**
 * @param a
 * @return array of singular values
 */
public static double[] calcSingularValues(Dataset a) {
	SingularValueDecomposition svd = new SingularValueDecomposition(createRealMatrix(a));
	return svd.getSingularValues();
}
 

/**
 * Calculate singular value decomposition {@code A = U S V^T}
 * @param a
 * @return array of U - orthogonal matrix, s - singular values vector, V - orthogonal matrix
 */
public static Dataset[] calcSingularValueDecomposition(Dataset a) {
	SingularValueDecomposition svd = new SingularValueDecomposition(createRealMatrix(a));
	return new Dataset[] {createDataset(svd.getU()), DatasetFactory.createFromObject(svd.getSingularValues()),
			createDataset(svd.getV())};
}
 

/**
 * Calculate matrix rank by singular value decomposition method
 * @param a
 * @return effective numerical rank of matrix
 */
public static int calcMatrixRank(Dataset a) {
	SingularValueDecomposition svd = new SingularValueDecomposition(createRealMatrix(a));
	return svd.getRank();
}
 

/**
 * Find a solution of the linear equation A X = B in the least square sense where
 * <ul>
 * <li>A is an n x m - matrix given as double[n][m]</li>
 * <li>B is an m x k - matrix given as double[m][k],</li>
 * <li>X is an n x k - matrix given as double[n][k],</li>
 * </ul>
 *
 * @param matrix The matrix A (left hand side of the linear equation).
 * @param rhs The matrix B (right hand of the linear equation).
 * @return A solution X to A X = B.
 */
public static double[][] solveLinearEquationLeastSquare(final double[][] matrix, final double[][] rhs) {
	// We use the linear algebra package apache commons math
	final DecompositionSolver solver = new SingularValueDecomposition(new Array2DRowRealMatrix(matrix, false)).getSolver();
	return solver.solve(new Array2DRowRealMatrix(rhs)).getData();
}