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

Example 11 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());
}

Example 12 with DecompositionSolver

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

the class LibCommonsMath method computeSolve.

/**
 * Function to solve a given system of equations.
 *
 * @param in1 matrix object 1
 * @param in2 matrix object 2
 * @return matrix block
 */
private static MatrixBlock computeSolve(MatrixObject in1, MatrixObject in2) {
    Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in1);
    Array2DRowRealMatrix vectorInput = DataConverter.convertToArray2DRowRealMatrix(in2);
    /*LUDecompositionImpl ludecompose = new LUDecompositionImpl(matrixInput);
		DecompositionSolver lusolver = ludecompose.getSolver();
		RealMatrix solutionMatrix = lusolver.solve(vectorInput);*/
    // Setup a solver based on QR Decomposition
    QRDecomposition qrdecompose = new QRDecomposition(matrixInput);
    DecompositionSolver solver = qrdecompose.getSolver();
    // Invoke solve
    RealMatrix solutionMatrix = solver.solve(vectorInput);
    return DataConverter.convertToMatrixBlock(solutionMatrix.getData());
}

Example 13 with DecompositionSolver

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

the class IrlsLearner method irlsRls.

/**
 * Do an 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 TrainingData<ClassificationTrainingRow> data, final RealMatrix beta, final int rC, final int tcC, final ExecutionMonitor exec) throws CanceledExecutionException {
    long rowCount = 0;
    int dim = (rC + 1) * (tcC - 1);
    RealMatrix xTwx = MatrixUtils.createRealMatrix(dim, dim);
    RealMatrix xTyu = MatrixUtils.createRealMatrix(dim, 1);
    double[] eBetaTx = new double[tcC - 1];
    double[] pi = new double[tcC - 1];
    final long totalRowCount = data.getRowCount();
    for (ClassificationTrainingRow row : data) {
        rowCount++;
        exec.checkCanceled();
        exec.setProgress(rowCount / (double) totalRowCount, "Row " + rowCount + "/" + totalRowCount);
        for (int k = 0; k < tcC - 1; k++) {
            double z = 0.0;
            for (FeatureIterator iter = row.getFeatureIterator(); iter.next(); ) {
                double featureVal = iter.getFeatureValue();
                int featureIdx = iter.getFeatureIndex();
                z += featureVal * beta.getEntry(0, k * (rC + 1) + featureIdx);
            }
            eBetaTx[k] = Math.exp(z);
        }
        double sumEBetaTx = 0;
        for (int k = 0; k < tcC - 1; k++) {
            sumEBetaTx += eBetaTx[k];
        }
        for (int k = 0; k < tcC - 1; k++) {
            double pik = eBetaTx[k] / (1 + sumEBetaTx);
            pi[k] = pik;
        }
        // fill xTwx (aka the hessian of the loglikelihood)
        for (FeatureIterator outer = row.getFeatureIterator(); outer.next(); ) {
            int i = outer.getFeatureIndex();
            double outerVal = outer.getFeatureValue();
            for (FeatureIterator inner = outer.spawn(); inner.next(); ) {
                int ii = inner.getFeatureIndex();
                double innerVal = inner.getFeatureValue();
                for (int k = 0; k < tcC - 1; k++) {
                    for (int kk = k; kk < tcC - 1; kk++) {
                        int o1 = k * (rC + 1);
                        int o2 = kk * (rC + 1);
                        double v = xTwx.getEntry(o1 + i, o2 + ii);
                        if (k == kk) {
                            double w = pi[k] * (1 - pi[k]);
                            v += outerVal * w * innerVal;
                            assert o1 == o2;
                        } else {
                            double w = -pi[k] * pi[kk];
                            v += outerVal * w * innerVal;
                        }
                        xTwx.setEntry(o1 + i, o2 + ii, v);
                        xTwx.setEntry(o1 + ii, o2 + i, v);
                        if (k != kk) {
                            xTwx.setEntry(o2 + ii, o1 + i, v);
                            xTwx.setEntry(o2 + i, o1 + ii, v);
                        }
                    }
                }
            }
        }
        int g = row.getCategory();
        // fill matrix xTyu
        for (FeatureIterator iter = row.getFeatureIterator(); iter.next(); ) {
            int idx = iter.getFeatureIndex();
            double val = iter.getFeatureValue();
            for (int k = 0; k < tcC - 1; k++) {
                int o = k * (rC + 1);
                double v = xTyu.getEntry(o + idx, 0);
                double y = k == g ? 1 : 0;
                v += (y - pi[k]) * val;
                xTyu.setEntry(o + idx, 0, v);
            }
        }
    }
    // currently not used but could become interesting in the future
    // 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()) {
        // but it's important to ensure this property
        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();
    RealMatrix betaNew = solver.solve(b);
    beta.setSubMatrix(betaNew.transpose().getData(), 0, 0);
}

Example 14 with DecompositionSolver

use of org.apache.commons.math3.linear.DecompositionSolver in project imagingbook-common by imagingbook.

the class AffineFit method fit.

@Override
public void fit(List<double[]> X, List<double[]> Y) {
    // fits n-dimensional data sets with affine model
    if (X.size() != Y.size())
        throw new IllegalArgumentException("point sequences X, Y must have same length");
    this.m = X.size();
    this.n = X.get(0).length;
    RealMatrix M = MatrixUtils.createRealMatrix(2 * m, 2 * (n + 1));
    RealVector b = new ArrayRealVector(2 * m);
    // mount matrix M:
    int row = 0;
    for (double[] x : X) {
        for (int j = 0; j < n; j++) {
            M.setEntry(row, j, x[j]);
            M.setEntry(row, n, 1);
            row++;
        }
        for (int j = 0; j < n; j++) {
            M.setEntry(row, j + n + 1, x[j]);
            M.setEntry(row, 2 * n + 1, 1);
            row++;
        }
    }
    // mount vector b
    row = 0;
    for (double[] y : Y) {
        for (int j = 0; j < n; j++) {
            b.setEntry(row, y[j]);
            row++;
        }
    }
    SingularValueDecomposition svd = new SingularValueDecomposition(M);
    DecompositionSolver solver = svd.getSolver();
    RealVector a = solver.solve(b);
    A = makeTransformationMatrix(a);
}

Example 15 with DecompositionSolver

use of org.apache.commons.math3.linear.DecompositionSolver in project imagingbook-common by imagingbook.

the class Matrix method solve.

// use general method, i.e. double[][] inverse(double[][] A)
// @Deprecated
// public static double[][] inverse3x3(final double[][] A) {
// double[][] B = duplicate(A);
// final double det = determinant3x3(B);
// if (Math.abs(det) < Arithmetic.EPSILON_DOUBLE)
// return null;
// else {
// final double a00 = B[0][0];
// final double a01 = B[0][1];
// final double a02 = B[0][2];
// final double a10 = B[1][0];
// final double a11 = B[1][1];
// final double a12 = B[1][2];
// final double a20 = B[2][0];
// final double a21 = B[2][1];
// final double a22 = B[2][2];
// B[0][0] = (a11 * a22 - a12 * a21) / det;
// B[0][1] = (a02 * a21 - a01 * a22) / det;
// B[0][2] = (a01 * a12 - a02 * a11) / det;
// 
// B[1][0] = (a12 * a20 - a10 * a22) / det;
// B[1][1] = (a00 * a22 - a02 * a20) / det;
// B[1][2] = (a02 * a10 - a00 * a12) / det;
// 
// B[2][0] = (a10 * a21 - a11 * a20) / det;
// B[2][1] = (a01 * a20 - a00 * a21) / det;
// B[2][2] = (a00 * a11 - a01 * a10) / det;
// return B;
// }
// }
// numerically stable? should be replaced by standard inversion
// @Deprecated
// public static float[][] inverse3x3(final float[][] A) {
// final float[][] B = duplicate(A);
// final double det = determinant3x3(B);
// // IJ.log("   determinant = " + det);
// if (Math.abs(det) < Arithmetic.EPSILON_DOUBLE)
// return null;
// else {
// final double a00 = B[0][0];
// final double a01 = B[0][1];
// final double a02 = B[0][2];
// final double a10 = B[1][0];
// final double a11 = B[1][1];
// final double a12 = B[1][2];
// final double a20 = B[2][0];
// final double a21 = B[2][1];
// final double a22 = B[2][2];
// B[0][0] = (float) ((a11 * a22 - a12 * a21) / det);
// B[0][1] = (float) ((a02 * a21 - a01 * a22) / det);
// B[0][2] = (float) ((a01 * a12 - a02 * a11) / det);
// 
// B[1][0] = (float) ((a12 * a20 - a10 * a22) / det);
// B[1][1] = (float) ((a00 * a22 - a02 * a20) / det);
// B[1][2] = (float) ((a02 * a10 - a00 * a12) / det);
// 
// B[2][0] = (float) ((a10 * a21 - a11 * a20) / det);
// B[2][1] = (float) ((a01 * a20 - a00 * a21) / det);
// B[2][2] = (float) ((a00 * a11 - a01 * a10) / det);
// return B;
// }
// }
// ------------------------------------------------------------------------
// Finds the EXACT solution x for A.x = b
public static double[] solve(final double[][] A, double[] b) {
    RealMatrix AA = MatrixUtils.createRealMatrix(A);
    RealVector bb = MatrixUtils.createRealVector(b);
    DecompositionSolver solver = new LUDecomposition(AA).getSolver();
    double[] x = null;
    try {
        x = solver.solve(bb).toArray();
    } catch (SingularMatrixException e) {
    }
    return x;
}