Java Code Examples for org.apache.commons.math3.linear.SingularValueDecomposition#getU()

/**
 * 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];
}
 

/**
 * 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};
}
 

/**
 * 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 };
}
 

/**
 * 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 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;
	}
}