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

Introduction

In this page you can find the example usage for org.apache.commons.math3.linear SingularValueDecomposition 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:com.idylwood.utils.OptimizationUtils.java

/**
 * Solves for Markowitz optimal portfolio by means of lagrange multiplier.
 * Returns null if it does not exist (the matrix is singular)
 * Otherwise it attempts to find the lowest variance portfolio for
 * the given <code>portfolio_return</code>
 * @param covariance Precalculated covariance matrix
 * @param returns Precalculated vector of returns
 * @param portfolio_return Return to optimize risk for
 * @author Harry C Kim/*from   ww w .ja va 2s .c  om*/
 * @return
 */
static final double[] MarkowitzSolve(final double[][] covariance, final double[] returns,
        final double portfolio_return) {
    if (covariance.length != covariance[0].length)
        throw new IllegalArgumentException("Covariance needs to be square matrix");
    if (returns.length != covariance.length)
        throw new IllegalArgumentException("Returns must be same length as covariance");

    /*
    for (int i = 0; i < covariance.length; i++)
       MathUtils.printArray(covariance[i]);
    System.out.println();
    MathUtils.printArray(returns);
    System.out.println();
    */

    final int timePoints = covariance.length;
    final double[][] lagrangeMatrix = new double[timePoints + 2][timePoints + 2];
    //b as in Ax = b
    final double[] b = new double[timePoints + 2];
    for (int i = 0; i < timePoints; i++) {
        for (int j = 0; j < timePoints; j++) {
            lagrangeMatrix[i][j] = 2 * covariance[i][j];
            b[i] = 0;
            // this is like riskTolerance*returns[i]; but since
            // returns[i]*weights[i] = portfolio_return it will go away in the derivative
        }
    }

    for (int j = 0; j < timePoints; j++) {
        lagrangeMatrix[timePoints][j] = returns[j];
        lagrangeMatrix[timePoints + 1][j] = 1;
        lagrangeMatrix[j][timePoints] = returns[j];
        lagrangeMatrix[j][timePoints + 1] = 1;
    }
    b[timePoints] = portfolio_return; //**** what is the constraint on total expected return?
    b[timePoints + 1] = 1;

    /*
    // Print out lagrangeMatrix augmented with b vector
    for(int i=0; i<timePoints+2; i++)
    {
       for(int j=0; j<timePoints+2;j++)
       {
    System.out.print(lagrangeMatrix[i][j] + " ");
       }
       System.out.println(b[i]);
    }
    */
    // TODO use Gaussian elimination solver, may be faster
    // TODO maybe refactor to use idylblas
    RealMatrix lagrangeReal = MatrixUtils.createRealMatrix(lagrangeMatrix);
    RealVector bReal = MatrixUtils.createRealVector(b);
    SingularValueDecomposition svd = new SingularValueDecomposition(lagrangeReal);

    if (!svd.getSolver().isNonSingular())
        return null;

    RealVector solution = svd.getSolver().solve(bReal);

    final double weights[] = new double[timePoints];

    // last two elements of solution are just lagrange multipliers
    for (int i = 0; i < weights.length; i++)
        weights[i] = solution.getEntry(i);

    // put these in some test class
    if (!MathUtils.fuzzyEquals(1, MathUtils.sum(weights)))
        throw new RuntimeException();
    if (!MathUtils.fuzzyEquals(portfolio_return, MathUtils.linearCombination(returns, weights)))
        throw new RuntimeException();
    //The following calculates the risk(variance) for the weights found
    // final double risk = MathUtils.linearCombination(MathUtils.matrixMultiply(covariance, weights),weights);
    return weights;
}

---

From source file:hivemall.utils.math.StatsUtils.java

/**
 * pdf(x, x_hat) = exp(-0.5 * (x-x_hat) * inv() * (x-x_hat)T) / ( 2^0.5d * det()^0.5)
 * //from w w  w .j av a 2s .c  om
 * @return value of probabilistic density function
 * @link https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Density_function
 */
public static double pdf(@Nonnull final RealVector x, @Nonnull final RealVector x_hat,
        @Nonnull final RealMatrix sigma) {
    final int dim = x.getDimension();
    Preconditions.checkArgument(x_hat.getDimension() == dim,
            "|x| != |x_hat|, |x|=" + dim + ", |x_hat|=" + x_hat.getDimension());
    Preconditions.checkArgument(sigma.getRowDimension() == dim,
            "|x| != |sigma|, |x|=" + dim + ", |sigma|=" + sigma.getRowDimension());
    Preconditions.checkArgument(sigma.isSquare(), "Sigma is not square matrix");

    LUDecomposition LU = new LUDecomposition(sigma);
    final double detSigma = LU.getDeterminant();
    double denominator = Math.pow(2.d * Math.PI, 0.5d * dim) * Math.pow(detSigma, 0.5d);
    if (denominator == 0.d) { // avoid divide by zero
        return 0.d;
    }

    final RealMatrix invSigma;
    DecompositionSolver solver = LU.getSolver();
    if (solver.isNonSingular() == false) {
        SingularValueDecomposition svd = new SingularValueDecomposition(sigma);
        invSigma = svd.getSolver().getInverse(); // least square solution
    } else {
        invSigma = solver.getInverse();
    }
    //EigenDecomposition eigen = new EigenDecomposition(sigma);
    //double detSigma = eigen.getDeterminant();
    //RealMatrix invSigma = eigen.getSolver().getInverse();

    RealVector diff = x.subtract(x_hat);
    RealVector premultiplied = invSigma.preMultiply(diff);
    double sum = premultiplied.dotProduct(diff);
    double numerator = Math.exp(-0.5d * sum);

    return numerator / denominator;
}

---

From source file:hivemall.utils.math.MatrixUtils.java

/**
 * L = A x R//from  ww w  . j av a2s.  com
 *
 * @return a matrix A that minimizes A x R - L
 */
@Nonnull
public static RealMatrix solve(@Nonnull final RealMatrix L, @Nonnull final RealMatrix R, final boolean exact)
        throws SingularMatrixException {
    LUDecomposition LU = new LUDecomposition(R);
    DecompositionSolver solver = LU.getSolver();
    final RealMatrix A;
    if (exact || solver.isNonSingular()) {
        A = LU.getSolver().solve(L);
    } else {
        SingularValueDecomposition SVD = new SingularValueDecomposition(R);
        A = SVD.getSolver().solve(L);
    }
    return A;
}

---

From source file:hivemall.utils.math.MatrixUtils.java

@Nonnull
public static RealMatrix inverse(@Nonnull final RealMatrix m, final boolean exact)
        throws SingularMatrixException {
    LUDecomposition LU = new LUDecomposition(m);
    DecompositionSolver solver = LU.getSolver();
    final RealMatrix inv;
    if (exact || solver.isNonSingular()) {
        inv = solver.getInverse();/*from w ww. j a v a  2s. co m*/
    } else {
        SingularValueDecomposition SVD = new SingularValueDecomposition(m);
        inv = SVD.getSolver().getInverse();
    }
    return inv;
}

---

From source file:com.opengamma.strata.math.impl.matrix.CommonsMatrixAlgebra.java

@Override
public DoubleMatrix getInverse(Matrix m) {
    ArgChecker.notNull(m, "matrix was null");
    if (m instanceof DoubleMatrix) {
        RealMatrix temp = CommonsMathWrapper.wrap((DoubleMatrix) m);
        SingularValueDecomposition sv = new SingularValueDecomposition(temp);
        RealMatrix inv = sv.getSolver().getInverse();
        return CommonsMathWrapper.unwrap(inv);
    }/*from  ww w. j  a v a 2s . co m*/
    throw new IllegalArgumentException("Can only find inverse of DoubleMatrix; have " + m.getClass());
}

---

From source file:com.github.kingtim1.jmdp.discounted.MatrixInversePolicyEvaluation.java

@Override
public DiscountedVFunction<S> eval(StationaryPolicy<S, A> policy) {
    int n = _smdp.numberOfStates();
    List<S> states = new ArrayList<S>(n);
    Iterable<S> istates = _smdp.states();
    for (S state : istates) {
        states.add(state);/*w w w.j a va2s  . co m*/
    }

    // Construct matrix A and vector b
    RealMatrix id = MatrixUtils.createRealIdentityMatrix(n);
    RealMatrix gpp = gammaPPi(policy, states);
    RealMatrix A = id.subtract(gpp);
    RealVector b = rPi(policy, states);

    // Solve for V^{\pi}
    SingularValueDecomposition decomp = new SingularValueDecomposition(A);
    DecompositionSolver dsolver = decomp.getSolver();
    RealVector vpi = dsolver.solve(b);

    // Construct the value function
    Map<S, Double> valueMap = new HashMap<S, Double>();
    for (int i = 0; i < states.size(); i++) {
        S state = states.get(i);
        double val = vpi.getEntry(i);
        valueMap.put(state, val);
    }

    return new MapVFunction<S>(valueMap, 0);
}

---

From source file:gamlss.utilities.WLSMultipleLinearRegression.java

/**
 * Calculates the fitted values without SAVING the 'beta' for later use
 * @param isSVD whether to use SVD or QR for the design matrix X
 * Note that SVD is more accurate but slower. 
 * @return the fitted values or null if not applicable
 *//*ww w .  jav  a2 s. c  o  m*/
public RealVector calculateFittedValues(boolean isSVD) {
    if (this.copyOriginal) {
        if (!isSVD) {
            return this.X.operate(this.calculateBeta());

        } else {
            SingularValueDecomposition svg = new SingularValueDecomposition(getX());

            return this.X.operate(svg.getSolver().solve(getY()));
        }
    }
    return null;
}

---

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

/**
 * @param svd The result of the SV decomposition, not null
 *///from  w  w  w. j av a  2 s  .  co  m
public SVDecompositionCommonsResult(SingularValueDecomposition svd) {
    ArgChecker.notNull(svd, "svd");
    _condition = svd.getConditionNumber();
    _norm = svd.getNorm();
    _rank = svd.getRank();
    _s = CommonsMathWrapper.unwrap(svd.getS());
    _singularValues = svd.getSingularValues();
    _u = CommonsMathWrapper.unwrap(svd.getU());
    _uTranspose = CommonsMathWrapper.unwrap(svd.getUT());
    _v = CommonsMathWrapper.unwrap(svd.getV());
    _vTranspose = CommonsMathWrapper.unwrap(svd.getVT());
    _solver = svd.getSolver();
}

---

From source file:edu.cmu.tetrad.search.IndTestMultiFisherZ.java

public boolean isIndependent(Node x, Node y, List<Node> z) {
    if (verbose) {
        System.out.println("\n" + x + " _||_ " + y + " | " + z);
    }/*ww  w  .  j a  va2  s  .com*/

    List<Node> aa = nodeMap.get(x);
    List<Node> bb = nodeMap.get(y);

    List<Node> cc = new ArrayList<Node>();

    for (Node _z : z) {
        cc.addAll(nodeMap.get(_z));
    }

    TetradMatrix submatrix = subMatrix(cov, aa, bb, cc);

    TetradMatrix inverse;
    int rank;

    try {
        inverse = submatrix.inverse();
        rank = inverse.columns();
    } catch (Exception e) {
        SingularValueDecomposition svd = new SingularValueDecomposition(submatrix.getRealMatrix());
        RealMatrix _inverse = svd.getSolver().getInverse();
        inverse = new TetradMatrix(_inverse, _inverse.getRowDimension(), _inverse.getColumnDimension());
        rank = svd.getRank();
    }

    final List<Double> pValues = new ArrayList<Double>();
    List<Integer> _i = new ArrayList<Integer>();
    List<Integer> _m = new ArrayList<Integer>();

    TetradMatrix stdSubmatrix = TestHippocampusUtils.subset(stdData, aa, bb, cc);

    for (int i = 0; i < aa.size(); i++) {
        for (int m = 0; m < bb.size(); m++) {
            int j = aa.size() + m;
            double a = -1.0 * inverse.get(i, j);
            double v0 = inverse.get(i, i);
            double v1 = inverse.get(j, j);
            double b = sqrt(v0 * v1);

            double r = a / b;

            int dof = cov.getSampleSize() - 1 - rank;

            if (dof < 0) {
                //                    out.println("Negative dof: " + dof + " n = " + cov.getSampleSize() + " cols = " + inverse.columns());
                dof = 0;
            }

            double t2 = moment22(stdSubmatrix, i, j);

            double z_ = sqrt(dof) * 0.5 * (Math.log(1.0 + r) - Math.log(1.0 - r));
            double p = 2.0 * (1.0 - RandomUtil.getInstance().normalCdf(0, sqrt(t2), abs(z_)));

            //                if (p < alpha) p = 0;

            //                System.out.println("z = " + z_ + " p = " + p);

            pValues.add(p);
            _i.add(i);
            _m.add(m);
        }
    }

    List<Integer> indices = new ArrayList<Integer>();
    for (int i = 0; i < pValues.size(); i++) {
        indices.add(i);
    }

    Collections.sort(indices, new Comparator<Integer>() {
        @Override
        public int compare(Integer o1, Integer o2) {
            return pValues.get(o1).compareTo(pValues.get(o2));
        }
    });

    List<Double> pValues2 = new ArrayList<Double>();
    List<Integer> _iIndices = new ArrayList<Integer>();
    List<Integer> _mIndices = new ArrayList<Integer>();

    for (int _y = 0; _y < indices.size(); _y++) {
        pValues2.add(pValues.get(indices.get(_y)));
        _iIndices.add(_i.get(indices.get(_y)));
        _mIndices.add(_m.get(indices.get(_y)));
    }

    int k = StatUtils.fdr(alpha, pValues2, false);
    double cutoff = StatUtils.fdrCutoff(alpha, pValues2, false);
    System.out.println("** cutoff = " + cutoff);

    int nonzero = -1;

    for (int i = 0; i < pValues2.size(); i++) {
        if (pValues2.get(i) != 0) {
            nonzero = i;
            break;
        }
    }

    double q = StatUtils.fdrQ(pValues2, nonzero);
    System.out.println("Q = " + q);

    //        boolean independent = q >= alpha;

    ////        boolean independent = k > nonzero;
    //
    System.out.println("k = " + k);
    System.out.println("nonzero = " + nonzero);
    System.out.println("fisher = " + fisherMethodP(pValues2));
    //
    //        boolean independent = k <= nonzero;

    double ratio = ((double) k) / pValues2.size();
    System.out.println("ratio = " + ratio);
    boolean independent = ratio < 0.1;
    //
    System.out.println(independent ? "Independent" : "Dependent");

    return independent;
}

---

From source file:edu.cmu.tetrad.search.IndTestMultiCci.java

public boolean isIndependent(Node x, Node y, List<Node> z) {
    if (verbose) {
        System.out.println("\n" + x + " _||_ " + y + " | " + z);
        out.println("\n" + x + " _||_ " + y + " | " + z);
    }/*from   w w w  .ja v  a2 s .c  o m*/

    List<Node> aa = nodeMap.get(x);
    List<Node> bb = nodeMap.get(y);

    //        int[][] twod = new int[aa.size()][bb.size()];

    List<Node> cc = new ArrayList<Node>();

    for (Node _z : z) {
        cc.addAll(nodeMap.get(_z));
    }

    TetradMatrix submatrix = TestHippocampusUtils.subMatrix(cov, aa, bb, cc);

    TetradMatrix inverse;
    int rank;

    try {
        inverse = submatrix.inverse();
        rank = inverse.columns();
    } catch (Exception e) {
        SingularValueDecomposition svd = new SingularValueDecomposition(submatrix.getRealMatrix());
        RealMatrix _inverse = svd.getSolver().getInverse();
        inverse = new TetradMatrix(_inverse, _inverse.getRowDimension(), _inverse.getColumnDimension());
        rank = svd.getRank();
    }

    final List<Double> pValues = new ArrayList<Double>();
    List<Integer> _i = new ArrayList<Integer>();
    List<Integer> _m = new ArrayList<Integer>();

    System.out.println("# voxels for " + x.getName() + " = " + aa.size());
    System.out.println("# voxels for " + y.getName() + " = " + bb.size());
    System.out.println("# p values = " + aa.size() * bb.size());

    IndTestConditionalCorrelation cciTest = new IndTestConditionalCorrelation(dataSet, alpha);

    for (int i = 0; i < aa.size(); i++) {
        for (int m = 0; m < bb.size(); m++) {
            int j = aa.size() + m;
            double a = -1.0 * inverse.get(i, j);
            double v0 = inverse.get(i, i);
            double v1 = inverse.get(j, j);
            double b = Math.sqrt(v0 * v1);

            double r = a / b;

            int dof = cov.getSampleSize() - 1 - rank;

            if (dof < 0) {
                out.println("Negative dof: " + dof + " n = " + cov.getSampleSize() + " cols = "
                        + inverse.columns());
                dof = 0;
            }

            double z_ = Math.sqrt(dof) * 0.5 * (Math.log(1.0 + r) - Math.log(1.0 - r));
            double p = 2.0 * (1.0 - RandomUtil.getInstance().normalCdf(0, 1, abs(z_)));

            cciTest.isIndependent(aa.get(i), bb.get(m), cc);
            pValues.add(p);

            if (m == 0) {
                System.out.println("i = " + i + " m = " + m + " p = " + p);
            }

            _i.add(i);
            _m.add(m);
        }
    }

    List<Integer> indices = new ArrayList<Integer>();
    for (int i = 0; i < pValues.size(); i++) {
        indices.add(i);
    }

    Collections.sort(indices, new Comparator<Integer>() {
        @Override
        public int compare(Integer o1, Integer o2) {
            return pValues.get(o1).compareTo(pValues.get(o2));
        }
    });

    List<Double> pValues2 = new ArrayList<Double>();
    List<Integer> _iIndices = new ArrayList<Integer>();
    List<Integer> _mIndices = new ArrayList<Integer>();

    for (int _y = 0; _y < indices.size(); _y++) {
        pValues2.add(pValues.get(indices.get(_y)));
        _iIndices.add(_i.get(indices.get(_y)));
        _mIndices.add(_m.get(indices.get(_y)));
    }

    int k = StatUtils.fdr(alpha, pValues2, false);
    double cutoff = StatUtils.fdrCutoff(alpha, pValues2, false);

    int nonzero = -1;

    for (int i = 0; i < pValues2.size(); i++) {
        if (pValues2.get(i) != 0) {
            nonzero = i;
            break;
        }
    }

    boolean dependent = k > nonzero;
    boolean independent = !dependent;

    TetradMatrix xCoords = coords.get(x);
    TetradMatrix yCoords = coords.get(y);

    int[][][] X = TestHippocampusUtils.threeDView(_iIndices, k, xCoords, independent, nonzero);
    int[][][] Y = TestHippocampusUtils.threeDView(_mIndices, k, yCoords, independent, nonzero);

    all3D.add(X);
    all3D.add(Y);

    String fact = SearchLogUtils.independenceFact(x, y, z);

    // Printing out stuff for Ruben. First print out dependent voxels.

    //        2.Second file, contain a list of all the dependencies between voxels.
    //
    //        10 -- 50
    //        30 -- 2

    int[][][] Cx = getC(X);
    int[][][] Cy = getC(Y);

    out.println("\n\n" + fact);

    for (int g = independent ? nonzero : 0; g < k; g++) {
        int i = getIndex(_iIndices, Cx, g, coords.get(x));
        int j = getIndex(_mIndices, Cy, g, coords.get(y));

        if (i == -1 || j == -1)
            throw new IllegalArgumentException();

        out.println(i + " -- " + j);
    }

    out.println();

    //        1. First file, containing info of both ROIs and all their voxels.
    //        Example:
    //
    //        ROI_LABEL  voxel_LABEL  COORDINATES  #Dependencies
    //        ENT          10         -80 50 38     6
    //        CA1          50         -70 15 90     2

    printDependencies(x, fact, X, Cx);
    printDependencies(y, fact, Y, Cy);

    // OK back to work.
    int xCount = countAboveThreshold(X, threshold);
    int yCount = countAboveThreshold(Y, threshold);

    System.out.println("Total above threshold count = " + (xCount + yCount));
    out.println("Total above threshold count = " + (xCount + yCount));

    boolean thresholdIndep = !(xCount > 0 && yCount > 0);

    String projection;

    projection = "Axial";

    TestHippocampusUtils.printChart(X, xCoords, 0, 1, x.getName(), projection, fact, false, false, out,
            threshold);
    TestHippocampusUtils.printChart(Y, yCoords, 0, 1, y.getName(), projection, fact, false, false, out,
            threshold);

    projection = "Coronal";

    TestHippocampusUtils.printChart(X, xCoords, 0, 2, x.getName(), projection, fact, true, false, out,
            threshold);
    TestHippocampusUtils.printChart(Y, yCoords, 0, 2, y.getName(), projection, fact, true, false, out,
            threshold);

    projection = "Saggital";

    TestHippocampusUtils.printChart(X, xCoords, 1, 2, x.getName(), projection, fact, true, false, out,
            threshold);
    TestHippocampusUtils.printChart(Y, yCoords, 1, 2, y.getName(), projection, fact, true, false, out,
            threshold);

    if (thresholdIndep) {
        //        if (independent) {
        if (verbose) {
            System.out.println("Independent");
            out.println("Independent");
        }

        out.flush();
        return true;
    } else {
        if (verbose) {
            System.out.println("Dependent\n");
            out.println("Dependent\n");
        }

        out.flush();
        return false;
    }
}

---

• «
• 1
• 2
• »