List of usage examples for org.apache.commons.math3.linear RealMatrix subtract
RealMatrix subtract(RealMatrix m) throws MatrixDimensionMismatchException;
From source file:hivemall.utils.math.MatrixUtils.java
/** * QR decomposition for a tridiagonal matrix T. * * @see https://gist.github.com/lightcatcher/8118181 * @see http://www.ericmart.in/blog/optimizing_julia_tridiag_qr * @param T target tridiagonal matrix/*from w w w.j ava2s .co m*/ * @param R output matrix for R which is the same shape as T * @param Qt output matrix for Q.T which is the same shape an T */ public static void tridiagonalQR(@Nonnull final RealMatrix T, @Nonnull final RealMatrix R, @Nonnull final RealMatrix Qt) { int n = T.getRowDimension(); Preconditions.checkArgument(n == R.getRowDimension() && n == R.getColumnDimension(), "T and R must be the same shape"); Preconditions.checkArgument(n == Qt.getRowDimension() && n == Qt.getColumnDimension(), "T and Qt must be the same shape"); // initial R = T R.setSubMatrix(T.getData(), 0, 0); // initial Qt = identity Qt.setSubMatrix(eye(n), 0, 0); for (int i = 0; i < n - 1; i++) { // Householder projection for a vector x // https://en.wikipedia.org/wiki/Householder_transformation RealVector x = T.getSubMatrix(i, i + 1, i, i).getColumnVector(0); x = unitL2norm(x); RealMatrix subR = R.getSubMatrix(i, i + 1, 0, n - 1); R.setSubMatrix(subR.subtract(x.outerProduct(subR.preMultiply(x)).scalarMultiply(2)).getData(), i, 0); RealMatrix subQt = Qt.getSubMatrix(i, i + 1, 0, n - 1); Qt.setSubMatrix(subQt.subtract(x.outerProduct(subQt.preMultiply(x)).scalarMultiply(2)).getData(), i, 0); } }
From source file:com.joptimizer.algebra.CholeskyDecompositionTest.java
/** * good decomposition./*from ww w . j a v a 2 s. c om*/ */ public void testDecomposition1() throws Exception { log.debug("testDecomposition1"); RealMatrix P1 = new Array2DRowRealMatrix(new double[][] { { 8.08073550734687, 1.59028724315583 }, { 1.59028724315583, 0.3250861184011492 } }); CholeskyDecomposition cFact1 = new CholeskyDecomposition(P1); log.debug("L: " + cFact1.getL()); log.debug("LT: " + cFact1.getLT()); // check L.LT-Q=0 RealMatrix P1Inv = cFact1.getL().multiply(cFact1.getLT()); double norm1 = P1Inv.subtract(P1).getNorm(); log.debug("norm1: " + norm1); assertTrue(norm1 < 1.E-12); }
From source file:com.itemanalysis.psychometrics.cfa.MaximumLikelihoodEstimation.java
public double gfi() { double fit = 0.0; double q = Double.valueOf(model.getNumberOfItems()).doubleValue(); RealMatrix I = new IdentityMatrix(nItems); LUDecomposition SLUD = new LUDecomposition(SIGMA); RealMatrix Sinv = SLUD.getSolver().getInverse(); RealMatrix P1 = Sinv.multiply(varcov); RealMatrix D = P1.subtract(I); double numerator = D.multiply(D).getTrace(); RealMatrix P2 = Sinv.multiply(varcov); double denom = P2.multiply(P2).getTrace(); fit = 1.0 - numerator / denom;/*from w w w . j ava2 s . c o m*/ return fit; }
From source file:com.itemanalysis.psychometrics.cfa.GeneralizedLeastSquares.java
public double gfi() { double fit = 0.0; double q = Double.valueOf(model.getNumberOfItems()).doubleValue(); RealMatrix I = new IdentityMatrix(nItems); RealMatrix P = SIGMA.multiply(VCinv); RealMatrix D = I.subtract(P); RealMatrix D2 = D.multiply(D);/*w w w.j a v a 2s . c o m*/ fit = 1.0 - D2.getTrace() / q; return fit; }
From source file:lirmm.inria.fr.math.TestUtils.java
/** verifies that two matrices are close (1-norm) */ public static void assertEquals(String msg, RealMatrix expected, RealMatrix observed, double tolerance) { Assert.assertNotNull(msg + "\nObserved should not be null", observed); if (expected.getColumnDimension() != observed.getColumnDimension() || expected.getRowDimension() != observed.getRowDimension()) { StringBuilder messageBuffer = new StringBuilder(msg); messageBuffer.append("\nObserved has incorrect dimensions."); messageBuffer/*from w w w.j a va 2 s.c om*/ .append("\nobserved is " + observed.getRowDimension() + " x " + observed.getColumnDimension()); messageBuffer .append("\nexpected " + expected.getRowDimension() + " x " + expected.getColumnDimension()); Assert.fail(messageBuffer.toString()); } RealMatrix delta = expected.subtract(observed); if (delta.getNorm() >= tolerance) { StringBuilder messageBuffer = new StringBuilder(msg); messageBuffer.append("\nExpected: " + expected); messageBuffer.append("\nObserved: " + observed); messageBuffer.append("\nexpected - observed: " + delta); Assert.fail(messageBuffer.toString()); } }
From source file:edu.ucdenver.bios.designcalculator.DesignCalculator.java
/** * Calculate the sum of squares hypothesis matrix (the H matrix) * @param params matrices input by user/*from w w w . j av a2 s . c om*/ * @return H matrix */ public RealMatrix getHypothesisSumOfSquares(RealMatrix C, RealMatrix beta, RealMatrix U, RealMatrix thetaNull, RealMatrix XtXInverse) { // M = C(X'X)^-1C' RealMatrix M = C.multiply(XtXInverse.multiply(C.transpose())); // thetaHat = C * Beta * U RealMatrix thetaHat = C.multiply(beta.multiply(U)); // thetaHat - thetaNull. Multiple by negative one to do subtraction RealMatrix thetaDiff = thetaHat.subtract(thetaNull); // calculate the hypothesis sum of squares: (thetaHat - thetaNull)'[C(X'X)-1C'](thetaHat - thetaNull) RealMatrix hss = thetaDiff.transpose().multiply(M.multiply(thetaDiff)); return hss; }
From source file:com.github.kingtim1.jmdp.discounted.MatrixInversePolicyEvaluation.java
@Override public DiscountedVFunction<S> eval(StationaryPolicy<S, A> policy) { int n = _smdp.numberOfStates(); List<S> states = new ArrayList<S>(n); Iterable<S> istates = _smdp.states(); for (S state : istates) { states.add(state);/* ww w .j av a 2 s . c om*/ } // Construct matrix A and vector b RealMatrix id = MatrixUtils.createRealIdentityMatrix(n); RealMatrix gpp = gammaPPi(policy, states); RealMatrix A = id.subtract(gpp); RealVector b = rPi(policy, states); // Solve for V^{\pi} SingularValueDecomposition decomp = new SingularValueDecomposition(A); DecompositionSolver dsolver = decomp.getSolver(); RealVector vpi = dsolver.solve(b); // Construct the value function Map<S, Double> valueMap = new HashMap<S, Double>(); for (int i = 0; i < states.size(); i++) { S state = states.get(i); double val = vpi.getEntry(i); valueMap.put(state, val); } return new MapVFunction<S>(valueMap, 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//w w w . j ava2s . c o 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:edu.oregonstate.eecs.mcplan.ml.GaussianMixtureModel.java
private boolean hasConverged(final int c, final double pi, final RealVector mu, final RealMatrix Sigma) { return (pi - pi_[c]) < epsilon_ && mu.getDistance(mu_[c]) < epsilon_ && Sigma.subtract(Sigma_[c]).getFrobeniusNorm() < epsilon_; }
From source file:com.analog.lyric.dimple.test.solvers.sumproduct.TestSampledFactors.java
/** * Adapted from MATLAB test4 in tests/algoGaussian/testSampledFactors.m *//*from w ww . ja v a 2 s. co m*/ @Test public void sampledComplexProduct() { // NOTE: test may fail if seed is changed! We keep the number of samples down so that the test doesn't // take too long. Increasing the samples produces better results. testRand.setSeed(42); try (CurrentModel cur = using(new FactorGraph())) { final Complex a = complex("a"); final Complex b = complex("b"); final Complex c = product(a, b); double[] aMean = new double[] { 10, 10 }; RealMatrix aCovariance = randCovariance(2); a.setPrior(new MultivariateNormal(aMean, aCovariance.getData())); double[] bMean = new double[] { -20, 20 }; RealMatrix bCovariance = randCovariance(2); b.setPrior(new MultivariateNormalParameters(bMean, bCovariance.getData())); GaussianRandomGenerator normalGenerator = new GaussianRandomGenerator(testRand); CorrelatedRandomVectorGenerator aGenerator = new CorrelatedRandomVectorGenerator(aMean, aCovariance, 1e-12, normalGenerator); CorrelatedRandomVectorGenerator bGenerator = new CorrelatedRandomVectorGenerator(bMean, bCovariance, 1e-12, normalGenerator); StorelessCovariance expectedCov = new StorelessCovariance(2); final int nSamples = 10000; RealVector expectedMean = MatrixUtils.createRealVector(new double[2]); double[] cSample = new double[2]; for (int i = 0; i < nSamples; ++i) { double[] aSample = aGenerator.nextVector(); double[] bSample = bGenerator.nextVector(); // Compute complex product cSample[0] = aSample[0] * bSample[0] - aSample[1] * bSample[1]; cSample[1] = aSample[0] * bSample[1] + aSample[1] * bSample[0]; expectedMean.addToEntry(0, cSample[0]); expectedMean.addToEntry(1, cSample[1]); expectedCov.increment(cSample); } expectedMean.mapDivideToSelf(nSamples); // normalize SumProductSolverGraph sfg = requireNonNull(cur.graph.setSolverFactory(new SumProductSolver())); sfg.setOption(GibbsOptions.numSamples, nSamples); sfg.solve(); MultivariateNormalParameters cBelief = requireNonNull(c.getBelief()); RealVector observedMean = MatrixUtils.createRealVector(cBelief.getMean()); double scaledMeanDistance = expectedMean.getDistance(observedMean) / expectedMean.getNorm(); // System.out.format("expectedMean = %s\n", expectedMean); // System.out.format("observedMean = %s\n", observedMean); // System.out.println(scaledMeanDistance); assertEquals(0.0, scaledMeanDistance, .02); RealMatrix expectedCovariance = expectedCov.getCovarianceMatrix(); RealMatrix observedCovariance = MatrixUtils.createRealMatrix(cBelief.getCovariance()); RealMatrix diffCovariance = expectedCovariance.subtract(observedCovariance); double scaledCovarianceDistance = diffCovariance.getNorm() / expectedCovariance.getNorm(); // System.out.println(expectedCovariance); // System.out.println(expectedCovariance.getNorm()); // System.out.println(diffCovariance); // System.out.println(diffCovariance.getNorm()); // System.out.println(diffCovariance.getNorm() / expectedCovariance.getNorm()); assertEquals(0.0, scaledCovarianceDistance, .2); } }