Java CholeskyDecomposition类代码示例

import org.apache.commons.math3.linear.CholeskyDecomposition; //导入依赖的package包/类
/**
 * Calculates the Cholesky decomposition of a matrix. The Cholesky
 * decomposition of a real symmetric positive-definite matrix A consists of
 * a lower triangular matrix L with same size such that: A = LLT. In a
 * sense, this is the square root of A.
 *
 * @param a The given matrix.
 * @return Result array.
 */
public static Array cholesky(Array a) {
    Array r = Array.factory(DataType.DOUBLE, a.getShape());
    double[][] aa = (double[][]) ArrayUtil.copyToNDJavaArray(a);
    RealMatrix matrix = new Array2DRowRealMatrix(aa, false);
    CholeskyDecomposition decomposition = new CholeskyDecomposition(matrix);
    RealMatrix L = decomposition.getL();
    int n = L.getColumnDimension();
    int m = L.getRowDimension();
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            r.setDouble(i * n + j, L.getEntry(i, j));
        }
    }

    return r;
}
 

import org.apache.commons.math3.linear.CholeskyDecomposition; //导入依赖的package包/类
@Override
public SlopeCoefficients estimateCoefficients(final DerivationEquation eq)
    throws EstimationException {
  final double[][] sourceTriangleMatrix = eq.getCovarianceLowerTriangularMatrix();
  // Copy matrix and enhance it to a full matrix as expected by CholeskyDecomposition
  // FIXME: Avoid copy job to speed-up the solving process e.g. by extending the CholeskyDecomposition constructor
  final int length = sourceTriangleMatrix.length;
  final double[][] matrix = new double[length][];
  for (int i = 0; i < length; i++) {
    matrix[i] = new double[length];
    final double[] s = sourceTriangleMatrix[i];
    final double[] t = matrix[i];
    for (int j = 0; j <= i; j++) {
      t[j] = s[j];
    }
    for (int j = i + 1; j < length; j++) {
      t[j] = sourceTriangleMatrix[j][i];
    }
  }
  final RealMatrix coefficients =
      new Array2DRowRealMatrix(matrix, false);
  try {
    final DecompositionSolver solver = new CholeskyDecomposition(coefficients).getSolver();
    final RealVector constants = new ArrayRealVector(eq.getConstraints(), true);
    final RealVector solution = solver.solve(constants);
    return new DefaultSlopeCoefficients(solution.toArray());
  } catch (final NonPositiveDefiniteMatrixException e) {
    throw new EstimationException("Matrix inversion error due to data is linearly dependent", e);
  }
}
 

import org.apache.commons.math3.linear.CholeskyDecomposition; //导入依赖的package包/类
/**
 * Returns the 'root' (U) of the inverse covariance matrix S^{-1},
 * such that S^{-1} = U^T . U
 * This matrix can be used to pre-transform the original sample
 * vectors X (by X → U . X) to a space where distance measurement (in the Mahalanobis
 * sense) can be calculated with the usual Euclidean norm.
 * The matrix U is invertible in case the reverse mapping is required.
 * 
 * @return The matrix for pre-transforming the original sample vectors.
 */
public double[][] getWhiteningTransformation() {
	IJ.log("in whitening");
	double relativeSymmetryThreshold = 1.0E-6; 		// CholeskyDecomposition.DEFAULT_RELATIVE_SYMMETRY_THRESHOLD == 1.0E-15; too small!
       double absolutePositivityThreshold = 1.0E-10;	// CholeskyDecomposition.DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD == 1.0E-10;
	CholeskyDecomposition cd = 
			new CholeskyDecomposition(MatrixUtils.createRealMatrix(iCov),
					relativeSymmetryThreshold, absolutePositivityThreshold);
	RealMatrix U = cd.getLT();
	return U.getData();
}
 

import org.apache.commons.math3.linear.CholeskyDecomposition; //导入依赖的package包/类
@Override
protected RealVector solve(final RealMatrix jacobian,
                           final RealVector residuals) {
    try {
        final Pair<RealMatrix, RealVector> normalEquation =
                computeNormalMatrix(jacobian, residuals);
        final RealMatrix normal = normalEquation.getFirst();
        final RealVector jTr = normalEquation.getSecond();
        return new CholeskyDecomposition(
                normal, SINGULARITY_THRESHOLD, SINGULARITY_THRESHOLD)
                .getSolver()
                .solve(jTr);
    } catch (NonPositiveDefiniteMatrixException e) {
        throw new ConvergenceException(LocalizedFormats.UNABLE_TO_SOLVE_SINGULAR_PROBLEM, e);
    }
}
 

import org.apache.commons.math3.linear.CholeskyDecomposition; //导入依赖的package包/类
/**
 * Correct the current state estimate with an actual measurement.
 *
 * @param z the measurement vector
 * @throws DimensionMismatchException if the dimension of the
 * measurement vector does not fit
 * @throws org.apache.commons.math3.linear.SingularMatrixException
 * if the covariance matrix could not be inverted
 */
public void correct(final RealVector z) {
    // sanity checks
    MathUtils.checkNotNull(z);
    if (z.getDimension() != measurementMatrix.getRowDimension()) {
        throw new DimensionMismatchException(z.getDimension(),
                                             measurementMatrix.getRowDimension());
    }

    // S = H * P(k) - * H' + R
    RealMatrix s = measurementMatrix.multiply(errorCovariance)
        .multiply(measurementMatrixT)
        .add(measurementModel.getMeasurementNoise());

    // invert S
    // as the error covariance matrix is a symmetric positive
    // semi-definite matrix, we can use the cholesky decomposition
    DecompositionSolver solver = new CholeskyDecomposition(s).getSolver();
    RealMatrix invertedS = solver.getInverse();

    // Inn = z(k) - H * xHat(k)-
    RealVector innovation = z.subtract(measurementMatrix.operate(stateEstimation));

    // calculate gain matrix
    // K(k) = P(k)- * H' * (H * P(k)- * H' + R)^-1
    // K(k) = P(k)- * H' * S^-1
    RealMatrix kalmanGain = errorCovariance.multiply(measurementMatrixT).multiply(invertedS);

    // update estimate with measurement z(k)
    // xHat(k) = xHat(k)- + K * Inn
    stateEstimation = stateEstimation.add(kalmanGain.operate(innovation));

    // update covariance of prediction error
    // P(k) = (I - K * H) * P(k)-
    RealMatrix identity = MatrixUtils.createRealIdentityMatrix(kalmanGain.getRowDimension());
    errorCovariance = identity.subtract(kalmanGain.multiply(measurementMatrix)).multiply(errorCovariance);
}
 

import org.apache.commons.math3.linear.CholeskyDecomposition; //导入依赖的package包/类
public Matrix invSPD() {
	try {
		return CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(new CholeskyDecomposition(matrix).getSolver()
				.getInverse());
	} catch (Exception e) {
		throw new RuntimeException(e);
	}
}
 

import org.apache.commons.math3.linear.CholeskyDecomposition; //导入依赖的package包/类
public Matrix solveSPD(Matrix b) {
	try {
		if (b instanceof AbstractCommonsMathDenseDoubleMatrix2D) {
			AbstractCommonsMathDenseDoubleMatrix2D b2 = (AbstractCommonsMathDenseDoubleMatrix2D) b;
			RealMatrix ret = new CholeskyDecomposition(matrix).getSolver().solve(b2.matrix);
			return CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(ret);
		} else {
			return super.solve(b);
		}
	} catch (Exception e) {
		throw new RuntimeException(e);
	}
}
 

import org.apache.commons.math3.linear.CholeskyDecomposition; //导入依赖的package包/类
public void fit(List<Location> locations){
	setLocations(locations);
	int n = locations.size();
	double[][] K = new double[n][n];

	for(int i=0; i<n; i++){
		Location loc1 = locations.get(i);
		double[] x1 = ModelAdaptUtils.locationToVec(loc1);
		for(int j=i; j<n; j++){
			Location loc2 = locations.get(j);
			double[] x2 = ModelAdaptUtils.locationToVec(loc2);
			double k =kernel.computeKernel(x1, x2);
			K[i][j] = k;
			K[j][i] = k;
		}
	}
	RealMatrix Kmat = MatrixUtils.createRealMatrix(K);
	RealMatrix lambdamat = MatrixUtils.createRealIdentityMatrix(n).scalarMultiply(sigma_n*sigma_n); //Tentative treatment

	RealMatrix Kymat = Kmat.add(lambdamat);

	CholeskyDecomposition chol = new CholeskyDecomposition(Kymat);
	RealMatrix Lmat = chol.getL();

	double[] normalRands = new double[n];
	for(int i=0; i<n; i++){
		normalRands[i] = rand.nextGaussian();
	}
	this.y = Lmat.operate(normalRands);

	RealMatrix invKymat = (new LUDecomposition(Kymat)).getSolver().getInverse();
	this.alphas = invKymat.operate(y);
}
 

import org.apache.commons.math3.linear.CholeskyDecomposition; //导入依赖的package包/类
/**
 * Function to compute Cholesky decomposition of the given input matrix. 
 * The input must be a real symmetric positive-definite matrix.
 * 
 * @param in commons-math3 Array2DRowRealMatrix
 * @return matrix block
 * @throws DMLRuntimeException if DMLRuntimeException occurs
 */
private static MatrixBlock computeCholesky(Array2DRowRealMatrix in) 
	throws DMLRuntimeException 
{	
	if ( !in.isSquare() )
		throw new DMLRuntimeException("Input to cholesky() must be square matrix -- given: a " + in.getRowDimension() + "x" + in.getColumnDimension() + " matrix.");

	CholeskyDecomposition cholesky = new CholeskyDecomposition(in, 1e-14,
		CholeskyDecomposition.DEFAULT_ABSOLUTE_POSITIVITY_THRESHOLD);
	RealMatrix rmL = cholesky.getL();
	
	return DataConverter.convertToMatrixBlock(rmL.getData());
}
 

import org.apache.commons.math3.linear.CholeskyDecomposition; //导入依赖的package包/类
/**
 * Correct the current state estimate with an actual measurement.
 *
 * @param z
 *            the measurement vector
 * @throws NullArgumentException
 *             if the measurement vector is {@code null}
 * @throws DimensionMismatchException
 *             if the dimension of the measurement vector does not fit
 * @throws SingularMatrixException
 *             if the covariance matrix could not be inverted
 */
public void correct(final RealVector z)
        throws NullArgumentException, DimensionMismatchException, SingularMatrixException {

    // sanity checks
    MathUtils.checkNotNull(z);
    if (z.getDimension() != measurementMatrix.getRowDimension()) {
        throw new DimensionMismatchException(z.getDimension(),
                                             measurementMatrix.getRowDimension());
    }

    // S = H * P(k) - * H' + R
    RealMatrix s = measurementMatrix.multiply(errorCovariance)
        .multiply(measurementMatrixT)
        .add(measurementModel.getMeasurementNoise());

    // invert S
    // as the error covariance matrix is a symmetric positive
    // semi-definite matrix, we can use the cholesky decomposition
    DecompositionSolver solver = new CholeskyDecomposition(s).getSolver();
    RealMatrix invertedS = solver.getInverse();

    // Inn = z(k) - H * xHat(k)-
    RealVector innovation = z.subtract(measurementMatrix.operate(stateEstimation));

    // calculate gain matrix
    // K(k) = P(k)- * H' * (H * P(k)- * H' + R)^-1
    // K(k) = P(k)- * H' * S^-1
    RealMatrix kalmanGain = errorCovariance.multiply(measurementMatrixT).multiply(invertedS);

    // update estimate with measurement z(k)
    // xHat(k) = xHat(k)- + K * Inn
    stateEstimation = stateEstimation.add(kalmanGain.operate(innovation));

    // update covariance of prediction error
    // P(k) = (I - K * H) * P(k)-
    RealMatrix identity = MatrixUtils.createRealIdentityMatrix(kalmanGain.getRowDimension());
    errorCovariance = identity.subtract(kalmanGain.multiply(measurementMatrix)).multiply(errorCovariance);
}
 

import org.apache.commons.math3.linear.CholeskyDecomposition; //导入依赖的package包/类
/**
 * Constructor.
 * @param ch The result of the Cholesky decomposition.
 */
public CholeskyDecompositionCommonsResult(CholeskyDecomposition ch) {
  ArgChecker.notNull(ch, "Cholesky decomposition");
  _determinant = ch.getDeterminant();
  _l = CommonsMathWrapper.unwrap(ch.getL());
  _lt = CommonsMathWrapper.unwrap(ch.getLT());
  _solver = ch.getSolver();
}
 

import org.apache.commons.math3.linear.CholeskyDecomposition; //导入依赖的package包/类
public static SampleLikelihood ellipticalSliceSampler(IntDoubleVector initial_pt, double initial_lnpdf, CholeskyDecomposition decomp, Function lnpdf) {	
	
	int D = lnpdf.getNumDimensions();
	double [] x = new double[D];
	for(int i=0; i<D; i++) x[i] = initial_pt.get(i);
	RealMatrix init_val = MatrixUtils.createColumnRealMatrix(x);
	RealMatrix L = decomp.getL();
	RealMatrix r = MatrixUtils.createColumnRealMatrix(getNormalVector(D));
	RealMatrix nu = L.multiply(r);
	double hh = Math.log(Prng.nextDouble()) + initial_lnpdf;
	
	// Set up the ellipse and the slice threshold
	double phi = Prng.nextDouble()*2.0*Math.PI;
	double phi_min = phi-2.0*Math.PI;
	double phi_max = phi;
	
	// Slice sampling loop
	while(true) {
		// Compute xx for proposed angle difference and check if it's on the slice
		RealMatrix xx_prop = init_val.scalarMultiply(phi).add(nu.scalarMultiply(Math.sin(phi)));
		double cur_lnpdf = lnpdf.getValue(new IntDoubleDenseVector(xx_prop.getColumn(0)));
		if(cur_lnpdf > hh) {
			return new SampleLikelihood(xx_prop.getColumnVector(0).toArray(), cur_lnpdf);
		}
		// Shrink slice to rejected point
		if(phi > 0) {
			phi_max = phi;
		} else if(phi < 0) {
			phi_min = phi;
		} else {
			throw new RuntimeException("Shrunk to current position and still not acceptable");
		}
		// Propose new angle difference
		phi = Prng.nextDouble()*(phi_max - phi_min) + phi_min;
	}
}
 

import org.apache.commons.math3.linear.CholeskyDecomposition; //导入依赖的package包/类
private void initialize() {
    final int dim = scaleMatrix.getColumnDimension();

    // Build gamma distributions for the diagonal
    gammas = new GammaDistribution[dim];
    for (int i = 0; i < dim; i++) {

        gammas[i] = new GammaDistribution(df-i-.99/2, 2);
    }

    // Build the cholesky decomposition
    cholesky = new CholeskyDecomposition(scaleMatrix);
}
 

import org.apache.commons.math3.linear.CholeskyDecomposition; //导入依赖的package包/类
/**
 * Correct the current state estimate with an actual measurement.
 *
 * @param z
 *            the measurement vector
 * @throws NullArgumentException
 *             if the measurement vector is {@code null}
 * @throws DimensionMismatchException
 *             if the dimension of the measurement vector does not fit
 * @throws SingularMatrixException
 *             if the covariance matrix could not be inverted
 */
public void correct(final RealVector z) throws NullArgumentException,
           DimensionMismatchException, SingularMatrixException
{

	// sanity checks
	MathUtils.checkNotNull(z);
	if (z.getDimension() != measurementMatrix.getRowDimension())
	{
		throw new DimensionMismatchException(z.getDimension(),
				measurementMatrix.getRowDimension());
	}

	// S = H * P(k) * H' + R
	RealMatrix s = measurementMatrix.multiply(errorCovariance)
			.multiply(measurementMatrixT)
			.add(measurementModel.getMeasurementNoise());

	// Inn = z(k) - H * xHat(k)-
	RealVector innovation = z.subtract(measurementMatrix
			.operate(stateEstimation));

	// calculate gain matrix
	// K(k) = P(k)- * H' * (H * P(k)- * H' + R)^-1
	// K(k) = P(k)- * H' * S^-1

	// instead of calculating the inverse of S we can rearrange the formula,
	// and then solve the linear equation A x X = B with A = S', X = K' and
	// B = (H * P)'

	// K(k) * S = P(k)- * H'
	// S' * K(k)' = H * P(k)-'
	RealMatrix kalmanGain = new CholeskyDecomposition(s).getSolver()
			.solve(measurementMatrix.multiply(errorCovariance.transpose()))
			.transpose();

	// update estimate with measurement z(k)
	// xHat(k) = xHat(k)- + K * Inn
	stateEstimation = stateEstimation.add(kalmanGain.operate(innovation));

	// update covariance of prediction error
	// P(k) = (I - K * H) * P(k)-
	RealMatrix identity = MatrixUtils.createRealIdentityMatrix(kalmanGain
			.getRowDimension());
	errorCovariance = identity.subtract(
			kalmanGain.multiply(measurementMatrix)).multiply(
			errorCovariance);
}
 

import org.apache.commons.math3.linear.CholeskyDecomposition; //导入依赖的package包/类
/**
 * Constructor
 * @param matrix    the input matrix
 */
XCD(RealMatrix matrix) {
    this.cd = new CholeskyDecomposition(matrix);
    this.l = LazyValue.of(() -> toDataFrame(cd.getL()));
}