List of usage examples for org.apache.commons.math3.linear DecompositionSolver getInverse
RealMatrix getInverse();
From source file:com.ibm.bi.dml.runtime.matrix.data.LibCommonsMath.java
/** * Function to compute matrix inverse via matrix decomposition. * /* ww w .j av a 2 s .co m*/ * @param in * @return * @throws DMLRuntimeException */ 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(); MatrixBlock inverse = DataConverter.convertToMatrixBlock(inverseMatrix.getData()); return inverse; }
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 ww . j a v a2 s . c o m*/ * @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
@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 {//from w ww. ja v a 2 s . co m SingularValueDecomposition SVD = new SingularValueDecomposition(m); inv = SVD.getSolver().getInverse(); } return inv; }
From source file:ellipsoidFit.FitPoints.java
/** * Find the center of the ellipsoid./*from w ww . j a v a2 s . c o m*/ * * @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); }
From source file:ellipsoidFit.FitPoints.java
/** * Solve the polynomial expression Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + * 2Gx + 2Hy + 2Iz from the provided points. * /* w ww. ja va 2s. co m*/ * @param points * the points that will be fit to the polynomial expression. * @return the solution vector to the polynomial expression. */ private RealVector solveSystem(ArrayList<ThreeSpacePoint> points) { // determine the number of points int numPoints = points.size(); // the design matrix // size: numPoints x 9 RealMatrix d = new Array2DRowRealMatrix(numPoints, 9); // Fit the ellipsoid in the form of // Ax^2 + By^2 + Cz^2 + 2Dxy + 2Exz + 2Fyz + 2Gx + 2Hy + 2Iz for (int i = 0; i < d.getRowDimension(); i++) { double xx = Math.pow(points.get(i).x, 2); double yy = Math.pow(points.get(i).y, 2); double zz = Math.pow(points.get(i).z, 2); double xy = 2 * (points.get(i).x * points.get(i).y); double xz = 2 * (points.get(i).x * points.get(i).z); double yz = 2 * (points.get(i).y * points.get(i).z); double x = 2 * points.get(i).x; double y = 2 * points.get(i).y; double z = 2 * points.get(i).z; d.setEntry(i, 0, xx); d.setEntry(i, 1, yy); d.setEntry(i, 2, zz); d.setEntry(i, 3, xy); d.setEntry(i, 4, xz); d.setEntry(i, 5, yz); d.setEntry(i, 6, x); d.setEntry(i, 7, y); d.setEntry(i, 8, z); } // solve the normal system of equations // v = (( d' * d )^-1) * ( d' * ones.mapAddToSelf(1)); // Multiply: d' * d RealMatrix dtd = d.transpose().multiply(d); // Create a vector of ones. RealVector ones = new ArrayRealVector(numPoints); ones.mapAddToSelf(1); // Multiply: d' * ones.mapAddToSelf(1) RealVector dtOnes = d.transpose().operate(ones); // Find ( d' * d )^-1 DecompositionSolver solver = new SingularValueDecomposition(dtd).getSolver(); RealMatrix dtdi = solver.getInverse(); // v = (( d' * d )^-1) * ( d' * ones.mapAddToSelf(1)); RealVector v = dtdi.operate(dtOnes); return v; }
From source file:jopensurf.FastHessian.java
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 ? " - " : "")); // }// w ww . ja va 2 s.c om // 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; }
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;//from w ww .j ava 2 s . co m 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:citation_prediction.CitationCore.java
/** * This function implements the algorithm designed by Josiah Neuberger and William Etcho used to solve for a WSB solution. * The general math for the Newton-Raphson method was provided by Dr. Allen Parks and can be found in the function 'getPartialsData'. * <br><br>/*ww w . ja v a 2s. c o m*/ * This algorithm was created under the direction of Supervising University of Mary Washington faculty, Dr. Melody Denhere along with * consultation from Dr. Jeff Solka and Computer Scientists Kristin Ash with the Dahlgren Naval Surface Warfare Center. * <br><br> * If your interested in more information about the research behind this algorithm please visit:<br> * http://josiahneuberger.github.io/citation_prediction/ * <br><br> * * @param data The citation data in days. * @param mu The initial mu guess to use in the Newton-Raphson method. * @param sigma The initial sigma guess to use in the Newton-Raphson method. * @param m The constant value, which is determined by the average number of references in each new paper for a journal. * @param t The last time value in the paper's citation history. * @param n The total number of citations for this paper. * @param l A list structure to store values for each iteration (should be null to start). * @param iteration The iteration (should be zero to start) * @param max_iteration The maximum number of iterations to try before stopping. * @param tolerance The tolerance level, which determines that a solution has been converged on. * @return A list containing the WSB solution of (lambda, mu, sigma, iterations). */ private LinkedHashMap<String, Double> newtonRaphson(double[][] data, double mu, double sigma, double m, double t, double n, LinkedHashMap<String, Double> l, int iteration, int max_iteration, double tolerance) { double lambda; LinkedHashMap<String, Double> r = new LinkedHashMap<String, Double>(); if (iteration > max_iteration) { System.out.println("Does not converge."); r.put("lambda", null); r.put("mu", null); r.put("sigma", null); r.put("iterations", null); return r; } else if (tolerance < 0.00000001) { System.out.println("Stopped due to tolerance."); r.put("lambda", getLambda(data, mu, sigma, m, t, n)); r.put("mu", mu); r.put("sigma", sigma); r.put("iterations", (double) iteration); return r; } else { l = getPartialsData(getIterationData(data, mu, sigma, m)); double[] array_xn = { mu, sigma }; double[] array_yn = { l.get("fn"), l.get("gn") }; RealMatrix xn = MatrixUtils.createColumnRealMatrix(array_xn); RealMatrix yn = MatrixUtils.createColumnRealMatrix(array_yn); //http://commons.apache.org/proper/commons-math/userguide/linear.html double[][] array_jacobian = { { l.get("df_dmu"), l.get("df_dsigma") }, { l.get("dg_dmu"), l.get("dg_dsigma") } }; RealMatrix jacobian = MatrixUtils.createRealMatrix(array_jacobian); LUDecomposition lud = new LUDecomposition(jacobian); DecompositionSolver decS = lud.getSolver(); l.put("iteration", (double) (iteration + 1)); //DEBUG: printList(l); if (!decS.isNonSingular()) { l.put("iteration", (double) (max_iteration + 1)); System.err.println("ERROR: Jacobian matrix was singular."); } else { RealMatrix solution = xn.subtract(decS.getInverse().multiply(yn)); RealMatrix xndiff = solution.subtract(xn); tolerance = Math.sqrt(Math.pow(xndiff.getEntry(0, 0), 2) + Math.pow(xndiff.getEntry(1, 0), 2)); //update values l.put("mu", solution.getEntry(0, 0)); l.put("sigma", solution.getEntry(1, 0)); //DEBUG: System.out.printf("\"%-20s=%25f\"\n", "NEW MU", l.get("mu")); //DEBUG: System.out.printf("\"%-20s=%25f\"\n", "NEW SIGMA", l.get("sigma")); } //DEBUG: System.out.println("****************************************"); return newtonRaphson(data, l.get("mu"), l.get("sigma"), m, l, (l.get("iteration").intValue()), tolerance); } }
From source file:org.apache.sysml.runtime.matrix.data.LibCommonsMath.java
/** * Function to compute matrix inverse via matrix decomposition. * /* w w w. ja va2s .c om*/ * @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()); }
From source file:org.orekit.orbits.Orbit.java
/** Create an inverse Jacobian. * @param type type of the position angle to use * @return inverse Jacobian//from w w w. j a va2 s .com */ private double[][] createInverseJacobian(final PositionAngle type) { // get the direct Jacobian final double[][] directJacobian = new double[6][6]; getJacobianWrtCartesian(type, directJacobian); // invert the direct Jacobian final RealMatrix matrix = MatrixUtils.createRealMatrix(directJacobian); final DecompositionSolver solver = new QRDecomposition(matrix).getSolver(); return solver.getInverse().getData(); }