List of usage examples for org.apache.commons.math3.linear CholeskyDecomposition CholeskyDecomposition
public CholeskyDecomposition(final RealMatrix matrix)
From source file:com.opengamma.strata.math.impl.linearalgebra.CholeskyDecompositionCommons.java
/** * {@inheritDoc}//from w ww. j av a 2s .co m */ @Override public CholeskyDecompositionResult apply(DoubleMatrix x) { ArgChecker.notNull(x, "x"); RealMatrix temp = CommonsMathWrapper.wrap(x); CholeskyDecomposition cholesky; try { cholesky = new CholeskyDecomposition(temp); } catch (Exception e) { throw new MathException(e.toString()); } return new CholeskyDecompositionCommonsResult(cholesky); }
From source file:edu.washington.gs.skyline.model.quantification.WeightedRegression.java
public static double[] weighted(double[][] x, double[] y, double[] weights, boolean intercept) { RealMatrix predictor;//from ww w .ja v a 2 s . co m if (intercept) { int nRows = x.length; int nCols = x[0].length + 1; predictor = MatrixUtils.createRealMatrix(nRows, nCols); for (int iRow = 0; iRow < nRows; iRow++) { predictor.setEntry(iRow, 0, 1.0); for (int iCol = 1; iCol < nCols; iCol++) { predictor.setEntry(iRow, iCol, x[iRow][iCol - 1]); } } } else { predictor = MatrixUtils.createRealMatrix(x); } RealVector responseVector = MatrixUtils.createRealVector(y); RealMatrix weightsMatrix = MatrixUtils.createRealDiagonalMatrix(weights); RealMatrix predictorTransposed = predictor.transpose(); RealMatrix predictorTransposedTimesWeights = predictorTransposed .multiply(weightsMatrix.multiply(predictor)); CholeskyDecomposition choleskyDecomposition = new CholeskyDecomposition(predictorTransposedTimesWeights); RealVector vectorToSolve = predictorTransposed.operate(weightsMatrix.operate(responseVector)); RealVector result = choleskyDecomposition.getSolver().solve(vectorToSolve); return result.toArray(); }
From source file:com.joptimizer.algebra.CholeskyDecompositionTest.java
/** * good decomposition.// w ww .j a va 2 s . co m */ public void testDecomposition1() throws Exception { log.debug("testDecomposition1"); RealMatrix P1 = new Array2DRowRealMatrix(new double[][] { { 8.08073550734687, 1.59028724315583 }, { 1.59028724315583, 0.3250861184011492 } }); CholeskyDecomposition cFact1 = new CholeskyDecomposition(P1); log.debug("L: " + cFact1.getL()); log.debug("LT: " + cFact1.getLT()); // check L.LT-Q=0 RealMatrix P1Inv = cFact1.getL().multiply(cFact1.getLT()); double norm1 = P1Inv.subtract(P1).getNorm(); log.debug("norm1: " + norm1); assertTrue(norm1 < 1.E-12); }
From source file:com.joptimizer.algebra.CholeskyDecompositionTest.java
/** * poor decomposition./*from w w w . j a v a 2s . c om*/ * rescaling can help in doing it better */ public void testDecomposition2() throws Exception { log.debug("testDecomposition2"); RealMatrix P1 = new Array2DRowRealMatrix(new double[][] { { 8.185301256666552E9, 1.5977225251367908E9 }, { 1.5977225251367908E9, 3.118660129093004E8 } }); CholeskyDecomposition cFact1 = new CholeskyDecomposition(P1); log.debug("L: " + cFact1.getL()); log.debug("LT: " + cFact1.getLT()); // check L.LT-Q=0 double norm1 = cFact1.getL().multiply(cFact1.getLT()).subtract(P1).getNorm(); log.debug("norm1: " + norm1); assertTrue(norm1 < 1.E-5); //poor precision, try to make it better //geometric eigenvalues mean DescriptiveStatistics ds = new DescriptiveStatistics(new double[] { 8.5E9, 0.00572 }); RealMatrix P2 = P1.scalarMultiply(1. / ds.getGeometricMean()); CholeskyDecomposition cFact2 = new CholeskyDecomposition(P2); log.debug("L: " + cFact2.getL()); log.debug("LT: " + cFact2.getLT()); // check L.LT-Q=0 double norm2 = cFact2.getL().multiply(cFact2.getLT()).subtract(P2).getNorm(); log.debug("norm2: " + norm2); assertTrue(norm2 < Utils.getDoubleMachineEpsilon()); }
From source file:com.joptimizer.solvers.CholeskyTest.java
/** * poor decomposition./*from ww w . j a va 2 s . co m*/ * rescaling can help in doing it better */ public void testDecomposition2() throws Exception { log.debug("testDecomposition2"); RealMatrix P1 = new Array2DRowRealMatrix(new double[][] { { 8.185301256666552E9, 1.5977225251367908E9 }, { 1.5977225251367908E9, 3.118660129093004E8 } }); CholeskyDecomposition cFact1 = new CholeskyDecomposition(P1); log.debug("L: " + cFact1.getL()); log.debug("LT: " + cFact1.getLT()); // check L.LT-Q=0 double norm1 = cFact1.getL().multiply(cFact1.getLT()).subtract(P1).getNorm(); log.debug("norm1: " + norm1); assertTrue(norm1 < 1.E-5); //poor precision, try to make it better //geometric eigenvalues mean DescriptiveStatistics ds = new DescriptiveStatistics(new double[] { 8.5E9, 0.00572 }); RealMatrix P2 = P1.scalarMultiply(1. / ds.getGeometricMean()); CholeskyDecomposition cFact2 = new CholeskyDecomposition(P2); log.debug("L: " + cFact2.getL()); log.debug("LT: " + cFact2.getLT()); // check L.LT-Q=0 double norm2 = cFact2.getL().multiply(cFact2.getLT()).subtract(P2).getNorm(); log.debug("norm2: " + norm2); assertTrue(norm2 < 1.E-9); }
From source file:com.joptimizer.functions.SDPLogarithmicBarrier.java
/** * @see "S.Boyd and L.Vandenberghe, Convex Optimization, p. 618" *///w w w .ja va 2 s . c om public double value(double[] X) { RealMatrix S = buildS(X); try { CholeskyDecomposition cFact = new CholeskyDecomposition(S); double detS = cFact.getDeterminant(); return -Math.log(detS); } catch (NonPositiveDefiniteMatrixException e) { return Double.NaN; } }
From source file:com.joptimizer.optimizers.NewtonUnconstrainedTest.java
/** * Quadratic objective./* ww w. j a v a 2 s . c o m*/ */ public void testOptimize() throws Exception { log.debug("testOptimize"); // START SNIPPET: newtonUnconstrained-1 RealMatrix PMatrix = new Array2DRowRealMatrix( new double[][] { { 1.68, 0.34, 0.38 }, { 0.34, 3.09, -1.59 }, { 0.38, -1.59, 1.54 } }); RealVector qVector = new ArrayRealVector(new double[] { 0.018, 0.025, 0.01 }); // Objective function. double theta = 0.01522; RealMatrix P = PMatrix.scalarMultiply(theta); RealVector q = qVector.mapMultiply(-1); PDQuadraticMultivariateRealFunction objectiveFunction = new PDQuadraticMultivariateRealFunction(P.getData(), q.toArray(), 0); OptimizationRequest or = new OptimizationRequest(); or.setF0(objectiveFunction); or.setInitialPoint(new double[] { 0.04, 0.50, 0.46 }); or.setTolerance(1.e-8); //optimization NewtonUnconstrained opt = new NewtonUnconstrained(); opt.setOptimizationRequest(or); int returnCode = opt.optimize(); // END SNIPPET: newtonUnconstrained-1 if (returnCode == OptimizationResponse.FAILED) { fail(); } OptimizationResponse response = opt.getOptimizationResponse(); double[] sol = response.getSolution(); log.debug("sol : " + ArrayUtils.toString(sol)); log.debug("value : " + objectiveFunction.value(sol)); // we know the analytic solution of the problem // sol = -PInv * q CholeskyDecomposition cFact = new CholeskyDecomposition(P); RealVector benchSol = cFact.getSolver().solve(q).mapMultiply(-1); log.debug("benchSol : " + ArrayUtils.toString(benchSol.toArray())); log.debug("benchValue : " + objectiveFunction.value(benchSol.toArray())); assertEquals(benchSol.getEntry(0), sol[0], 0.00000000000001); assertEquals(benchSol.getEntry(1), sol[1], 0.00000000000001); assertEquals(benchSol.getEntry(2), sol[2], 0.00000000000001); }
From source file:com.joptimizer.solvers.KKTSolver.java
/** * Solve the KKT system as a whole./*ww w . j a v a 2 s . c om*/ * Useful only if A not null. * @see "S.Boyd and L.Vandenberghe, Convex Optimization, p. 547" */ protected double[][] solveFullKKT(KKTSolver kktSolver) throws Exception { Log.d(MainActivity.JOPTIMIZER_LOGTAG, "solveFullKKT"); //if the KKT matrix is nonsingular, then H + [A]T.A > 0. RealMatrix HATA = H.add(AT.multiply(A)); try { CholeskyDecomposition cFact = new CholeskyDecomposition(HATA); cFact.getSolver().getInverse(); } catch (Exception e) { throw new Exception("singular KKT system"); } kktSolver.setHMatrix(HATA.getData());//this is positive kktSolver.setAMatrix(A.getData()); kktSolver.setATMatrix(AT.getData()); kktSolver.setGVector(g.toArray()); if (h != null) { RealVector ATQh = AT.operate(MatrixUtils.createRealIdentityMatrix(A.getRowDimension()).operate(h)); RealVector gATQh = g.add(ATQh); kktSolver.setGVector(gATQh.toArray()); kktSolver.setHVector(h.toArray()); } return kktSolver.solve(); }
From source file:com.joptimizer.functions.SDPLogarithmicBarrier.java
/** * @see "S.Boyd and L.Vandenberghe, Convex Optimization, p. 618" *//* w w w . j a v a 2s . c o m*/ public double[] gradient(double[] X) { double[] ret = new double[dim]; RealMatrix S = buildS(X); CholeskyDecomposition cFact = new CholeskyDecomposition(S); RealMatrix SInv = cFact.getSolver().getInverse(); for (int i = 0; i < dim; i++) { ret[i] = SInv.multiply(this.Fi[i]).getTrace(); } return ret; }
From source file:com.joptimizer.functions.SDPLogarithmicBarrier.java
/** * @see "S.Boyd and L.Vandenberghe, Convex Optimization, p. 618" */// w w w. j a v a2 s. c om public double[][] hessian(double[] X) { RealMatrix S = buildS(X); CholeskyDecomposition cFact = new CholeskyDecomposition(S); RealMatrix SInv = cFact.getSolver().getInverse(); double[][] ret = new double[dim][dim]; for (int i = 0; i < dim; i++) { for (int j = i; j < dim; j++) { double h = SInv.multiply(this.Fi[i]).multiply(SInv.multiply(this.Fi[j])).getTrace(); ret[i][j] = h; ret[j][i] = h; } } return ret; }