/** * 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; }
@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); } }
/** * 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(); }
@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); } }
/** * 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); }
public Matrix invSPD() { try { return CommonsMathDenseDoubleMatrix2DFactory.INSTANCE.dense(new CholeskyDecomposition(matrix).getSolver() .getInverse()); } catch (Exception e) { throw new RuntimeException(e); } }
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); } }
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); }
/** * 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()); }
/** * 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); }
/** * 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(); }
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; } }
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); }
/** * 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); }
/** * Constructor * @param matrix the input matrix */ XCD(RealMatrix matrix) { this.cd = new CholeskyDecomposition(matrix); this.l = LazyValue.of(() -> toDataFrame(cd.getL())); }
/** * 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); }
@Override public Solution[] evolve(Solution[] parents) { int k = parents.length; int n = parents[0].getNumberOfVariables(); RealMatrix x = new Array2DRowRealMatrix(k, n); for (int i=0; i<k; i++) { x.setRow(i, EncodingUtils.getReal(parents[i])); } try { //perform Cholesky factorization and get the upper triangular matrix double jumpRate = Math.pow(jumpRateCoefficient / Math.sqrt(n), 2.0); RealMatrix chol = new CholeskyDecomposition( new Covariance(x.scalarMultiply(jumpRate)) .getCovarianceMatrix()).getLT(); //produce the offspring Solution[] offspring = new Solution[numberOfOffspring]; for (int i=0; i<numberOfOffspring; i++) { Solution child = parents[PRNG.nextInt(parents.length)].copy(); //apply adaptive metropolis step to solution RealVector muC = new ArrayRealVector( EncodingUtils.getReal(child)); RealVector ru = new ArrayRealVector(n); for (int j=0; j<n; j++) { ru.setEntry(j, PRNG.nextGaussian()); } double[] variables = muC.add(chol.preMultiply(ru)).toArray(); //assign variables back to solution for (int j=0; j<n; j++) { RealVariable variable = (RealVariable)child.getVariable(j); double value = variables[j]; if (value < variable.getLowerBound()) { value = variable.getLowerBound(); } else if (value > variable.getUpperBound()) { value = variable.getUpperBound(); } variable.setValue(value); } offspring[i] = child; } return offspring; } catch (Exception e) { return new Solution[0]; } }
public static SampleLikelihood ellipticalSliceSampler(IntDoubleVector initial_pt, CholeskyDecomposition decomp, Function lnpdf) { double initial_lnpdf = lnpdf.getValue(initial_pt); // FIXME: what's the point of the line above? return ellipticalSliceSampler(initial_pt, decomp, lnpdf); }
/** * Calculate Cholesky decomposition A = L L^T * @param a * @return L */ public static Dataset calcCholeskyDecomposition(Dataset a) { CholeskyDecomposition cd = new CholeskyDecomposition(createRealMatrix(a)); return createDataset(cd.getL()); }
/** * This will return the cholesky decomposition of * the given matrix * @param m the matrix to convert * @return the cholesky decomposition of the given * matrix. * See: * http://en.wikipedia.org/wiki/Cholesky_decomposition * @throws NonSquareMatrixException */ public CholeskyDecomposition choleskyFromMatrix(RealMatrix m) throws Exception { return new CholeskyDecomposition(m); }
