List of usage examples for org.apache.commons.math3.linear RealMatrix multiply
RealMatrix multiply(RealMatrix m) throws DimensionMismatchException;
From source file:com.analog.lyric.dimple.solvers.sumproduct.customFactors.MutlivariateGaussianMatrixProduct.java
public MutlivariateGaussianMatrixProduct(double[][] A) { int i; //,m; M = A.length;/*from w w w .j av a 2s . c o m*/ N = A[0].length; /* Here we precompute and store matrices for future message computations. * First, compute an SVD of the matrix A using EigenDecompositions of A*A^T and A^T*A * This way, we get nullspaces for free along with regularized inverse. */ RealMatrix Amat = wrapRealMatrix(A); SingularValueDecomposition svd = new SingularValueDecomposition(Amat); RealMatrix tmp = svd.getVT(); tmp = svd.getS().multiply(tmp); tmp = svd.getU().multiply(tmp); A_clean = matrixGetDataRef(tmp); RealMatrix ST = svd.getS().transpose(); int numS = Math.min(ST.getColumnDimension(), ST.getRowDimension()); for (i = 0; i < numS; i++) { double d = ST.getEntry(i, i); if (d < eps) d = eps; else if (d > 1 / eps) d = 1 / eps; ST.setEntry(i, i, 1.0 / d); } A_pinv = matrixGetDataRef(svd.getV().multiply(ST.multiply(svd.getUT()))); }
From source file:com.opengamma.strata.math.impl.matrix.CommonsMatrixAlgebra.java
/** * {@inheritDoc}//w w w .j a v a2s . c om * The following combinations of input matrices m1 and m2 are allowed: * <ul> * <li> m1 = 2-D matrix, m2 = 2-D matrix, returns $\mathbf{C} = \mathbf{AB}$ * <li> m1 = 2-D matrix, m2 = 1-D matrix, returns $\mathbf{C} = \mathbf{A}b$ * </ul> */ @Override public Matrix multiply(Matrix m1, Matrix m2) { ArgChecker.notNull(m1, "m1"); ArgChecker.notNull(m2, "m2"); ArgChecker.isTrue(!(m1 instanceof DoubleArray), "Cannot have 1D matrix as first argument"); if (m1 instanceof DoubleMatrix) { RealMatrix t1 = CommonsMathWrapper.wrap((DoubleMatrix) m1); RealMatrix t2; if (m2 instanceof DoubleArray) { t2 = CommonsMathWrapper.wrapAsMatrix((DoubleArray) m2); } else if (m2 instanceof DoubleMatrix) { t2 = CommonsMathWrapper.wrap((DoubleMatrix) m2); } else { throw new IllegalArgumentException("Can only have 1D or 2D matrix as second argument"); } return CommonsMathWrapper.unwrap(t1.multiply(t2)); } throw new IllegalArgumentException("Can only multiply 2D and 1D matrices"); }
From source file:hulo.localization.models.obs.GaussianProcess.java
public double looMSE() { double[][] Y = getY(); double[][] dY = getdY(); int ns = X.length; int ny = Y[0].length; RealMatrix invKy = MatrixUtils.createRealMatrix(this.invKy); RealMatrix K = MatrixUtils.createRealMatrix(this.K); RealMatrix Hmat = invKy.multiply(K); double[][] H = Hmat.getData(); double sum = 0; double count = 0; for (int j = 0; j < ny; j++) { for (int i = 0; i < ns; i++) { double preddY = 0; for (int k = 0; k < ns; k++) { preddY += H[i][k] * dY[k][j]; }//w ww.j a v a2 s . c om double diff = (dY[i][j] - preddY) / (1.0 - H[i][i]); sum += (diff * diff); count += 1; } } sum /= count; return sum; }
From source file:com.github.thorbenlindhauer.factor.CanonicalGaussianFactor.java
@Override public GaussianFactor marginal(Scope marginalizationScope) { // the following assumes that the precision matrix (and mean vector) can be restructured as follows: // ( SUBMATRIX_XX SUBMATRIX_XY ) // ( SUBMATRIX_YX SUBMATRIX_YY ) // where X indicates the entries for the variables that are kept (i.e. in the scope argument) and // Y the variables that are marginalized out if (marginalizationScope.contains(scope)) { return this; }/*w w w . ja v a 2 s. c om*/ Scope newScope = scope.intersect(marginalizationScope); Scope scopeToMarginalize = scope.reduceBy(newScope); int[] xMapping = newScope.createContinuousVariableMapping(scope); RealMatrix xxMatrix = precisionMatrix.getSubMatrix(xMapping, xMapping); int[] yMapping = scopeToMarginalize.createContinuousVariableMapping(scope); RealMatrix yyMatrix = precisionMatrix.getSubMatrix(yMapping, yMapping); RealMatrix xyMatrix = precisionMatrix.getSubMatrix(xMapping, yMapping); RealMatrix yxMatrix = precisionMatrix.getSubMatrix(yMapping, xMapping); MathUtil yyUtil = new MathUtil(yyMatrix); RealMatrix yyInverse = yyUtil.invert(); RealMatrix newPrecisionMatrix = xxMatrix.subtract(xyMatrix.multiply(yyInverse.multiply(yxMatrix))); RealVector xVector = getSubVector(scaledMeanVector, xMapping); RealVector yVector = getSubVector(scaledMeanVector, yMapping); RealVector newScaledMeanVector = xVector.subtract(xyMatrix.operate(yyInverse.operate(yVector))); MathUtil inverseUtil = new MathUtil(yyInverse.scalarMultiply(2.0d * Math.PI)); double newNormalizationConstant = normalizationConstant + 0.5d * (Math.log(inverseUtil.determinant()) + yVector.dotProduct(yyInverse.operate(yVector))); return new CanonicalGaussianFactor(newScope, newPrecisionMatrix, newScaledMeanVector, newNormalizationConstant); }
From source file:io.yields.math.concepts.operator.Smoothness.java
private RealMatrix computeDistance(List<Tuple> normalizedData) { RealMatrix a = new Array2DRowRealMatrix(normalizedData.size(), order + 1); RealMatrix b = new Array2DRowRealMatrix(normalizedData.size(), 1); int row = 0;//from ww w. j a v a 2s . c om for (Tuple tuple : normalizedData) { final int currentRow = row; IntStream.rangeClosed(0, order) .forEach(power -> a.setEntry(currentRow, power, Math.pow(tuple.getX(), power))); b.setEntry(currentRow, 0, tuple.getY()); row++; } DecompositionSolver solver = new QRDecomposition(a, EPS).getSolver(); if (solver.isNonSingular() && !isConstant(b)) { RealMatrix solution = solver.solve(b); return a.multiply(solution).subtract(b); } else { return new Array2DRowRealMatrix(normalizedData.size(), 1); } }
From source file:com.joptimizer.functions.SOCPLogarithmicBarrier.java
public double[][] hessian(double[] X) { RealVector x = new ArrayRealVector(X); RealMatrix ret = new Array2DRowRealMatrix(dim, dim); for (int i = 0; i < socpConstraintParametersList.size(); i++) { SOCPConstraintParameters param = socpConstraintParametersList.get(i); double t = this.buildT(param, x); RealVector u = this.buildU(param, x); double t2uu = t * t - u.dotProduct(u); RealVector t2u = u.mapMultiply(-2 * t); RealMatrix Jacob = this.buildJ(param, x); int k = u.getDimension(); RealMatrix H = new Array2DRowRealMatrix(k + 1, k + 1); RealMatrix ID = MatrixUtils.createRealIdentityMatrix(k); H.setSubMatrix(ID.scalarMultiply(t2uu).add(u.outerProduct(u).scalarMultiply(2)).getData(), 0, 0); H.setSubMatrix(new double[][] { t2u.toArray() }, k, 0); for (int j = 0; j < k; j++) { H.setEntry(j, k, t2u.getEntry(j)); }//w w w . j av a 2 s .c om H.setEntry(k, k, t * t + u.dotProduct(u)); RealMatrix ret_i = Jacob.multiply(H).multiply(Jacob.transpose()).scalarMultiply(2. / Math.pow(t2uu, 2)); ret = ret.add(ret_i); } return ret.getData(); }
From source file:io.github.malapert.jwcs.JWcs.java
/** * Transforms the sky position given by (longitude, latitude) in a pixel * position./*from ww w.jav a 2 s . c o m*/ * * @param longitude longitude of the sky position * @param latitude latitude of the sky position * @return the sky position in the pixel grid. * @throws io.github.malapert.jwcs.proj.exception.ProjectionException when there is a projection error */ @Override public double[] wcs2pix(double longitude, double latitude) throws ProjectionException { checkLongitudeLatitude(longitude, latitude); double[] coordVal = this.getProj().wcs2projectionPlane(longitude, latitude); double[][] coord = { { coordVal[0], coordVal[1] } }; RealMatrix coordM = MatrixUtils.createRealMatrix(coord); RealMatrix matrix = coordM.multiply(getCdInverse()); return new double[] { matrix.getEntry(0, 0) + crpix(1), matrix.getEntry(0, 1) + crpix(2) }; }
From source file:com.datumbox.framework.core.machinelearning.featureselection.PCA.java
/** {@inheritDoc} */ @Override/*from www . j a va2s . co m*/ protected void _fit(Dataframe trainingData) { ModelParameters modelParameters = knowledgeBase.getModelParameters(); int n = trainingData.size(); int d = trainingData.xColumnSize(); //convert data into matrix Map<Object, Integer> featureIds = modelParameters.getFeatureIds(); DataframeMatrix matrixDataset = DataframeMatrix.newInstance(trainingData, false, null, featureIds); RealMatrix X = matrixDataset.getX(); //calculate means and subtract them from data RealVector meanValues = new OpenMapRealVector(d); for (Integer columnId : featureIds.values()) { double mean = 0.0; for (int row = 0; row < n; row++) { mean += X.getEntry(row, columnId); } mean /= n; for (int row = 0; row < n; row++) { X.addToEntry(row, columnId, -mean); } meanValues.setEntry(columnId, mean); } modelParameters.setMean(meanValues); //dxd matrix RealMatrix covarianceDD = (X.transpose().multiply(X)).scalarMultiply(1.0 / (n - 1.0)); EigenDecomposition decomposition = new EigenDecomposition(covarianceDD); RealVector eigenValues = new ArrayRealVector(decomposition.getRealEigenvalues(), false); RealMatrix components = decomposition.getV(); //Whiten Components W = U*L^0.5; To whiten them we multiply with L^0.5. if (knowledgeBase.getTrainingParameters().isWhitened()) { RealMatrix sqrtEigenValues = new DiagonalMatrix(d); for (int i = 0; i < d; i++) { sqrtEigenValues.setEntry(i, i, FastMath.sqrt(eigenValues.getEntry(i))); } components = components.multiply(sqrtEigenValues); } //the eigenvalues and their components are sorted by descending order no need to resort them Integer maxDimensions = knowledgeBase.getTrainingParameters().getMaxDimensions(); Double variancePercentageThreshold = knowledgeBase.getTrainingParameters().getVariancePercentageThreshold(); if (variancePercentageThreshold != null && variancePercentageThreshold <= 1) { double totalVariance = 0.0; for (int i = 0; i < d; i++) { totalVariance += eigenValues.getEntry(i); } double sum = 0.0; int varCounter = 0; for (int i = 0; i < d; i++) { sum += eigenValues.getEntry(i) / totalVariance; varCounter++; if (sum >= variancePercentageThreshold) { break; } } if (maxDimensions == null || maxDimensions > varCounter) { maxDimensions = varCounter; } } if (maxDimensions != null && maxDimensions < d) { //keep only the maximum selected eigenvalues eigenValues = eigenValues.getSubVector(0, maxDimensions); //keep only the maximum selected eigenvectors components = components.getSubMatrix(0, components.getRowDimension() - 1, 0, maxDimensions - 1); } modelParameters.setEigenValues(eigenValues); modelParameters.setComponents(components); }
From source file:edu.cmu.tetrad.util.TetradMatrix1.java
public TetradMatrix1 sqrt() { SingularValueDecomposition svd = new SingularValueDecomposition(getRealMatrix()); RealMatrix U = svd.getU(); RealMatrix V = svd.getV();//from w ww . j a v a 2 s . c o m double[] s = svd.getSingularValues(); for (int i = 0; i < s.length; i++) s[i] = 1.0 / s[i]; RealMatrix S = new BlockRealMatrix(s.length, s.length); for (int i = 0; i < s.length; i++) S.setEntry(i, i, s[i]); RealMatrix sqrt = U.multiply(S).multiply(V); return new TetradMatrix1(sqrt); }
From source file:hulo.localization.models.obs.GaussianProcess.java
public double looPredLogLikelihood() { double[][] Y = getY(); double[][] dY = getdY(); int ns = X.length; int ny = Y[0].length; RealMatrix Ky = MatrixUtils.createRealMatrix(this.Ky); RealMatrix invKy = new LUDecomposition(Ky).getSolver().getInverse(); RealMatrix dYmat = MatrixUtils.createRealMatrix(dY); double[] LOOPredLL = new double[ny]; for (int j = 0; j < ny; j++) { RealMatrix dy = MatrixUtils.createColumnRealMatrix(dYmat.getColumn(j)); RealMatrix invKdy = invKy.multiply(dy); double sum = 0; for (int i = 0; i < ns; i++) { double mu_i = dYmat.getEntry(i, j) - invKdy.getEntry(i, 0) / invKy.getEntry(i, i); double sigma_i2 = 1 / invKy.getEntry(i, i); double logLL = StatUtils.logProbaNormal(dYmat.getEntry(i, j), mu_i, Math.sqrt(sigma_i2)); sum += logLL;// w ww . ja v a 2 s. com } LOOPredLL[j] = sum; } double sumLOOPredLL = 0; for (int j = 0; j < ny; j++) { sumLOOPredLL += LOOPredLL[j]; } return sumLOOPredLL; }