org.ejml.factory.DecompositionFactory Java Examples

private static Complex64F[] findRoots(double... coefficients) {
    int N = coefficients.length - 1;

    // Construct the companion matrix. This is a square N x N matrix.
    final DenseMatrix64F c = new DenseMatrix64F(N, N);

    double a = coefficients[N];
    for (int i = 0; i < N; i++) {
        c.set(i, N - 1, -coefficients[i] / a);
    }
    for (int i = 1; i < N; i++) {
        c.set(i, i - 1, 1);
    }

    // Use generalized eigenvalue decomposition to find the roots.
    EigenDecomposition<DenseMatrix64F> evd = DecompositionFactory.eig(N, false);

    evd.decompose(c);

    final Complex64F[] roots = new Complex64F[N];

    for (int i = 0; i < N; i++) {
        roots[i] = evd.getEigenvalue(i);
    }

    return roots;
}
 

private static void rotateNd(double[] x, int dim) {

        // Get first `dim` locations
        DenseMatrix64F matrix = new DenseMatrix64F(dim, dim);
        for (int row = 0; row < dim; ++row) {
            for (int col = 0; col < dim; ++col) {
                matrix.set(row, col, x[col * dim + row]);
            }
        }

        // Do a QR decomposition
        QRDecomposition<DenseMatrix64F> qr = DecompositionFactory.qr(dim, dim);
        qr.decompose(matrix);
        DenseMatrix64F qm = qr.getQ(null, true);
        DenseMatrix64F rm = qr.getR(null, true);

        // Reflection invariance
        if (rm.get(0,0) < 0) {
            CommonOps.scale(-1, rm);
            CommonOps.scale(-1, qm);
        }

        // Compute Q^{-1}
        DenseMatrix64F qInv = new DenseMatrix64F(dim, dim);
        CommonOps.transpose(qm, qInv);

        // Apply to each location
        for (int location = 0; location < x.length / dim; ++location) {
            WrappedVector locationVector = new WrappedVector.Raw(x, location * dim, dim);
            MissingOps.matrixVectorMultiple(qInv, locationVector, locationVector, dim);
        }
    }
 

@Override
public EigenDecomposition decomposeStrenghtOfSelection(MatrixParameterInterface AParam, int dim, boolean isSymmetric) {
    DenseMatrix64F A = MissingOps.wrap(AParam);
    org.ejml.interfaces.decomposition.EigenDecomposition eigA = DecompositionFactory.eig(dim, true, isSymmetric);
    if (!eigA.decompose(A)) throw new RuntimeException("Eigen decomposition failed.");
    return new EigenDecomposition(eigenVectorsMatrix(eigA),
            null,
            eigenValuesMatrix(eigA, dim));
}
 

static DenseMatrix64F getCholeskyOfVariance(DenseMatrix64F variance, final int dim) {

            org.ejml.interfaces.decomposition.CholeskyDecomposition<DenseMatrix64F> engine =
                    DecompositionFactory.chol(dim, true);
            engine.decompose(variance);

            return engine.getT(null);
        }
 

private EigenDecomposition decomposeStrenghtOfSelection(double[] Aparam) {
    int n = getDimTrait();
    DenseMatrix64F A = MissingOps.wrap(Aparam, 0, n, n);
    // Decomposition
    EigenDecomposition eigA = DecompositionFactory.eig(n, true, false);
    if (!eigA.decompose(A)) throw new RuntimeException("Eigen decomposition failed.");
    return eigA;
}
 

private EigenDecomposition decomposeStrenghtOfSelection(double[] Aparam) {
    int n = getDimTrait();
    DenseMatrix64F A = MissingOps.wrap(Aparam, 0, n, n);
    // Decomposition
    EigenDecomposition eigA = DecompositionFactory.eig(n, true, false);
    if (!eigA.decompose(A)) throw new RuntimeException("Eigen decomposition failed.");
    return eigA;
}
 

public static InversionResult safeDeterminant(DenseMatrix64F source, boolean invert) {
        final int finiteCount = countFiniteNonZeroDiagonals(source);

        InversionResult result;

        if (finiteCount == 0) {
            result = new InversionResult(NOT_OBSERVED, 0, Double.NEGATIVE_INFINITY, true);
        } else {
//            LinearSolver<DenseMatrix64F> solver = LinearSolverFactory.pseudoInverse(true);
//            solver.setA(source);
//
//            SingularValueDecomposition<DenseMatrix64F> svd = solver.getDecomposition();
//            double[] values = svd.getSingularValues();
//
//            if (values == null) {
//                throw new RuntimeException("Unable to perform SVD");
//            }

            SingularValueDecomposition<DenseMatrix64F> svd = DecompositionFactory.svd(source.getNumRows(), source.getNumCols(), false, false, false);
            if (!svd.decompose(source)) {
                if (SingularOps.rank(svd) == 0)
                    return new InversionResult(NOT_OBSERVED, 0, Double.NEGATIVE_INFINITY, true);
                throw new RuntimeException("SVD decomposition failed");
            }
            double[] values = svd.getSingularValues();
            double tol = SingularOps.singularThreshold(svd);
//            double tol = 0.0;

            int dim = 0;
            double logDet = 0;
            for (int i = 0; i < values.length; i++) {
                final double lambda = values[i];
                if (lambda > tol) {
                    logDet += Math.log(lambda);
                    ++dim;
                }
            }

            if (invert) {
                logDet = -logDet;
            }

            result = new InversionResult(dim == source.getNumCols() ? FULLY_OBSERVED : PARTIALLY_OBSERVED, dim, logDet, true);
        }

        return result;
    }
 

/**
 * Computes a basis (the principle components) from the most dominant eigenvectors.
 */
public void computeBasis() {
	if (sampleIndex != numSamples)
		throw new IllegalArgumentException("Not all the data has been added");
	if (numComponents > numSamples)
		throw new IllegalArgumentException(
				"More data needed to compute the desired number of components");

	means = new DenseMatrix64F(sampleSize, 1);
	// compute the mean of all the samples
	for (int i = 0; i < numSamples; i++) {
		for (int j = 0; j < sampleSize; j++) {
			double val = means.get(j);
			means.set(j, val + A.get(i, j));
		}
	}
	for (int j = 0; j < sampleSize; j++) {
		double avg = means.get(j) / numSamples;
		means.set(j, avg);
	}

	// subtract the mean from the original data
	for (int i = 0; i < numSamples; i++) {
		for (int j = 0; j < sampleSize; j++) {
			A.set(i, j, A.get(i, j) - means.get(j));
		}
	}

	// compute SVD and save time by not computing U
	SingularValueDecomposition<DenseMatrix64F> svd = DecompositionFactory.svd(numSamples, sampleSize,
			false, true, compact);
	if (!svd.decompose(A))
		throw new RuntimeException("SVD failed");

	V_t = svd.getV(null, true);
	W = svd.getW(null);

	// singular values are in an arbitrary order initially and need to be sorted in descending order
	SingularOps.descendingOrder(null, false, W, V_t, true);

	// strip off unneeded components and find the basis
	V_t.reshape(numComponents, sampleSize, true);

}