List of usage examples for org.apache.commons.math3.linear RealMatrix scalarMultiply
RealMatrix scalarMultiply(double d);
From source file:edu.cudenver.bios.matrix.MatrixUtils.java
/** * Calculate the horizontal direct product of two matrices * * @param matrix1 first matrix/*ww w . j a va2 s .c o m*/ * @param matrix2 second matrix * @return horizontal direct product of matrix 1 and matrix 2 */ public static RealMatrix getHorizontalDirectProduct(RealMatrix matrix1, RealMatrix matrix2) throws IllegalArgumentException { if (matrix1 == null || matrix2 == null) throw new IllegalArgumentException("null input matrix"); if (matrix1.getRowDimension() != matrix2.getRowDimension()) throw new IllegalArgumentException("input matrices must have equal row dimension"); int mRows = matrix1.getRowDimension(); int m1Cols = matrix1.getColumnDimension(); int m2Cols = matrix2.getColumnDimension(); double[][] productData = new double[mRows][m1Cols * m2Cols]; RealMatrix productMatrix = new Array2DRowRealMatrix(productData); for (int col = 0; col < m1Cols; col++) { for (int row = 0; row < mRows; row++) { RealMatrix m2Row = matrix2.getRowMatrix(row); productMatrix.setSubMatrix((m2Row.scalarMultiply(matrix1.getEntry(row, col))).getData(), row, col * m2Cols); } } return productMatrix; }
From source file:com.github.thorbenlindhauer.factor.CanonicalGaussianFactor.java
/** * Transforms a conditional linear gaussian (i.e. a Gaussian of the form N(x; a + B^T * Y, C) into canonical form. * * @param meanVector a/*from w w w . j a va 2 s .co m*/ * @param weightMatrix B * @param covarianceMatrix C */ public static CanonicalGaussianFactor fromConditionalForm(Scope scope, Scope conditioningScope, RealVector meanVector, RealMatrix covarianceMatrix, RealMatrix weightMatrix) { // TODO: perform cardinality checks etc. // the following assumes that the resulting precision matrix (and mean vector) can be restructured as follows: // ( SUBMATRIX_XX SUBMATRIX_XY ) // ( SUBMATRIX_YX SUBMATRIX_YY ) // where X indicates variables that are part of the prediction scope and Y are variables being part of the conditioning scope // assuming // meanVector: a; |x| vector // covarianceMatrix: C; |x| cross |x| matrix // weightMatrix: B; |x| cross |y| matrix // XX = C^-1 // XY = -C^-1 * B // YX = -B^T * C^-1 // YY = B^T * C^-1 * B^T MathUtil mathUtil = new MathUtil(covarianceMatrix); // C^(-1) RealMatrix xxMatrix = null; xxMatrix = mathUtil.invert(); // if (!mathUtil.isZeroMatrix()) { // xxMatrix = mathUtil.invert(); // } else { // // // this is a special case for convolution in which the "summing" variable has no variance itself // // although a 0 variance is not valid in general // xxMatrix = covarianceMatrix; // } // B^T * C^(-1) RealMatrix weightedXXMatrix = weightMatrix.transpose().multiply(xxMatrix); // -B^T * C^(-1) RealMatrix yxMatrix = weightedXXMatrix.scalarMultiply(-1.0d); // -C^(-1)^T * B RealMatrix xyMatrix = xxMatrix.transpose().multiply(weightMatrix).scalarMultiply(-1.0d); // B^T * C^(-1) * B RealMatrix yyMatrix = weightedXXMatrix.multiply(weightMatrix); // K RealMatrix conditionedPrecisionMatrix = new Array2DRowRealMatrix(scope.size(), scope.size()); // Matrix to generate h RealMatrix conditionedMeanTransformationMatrix = new Array2DRowRealMatrix(scope.size(), scope.size()); Scope predictionScope = scope.reduceBy(conditioningScope); int[] predictionMapping = scope.createContinuousVariableMapping(predictionScope); int[] conditioningMapping = scope.createContinuousVariableMapping(conditioningScope); for (int i = 0; i < scope.size(); i++) { RealVector precisionColumn = conditionedPrecisionMatrix.getColumnVector(i); if (predictionMapping[i] >= 0) { precisionColumn = precisionColumn.add( padVector(xxMatrix.getColumnVector(predictionMapping[i]), scope.size(), predictionMapping)); conditionedMeanTransformationMatrix.setColumnVector(i, precisionColumn); precisionColumn = precisionColumn.add(padVector(yxMatrix.getColumnVector(predictionMapping[i]), scope.size(), conditioningMapping)); conditionedPrecisionMatrix.setColumnVector(i, precisionColumn); } if (conditioningMapping[i] >= 0) { precisionColumn = precisionColumn.add(padVector(xyMatrix.getColumnVector(conditioningMapping[i]), scope.size(), predictionMapping)); conditionedMeanTransformationMatrix.setColumnVector(i, precisionColumn); precisionColumn = precisionColumn.add(padVector(yyMatrix.getColumnVector(conditioningMapping[i]), scope.size(), conditioningMapping)); conditionedPrecisionMatrix.setColumnVector(i, precisionColumn); } } // h = (a, 0)^T * (XX, XY; 0, 0) RealMatrix scaledMeanMatrix = new Array2DRowRealMatrix(1, scope.size()); scaledMeanMatrix.setRowVector(0, padVector(meanVector, scope.size(), predictionMapping)); scaledMeanMatrix = scaledMeanMatrix.multiply(conditionedMeanTransformationMatrix); RealVector scaledMeanVector = scaledMeanMatrix.getRowVector(0); // g = a^T * C^-1 * a - log((2 * PI) ^ m/2 * det(C)^0.5) where m is the size of the prediction scope RealMatrix meanMatrix = new Array2DRowRealMatrix(predictionScope.size(), 1); meanMatrix.setColumnVector(0, meanVector); double normalizationConstant = -0.5d * meanVector.dotProduct(xxMatrix.operate(meanVector)) - Math.log( Math.pow(2 * Math.PI, (double) predictionScope.size() / 2.0d) * Math.sqrt(mathUtil.determinant())); return new CanonicalGaussianFactor(scope, conditionedPrecisionMatrix, scaledMeanVector, normalizationConstant); }
From source file:edu.cudenver.bios.matrix.MatrixUtils.java
/** * Calculate the Kronecker product of two matrices * * @param matrix1 first matrix/*from ww w. jav a 2s . c o m*/ * @param matrix2 second matrix * @return Kronecker product of matrix 1 and matrix 2 */ public static RealMatrix getKroneckerProduct(RealMatrix matrix1, RealMatrix matrix2) { if (matrix1 == null || matrix2 == null) throw new IllegalArgumentException("null input matrix"); int m1Rows = matrix1.getRowDimension(); int m1Cols = matrix1.getColumnDimension(); int m2Rows = matrix2.getRowDimension(); int m2Cols = matrix2.getColumnDimension(); double[][] productData = new double[m1Rows * m2Rows][m1Cols * m2Cols]; RealMatrix productMatrix = new Array2DRowRealMatrix(productData); for (int col = 0; col < m1Cols; col++) { for (int row = 0; row < m1Rows; row++) { productMatrix.setSubMatrix((matrix2.scalarMultiply(matrix1.getEntry(row, col))).getData(), row * m2Rows, col * m2Cols); } } return productMatrix; }
From source file:edu.washington.gs.skyline.model.quantification.LinearModel.java
public List<LinearFitResult> fit(double[] observations) { if (observations.length != designMatrix.getRowDimension()) { throw new IllegalArgumentException("Wrong number of rows"); }/*w w w . j a v a2s. c o m*/ RealVector coefficients = qrDecomposition.getSolver().solve(new ArrayRealVector(observations)); RealVector fittedValues = new ArrayRealVector(observations.length); RealVector residuals = new ArrayRealVector(observations.length); for (int iRow = 0; iRow < observations.length; iRow++) { RealVector designRow = designMatrix.getRowVector(iRow); fittedValues.setEntry(iRow, designRow.dotProduct(coefficients)); residuals.setEntry(iRow, observations[iRow] - fittedValues.getEntry(iRow)); } double rss = residuals.dotProduct(residuals); int degreesOfFreedom = observations.length - qrDecomposition.getR().getColumnDimension(); double resVar = rss / degreesOfFreedom; double sigma = Math.sqrt(resVar); RealMatrix covarianceUnscaled = matrixCrossproductInverse; RealMatrix scaledCovariance = covarianceUnscaled.scalarMultiply(sigma * sigma); List<LinearFitResult> results = new ArrayList<>(); for (int iContrastRow = 0; iContrastRow < contrastValues.getRowDimension(); iContrastRow++) { RealVector contrastRow = contrastValues.getRowVector(iContrastRow); double standardError = 0; for (int iRow = 0; iRow < independentColumnIndices.length; iRow++) { for (int iCol = 0; iCol < independentColumnIndices.length; iCol++) { standardError = contrastRow.getEntry(independentColumnIndices[iRow]) * scaledCovariance.getEntry(iRow, iCol) * contrastRow.getEntry(independentColumnIndices[iCol]); } } standardError = Math.sqrt(standardError); double foldChange = coefficients.dotProduct(contrastRow); LinearFitResult linearFitResult = new LinearFitResult(foldChange); double tValue = foldChange / standardError; linearFitResult.setTValue(tValue); linearFitResult.setStandardError(standardError); linearFitResult.setDegreesOfFreedom(degreesOfFreedom); if (0 == degreesOfFreedom) { linearFitResult.setPValue(1.0); } else { TDistribution tDistribution = new TDistribution(degreesOfFreedom); double pValue = (1 - tDistribution.cumulativeProbability(Math.abs(tValue))) * 2; linearFitResult.setPValue(pValue); } results.add(linearFitResult); } return results; }
From source file:com.itemanalysis.psychometrics.cfa.ParallelModel.java
public RealMatrix getImpliedCovariance(double[] argument) { setParameters(argument);//from w w w . j av a 2s.co m RealMatrix I = new IdentityMatrix(nItems); B = Ivec.scalarMultiply(argument[0] * argument[0]); RealMatrix THETA = I.scalarMultiply(argument[1]); SIGMA = B.multiply(Ivec.transpose()).add(THETA); return SIGMA; }
From source file:com.analog.lyric.dimple.test.solvers.sumproduct.TestSampledFactors.java
/** * Generates a random covariance matrix with given dimension. *//*from ww w . java2 s. co m*/ RealMatrix randCovariance(int n) { RealMatrix A = MatrixUtils.createRealMatrix(n, n); // Randomize A.walkInOptimizedOrder(new DefaultRealMatrixChangingVisitor() { @Override public double visit(int row, int column, double value) { return testRand.nextDouble(); } }); RealMatrix B = A.add(A.transpose()); // B is symmetric double minEig = Doubles.min(new EigenDecomposition(B).getRealEigenvalues()); double r = testRand.nextGaussian(); r *= r; r += Math.ulp(1.0); RealMatrix I = MatrixUtils.createRealIdentityMatrix(n); RealMatrix C = B.add(I.scalarMultiply(r - minEig)); return C; }
From source file:edu.oregonstate.eecs.mcplan.ml.SequentialProjectionHashLearner.java
@Override public void run() { final RealMatrix cov_reg = MatrixUtils.createRealIdentityMatrix(X.getRowDimension()) .scalarMultiply(shrinkage);/* w w w . java2 s . c o m*/ for (int k = 0; k < K; ++k) { System.out.println("k = " + k); System.out.println("\tCovariance"); final RealMatrix cov = Xi_.multiply(Xi_.transpose()).add(cov_reg); // System.out.println( cov ); System.out.println("\tM"); final RealMatrix M = cov; // XL.multiply( Si_ ).multiply( XLt ).add( cov.scalarMultiply( eta ) ); // TODO: You only need the largest eigenvalue, so the full // decomposition is a ton of unnecessary work. System.out.println("\tM[" + M.getRowDimension() + "x" + M.getColumnDimension() + "]"); final EigenDecomposition evd = new EigenDecomposition(M); final RealVector w = evd.getEigenvector(0); w.mapMultiplyToSelf(1.0 / w.getNorm()); // if( Math.abs( w.getNorm() - 1.0 ) > 1e-2 ) { // System.out.println( "! Non-unit eigenvector: ||w|| = " + w.getNorm() ); // System.out.println( Math.abs( w.getNorm() - 1.0 ) ); // assert( false ); // } W.add(w); final RealMatrix w_wt = w.outerProduct(w); final RealMatrix S_tilde = XLt.multiply(w_wt).multiply(XL); T(S_tilde, Si_); Si_ = Si_.subtract(S_tilde.scalarMultiply(alpha)); Xi_ = Xi_.subtract(w_wt.multiply(Xi_)); } }
From source file:com.joptimizer.functions.SDPLogarithmicBarrier.java
/** * Calculates the initial value for the s parameter in Phase I. * Return s so that F(x)-s.I is negative definite * @see "S.Boyd and L.Vandenberghe, Convex Optimization, 11.6.2" * @see "S.Boyd and L.Vandenberghe, Semidefinite programming, 6.1" *///from ww w . j av a2 s . c o m public double calculatePhase1InitialFeasiblePoint(double[] originalNotFeasiblePoint, double tolerance) { RealMatrix F = this.buildS(originalNotFeasiblePoint).scalarMultiply(-1); RealMatrix S = F.scalarMultiply(-1); try { new CholeskyDecomposition(S); //already feasible return -1; } catch (NonPositiveDefiniteMatrixException ee) { //it does NOT mean that F is negative, it can be not definite EigenDecomposition eFact = new EigenDecomposition(F, Double.NaN); double[] eValues = eFact.getRealEigenvalues(); double minEigenValue = eValues[Utils.getMinIndex(eValues)]; return -Math.min(minEigenValue * Math.pow(tolerance, -0.5), 0.); } }
From source file:edu.cudenver.bios.matrix.GramSchmidtOrthonormalization.java
/** * Perform Gram Schmidt Orthonormalization on the specified * matrix. The matrix A (mxn) is decomposed into two matrices * Q (mxn), R (nxn) such that//from w w w . j a v a2s.co m * <ul> * <li>A = QR * <li>Q'Q = Identity * <li>R is upper triangular * </ul> * The resulting Q, R matrices can be retrieved with the getQ() * and getR() functions. * * @param matrix */ public GramSchmidtOrthonormalization(RealMatrix matrix) { if (matrix == null) throw new IllegalArgumentException("Null matrix"); // create the Q, R matrices int m = matrix.getRowDimension(); int n = matrix.getColumnDimension(); Q = MatrixUtils.createRealMatrix(m, n); R = MatrixUtils.createRealMatrix(n, n); // perform Gram Schmidt process using the following algorithm // let w<n> be the resulting orthonormal column vectors // let v<n> be the columns of the incoming matrix // w1 = (1/norm(v1))*v1 // ... // wj = 1/norm(vj - projectionVj-1Vj)*[vj - projectionVj-1Vj] // where projectionVj-1Vj = (w1 * vj) * w1 + (w2 * vj) * w2 + ... + (wj-1 * vj) * wj-1 // for (int i = 0; i < n; i++) { RealMatrix v = matrix.getColumnMatrix(i); for (int j = 0; j < i; j++) { RealMatrix Qj = Q.getColumnMatrix(j); double value = Qj.transpose().multiply(v).getEntry(0, 0); R.setEntry(j, i, value); v = v.subtract(Qj.scalarMultiply(value)); } double norm = v.getFrobeniusNorm(); R.setEntry(i, i, norm); Q.setColumnMatrix(i, v.scalarMultiply(1 / norm)); } }
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)); }// ww w .j a v 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(); }