Examples with DecompositionSolver org.apache.commons.math3.linear.DecompositionSolver used on opensource projects

Example 1 with DecompositionSolver

use of org.apache.commons.math3.linear.DecompositionSolver in project systemml by apache.

the class LibCommonsMath method computeMatrixInverse.

/**
 * Function to compute matrix inverse via matrix decomposition.
 *
 * @param in commons-math3 Array2DRowRealMatrix
 * @return matrix block
 */
private static MatrixBlock computeMatrixInverse(Array2DRowRealMatrix in) {
    if (!in.isSquare())
        throw new DMLRuntimeException("Input to inv() must be square matrix -- given: a " + in.getRowDimension() + "x" + in.getColumnDimension() + " matrix.");
    QRDecomposition qrdecompose = new QRDecomposition(in);
    DecompositionSolver solver = qrdecompose.getSolver();
    RealMatrix inverseMatrix = solver.getInverse();
    return DataConverter.convertToMatrixBlock(inverseMatrix.getData());
}

Example 2 with DecompositionSolver

use of org.apache.commons.math3.linear.DecompositionSolver in project FSensor by KalebKE.

the class FitPoints method findCenter.

/**
 * Find the offset of the ellipsoid.
 *
 * @param a the algebraic from of the polynomial.
 * @return a vector containing the offset of the ellipsoid.
 */
private RealVector findCenter(RealMatrix a) {
    RealMatrix subA = a.getSubMatrix(0, 2, 0, 2);
    for (int q = 0; q < subA.getRowDimension(); q++) {
        for (int s = 0; s < subA.getColumnDimension(); s++) {
            subA.multiplyEntry(q, s, -1.0);
        }
    }
    RealVector subV = a.getRowVector(3).getSubVector(0, 3);
    // inv (dtd)
    DecompositionSolver solver = new SingularValueDecomposition(subA).getSolver();
    RealMatrix subAi = solver.getInverse();
    return subAi.operate(subV);
}

Example 3 with DecompositionSolver

use of org.apache.commons.math3.linear.DecompositionSolver in project knime-core by knime.

the class Learner method irlsRls.

/**
 * Do a irls step. The result is stored in beta.
 *
 * @param data over trainings data.
 * @param beta parameter vector
 * @param rC regressors count
 * @param tcC target category count
 * @throws CanceledExecutionException when method is cancelled
 */
private void irlsRls(final RegressionTrainingData data, final RealMatrix beta, final int rC, final int tcC, final ExecutionMonitor exec) throws CanceledExecutionException {
    Iterator<RegressionTrainingRow> iter = data.iterator();
    long rowCount = 0;
    int dim = (rC + 1) * (tcC - 1);
    RealMatrix xTwx = new Array2DRowRealMatrix(dim, dim);
    RealMatrix xTyu = new Array2DRowRealMatrix(dim, 1);
    RealMatrix x = new Array2DRowRealMatrix(1, rC + 1);
    RealMatrix eBetaTx = new Array2DRowRealMatrix(1, tcC - 1);
    RealMatrix pi = new Array2DRowRealMatrix(1, tcC - 1);
    final long totalRowCount = data.getRowCount();
    while (iter.hasNext()) {
        rowCount++;
        RegressionTrainingRow row = iter.next();
        exec.checkCanceled();
        exec.setProgress(rowCount / (double) totalRowCount, "Row " + rowCount + "/" + totalRowCount);
        x.setEntry(0, 0, 1);
        x.setSubMatrix(row.getParameter().getData(), 0, 1);
        for (int k = 0; k < tcC - 1; k++) {
            RealMatrix betaITx = x.multiply(beta.getSubMatrix(0, 0, k * (rC + 1), (k + 1) * (rC + 1) - 1).transpose());
            eBetaTx.setEntry(0, k, Math.exp(betaITx.getEntry(0, 0)));
        }
        double sumEBetaTx = 0;
        for (int k = 0; k < tcC - 1; k++) {
            sumEBetaTx += eBetaTx.getEntry(0, k);
        }
        for (int k = 0; k < tcC - 1; k++) {
            double pik = eBetaTx.getEntry(0, k) / (1 + sumEBetaTx);
            pi.setEntry(0, k, pik);
        }
        // fill the diagonal blocks of matrix xTwx (k = k')
        for (int k = 0; k < tcC - 1; k++) {
            for (int i = 0; i < rC + 1; i++) {
                for (int ii = i; ii < rC + 1; ii++) {
                    int o = k * (rC + 1);
                    double v = xTwx.getEntry(o + i, o + ii);
                    double w = pi.getEntry(0, k) * (1 - pi.getEntry(0, k));
                    v += x.getEntry(0, i) * w * x.getEntry(0, ii);
                    xTwx.setEntry(o + i, o + ii, v);
                    xTwx.setEntry(o + ii, o + i, v);
                }
            }
        }
        // fill the rest of xTwx (k != k')
        for (int k = 0; k < tcC - 1; k++) {
            for (int kk = k + 1; kk < tcC - 1; kk++) {
                for (int i = 0; i < rC + 1; i++) {
                    for (int ii = i; ii < rC + 1; ii++) {
                        int o1 = k * (rC + 1);
                        int o2 = kk * (rC + 1);
                        double v = xTwx.getEntry(o1 + i, o2 + ii);
                        double w = -pi.getEntry(0, k) * pi.getEntry(0, kk);
                        v += x.getEntry(0, i) * w * x.getEntry(0, ii);
                        xTwx.setEntry(o1 + i, o2 + ii, v);
                        xTwx.setEntry(o1 + ii, o2 + i, v);
                        xTwx.setEntry(o2 + ii, o1 + i, v);
                        xTwx.setEntry(o2 + i, o1 + ii, v);
                    }
                }
            }
        }
        int g = (int) row.getTarget();
        // fill matrix xTyu
        for (int k = 0; k < tcC - 1; k++) {
            for (int i = 0; i < rC + 1; i++) {
                int o = k * (rC + 1);
                double v = xTyu.getEntry(o + i, 0);
                double y = k == g ? 1 : 0;
                v += (y - pi.getEntry(0, k)) * x.getEntry(0, i);
                xTyu.setEntry(o + i, 0, v);
            }
        }
    }
    if (m_penaltyTerm > 0.0) {
        RealMatrix stdError = getStdErrorMatrix(xTwx);
        // do not penalize the constant terms
        for (int i = 0; i < tcC - 1; i++) {
            stdError.setEntry(i * (rC + 1), i * (rC + 1), 0);
        }
        xTwx = xTwx.add(stdError.scalarMultiply(-0.00001));
    }
    exec.checkCanceled();
    b = xTwx.multiply(beta.transpose()).add(xTyu);
    A = xTwx;
    if (rowCount < A.getColumnDimension()) {
        throw new IllegalStateException("The dataset must have at least " + A.getColumnDimension() + " rows, but it has only " + rowCount + " rows. It is recommended to use a " + "larger dataset in order to increase accuracy.");
    }
    DecompositionSolver solver = new QRDecomposition(A).getSolver();
    boolean isNonSingular = solver.isNonSingular();
    if (isNonSingular) {
        RealMatrix betaNew = solver.solve(b);
        beta.setSubMatrix(betaNew.transpose().getData(), 0, 0);
    } else {
        throw new RuntimeException(FAILING_MSG);
    }
}

Example 4 with DecompositionSolver

use of org.apache.commons.math3.linear.DecompositionSolver in project knime-core by knime.

the class Learner method irlsRls.

/**
 * Do a irls step. The result is stored in beta.
 *
 * @param data over trainings data.
 * @param beta parameter vector
 * @param rC regressors count
 * @param tcC target category count
 * @throws CanceledExecutionException when method is cancelled
 */
private void irlsRls(final RegressionTrainingData data, final RealMatrix beta, final int rC, final int tcC, final ExecutionMonitor exec) throws CanceledExecutionException {
    Iterator<RegressionTrainingRow> iter = data.iterator();
    long rowCount = 0;
    int dim = (rC + 1) * (tcC - 1);
    RealMatrix xTwx = new Array2DRowRealMatrix(dim, dim);
    RealMatrix xTyu = new Array2DRowRealMatrix(dim, 1);
    RealMatrix x = new Array2DRowRealMatrix(1, rC + 1);
    RealMatrix eBetaTx = new Array2DRowRealMatrix(1, tcC - 1);
    RealMatrix pi = new Array2DRowRealMatrix(1, tcC - 1);
    final long totalRowCount = data.getRowCount();
    while (iter.hasNext()) {
        rowCount++;
        RegressionTrainingRow row = iter.next();
        exec.checkCanceled();
        exec.setProgress(rowCount / (double) totalRowCount, "Row " + rowCount + "/" + totalRowCount);
        x.setEntry(0, 0, 1);
        x.setSubMatrix(row.getParameter().getData(), 0, 1);
        for (int k = 0; k < tcC - 1; k++) {
            RealMatrix betaITx = x.multiply(beta.getSubMatrix(0, 0, k * (rC + 1), (k + 1) * (rC + 1) - 1).transpose());
            eBetaTx.setEntry(0, k, Math.exp(betaITx.getEntry(0, 0)));
        }
        double sumEBetaTx = 0;
        for (int k = 0; k < tcC - 1; k++) {
            sumEBetaTx += eBetaTx.getEntry(0, k);
        }
        for (int k = 0; k < tcC - 1; k++) {
            double pik = eBetaTx.getEntry(0, k) / (1 + sumEBetaTx);
            pi.setEntry(0, k, pik);
        }
        // fill the diagonal blocks of matrix xTwx (k = k')
        for (int k = 0; k < tcC - 1; k++) {
            for (int i = 0; i < rC + 1; i++) {
                for (int ii = i; ii < rC + 1; ii++) {
                    int o = k * (rC + 1);
                    double v = xTwx.getEntry(o + i, o + ii);
                    double w = pi.getEntry(0, k) * (1 - pi.getEntry(0, k));
                    v += x.getEntry(0, i) * w * x.getEntry(0, ii);
                    xTwx.setEntry(o + i, o + ii, v);
                    xTwx.setEntry(o + ii, o + i, v);
                }
            }
        }
        // fill the rest of xTwx (k != k')
        for (int k = 0; k < tcC - 1; k++) {
            for (int kk = k + 1; kk < tcC - 1; kk++) {
                for (int i = 0; i < rC + 1; i++) {
                    for (int ii = i; ii < rC + 1; ii++) {
                        int o1 = k * (rC + 1);
                        int o2 = kk * (rC + 1);
                        double v = xTwx.getEntry(o1 + i, o2 + ii);
                        double w = -pi.getEntry(0, k) * pi.getEntry(0, kk);
                        v += x.getEntry(0, i) * w * x.getEntry(0, ii);
                        xTwx.setEntry(o1 + i, o2 + ii, v);
                        xTwx.setEntry(o1 + ii, o2 + i, v);
                        xTwx.setEntry(o2 + ii, o1 + i, v);
                        xTwx.setEntry(o2 + i, o1 + ii, v);
                    }
                }
            }
        }
        int g = (int) row.getTarget();
        // fill matrix xTyu
        for (int k = 0; k < tcC - 1; k++) {
            for (int i = 0; i < rC + 1; i++) {
                int o = k * (rC + 1);
                double v = xTyu.getEntry(o + i, 0);
                double y = k == g ? 1 : 0;
                v += (y - pi.getEntry(0, k)) * x.getEntry(0, i);
                xTyu.setEntry(o + i, 0, v);
            }
        }
    }
    if (m_penaltyTerm > 0.0) {
        RealMatrix stdError = getStdErrorMatrix(xTwx);
        // do not penalize the constant terms
        for (int i = 0; i < tcC - 1; i++) {
            stdError.setEntry(i * (rC + 1), i * (rC + 1), 0);
        }
        xTwx = xTwx.add(stdError.scalarMultiply(-0.00001));
    }
    exec.checkCanceled();
    b = xTwx.multiply(beta.transpose()).add(xTyu);
    A = xTwx;
    if (rowCount < A.getColumnDimension()) {
        throw new IllegalStateException("The dataset must have at least " + A.getColumnDimension() + " rows, but it has only " + rowCount + " rows. It is recommended to use a " + "larger dataset in order to increase accuracy.");
    }
    DecompositionSolver solver = new SingularValueDecomposition(A).getSolver();
    // boolean isNonSingular = solver.isNonSingular();
    // if (isNonSingular) {
    RealMatrix betaNew = solver.solve(b);
    beta.setSubMatrix(betaNew.transpose().getData(), 0, 0);
// } else {
// throw new RuntimeException(FAILING_MSG);
// }
}

Example 5 with DecompositionSolver

use of org.apache.commons.math3.linear.DecompositionSolver in project incubator-systemml by apache.

the class LibCommonsMath method computeMatrixInverse.

/**
 * Function to compute matrix inverse via matrix decomposition.
 *
 * @param in commons-math3 Array2DRowRealMatrix
 * @return matrix block
 */
private static MatrixBlock computeMatrixInverse(Array2DRowRealMatrix in) {
    if (!in.isSquare())
        throw new DMLRuntimeException("Input to inv() must be square matrix -- given: a " + in.getRowDimension() + "x" + in.getColumnDimension() + " matrix.");
    QRDecomposition qrdecompose = new QRDecomposition(in);
    DecompositionSolver solver = qrdecompose.getSolver();
    RealMatrix inverseMatrix = solver.getInverse();
    return DataConverter.convertToMatrixBlock(inverseMatrix.getData());
}