Java Code Examples for org.apache.commons.math3.linear.RealMatrix#multiply()

@Test
public void testTridiagonalEigen() {
    // Tridiagonal Matrix
    RealMatrix T = new Array2DRowRealMatrix(
        new double[][] {new double[] {40, 60, 0, 0}, new double[] {60, 10, 120, 0},
                new double[] {0, 120, 10, 120}, new double[] {0, 0, 120, 10}});

    double[] eigvals = new double[4];
    RealMatrix eigvecs = new Array2DRowRealMatrix(new double[4][4]);

    MatrixUtils.tridiagonalEigen(T, 2, eigvals, eigvecs);

    RealMatrix actual = eigvecs.multiply(eigvecs.transpose());

    RealMatrix expected = new Array2DRowRealMatrix(new double[4][4]);
    for (int i = 0; i < 4; i++) {
        expected.setEntry(i, i, 1);
    }

    for (int i = 0; i < 4; i++) {
        for (int j = 0; j < 4; j++) {
            Assert.assertEquals(expected.getEntry(i, j), actual.getEntry(i, j), 0.001d);
        }
    }
}
 

public CorrelatedRandomVectorGeneratorTest() {
    mean = new double[] { 0.0, 1.0, -3.0, 2.3 };

    RealMatrix b = MatrixUtils.createRealMatrix(4, 3);
    int counter = 0;
    for (int i = 0; i < b.getRowDimension(); ++i) {
        for (int j = 0; j < b.getColumnDimension(); ++j) {
            b.setEntry(i, j, 1.0 + 0.1 * ++counter);
        }
    }
    RealMatrix bbt = b.multiply(b.transpose());
    covariance = MatrixUtils.createRealMatrix(mean.length, mean.length);
    for (int i = 0; i < covariance.getRowDimension(); ++i) {
        covariance.setEntry(i, i, bbt.getEntry(i, i));
        for (int j = 0; j < covariance.getColumnDimension(); ++j) {
            double s = bbt.getEntry(i, j);
            covariance.setEntry(i, j, s);
            covariance.setEntry(j, i, s);
        }
    }

    RandomGenerator rg = new JDKRandomGenerator();
    rg.setSeed(17399225432l);
    GaussianRandomGenerator rawGenerator = new GaussianRandomGenerator(rg);
    generator = new CorrelatedRandomVectorGenerator(mean,
                                                    covariance,
                                                    1.0e-12 * covariance.getNorm(),
                                                    rawGenerator);
}
 

public CorrelatedRandomVectorGeneratorTest() {
    mean = new double[] { 0.0, 1.0, -3.0, 2.3 };

    RealMatrix b = MatrixUtils.createRealMatrix(4, 3);
    int counter = 0;
    for (int i = 0; i < b.getRowDimension(); ++i) {
        for (int j = 0; j < b.getColumnDimension(); ++j) {
            b.setEntry(i, j, 1.0 + 0.1 * ++counter);
        }
    }
    RealMatrix bbt = b.multiply(b.transpose());
    covariance = MatrixUtils.createRealMatrix(mean.length, mean.length);
    for (int i = 0; i < covariance.getRowDimension(); ++i) {
        covariance.setEntry(i, i, bbt.getEntry(i, i));
        for (int j = 0; j < covariance.getColumnDimension(); ++j) {
            double s = bbt.getEntry(i, j);
            covariance.setEntry(i, j, s);
            covariance.setEntry(j, i, s);
        }
    }

    RandomGenerator rg = new JDKRandomGenerator();
    rg.setSeed(17399225432l);
    GaussianRandomGenerator rawGenerator = new GaussianRandomGenerator(rg);
    generator = new CorrelatedRandomVectorGenerator(mean,
                                                    covariance,
                                                    1.0e-12 * covariance.getNorm(),
                                                    rawGenerator);
}
 

@Test
void givenTwoMatrices_whenMultiply_thenMultiplicatedMatrix() {
    RealMatrix firstMatrix = new Array2DRowRealMatrix(
      new double[][] {
        new double[] {1d, 5d},
        new double[] {2d, 3d},
        new double[] {1d ,7d}
      }
    );

    RealMatrix secondMatrix = new Array2DRowRealMatrix(
      new double[][] {
        new double[] {1d, 2d, 3d, 7d},
        new double[] {5d, 2d, 8d, 1d}
      }
    );

    RealMatrix expected = new Array2DRowRealMatrix(
      new double[][] {
        new double[] {26d, 12d, 43d, 12d},
        new double[] {17d, 10d, 30d, 17d},
        new double[] {36d, 16d, 59d, 14d}
      }
    );

    RealMatrix actual = firstMatrix.multiply(secondMatrix);

    assertThat(actual).isEqualTo(expected);
}
 

@Test
public void testRootMatrix() {
    RealMatrix b = generator.getRootMatrix();
    RealMatrix bbt = b.multiply(b.transpose());
    for (int i = 0; i < covariance.getRowDimension(); ++i) {
        for (int j = 0; j < covariance.getColumnDimension(); ++j) {
            Assert.assertEquals(covariance.getEntry(i, j), bbt.getEntry(i, j), 1.0e-12);
        }
    }
}
 

@Test
public void testRootMatrix() {
    RealMatrix b = generator.getRootMatrix();
    RealMatrix bbt = b.multiply(b.transpose());
    for (int i = 0; i < covariance.getRowDimension(); ++i) {
        for (int j = 0; j < covariance.getColumnDimension(); ++j) {
            Assert.assertEquals(covariance.getEntry(i, j), bbt.getEntry(i, j), 1.0e-12);
        }
    }
}
 

public CorrelatedRandomVectorGeneratorTest() {
    mean = new double[] { 0.0, 1.0, -3.0, 2.3 };

    RealMatrix b = MatrixUtils.createRealMatrix(4, 3);
    int counter = 0;
    for (int i = 0; i < b.getRowDimension(); ++i) {
        for (int j = 0; j < b.getColumnDimension(); ++j) {
            b.setEntry(i, j, 1.0 + 0.1 * ++counter);
        }
    }
    RealMatrix bbt = b.multiply(b.transpose());
    covariance = MatrixUtils.createRealMatrix(mean.length, mean.length);
    for (int i = 0; i < covariance.getRowDimension(); ++i) {
        covariance.setEntry(i, i, bbt.getEntry(i, i));
        for (int j = 0; j < covariance.getColumnDimension(); ++j) {
            double s = bbt.getEntry(i, j);
            covariance.setEntry(i, j, s);
            covariance.setEntry(j, i, s);
        }
    }

    RandomGenerator rg = new JDKRandomGenerator();
    rg.setSeed(17399225432l);
    GaussianRandomGenerator rawGenerator = new GaussianRandomGenerator(rg);
    generator = new CorrelatedRandomVectorGenerator(mean,
                                                    covariance,
                                                    1.0e-12 * covariance.getNorm(),
                                                    rawGenerator);
}
 

/**
 * Asserts that a collection of beta-hats corresponds to the expected value given
 * the input pre-tangent normalization matrix and the PoN file.
 */
private void assertBetaHats(final ReadCountCollection preTangentNormalized,
                            final RealMatrix actual, final File ponFile) {
    Assert.assertEquals(actual.getColumnDimension(), preTangentNormalized.columnNames().size());
    final double epsilon = PCATangentNormalizationUtils.EPSILON;

    try (final HDF5File ponReader = new HDF5File(ponFile)) {
        final PCACoveragePoN pon = new HDF5PCACoveragePoN(ponReader);
        final List<String> ponTargets = pon.getPanelTargetNames();
        final RealMatrix inCounts = reorderTargetsToPoNOrder(preTangentNormalized, ponTargets);

        // obtain subset of relevant targets to calculate the beta-hats;
        final int[][] betaHatTargets = new int[inCounts.getColumnDimension()][];
        for (int i = 0; i < inCounts.getColumnDimension(); i++) {
            final List<Integer> relevantTargets = new ArrayList<>();
            for (int j = 0; j < inCounts.getRowDimension(); j++) {
                if (inCounts.getEntry(j, i) > 1  +  (Math.log(epsilon) / Math.log(2))) {
                    relevantTargets.add(j);
                }
            }
            betaHatTargets[i] = relevantTargets.stream().mapToInt(Integer::intValue).toArray();
        }
        // calculate beta-hats per column and check with actual values.
        final RealMatrix normalsInv = pon.getReducedPanelPInverseCounts();
        Assert.assertEquals(actual.getRowDimension(), normalsInv.getRowDimension());
        final RealMatrix normalsInvT = normalsInv.transpose();
        for (int i = 0; i < inCounts.getColumnDimension(); i++) {
            final RealMatrix inValues = inCounts.getColumnMatrix(i).transpose().getSubMatrix(new int[] { 0 }, betaHatTargets[i]);
            final RealMatrix normalValues = normalsInvT.getSubMatrix(betaHatTargets[i], IntStream.range(0, normalsInvT.getColumnDimension()).toArray());
            final RealMatrix betaHats = inValues.multiply(normalValues);
            for (int j = 0; j < actual.getRowDimension(); j++) {
                Assert.assertEquals(actual.getEntry(j, i), betaHats.getEntry(0, j),0.000001,"Col " + i + " row " + j);
            }
        }
    }
}
 

@Test
public void testRootMatrix() {
    RealMatrix b = generator.getRootMatrix();
    RealMatrix bbt = b.multiply(b.transpose());
    for (int i = 0; i < covariance.getRowDimension(); ++i) {
        for (int j = 0; j < covariance.getColumnDimension(); ++j) {
            Assert.assertEquals(covariance.getEntry(i, j), bbt.getEntry(i, j), 1.0e-12);
        }
    }
}
 

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];
				}
				double diff = (dY[i][j] - preddY)/(1.0 - H[i][i]);
				sum += (diff*diff);
				count += 1;
			}
		}
		sum/=count;
		return sum;
	}
 

@Override
public double looPredLogLikelihood(){

	double[][] Y = getY();
	double[][] dY = getdY();
	double[][] mask = getMask();

	int ns = Y.length;
	int ny = Y[0].length;

	RealMatrix Ky = MatrixUtils.createRealMatrix(this.Ky);
	RealMatrix invK = 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 = invK.multiply(dy);
		double sum=0;
		for(int i=0; i<ns; i++){
			double mu_i = dYmat.getEntry(i, j) - invKdy.getEntry(i, 0)/invK.getEntry(i, i);
			double sigma_i2 = 1/invK.getEntry(i, i);
			double logLL = StatUtils.logProbaNormal(dYmat.getEntry(i, j), mu_i, Math.sqrt(sigma_i2));
			sum += logLL * mask[i][j];
		}
		LOOPredLL[j] = sum;
	}

	double sumLOOPredLL=0;
	for(int j=0;j<ny; j++){
		sumLOOPredLL += LOOPredLL[j];
	}

	return sumLOOPredLL;
}
 

@Override
public double looMSE(){

	double[][] Y = getY();
	double[][] dY = getdY();
	double[][] mask = getMask();

	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(); // H = (K+sI)^(-1) K = invKy K

	double sum=0;
	double count=0;
	for(int j=0;j<ny; j++){
		for(int i=0; i<ns; i++){
			if(mask[i][j]>0.0){
				double preddY=0;
				for(int k=0; k<ns; k++){
					preddY += H[i][k]*dY[k][j];
				}
				double diff = (dY[i][j] - preddY)/(1.0 - H[i][i]);
				sum += (diff*diff) * mask[i][j];
				count += mask[i][j];
			}
		}
	}
	sum/=count;
	return sum;
}
 

@Test
public void testRootMatrix() {
    RealMatrix b = generator.getRootMatrix();
    RealMatrix bbt = b.multiply(b.transpose());
    for (int i = 0; i < covariance.getRowDimension(); ++i) {
        for (int j = 0; j < covariance.getColumnDimension(); ++j) {
            Assert.assertEquals(covariance.getEntry(i, j), bbt.getEntry(i, j), 1.0e-12);
        }
    }
}
 

/**
 * Find eigenvalues and eigenvectors of given tridiagonal matrix T.
 *
 * @see http://web.csulb.edu/~tgao/math423/s94.pdf
 * @see http://stats.stackexchange.com/questions/20643/finding-matrix-eigenvectors-using-qr-
 *      decomposition
 * @param T target tridiagonal matrix
 * @param nIter number of iterations for the QR method
 * @param eigvals eigenvalues are stored here
 * @param eigvecs eigenvectors are stored here
 */
public static void tridiagonalEigen(@Nonnull final RealMatrix T, @Nonnull final int nIter,
        @Nonnull final double[] eigvals, @Nonnull final RealMatrix eigvecs) {
    Preconditions.checkArgument(Arrays.deepEquals(T.getData(), T.transpose().getData()),
        "Target matrix T must be a symmetric (tridiagonal) matrix");
    Preconditions.checkArgument(eigvecs.getRowDimension() == eigvecs.getColumnDimension(),
        "eigvecs must be a square matrix");
    Preconditions.checkArgument(T.getRowDimension() == eigvecs.getRowDimension(),
        "T and eigvecs must be the same shape");
    Preconditions.checkArgument(eigvals.length == eigvecs.getRowDimension(),
        "Number of eigenvalues and eigenvectors must be same");

    int nEig = eigvals.length;

    // initialize eigvecs as an identity matrix
    eigvecs.setSubMatrix(eye(nEig), 0, 0);

    RealMatrix T_ = T.copy();

    for (int i = 0; i < nIter; i++) {
        // QR decomposition for the tridiagonal matrix T
        RealMatrix R = new Array2DRowRealMatrix(nEig, nEig);
        RealMatrix Qt = new Array2DRowRealMatrix(nEig, nEig);
        tridiagonalQR(T_, R, Qt);

        RealMatrix Q = Qt.transpose();
        T_ = R.multiply(Q);
        eigvecs.setSubMatrix(eigvecs.multiply(Q).getData(), 0, 0);
    }

    // diagonal elements correspond to the eigenvalues
    for (int i = 0; i < nEig; i++) {
        eigvals[i] = T_.getEntry(i, i);
    }
}
 

@Test
public void testDerivatives() {

    double eps      = 5.0e-16;
    double kx       = 2;
    double ky       = -3;
    double kz       = 5;
    double n2       = kx * kx + ky * ky + kz * kz;
    double n        = FastMath.sqrt(n2);
    double theta    = 1.7;
    double cosTheta = FastMath.cos(theta);
    double sinTheta = FastMath.sin(theta);
    FieldRotation<DerivativeStructure> r    = new FieldRotation<DerivativeStructure>(createAxis(kx, ky, kz), createAngle(theta));
    Vector3D a      = new Vector3D(kx / n, ky / n, kz / n);

    // Jacobian of the normalized rotation axis a with respect to the Cartesian vector k
    RealMatrix dadk = MatrixUtils.createRealMatrix(new double[][] {
        { (ky * ky + kz * kz) / ( n * n2),            -kx * ky / ( n * n2),            -kx * kz / ( n * n2) },
        {            -kx * ky / ( n * n2), (kx * kx + kz * kz) / ( n * n2),            -ky * kz / ( n * n2) },
        {            -kx * kz / ( n * n2),            -ky * kz / ( n * n2), (kx * kx + ky * ky) / ( n * n2) }
    });

    for (double x = -0.9; x < 0.9; x += 0.2) {
        for (double y = -0.9; y < 0.9; y += 0.2) {
            for (double z = -0.9; z < 0.9; z += 0.2) {
                Vector3D   u = new Vector3D(x, y, z);
                FieldVector3D<DerivativeStructure> v = r.applyTo(createVector(x, y, z));

                // explicit formula for rotation of vector u around axis a with angle theta
                double dot     = Vector3D.dotProduct(u, a);
                Vector3D cross = Vector3D.crossProduct(a, u);
                double c1      = 1 - cosTheta;
                double c2      = c1 * dot;
                Vector3D rt    = new Vector3D(cosTheta, u, c2, a, sinTheta, cross);
                Assert.assertEquals(rt.getX(), v.getX().getReal(), eps);
                Assert.assertEquals(rt.getY(), v.getY().getReal(), eps);
                Assert.assertEquals(rt.getZ(), v.getZ().getReal(), eps);

                // Jacobian of the image v = r(u) with respect to rotation axis a
                // (analytical differentiation of the explicit formula)
                RealMatrix dvda = MatrixUtils.createRealMatrix(new double[][] {
                    { c1 * x * a.getX() + c2,           c1 * y * a.getX() + sinTheta * z, c1 * z * a.getX() - sinTheta * y },
                    { c1 * x * a.getY() - sinTheta * z, c1 * y * a.getY() + c2,           c1 * z * a.getY() + sinTheta * x },
                    { c1 * x * a.getZ() + sinTheta * y, c1 * y * a.getZ() - sinTheta * x, c1 * z * a.getZ() + c2           }
                });

                // compose Jacobians
                RealMatrix dvdk = dvda.multiply(dadk);

                // derivatives with respect to un-normalized axis
                Assert.assertEquals(dvdk.getEntry(0, 0), v.getX().getPartialDerivative(1, 0, 0, 0), eps);
                Assert.assertEquals(dvdk.getEntry(0, 1), v.getX().getPartialDerivative(0, 1, 0, 0), eps);
                Assert.assertEquals(dvdk.getEntry(0, 2), v.getX().getPartialDerivative(0, 0, 1, 0), eps);
                Assert.assertEquals(dvdk.getEntry(1, 0), v.getY().getPartialDerivative(1, 0, 0, 0), eps);
                Assert.assertEquals(dvdk.getEntry(1, 1), v.getY().getPartialDerivative(0, 1, 0, 0), eps);
                Assert.assertEquals(dvdk.getEntry(1, 2), v.getY().getPartialDerivative(0, 0, 1, 0), eps);
                Assert.assertEquals(dvdk.getEntry(2, 0), v.getZ().getPartialDerivative(1, 0, 0, 0), eps);
                Assert.assertEquals(dvdk.getEntry(2, 1), v.getZ().getPartialDerivative(0, 1, 0, 0), eps);
                Assert.assertEquals(dvdk.getEntry(2, 2), v.getZ().getPartialDerivative(0, 0, 1, 0), eps);

                // derivative with respect to rotation angle
                // (analytical differentiation of the explicit formula)
                Vector3D dvdTheta =
                        new Vector3D(-sinTheta, u, sinTheta * dot, a, cosTheta, cross);
                Assert.assertEquals(dvdTheta.getX(), v.getX().getPartialDerivative(0, 0, 0, 1), eps);
                Assert.assertEquals(dvdTheta.getY(), v.getY().getPartialDerivative(0, 0, 0, 1), eps);
                Assert.assertEquals(dvdTheta.getZ(), v.getZ().getPartialDerivative(0, 0, 0, 1), eps);

            }
        }
    }
 }
 

public static double[][] multiply(double[][] X1, double[][] X2){
	RealMatrix M1 = MatrixUtils.createRealMatrix(X1);
	RealMatrix M2 = MatrixUtils.createRealMatrix(X2);
	RealMatrix M3 = M1.multiply(M2);
	return M3.getData();
}
 

/**
 * <p>Calculates the variance-covariance matrix of the regression parameters.
 * </p>
 * <p>Var(b) = (X<sup>T</sup>X)<sup>-1</sup>
 * </p>
 * <p>Uses QR decomposition to reduce (X<sup>T</sup>X)<sup>-1</sup>
 * to (R<sup>T</sup>R)<sup>-1</sup>, with only the top p rows of
 * R included, where p = the length of the beta vector.</p>
 *
 * <p>Data for the model must have been successfully loaded using one of
 * the {@code newSampleData} methods before invoking this method; otherwise
 * a {@code NullPointerException} will be thrown.</p>
 *
 * @return The beta variance-covariance matrix
 */
@Override
protected RealMatrix calculateBetaVariance() {
    int p = getX().getColumnDimension();
    RealMatrix Raug = qr.getR().getSubMatrix(0, p - 1 , 0, p - 1);
    RealMatrix Rinv = new LUDecomposition(Raug).getSolver().getInverse();
    return Rinv.multiply(Rinv.transpose());
}
 

/**
 * <p>Calculates the variance-covariance matrix of the regression parameters.
 * </p>
 * <p>Var(b) = (X<sup>T</sup>X)<sup>-1</sup>
 * </p>
 * <p>Uses QR decomposition to reduce (X<sup>T</sup>X)<sup>-1</sup>
 * to (R<sup>T</sup>R)<sup>-1</sup>, with only the top p rows of
 * R included, where p = the length of the beta vector.</p>
 *
 * <p>Data for the model must have been successfully loaded using one of
 * the {@code newSampleData} methods before invoking this method; otherwise
 * a {@code NullPointerException} will be thrown.</p>
 *
 * @return The beta variance-covariance matrix
 * @throws org.apache.commons.math3.linear.SingularMatrixException if the design matrix is singular
 * @throws NullPointerException if the data for the model have not been loaded
 */
@Override
protected RealMatrix calculateBetaVariance() {
    int p = getX().getColumnDimension();
    RealMatrix Raug = qr.getR().getSubMatrix(0, p - 1 , 0, p - 1);
    RealMatrix Rinv = new LUDecomposition(Raug).getSolver().getInverse();
    return Rinv.multiply(Rinv.transpose());
}
 

/**
 * Applies tangent normalization.
 * <p>
 *     The input row order should match the panel's target order.
 * </p>
 *
 * @param normals the log-normalized or reduced-panel counts from a panel of normals
 * @param input the input counts to normalize. This matrix is TxS where T is the number of targets
 *              and S the number of count columns.
 * @param betaHats the beta-hats for the projection to use for the normalization. This matrix
 *                  is NxS where N is the number of samples in the panel of choice and S is the number of count columns.
 * @return never {@code null}.
 */
private static RealMatrix tangentNormalize(final RealMatrix normals,
                                           final RealMatrix input,
                                           final RealMatrix betaHats) {
    Utils.validateArg(input.getColumnDimension() == betaHats.getColumnDimension(),
            String.format("the input count column count (%d) does not match the number of columns in the beta-hats (%d)",
                    input.getColumnDimension(), betaHats.getColumnDimension()));
    Utils.validateArg(normals.getColumnDimension() == betaHats.getRowDimension(),
            String.format("beta-hats component count (%d) does not match the number of samples in the PoN (%d)",
                    normals.getRowDimension(), normals.getColumnDimension()));
    final RealMatrix projection = normals.multiply(betaHats);
    return input.subtract(projection);
}
 

/**
 * <p>Calculates the variance-covariance matrix of the regression parameters.
 * </p>
 * <p>Var(b) = (X<sup>T</sup>X)<sup>-1</sup>
 * </p>
 * <p>Uses QR decomposition to reduce (X<sup>T</sup>X)<sup>-1</sup>
 * to (R<sup>T</sup>R)<sup>-1</sup>, with only the top p rows of
 * R included, where p = the length of the beta vector.</p>
 *
 * <p>Data for the model must have been successfully loaded using one of
 * the {@code newSampleData} methods before invoking this method; otherwise
 * a {@code NullPointerException} will be thrown.</p>
 *
 * @return The beta variance-covariance matrix
 */
@Override
protected RealMatrix calculateBetaVariance() {
    int p = getX().getColumnDimension();
    RealMatrix Raug = qr.getR().getSubMatrix(0, p - 1 , 0, p - 1);
    RealMatrix Rinv = new LUDecomposition(Raug).getSolver().getInverse();
    return Rinv.multiply(Rinv.transpose());
}