Example usage for org.apache.commons.math3.linear RealMatrix preMultiply

List of usage examples for org.apache.commons.math3.linear RealMatrix preMultiply

  1. HOME
  2. Java
  3. org.apache.commons
  4. org.apache.commons.math3.linear.*
  5. RealMatrix
  6. preMultiply

Introduction

In this page you can find the example usage for org.apache.commons.math3.linear RealMatrix preMultiply.

Prototype

RealVector preMultiply(RealVector v) throws DimensionMismatchException;

Source Link

Document

Returns the (row) vector result of premultiplying this by the vector v .

Usage

From source file:edu.byu.nlp.stats.DirichletMLEOptimizableTest.java

/**
 * Computes a Newton-Raphson update in-place to alpha.
 *//*w  ww. j  ava2s .  c  o  m*/
private static RealVector newtonRaphsonUpdate(final double[][] data, double[] alpha) {
    // We'll compute the gold-standard value the "long" way (taking the inverse of the Hessian)
    RealMatrix hessian = new Array2DRowRealMatrix(alpha.length, alpha.length);
    for (int r = 0; r < hessian.getRowDimension(); r++) {
        for (int c = 0; c < hessian.getColumnDimension(); c++) {
            hessian.addToEntry(r, c, data.length * Gamma.trigamma(DoubleArrays.sum(alpha)));
            if (r == c) {
                hessian.addToEntry(r, c, -data.length * Gamma.trigamma(alpha[r]));
            }
        }
    }
    RealVector derivative = new ArrayRealVector(alpha.length);
    for (int k = 0; k < alpha.length; k++) {
        derivative.setEntry(k,
                data.length * (Gamma.digamma(DoubleArrays.sum(alpha)) - Gamma.digamma(alpha[k])));
        for (double[] theta : data) {
            derivative.addToEntry(k, theta[k]);
        }
    }

    RealMatrix hessianInverse = new LUDecomposition(hessian).getSolver().getInverse();
    RealVector negDiff = hessianInverse.preMultiply(derivative);
    negDiff.mapMultiplyToSelf(-1.0);

    RealVector expected = new ArrayRealVector(alpha, true);
    return expected.add(negDiff);
}

From source file:hivemall.utils.math.StatsUtils.java

/**
 * @param mu1 mean vector of the first normal distribution
 * @param sigma1 covariance matrix of the first normal distribution
 * @param mu2 mean vector of the second normal distribution
 * @param sigma2 covariance matrix of the second normal distribution
 * @return the Hellinger distance between two multivariate normal distributions
 * @link https://en.wikipedia.org/wiki/Hellinger_distance#Examples
 *///from   www  .  j  a  v a  2s  .c o  m
public static double hellingerDistance(@Nonnull final RealVector mu1, @Nonnull final RealMatrix sigma1,
        @Nonnull final RealVector mu2, @Nonnull final RealMatrix sigma2) {
    RealVector muSub = mu1.subtract(mu2);
    RealMatrix sigmaMean = sigma1.add(sigma2).scalarMultiply(0.5d);
    LUDecomposition LUsigmaMean = new LUDecomposition(sigmaMean);
    double denominator = Math.sqrt(LUsigmaMean.getDeterminant());
    if (denominator == 0.d) {
        return 1.d; // avoid divide by zero
    }
    RealMatrix sigmaMeanInv = LUsigmaMean.getSolver().getInverse(); // has inverse iff det != 0
    double sigma1Det = MatrixUtils.det(sigma1);
    double sigma2Det = MatrixUtils.det(sigma2);

    double numerator = Math.pow(sigma1Det, 0.25d) * Math.pow(sigma2Det, 0.25d)
            * Math.exp(-0.125d * sigmaMeanInv.preMultiply(muSub).dotProduct(muSub));
    return 1.d - numerator / denominator;
}

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  ww  w.j  a  v a2  s . c o 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: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 w  w.  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:de.andreasschoknecht.LS3.Query.java

/**
 * Calculate pseudo document by applying the formula p^T * Uk * Sk^-1.
 *
 * @param Uk The matrix of term vectors Uk
 * @param Sk The matrix of singular values Sk
 *//*w  w w  . j a v  a2s. c o  m*/
void calculatePseudoDocument(RealMatrix Uk, RealMatrix Sk) {
    // temp = q^T * Uk: calculate transpose of q times Uk
    pseudoDocument = Uk.preMultiply(weightedTermFrequencies);

    // calculate inverse of Sk
    RealMatrix inverseSk = MatrixUtils.inverse(Sk);

    // pseudoDocument = temp * inverseSk: calculate pseudoDocument by multiplying temp and inverseSk.
    pseudoDocument = inverseSk.preMultiply(pseudoDocument);
}

From source file:com.opengamma.strata.math.impl.regression.WeightedLeastSquaresRegression.java

public LeastSquaresRegressionResult regress(double[][] x, double[] weights, double[] y, boolean useIntercept) {
    if (weights == null) {
        throw new IllegalArgumentException("Cannot perform WLS regression without an array of weights");
    }//from  w  w w  .j av  a 2s .com
    checkData(x, weights, y);
    double[][] dep = addInterceptVariable(x, useIntercept);
    double[] w = new double[weights.length];
    for (int i = 0; i < y.length; i++) {
        w[i] = weights[i];
    }
    DoubleMatrix matrix = DoubleMatrix.copyOf(dep);
    DoubleArray vector = DoubleArray.copyOf(y);
    RealMatrix wDiag = new DiagonalMatrix(w);
    DoubleMatrix transpose = s_algebra.getTranspose(matrix);

    DoubleMatrix wDiagTimesMatrix = DoubleMatrix
            .ofUnsafe(wDiag.multiply(new Array2DRowRealMatrix(matrix.toArrayUnsafe())).getData());
    DoubleMatrix tmp = (DoubleMatrix) s_algebra
            .multiply(s_algebra.getInverse(s_algebra.multiply(transpose, wDiagTimesMatrix)), transpose);

    DoubleMatrix wTmpTimesDiag = DoubleMatrix
            .copyOf(wDiag.preMultiply(new Array2DRowRealMatrix(tmp.toArrayUnsafe())).getData());
    DoubleMatrix betasVector = (DoubleMatrix) s_algebra.multiply(wTmpTimesDiag, vector);
    double[] yModel = super.writeArrayAsVector(
            ((DoubleMatrix) s_algebra.multiply(matrix, betasVector)).toArray());
    double[] betas = super.writeArrayAsVector(betasVector.toArray());
    return getResultWithStatistics(x, convertArray(wDiag.getData()), y, betas, yModel, transpose, matrix,
            useIntercept);
}

From source file:iDynoOptimizer.MOEAFramework26.src.org.moeaframework.core.operator.real.AdaptiveMetropolis.java

@Override
public Solution[] evolve(Solution[] parents) {
    int k = parents.length;
    int n = parents[0].getNumberOfVariables();
    RealMatrix x = new Array2DRowRealMatrix(k, n);

    for (int i = 0; i < k; i++) {
        x.setRow(i, EncodingUtils.getReal(parents[i]));
    }/* w  w  w  . j av  a2s .  c  o  m*/

    try {
        //perform Cholesky factorization and get the upper triangular matrix
        double jumpRate = Math.pow(jumpRateCoefficient / Math.sqrt(n), 2.0);

        RealMatrix chol = new CholeskyDecomposition(
                new Covariance(x.scalarMultiply(jumpRate)).getCovarianceMatrix()).getLT();

        //produce the offspring
        Solution[] offspring = new Solution[numberOfOffspring];

        for (int i = 0; i < numberOfOffspring; i++) {
            Solution child = parents[PRNG.nextInt(parents.length)].copy();

            //apply adaptive metropolis step to solution
            RealVector muC = new ArrayRealVector(EncodingUtils.getReal(child));
            RealVector ru = new ArrayRealVector(n);

            for (int j = 0; j < n; j++) {
                ru.setEntry(j, PRNG.nextGaussian());
            }

            double[] variables = muC.add(chol.preMultiply(ru)).toArray();

            //assign variables back to solution
            for (int j = 0; j < n; j++) {
                RealVariable variable = (RealVariable) child.getVariable(j);
                double value = variables[j];

                if (value < variable.getLowerBound()) {
                    value = variable.getLowerBound();
                } else if (value > variable.getUpperBound()) {
                    value = variable.getUpperBound();
                }

                variable.setValue(value);
            }

            offspring[i] = child;
        }

        return offspring;
    } catch (Exception e) {
        return new Solution[0];
    }
}

From source file:lirmm.inria.fr.math.BigSparseRealMatrixTest.java

/**
 * test preMultiply by vector//from w w w  . j  a va  2  s  . c  om
 */
@Test
public void testPremultiplyVector() {
    RealMatrix m = new BigSparseRealMatrix(testData);
    TestUtils.assertEquals("premultiply", m.preMultiply(testVector), preMultTest, normTolerance);
    TestUtils.assertEquals("premultiply", m.preMultiply(new ArrayRealVector(testVector).toArray()), preMultTest,
            normTolerance);
    m = new BigSparseRealMatrix(bigSingular);
    try {
        m.preMultiply(testVector);
        Assert.fail("expecting MathIllegalArgumentException");
    } catch (MathIllegalArgumentException ex) {
        // ignored
    }
}

From source file:lirmm.inria.fr.math.BigSparseRealMatrixTest.java

@Test
public void testPremultiply() {
    RealMatrix m3 = new BigSparseRealMatrix(d3);
    RealMatrix m4 = new BigSparseRealMatrix(d4);
    RealMatrix m5 = new BigSparseRealMatrix(d5);
    TestUtils.assertEquals("m3*m4=m5", m4.preMultiply(m3), m5, entryTolerance);

    BigSparseRealMatrix m = new BigSparseRealMatrix(testData);
    BigSparseRealMatrix mInv = new BigSparseRealMatrix(testDataInv);
    BigSparseRealMatrix identity = new BigSparseRealMatrix(id);
    TestUtils.assertEquals("inverse multiply", m.preMultiply(mInv), identity, entryTolerance);
    TestUtils.assertEquals("inverse multiply", mInv.preMultiply(m), identity, entryTolerance);
    TestUtils.assertEquals("identity multiply", m.preMultiply(identity), m, entryTolerance);
    TestUtils.assertEquals("identity multiply", identity.preMultiply(mInv), mInv, entryTolerance);
    try {// w  w  w .  j  a va  2  s .com
        m.preMultiply(new BigSparseRealMatrix(bigSingular));
        Assert.fail("Expecting illegalArgumentException");
    } catch (MathIllegalArgumentException ex) {
        // ignored
    }
}

From source file:com.analog.lyric.dimple.solvers.core.parameterizedMessages.MultivariateNormalParameters.java

/**
 * {@inheritDoc}//  w w  w.j  a v a 2 s. com
 * <p>
 * For multivariate normal distributions, the formula is given by:
 * 
 * <blockquote>
 * &frac12; {
 * trace(&Sigma;<sub>Q</sub><sup>-1</sup>&Sigma;<sub>P</sub>) +
 * (&mu;<sub>Q</sub>-&mu;<sub>P</sub>)<sup>T</sup>&Sigma;<sub>Q</sub><sup>-1</sup>(&mu;<sub>Q</sub>-&mu;<sub>P</sub>)
 * -K - ln(det(&Sigma;<sub>P</sub>)/det(&Sigma;<sub>Q</sub>)))
 * }
 * </blockquote>
 * Note that this assumes that the determinants of the covariance matrices are non-zero.
 */
@Override
public double computeKLDivergence(IParameterizedMessage that) {
    if (that instanceof MultivariateNormalParameters) {
        // http://en.wikipedia.org/wiki/Multivariate_normal_distribution#Kullback.E2.80.93Leibler_divergence
        //
        // K: size of vectors, # rows/columns of matrices
        // up, uq: vectors of means for P and Q
        // CP, CQ: covariance matrices for P and Q
        // inv(x): inverse of x
        // det(x): determinant of x
        // tr(x): trace of x
        // x': transpose of x
        //
        // KL(P|Q) == .5 * ( tr(inv(CQ) * CP) + (uq - up)' * inv(CQ) * (uq - up) - K - ln(det(CP)/det(CQ)) )
        //

        final MultivariateNormalParameters P = this, Q = (MultivariateNormalParameters) that;
        final int K = P.getVectorLength();
        assertSameSize(Q.getVectorLength());

        if (P.isDiagonal() && Q.isDiagonal()) {
            // If both are diagonal, we can simply add up the KL for the univariate cases along the diagonal.

            final double[] Pmeans = P._mean, Qmeans = Q._mean;
            final double[] Pprecisions = P._precision, Qprecisions = Q._precision;

            double kl = 0.0;

            for (int i = 0; i < K; ++i) {
                kl += NormalParameters.computeKLDiverence(Pmeans[i], Pprecisions[i], Qmeans[i], Qprecisions[i]);
            }

            return kl;
        }

        // TODO - if we ever start storing the eigendecomposition in this object, we can simply use
        // the eigenvalues to efficiently compute the trace and determinants. Perhaps it would be worthwhile
        // to save the eigenvalues if nothing else.

        RealVector mP = wrapRealVector(P.getMean());
        RealVector mQ = wrapRealVector(Q.getMean());
        RealMatrix CP = wrapRealMatrix(P.getCovariance());
        RealMatrix CQinv = wrapRealMatrix(Q.getInformationMatrix());

        RealVector mdiff = mQ.subtract(mP);

        // FIXME: do we need to worry about singular covariance matrices?

        double divergence = -K;

        // trace of product of matrices is equivalent to the dot-product of the vectorized versions
        // of the matrices - this is much faster than doing the actual matrix product
        // divergence += CQinv.multiply(CP).trace();
        for (int i = 0; i < K; ++i)
            for (int j = 0; j < K; ++j)
                divergence += CQinv.getEntry(i, j) * CP.getEntry(i, j);

        divergence += CQinv.preMultiply(mdiff).dotProduct(mdiff);
        divergence -= Math.log(
                new EigenDecomposition(CP).getDeterminant() * new EigenDecomposition(CQinv).getDeterminant());
        return Math.abs(divergence / 2); // use abs to guard against precision errors causing this to go negative.
    }

    throw new IllegalArgumentException(
            String.format("Expected '%s' but got '%s'", getClass(), that.getClass()));
}