/** * @param data dense matrix represented in row-major form * @return solver for the system Ax = b */ static Solver getSolver(double[][] data) { if (data == null) { return null; } RealMatrix M = new Array2DRowRealMatrix(data, false); double infNorm = M.getNorm(); double singularityThreshold = infNorm * SINGULARITY_THRESHOLD_RATIO; RRQRDecomposition decomposition = new RRQRDecomposition(M, singularityThreshold); DecompositionSolver solver = decomposition.getSolver(); if (solver.isNonSingular()) { return new Solver(solver); } // Otherwise try to report apparent rank int apparentRank = decomposition.getRank(0.01); // Better value? log.warn("{} x {} matrix is near-singular (threshold {}). Add more data or decrease the " + "number of features, to <= about {}", M.getRowDimension(), M.getColumnDimension(), singularityThreshold, apparentRank); throw new SingularMatrixSolverException(apparentRank, "Apparent rank: " + apparentRank); }
/** * Find the center of the ellipsoid. * * @param a the algebraic from of the polynomial. * @return a vector containing the center 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); }
/** * Find the center of the ellipsoid. * * @param a * the algebraic from of the polynomial. * @return a vector containing the center 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); }
/** * Function to solve a given system of equations. * * @param in1 matrix object 1 * @param in2 matrix object 2 * @return matrix block * @throws DMLRuntimeException if DMLRuntimeException occurs */ private static MatrixBlock computeSolve(MatrixObject in1, MatrixObject in2) throws DMLRuntimeException { 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()); }
/** * Get the covariance matrix of the optimized parameters. * <br/> * Note that this operation involves the inversion of the * <code>J<sup>T</sup>J</code> matrix, where {@code J} is the * Jacobian matrix. * The {@code threshold} parameter is a way for the caller to specify * that the result of this computation should be considered meaningless, * and thus trigger an exception. * * @param threshold Singularity threshold. * @return the covariance matrix. * @throws org.apache.commons.math3.linear.SingularMatrixException * if the covariance matrix cannot be computed (singular problem). */ public double[][] getCovariances(double threshold) { // Set up the jacobian. updateJacobian(); // Compute transpose(J)J, without building intermediate matrices. double[][] jTj = new double[cols][cols]; for (int i = 0; i < cols; ++i) { for (int j = i; j < cols; ++j) { double sum = 0; for (int k = 0; k < rows; ++k) { sum += weightedResidualJacobian[k][i] * weightedResidualJacobian[k][j]; } jTj[i][j] = sum; jTj[j][i] = sum; } } // Compute the covariances matrix. final DecompositionSolver solver = new QRDecomposition(MatrixUtils.createRealMatrix(jTj), threshold).getSolver(); return solver.getInverse().getData(); }
/** * Find a solution of the linear equation A x = b where * <ul> * <li>A is an symmetric n x n - matrix given as double[n][n]</li> * <li>b is an n - vector given as double[n],</li> * <li>x is an n - vector given as double[n],</li> * </ul> * * @param matrix The matrix A (left hand side of the linear equation). * @param vector The vector b (right hand of the linear equation). * @return A solution x to A x = b. */ public static double[] solveLinearEquationSymmetric(double[][] matrix, double[] vector) { if(isSolverUseApacheCommonsMath) { DecompositionSolver solver = new LUDecomposition(new Array2DRowRealMatrix(matrix)).getSolver(); return solver.solve(new Array2DRowRealMatrix(vector)).getColumn(0); } else { return org.jblas.Solve.solveSymmetric(new org.jblas.DoubleMatrix(matrix), new org.jblas.DoubleMatrix(vector)).data; /* To use the linear algebra package colt from cern. cern.colt.matrix.linalg.Algebra linearAlgebra = new cern.colt.matrix.linalg.Algebra(); double[] x = linearAlgebra.solve(new DenseDoubleMatrix2D(A), linearAlgebra.transpose(new DenseDoubleMatrix2D(new double[][] { b }))).viewColumn(0).toArray(); return x; */ } }
@Override public Solver getSolver(RealMatrix M) { if (M == null) { return null; } RRQRDecomposition decomposition = new RRQRDecomposition(M, SINGULARITY_THRESHOLD); DecompositionSolver solver = decomposition.getSolver(); if (solver.isNonSingular()) { return new CommonsMathSolver(solver); } // Otherwise try to report apparent rank int apparentRank = decomposition.getRank(0.01); // Better value? log.warn("{} x {} matrix is near-singular (threshold {}). Add more data or decrease the value of model.features, " + "to <= about {}", M.getRowDimension(), M.getColumnDimension(), SINGULARITY_THRESHOLD, apparentRank); throw new SingularMatrixSolverException(apparentRank, "Apparent rank: " + apparentRank); }
@Override public SlopeCoefficients estimateCoefficients(final DerivationEquation eq) throws EstimationException { final double[][] sourceTriangleMatrix = eq.getCovarianceLowerTriangularMatrix(); // Copy matrix and enhance it to a full matrix as expected by CholeskyDecomposition // FIXME: Avoid copy job to speed-up the solving process e.g. by extending the CholeskyDecomposition constructor final int length = sourceTriangleMatrix.length; final double[][] matrix = new double[length][]; for (int i = 0; i < length; i++) { matrix[i] = new double[length]; final double[] s = sourceTriangleMatrix[i]; final double[] t = matrix[i]; for (int j = 0; j <= i; j++) { t[j] = s[j]; } for (int j = i + 1; j < length; j++) { t[j] = sourceTriangleMatrix[j][i]; } } final RealMatrix coefficients = new Array2DRowRealMatrix(matrix, false); try { final DecompositionSolver solver = new CholeskyDecomposition(coefficients).getSolver(); final RealVector constants = new ArrayRealVector(eq.getConstraints(), true); final RealVector solution = solver.solve(constants); return new DefaultSlopeCoefficients(solution.toArray()); } catch (final NonPositiveDefiniteMatrixException e) { throw new EstimationException("Matrix inversion error due to data is linearly dependent", e); } }
@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); }
/** * Constructor for more than 4 point pairs, finds a least-squares solution * for the homography parameters. * NOTE: this is UNFINISHED code! * @param P sequence of points (source) * @param Q sequence of points (target) * @param dummy unused (only to avoid duplicate signature) */ public ProjectiveMapping(Point2D[] P, Point2D[] Q, boolean dummy) { final int n = P.length; double[] ba = new double[2 * n]; double[][] Ma = new double[2 * n][]; for (int i = 0; i < n; i++) { double x = P[i].getX(); double y = P[i].getY(); double u = Q[i].getX(); double v = Q[i].getY(); ba[2 * i + 0] = u; ba[2 * i + 1] = v; Ma[2 * i + 0] = new double[] { x, y, 1, 0, 0, 0, -u * x, -u * y }; Ma[2 * i + 1] = new double[] { 0, 0, 0, x, y, 1, -v * x, -v * y }; } RealMatrix M = MatrixUtils.createRealMatrix(Ma); RealVector b = MatrixUtils.createRealVector(ba); DecompositionSolver solver = new SingularValueDecomposition(M).getSolver(); RealVector h = solver.solve(b); a00 = h.getEntry(0); a01 = h.getEntry(1); a02 = h.getEntry(2); a10 = h.getEntry(3); a11 = h.getEntry(4); a12 = h.getEntry(5); a20 = h.getEntry(6); a21 = h.getEntry(7); a22 = 1; }
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; }
/** * @see AbstractSolver#solve(AbstractMatrix, AbstractVector) */ @Override public AbstractVector solve(AbstractMatrix m, AbstractVector b) { if (m instanceof ApacheMatrix && b instanceof ApacheVector) { DecompositionSolver solver = new LUDecomposition(((ApacheMatrix) m).getMatrix()).getSolver(); RealVector bApache = ((ApacheVector) b).getVector(); RealVector solution = solver.solve(bApache); return new ApacheVector(solution); } throw new UnsupportedOperationException(); }
/** * pdf(x, x_hat) = exp(-0.5 * (x-x_hat) * inv(Σ) * (x-x_hat)T) / ( 2π^0.5d * det(Σ)^0.5) * * @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; }
@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(); } else { SingularValueDecomposition SVD = new SingularValueDecomposition(m); inv = SVD.getSolver().getInverse(); } return inv; }
/** * L = A x R * * @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; }
/** {@inheritDoc} */ public RealMatrix getCovariances(double threshold) { // Set up the Jacobian. final RealMatrix j = this.getJacobian(); // Compute transpose(J)J. final RealMatrix jTj = j.transpose().multiply(j); // Compute the covariances matrix. final DecompositionSolver solver = new QRDecomposition(jTj, threshold).getSolver(); return solver.getInverse(); }
/** * Correct the current state estimate with an actual measurement. * * @param z the measurement vector * @throws DimensionMismatchException if the dimension of the * measurement vector does not fit * @throws org.apache.commons.math3.linear.SingularMatrixException * if the covariance matrix could not be inverted */ public void correct(final RealVector z) { // sanity checks MathUtils.checkNotNull(z); if (z.getDimension() != measurementMatrix.getRowDimension()) { throw new DimensionMismatchException(z.getDimension(), measurementMatrix.getRowDimension()); } // S = H * P(k) - * H' + R RealMatrix s = measurementMatrix.multiply(errorCovariance) .multiply(measurementMatrixT) .add(measurementModel.getMeasurementNoise()); // invert S // as the error covariance matrix is a symmetric positive // semi-definite matrix, we can use the cholesky decomposition DecompositionSolver solver = new CholeskyDecomposition(s).getSolver(); RealMatrix invertedS = solver.getInverse(); // Inn = z(k) - H * xHat(k)- RealVector innovation = z.subtract(measurementMatrix.operate(stateEstimation)); // calculate gain matrix // K(k) = P(k)- * H' * (H * P(k)- * H' + R)^-1 // K(k) = P(k)- * H' * S^-1 RealMatrix kalmanGain = errorCovariance.multiply(measurementMatrixT).multiply(invertedS); // update estimate with measurement z(k) // xHat(k) = xHat(k)- + K * Inn stateEstimation = stateEstimation.add(kalmanGain.operate(innovation)); // update covariance of prediction error // P(k) = (I - K * H) * P(k)- RealMatrix identity = MatrixUtils.createRealIdentityMatrix(kalmanGain.getRowDimension()); errorCovariance = identity.subtract(kalmanGain.multiply(measurementMatrix)).multiply(errorCovariance); }
/** Guess a start line using the last four results. * @param n number of cached results available * @param target target ground point * @return guessed start line */ private double guessStartLine(final int n, final Vector3D target) { try { // assume a linear model of the form: l = ax + by + cz + d final RealMatrix m = new Array2DRowRealMatrix(n, 4); final RealVector v = new ArrayRealVector(n); for (int i = 0; i < n; ++i) { m.setEntry(i, 0, cachedResults[i].getTarget().getX()); m.setEntry(i, 1, cachedResults[i].getTarget().getY()); m.setEntry(i, 2, cachedResults[i].getTarget().getZ()); m.setEntry(i, 3, 1.0); v.setEntry(i, cachedResults[i].getLine()); } final DecompositionSolver solver = new QRDecomposition(m, Precision.SAFE_MIN).getSolver(); final RealVector coefficients = solver.solve(v); // apply the linear model return target.getX() * coefficients.getEntry(0) + target.getY() * coefficients.getEntry(1) + target.getZ() * coefficients.getEntry(2) + coefficients.getEntry(3); } catch (SingularMatrixException sme) { // the points are not independent return Double.NaN; } }
/** * Solve a linear matrix equation, or system of linear scalar equations. * * @param a Coefficient matrix. * @param b Ordinate or “dependent variable” values. * @return Solution to the system a x = b. Returned shape is identical to b. */ public static Array solve(Array a, Array b) { Array r = Array.factory(DataType.DOUBLE, b.getShape()); double[][] aa = (double[][]) ArrayUtil.copyToNDJavaArray(a); RealMatrix coefficients = new Array2DRowRealMatrix(aa, false); DecompositionSolver solver = new LUDecomposition(coefficients).getSolver(); double[] bb = (double[]) ArrayUtil.copyToNDJavaArray(b); RealVector constants = new ArrayRealVector(bb, false); RealVector solution = solver.solve(constants); for (int i = 0; i < r.getSize(); i++) { r.setDouble(i, solution.getEntry(i)); } return r; }
/** * Function to compute matrix inverse via matrix decomposition. * * @param in commons-math3 Array2DRowRealMatrix * @return matrix block * @throws DMLRuntimeException if DMLRuntimeException occurs */ private static MatrixBlock computeMatrixInverse(Array2DRowRealMatrix in) throws DMLRuntimeException { 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()); }
/** * Correct the current state estimate with an actual measurement. * * @param z * the measurement vector * @throws NullArgumentException * if the measurement vector is {@code null} * @throws DimensionMismatchException * if the dimension of the measurement vector does not fit * @throws SingularMatrixException * if the covariance matrix could not be inverted */ public void correct(final RealVector z) throws NullArgumentException, DimensionMismatchException, SingularMatrixException { // sanity checks MathUtils.checkNotNull(z); if (z.getDimension() != measurementMatrix.getRowDimension()) { throw new DimensionMismatchException(z.getDimension(), measurementMatrix.getRowDimension()); } // S = H * P(k) - * H' + R RealMatrix s = measurementMatrix.multiply(errorCovariance) .multiply(measurementMatrixT) .add(measurementModel.getMeasurementNoise()); // invert S // as the error covariance matrix is a symmetric positive // semi-definite matrix, we can use the cholesky decomposition DecompositionSolver solver = new CholeskyDecomposition(s).getSolver(); RealMatrix invertedS = solver.getInverse(); // Inn = z(k) - H * xHat(k)- RealVector innovation = z.subtract(measurementMatrix.operate(stateEstimation)); // calculate gain matrix // K(k) = P(k)- * H' * (H * P(k)- * H' + R)^-1 // K(k) = P(k)- * H' * S^-1 RealMatrix kalmanGain = errorCovariance.multiply(measurementMatrixT).multiply(invertedS); // update estimate with measurement z(k) // xHat(k) = xHat(k)- + K * Inn stateEstimation = stateEstimation.add(kalmanGain.operate(innovation)); // update covariance of prediction error // P(k) = (I - K * H) * P(k)- RealMatrix identity = MatrixUtils.createRealIdentityMatrix(kalmanGain.getRowDimension()); errorCovariance = identity.subtract(kalmanGain.multiply(measurementMatrix)).multiply(errorCovariance); }
/** * Calculates the parameters of a bivariate quadratic equation. * * @param elevationValues the window of points to use. * @return the parameters of the bivariate quadratic equation as [a, b, c, d, e, f] */ private static double[] calculateParameters( final double[][] elevationValues ) { int rows = elevationValues.length; int cols = elevationValues[0].length; int pointsNum = rows * cols; final double[][] xyMatrix = new double[pointsNum][6]; final double[] valueArray = new double[pointsNum]; // TODO check on resolution int index = 0; for( int y = 0; y < rows; y++ ) { for( int x = 0; x < cols; x++ ) { xyMatrix[index][0] = x * x; // x^2 xyMatrix[index][1] = y * y; // y^2 xyMatrix[index][2] = x * y; // xy xyMatrix[index][3] = x; // x xyMatrix[index][4] = y; // y xyMatrix[index][5] = 1; valueArray[index] = elevationValues[y][x]; index++; } } RealMatrix A = MatrixUtils.createRealMatrix(xyMatrix); RealVector z = MatrixUtils.createRealVector(valueArray); DecompositionSolver solver = new RRQRDecomposition(A).getSolver(); RealVector solution = solver.solve(z); // start values for a, b, c, d, e, f, all set to 0.0 final double[] parameters = solution.toArray(); return parameters; }
/** * Find a solution of the linear equation A x = b where * <ul> * <li>A is an n x m - matrix given as double[n][m]</li> * <li>b is an m - vector given as double[m],</li> * <li>x is an n - vector given as double[n],</li> * </ul> * * @param A The matrix (left hand side of the linear equation). * @param b The vector (right hand of the linear equation). * @return A solution x to A x = b. */ public static double[] solveLinearEquation(double[][] A, double[] b) { if(isSolverUseApacheCommonsMath) { Array2DRowRealMatrix matrix = new Array2DRowRealMatrix(A); DecompositionSolver solver; if(matrix.getColumnDimension() == matrix.getRowDimension()) { solver = new LUDecomposition(matrix).getSolver(); } else { solver = new QRDecomposition(new Array2DRowRealMatrix(A)).getSolver(); } // Using SVD - very slow // solver = new SingularValueDecomposition(new Array2DRowRealMatrix(A)).getSolver(); return solver.solve(new Array2DRowRealMatrix(b)).getColumn(0); } else { return org.jblas.Solve.solve(new org.jblas.DoubleMatrix(A), new org.jblas.DoubleMatrix(b)).data; // For use of colt: // cern.colt.matrix.linalg.Algebra linearAlgebra = new cern.colt.matrix.linalg.Algebra(); // return linearAlgebra.solve(new DenseDoubleMatrix2D(A), linearAlgebra.transpose(new DenseDoubleMatrix2D(new double[][] { b }))).viewColumn(0).toArray(); // For use of parallel colt: // return new cern.colt.matrix.tdouble.algo.decomposition.DenseDoubleLUDecomposition(new cern.colt.matrix.tdouble.impl.DenseDoubleMatrix2D(A)).solve(new cern.colt.matrix.tdouble.impl.DenseDoubleMatrix1D(b)).toArray(); } }
private double[] interpolateStep(int r, int c, ResponseLayer t, ResponseLayer m, ResponseLayer b){ double[] values = new double[3]; RealMatrix partialDerivs = getPartialDerivativeMatrix(r,c,t,m,b); RealMatrix hessian3D = getHessian3DMatrix(r,c,t,m,b); DecompositionSolver solver = new LUDecomposition(hessian3D).getSolver(); RealMatrix X = solver.getInverse().multiply(partialDerivs); // System.out.println("X = " + X.getColumnDimension() + " col x " + X.getRowDimension() + " rows"); // for ( int i = 0; i < X.getRowDimension(); i++ ){ // for ( int j = 0; j < X.getColumnDimension(); j++ ){ // System.out.print(X.getEntry(i,j) + (j != X.getColumnDimension()-1 ? " - " : "")); // } // System.out.println(); // } // System.out.println(); // //values of them are used //xi values[0] = -X.getEntry(2,0); //xr values[1] = -X.getEntry(1,0); //xc values[2] = -X.getEntry(0,0); return values; }
