Example usage for org.apache.commons.math3.linear CholeskyDecomposition getSolver

List of usage examples for org.apache.commons.math3.linear CholeskyDecomposition getSolver

  1. HOME
  2. Java
  3. org.apache.commons
  4. org.apache.commons.math3.linear.*
  5. CholeskyDecomposition
  6. getSolver

Introduction

In this page you can find the example usage for org.apache.commons.math3.linear CholeskyDecomposition getSolver.

Prototype

public DecompositionSolver getSolver()

Source Link

Document

Get a solver for finding the A × X = B solution in least square sense.

Usage

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;/*  ww  w  .ja v  a2s .  c o  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.util.Utils.java

/**
 * @see http://en.wikipedia.org/wiki/Machine_epsilon#Approximation_using_Java
 *//*from   ww  w.  j a v  a  2 s.c o m*/
//   public static final float calculateFloatMachineEpsilon() {
//      float eps = 1.0f;
//      do {
//         eps /= 2.0f;
//      } while ((float) (1.0 + (eps / 2.0)) != 1.0);
//      Log.d(MainActivity.JOPTIMIZER_LOGTAG,"Calculated float machine epsilon: " + eps);
//      return eps;
//   }

public static RealMatrix squareMatrixInverse(RealMatrix M) throws SingularMatrixException {
    if (!M.isSquare()) {
        throw new IllegalArgumentException("Not square matrix!");
    }

    // try commons-math cholesky
    try {
        CholeskyDecomposition cd = new CholeskyDecomposition(M);
        return cd.getSolver().getInverse();
    } catch (SingularMatrixException e) {
        throw e;
    } catch (Exception e) {
        Log.d(MainActivity.JOPTIMIZER_LOGTAG, e.getMessage());
    }

    // try joptimizer cholesky
    try {
        CholeskyFactorization cd = new CholeskyFactorization(M.getData());
        double[][] MInv = cd.getInverse();
        return MatrixUtils.createRealMatrix(MInv);
    } catch (Exception e) {
        Log.d(MainActivity.JOPTIMIZER_LOGTAG, e.getMessage());
    }

    // try LU
    try {
        LUDecomposition ld = new LUDecomposition(M);
        return ld.getSolver().getInverse();
    } catch (SingularMatrixException e) {
        throw e;
    } catch (Exception e) {
        Log.d(MainActivity.JOPTIMIZER_LOGTAG, e.getMessage());
    }

    // try QR
    try {
        QRDecomposition qr = new org.apache.commons.math3.linear.QRDecomposition(M);
        return qr.getSolver().getInverse();
    } catch (SingularMatrixException e) {
        throw e;
    } catch (Exception e) {
        Log.d(MainActivity.JOPTIMIZER_LOGTAG, e.getMessage());
    }

    return null;
}

From source file:com.joptimizer.solvers.KKTSolver.java

/**
 * Solve the KKT system as a whole./*from w w  w  .j  ava 2  s  .  c  o  m*/
 * 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.optimizers.NewtonUnconstrainedTest.java

/**
 * Quadratic objective.//w  w w. java  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.functions.SDPLogarithmicBarrier.java

/**
 * @see "S.Boyd and L.Vandenberghe, Convex Optimization, p. 618"
 *///from  w w  w  .j a v a 2 s  .  co  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"
 *//*from   w w w.  ja  v  a  2 s . c  o m*/
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;
}

From source file:com.opengamma.strata.math.impl.linearalgebra.CholeskyDecompositionCommonsResult.java

/**
 * Constructor.// w  w  w. ja  v a2 s  .  c  om
 * @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();
}

From source file:com.linkedin.mlease.regression.liblinearfunc.LibLinear.java

public void train(LibLinearDataset dataset, Map<String, Double> initParam, Map<String, Double> priorMean,
        Map<String, Double> priorVar, double defaultPriorMean, double defaultPriorVar, String option,
        boolean computePosteriorVar) throws Exception {
    if (!dataset.isFinished())
        throw new IOException("Cannot train a model using unfinished dataset");
    bias = dataset.bias;/* w ww.  ja  v a 2  s.c  om*/
    parseOption(option);

    // setup initial parameter vector
    param = new double[dataset.nFeatures()];
    initSetup(param, initParam, dataset, 0.0);

    // setup prior mean
    this.priorMean = new double[dataset.nFeatures()];
    initSetup(this.priorMean, priorMean, dataset, defaultPriorMean);

    // setup prior var
    this.priorVar = new double[dataset.nFeatures()];
    initSetup(this.priorVar, priorVar, dataset, defaultPriorVar);

    // setup initial posterior variance
    postVar = null;
    postVarMap = null;
    postVarMatrix = null;
    postVarMatrixMap = null;
    if (computePosteriorVar) {
        // initialize the diagonal posterior variance
        postVar = new double[this.priorVar.length];
        for (int i = 0; i < postVar.length; i++)
            postVar[i] = this.priorVar[i];

        if (computeFullPostVar) {
            // initialize the full posterior variance matrix
            postVarMatrix = new double[postVar.length][];
            for (int i = 0; i < postVarMatrix.length; i++) {
                postVarMatrix[i] = new double[postVarMatrix.length];
                for (int j = 0; j < postVarMatrix.length; j++) {
                    if (i == j)
                        postVarMatrix[i][j] = this.priorVar[i];
                    else
                        postVarMatrix[i][j] = 0;
                }
            }
        }
    }

    int pos = 0;
    for (int i = 0; i < dataset.nInstances(); i++)
        if (dataset.y[i] == 1)
            pos++;
    int neg = dataset.nInstances() - pos;

    if (type.equals(Logistic_L2_primal)) {
        double multiplier = 1;
        LibLinearFunction func;

        if (dataset instanceof LibLinearBinaryDataset) {
            func = new LogisticRegressionL2BinaryFeature((LibLinearBinaryDataset) dataset, this.priorMean,
                    this.priorVar, multiplier, positive_weight, 1);
        } else {
            func = new LogisticRegressionL2(dataset, this.priorMean, this.priorVar, multiplier, positive_weight,
                    1);
        }
        if (reporter != null) {
            reporter.setStatus("Start LibLinear of type " + type);
            func.setReporter(reporter, reportFrequency);
        }

        // Find the posterior mode
        Tron tron = new Tron(func, epsilon * Math.min(pos, neg) / dataset.nInstances(), max_iter);
        tron.tron(param);

        // Compute the posterior variance
        if (computePosteriorVar) {
            if (computeFullPostVar) {
                // Compute the full posterior variance matrix
                func.hessian(param, postVarMatrix);
                RealMatrix H = MatrixUtils.createRealMatrix(postVarMatrix);
                CholeskyDecomposition decomp = new CholeskyDecomposition(H);
                DecompositionSolver solver = decomp.getSolver();
                RealMatrix Var = solver.getInverse();
                postVarMatrix = Var.getData();
                for (int i = 0; i < postVar.length; i++)
                    postVar[i] = postVarMatrix[i][i];
            } else {
                // Compute the diagonal elements of the variance
                func.hessianDiagonal(param, postVar);
                for (int i = 0; i < postVar.length; i++)
                    postVar[i] = 1.0 / postVar[i];
            }
        }
    } else if (type.equals(Do_nothing)) {
        // do nothing
    } else
        throw new Exception("Unknown type: " + type);

    coeff = new HashMap<String, Double>();
    if (computePosteriorVar)
        postVarMap = new HashMap<String, Double>();
    for (int index = 1; index <= dataset.nFeatures(); index++) {
        String featureName = dataset.getFeatureName(index);
        coeff.put(featureName, param[index - 1]);
        if (computePosteriorVar)
            postVarMap.put(featureName, postVar[index - 1]);
    }

    if (computePosteriorVar && computeFullPostVar) {
        postVarMatrixMap = new HashMap<List<String>, Double>();
        for (int i = 1; i <= dataset.nFeatures(); i++) {
            String name_i = dataset.getFeatureName(i);
            for (int j = 1; j <= dataset.nFeatures(); j++) {
                double cov = postVarMatrix[i - 1][j - 1];
                if (cov != 0) {
                    String name_j = dataset.getFeatureName(j);
                    ArrayList<String> pair = new ArrayList<String>(2);
                    pair.add(name_i);
                    pair.add(name_j);
                    postVarMatrixMap.put(pair, cov);
                }
            }
        }
    }

    // check for features with non-zero prior that do not appear in the dataset
    if (priorMean != null) {
        for (String key : priorMean.keySet()) {
            if (!coeff.containsKey(key)) {
                coeff.put(key, priorMean.get(key));
            }
        }
    }
    if (priorVar != null && computePosteriorVar) {
        for (String key : priorVar.keySet()) {
            if (!postVarMap.containsKey(key))
                postVarMap.put(key, priorVar.get(key));
            if (computeFullPostVar) {
                ArrayList<String> pair = new ArrayList<String>(2);
                pair.add(key);
                pair.add(key);
                if (!postVarMatrixMap.containsKey(pair))
                    postVarMatrixMap.put(pair, priorVar.get(key));
            }
        }
    }
}

From source file:org.ojalgo.benchmark.contestant.ACM.java

@Override
public BenchmarkContestant<RealMatrix>.HermitianSolver getHermitianSolver() {
    return new HermitianSolver() {

        @Override/*from  w  ww. j  a v  a 2 s.c o  m*/
        public RealMatrix apply(final RealMatrix body, final RealMatrix rhs) {

            final CholeskyDecomposition tmpCholesky = new CholeskyDecomposition(body);

            return tmpCholesky.getSolver().solve(rhs);
        }

    };
}

From source file:org.ojalgo.benchmark.lab.library.ACM.java

@Override
public ProducingBinaryMatrixMatrixOperation<RealMatrix, Array2DRowRealMatrix> getOperationEquationSystemSolver(
        final int numbEquations, final int numbVariables, final int numbSolutions, final boolean spd) {
    if (numbEquations == numbVariables) {
        if (spd) {
            return (body, rhs) -> {
                final CholeskyDecomposition cholesky = new CholeskyDecomposition(body);
                return cholesky.getSolver().solve(rhs);
            };/*from w w  w .j  a  v  a 2  s  .com*/
        } else {
            return (body, rhs) -> {
                final LUDecomposition lu = new LUDecomposition(body);
                return lu.getSolver().solve(rhs);
            };
        }
    } else if (numbEquations > numbVariables) {
        return (body, rhs) -> {
            final QRDecomposition qr = new QRDecomposition(body);
            return qr.getSolver().solve(rhs);
        };
    } else {
        return null;
    }
}