Java Code Examples for org.apache.commons.math3.linear.MatrixUtils#createRealMatrix()

public static RealMatrix averageHankelMatrix(RealMatrix hankelMat, int windowSize) {

        int k = hankelMat.getRowDimension();
        int m = hankelMat.getColumnDimension() / windowSize;
        int n = k + windowSize - 1;

        RealMatrix result = MatrixUtils.createRealMatrix(n, m);

        for (int t = 0; t < n; ++t) {
            int i = (t < windowSize) ? 0 : (t - windowSize + 1);
            int j = (t < windowSize) ? t : (windowSize - 1);
            int counter = 0;

            for (; i < k && j >= 0; ++i, --j, ++counter) {
                for (int h = 0; h < m; ++h) {
                    result.addToEntry(t, h, hankelMat.getEntry(i, j * m + h));
                }
            }

            for (int h = 0; h < m; ++h) {
                result.setEntry(t, h, result.getEntry(t, h) / counter);
            }
        }

        return result;
    }
 

public RidgeRegression(double[][] x, double[] y) {
	this.X = MatrixUtils.createRealMatrix(x);
	this.X_svd = null;
	this.Y = y;
	this.l2penalty = 0;
	this.coefficients = null;
	
	this.fitted = new double[y.length];
	this.residuals = new double[y.length];
}
 

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;
			}
			LOOPredLL[j] = sum;
		}

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

		return sumLOOPredLL;
	}
 

/**
 * Calculates the parameters of a bivariate quadratic equation.
 * 
 * @param elevationValues the window of points to use.
 * @return the parameters of the bivariate quadratic equation as [a, b, c, d, e, f]
 */
private static double[] calculateParameters( final double[][] elevationValues ) {
    int rows = elevationValues.length;
    int cols = elevationValues[0].length;
    int pointsNum = rows * cols;

    final double[][] xyMatrix = new double[pointsNum][6];
    final double[] valueArray = new double[pointsNum];

    // TODO check on resolution
    int index = 0;
    for( int y = 0; y < rows; y++ ) {
        for( int x = 0; x < cols; x++ ) {
            xyMatrix[index][0] = x * x; // x^2
            xyMatrix[index][1] = y * y; // y^2
            xyMatrix[index][2] = x * y; // xy
            xyMatrix[index][3] = x; // x
            xyMatrix[index][4] = y; // y
            xyMatrix[index][5] = 1;
            valueArray[index] = elevationValues[y][x];
            index++;
        }
    }

    RealMatrix A = MatrixUtils.createRealMatrix(xyMatrix);
    RealVector z = MatrixUtils.createRealVector(valueArray);

    DecompositionSolver solver = new RRQRDecomposition(A).getSolver();
    RealVector solution = solver.solve(z);

    // start values for a, b, c, d, e, f, all set to 0.0
    final double[] parameters = solution.toArray();
    return parameters;
}
 

@Test
public void testMatricesValues() {
   BiDiagonalTransformer transformer =
        new BiDiagonalTransformer(MatrixUtils.createRealMatrix(testSquare));
   final double s17 = FastMath.sqrt(17.0);
    RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
            {  -8 / (5 * s17), 19 / (5 * s17) },
            { -19 / (5 * s17), -8 / (5 * s17) }
    });
    RealMatrix bRef = MatrixUtils.createRealMatrix(new double[][] {
            { -3 * s17 / 5, 32 * s17 / 85 },
            {      0.0,     -5 * s17 / 17 }
    });
    RealMatrix vRef = MatrixUtils.createRealMatrix(new double[][] {
            { 1.0,  0.0 },
            { 0.0, -1.0 }
    });

    // check values against known references
    RealMatrix u = transformer.getU();
    Assert.assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-14);
    RealMatrix b = transformer.getB();
    Assert.assertEquals(0, b.subtract(bRef).getNorm(), 1.0e-14);
    RealMatrix v = transformer.getV();
    Assert.assertEquals(0, v.subtract(vRef).getNorm(), 1.0e-14);

    // check the same cached instance is returned the second time
    Assert.assertTrue(u == transformer.getU());
    Assert.assertTrue(b == transformer.getB());
    Assert.assertTrue(v == transformer.getV());

}
 

@Test
public void f() {
    double[][] mat = { {1, 2}, {3, 2}, {2, 5}, {7, 8}, {3, 4}, {8, 9}, {3, 3}};
    RealMatrix data = MatrixUtils.createRealMatrix(mat);

    RealMatrix res = SpectralMethods.createHankelMatrix(data, 3);
    for (int i = 0; i < res.getRowDimension(); ++i) {
        System.out.println(Arrays.toString(res.getRow(i)));
    }

    RealMatrix data2 = SpectralMethods.averageHankelMatrix(res, 3);
    for (int i = 0; i < data2.getRowDimension(); ++i) {
        System.out.println(Arrays.toString(data2.getRow(i)));
    }
}
 

final public RealMatrix getMatrixVarByName(String name) {
    readLock.lock();
    try {
        int row = getVectorVarSizeByName(name);
        int column = getVectorVarDimensionByName(name);
        RealMatrix matrix = MatrixUtils.createRealMatrix(row, column);
        for (int i=0; i<row; i++) {
            matrix.setRowVector(i, vectorVars.get(name).get(i));
        }
        return matrix;
    } finally {
        readLock.unlock();
    }
}
 

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);
}
 

/**
 * Compute the normal matrix, J<sup>T</sup>J.
 *
 * @param jacobian  the m by n jacobian matrix, J. Input.
 * @param residuals the m by 1 residual vector, r. Input.
 * @return  the n by n normal matrix and  the n by 1 J<sup>Tr vector.
 */
private static Pair<RealMatrix, RealVector> computeNormalMatrix(final RealMatrix jacobian,
                                                                final RealVector residuals) {
    //since the normal matrix is symmetric, we only need to compute half of it.
    final int nR = jacobian.getRowDimension();
    final int nC = jacobian.getColumnDimension();
    //allocate space for return values
    final RealMatrix normal = MatrixUtils.createRealMatrix(nC, nC);
    final RealVector jTr = new ArrayRealVector(nC);
    //for each measurement
    for (int i = 0; i < nR; ++i) {
        //compute JTr for measurement i
        for (int j = 0; j < nC; j++) {
            jTr.setEntry(j, jTr.getEntry(j) +
                    residuals.getEntry(i) * jacobian.getEntry(i, j));
        }

        // add the the contribution to the normal matrix for measurement i
        for (int k = 0; k < nC; ++k) {
            //only compute the upper triangular part
            for (int l = k; l < nC; ++l) {
                normal.setEntry(k, l, normal.getEntry(k, l) +
                        jacobian.getEntry(i, k) * jacobian.getEntry(i, l));
            }
        }
    }
    //copy the upper triangular part to the lower triangular part.
    for (int i = 0; i < nC; i++) {
        for (int j = 0; j < i; j++) {
            normal.setEntry(i, j, normal.getEntry(j, i));
        }
    }
    return new Pair<RealMatrix, RealVector>(normal, jTr);
}
 

/**
 * Converts an Earth-Centered Earth-Fixed (ECEF) coordinate to ENU. 
 * 
 * @param cEcef the ECEF coordinate.
 * @return the ENU coordinate.
 * @throws MatrixException 
 */
public Coordinate ecefToEnu( Coordinate cEcef ) {
    double deltaX = cEcef.x - _ecefROriginX;
    double deltaY = cEcef.y - _ecefROriginY;
    double deltaZ = cEcef.z - _ecefROriginZ;

    double[][] deltas = new double[][]{{deltaX}, {deltaY}, {deltaZ}};
    RealMatrix deltaMatrix = MatrixUtils.createRealMatrix(deltas);

    RealMatrix enuMatrix = _rotationMatrix.multiply(deltaMatrix);
    double[] column = enuMatrix.getColumn(0);

    return new Coordinate(column[0], column[1], column[2]);
}
 

public static Map<Integer, double[]> spectralEmbedding(GraphReadMethods rg, final Set<Integer> includedVertices) {

        Map<Integer, double[]> vertexPositions = new HashMap<>();

        // Don't position anything if there are fewer than 3 vertices to embedd - this embedding shouldn't be used in these cases.
        if (includedVertices.size() <= 2) {
            return vertexPositions;
        }

        GraphMatrix l = GraphMatrix.adjacencyFromGraph(rg, includedVertices, new HashSet<>());

        final EigenDecomposition e = new EigenDecomposition(MatrixUtils.createRealMatrix(l.laplacianMatrix));
        final int numVectors = e.getRealEigenvalues().length;

        for (int i = 0; i < l.dimension; i++) {
            double xPos = 0;
            double yPos = 0;
            double norm = 0;
            for (int j = 0; j < numVectors - 1; j++) {
                xPos += e.getRealEigenvalue(j) * e.getEigenvector(j).getEntry(i);
                yPos += e.getRealEigenvalue(j) * e.getEigenvector(j).getEntry(i) * (j % 2 == 0 ? -1 : 1);
                norm += Math.abs(e.getRealEigenvalue(j)) * (Math.pow(e.getEigenvector(j).getEntry(i), 2));
            }
            norm = Math.sqrt(norm);
            vertexPositions.put(l.matrixPositionToID.get(i), new double[]{norm * xPos, norm * yPos});
        }

        return vertexPositions;

    }
 

private void initMatrices() {
	this.L = MatrixUtils.createRealMatrix(this.X.getRowDimension(), this.X.getColumnDimension());
	this.S = MatrixUtils.createRealMatrix(this.X.getRowDimension(), this.X.getColumnDimension());
	this.E = MatrixUtils.createRealMatrix(this.X.getRowDimension(), this.X.getColumnDimension());
}
 

/**
 * Generate an error covariance matrix and sample data representing models
 * with this error structure. Then verify that GLS estimated coefficients,
 * on average, perform better than OLS.
 */
@Test
public void testGLSEfficiency() {
    RandomGenerator rg = new JDKRandomGenerator();
    rg.setSeed(200);  // Seed has been selected to generate non-trivial covariance
    
    // Assume model has 16 observations (will use Longley data).  Start by generating
    // non-constant variances for the 16 error terms.
    final int nObs = 16;
    double[] sigma = new double[nObs];
    for (int i = 0; i < nObs; i++) {
        sigma[i] = 10 * rg.nextDouble();
    }
    
    // Now generate 1000 error vectors to use to estimate the covariance matrix
    // Columns are draws on N(0, sigma[col])
    final int numSeeds = 1000;
    RealMatrix errorSeeds = MatrixUtils.createRealMatrix(numSeeds, nObs);
    for (int i = 0; i < numSeeds; i++) {
        for (int j = 0; j < nObs; j++) {
            errorSeeds.setEntry(i, j, rg.nextGaussian() * sigma[j]);
        }
    }
    
    // Get covariance matrix for columns
    RealMatrix cov = (new Covariance(errorSeeds)).getCovarianceMatrix();
      
    // Create a CorrelatedRandomVectorGenerator to use to generate correlated errors
    GaussianRandomGenerator rawGenerator = new GaussianRandomGenerator(rg);
    double[] errorMeans = new double[nObs];  // Counting on init to 0 here
    CorrelatedRandomVectorGenerator gen = new CorrelatedRandomVectorGenerator(errorMeans, cov,
     1.0e-12 * cov.getNorm(), rawGenerator);
    
    // Now start generating models.  Use Longley X matrix on LHS
    // and Longley OLS beta vector as "true" beta.  Generate
    // Y values by XB + u where u is a CorrelatedRandomVector generated
    // from cov.
    OLSMultipleLinearRegression ols = new OLSMultipleLinearRegression();
    ols.newSampleData(longley, nObs, 6);
    final RealVector b = ols.calculateBeta().copy();
    final RealMatrix x = ols.getX().copy();
    
    // Create a GLS model to reuse
    GLSMultipleLinearRegression gls = new GLSMultipleLinearRegression();
    gls.newSampleData(longley, nObs, 6);
    gls.newCovarianceData(cov.getData());
    
    // Create aggregators for stats measuring model performance
    DescriptiveStatistics olsBetaStats = new DescriptiveStatistics();
    DescriptiveStatistics glsBetaStats = new DescriptiveStatistics();
    
    // Generate Y vectors for 10000 models, estimate GLS and OLS and
    // Verify that OLS estimates are better
    final int nModels = 10000;
    for (int i = 0; i < nModels; i++) {
        
        // Generate y = xb + u with u cov
        RealVector u = MatrixUtils.createRealVector(gen.nextVector());
        double[] y = u.add(x.operate(b)).toArray();
        
        // Estimate OLS parameters
        ols.newYSampleData(y);
        RealVector olsBeta = ols.calculateBeta();
        
        // Estimate GLS parameters
        gls.newYSampleData(y);
        RealVector glsBeta = gls.calculateBeta();
        
        // Record deviations from "true" beta
        double dist = olsBeta.getDistance(b);
        olsBetaStats.addValue(dist * dist);
        dist = glsBeta.getDistance(b);
        glsBetaStats.addValue(dist * dist);
        
    }
    
    // Verify that GLS is on average more efficient, lower variance
    assert(olsBetaStats.getMean() > 1.5 * glsBetaStats.getMean());
    assert(olsBetaStats.getStandardDeviation() > glsBetaStats.getStandardDeviation());  
}
 

private void computeADFStatistics() {
	double[] y = diff(ts);
	RealMatrix designMatrix = null;
	int k = lag+1;
	int n = ts.length - 1;
	
	RealMatrix z = MatrixUtils.createRealMatrix(laggedMatrix(y, k)); //has rows length(ts) - 1 - k + 1
	RealVector zcol1 = z.getColumnVector(0); //has length length(ts) - 1 - k + 1
	double[] xt1 = subsetArray(ts, k-1, n-1);  //ts[k:(length(ts) - 1)], has length length(ts) - 1 - k + 1
	double[] trend = sequence(k,n); //trend k:n, has length length(ts) - 1 - k + 1
	if (k > 1) {
		RealMatrix yt1 = z.getSubMatrix(0, ts.length - 1 - k, 1, k-1); //same as z but skips first column
		//build design matrix as cbind(xt1, 1, trend, yt1)
		designMatrix = MatrixUtils.createRealMatrix(ts.length - 1 - k + 1, 3 + k - 1);
		designMatrix.setColumn(0, xt1);
		designMatrix.setColumn(1, ones(ts.length - 1 - k + 1));
		designMatrix.setColumn(2, trend);
		designMatrix.setSubMatrix(yt1.getData(), 0, 3);
		
	} else {
		//build design matrix as cbind(xt1, 1, tt)
		designMatrix = MatrixUtils.createRealMatrix(ts.length - 1 - k + 1, 3);
		designMatrix.setColumn(0, xt1);
		designMatrix.setColumn(1, ones(ts.length - 1 - k + 1));
		designMatrix.setColumn(2, trend);
	}
	/*OLSMultipleLinearRegression regression = new OLSMultipleLinearRegression();
	regression.setNoIntercept(true);
	regression.newSampleData(zcol1.toArray(), designMatrix.getData());
	double[] beta = regression.estimateRegressionParameters();
	double[] sd = regression.estimateRegressionParametersStandardErrors();
	*/
	RidgeRegression regression = new RidgeRegression(designMatrix.getData(), zcol1.toArray());
	regression.updateCoefficients(.0001);
	double[] beta = regression.getCoefficients();
	double[] sd = regression.getStandarderrors();
	
	double t = beta[0] / sd[0];
	if (t <= PVALUE_THRESHOLD) {
		this.needsDiff = true;
	} else {
		this.needsDiff = false;
	}	
}
 

/**
 * Generate an error covariance matrix and sample data representing models
 * with this error structure. Then verify that GLS estimated coefficients,
 * on average, perform better than OLS.
 */
@Test
public void testGLSEfficiency() {
    RandomGenerator rg = new JDKRandomGenerator();
    rg.setSeed(200);  // Seed has been selected to generate non-trivial covariance
    
    // Assume model has 16 observations (will use Longley data).  Start by generating
    // non-constant variances for the 16 error terms.
    final int nObs = 16;
    double[] sigma = new double[nObs];
    for (int i = 0; i < nObs; i++) {
        sigma[i] = 10 * rg.nextDouble();
    }
    
    // Now generate 1000 error vectors to use to estimate the covariance matrix
    // Columns are draws on N(0, sigma[col])
    final int numSeeds = 1000;
    RealMatrix errorSeeds = MatrixUtils.createRealMatrix(numSeeds, nObs);
    for (int i = 0; i < numSeeds; i++) {
        for (int j = 0; j < nObs; j++) {
            errorSeeds.setEntry(i, j, rg.nextGaussian() * sigma[j]);
        }
    }
    
    // Get covariance matrix for columns
    RealMatrix cov = (new Covariance(errorSeeds)).getCovarianceMatrix();
      
    // Create a CorrelatedRandomVectorGenerator to use to generate correlated errors
    GaussianRandomGenerator rawGenerator = new GaussianRandomGenerator(rg);
    double[] errorMeans = new double[nObs];  // Counting on init to 0 here
    CorrelatedRandomVectorGenerator gen = new CorrelatedRandomVectorGenerator(errorMeans, cov,
     1.0e-12 * cov.getNorm(), rawGenerator);
    
    // Now start generating models.  Use Longley X matrix on LHS
    // and Longley OLS beta vector as "true" beta.  Generate
    // Y values by XB + u where u is a CorrelatedRandomVector generated
    // from cov.
    OLSMultipleLinearRegression ols = new OLSMultipleLinearRegression();
    ols.newSampleData(longley, nObs, 6);
    final RealVector b = ols.calculateBeta().copy();
    final RealMatrix x = ols.getX().copy();
    
    // Create a GLS model to reuse
    GLSMultipleLinearRegression gls = new GLSMultipleLinearRegression();
    gls.newSampleData(longley, nObs, 6);
    gls.newCovarianceData(cov.getData());
    
    // Create aggregators for stats measuring model performance
    DescriptiveStatistics olsBetaStats = new DescriptiveStatistics();
    DescriptiveStatistics glsBetaStats = new DescriptiveStatistics();
    
    // Generate Y vectors for 10000 models, estimate GLS and OLS and
    // Verify that OLS estimates are better
    final int nModels = 10000;
    for (int i = 0; i < nModels; i++) {
        
        // Generate y = xb + u with u cov
        RealVector u = MatrixUtils.createRealVector(gen.nextVector());
        double[] y = u.add(x.operate(b)).toArray();
        
        // Estimate OLS parameters
        ols.newYSampleData(y);
        RealVector olsBeta = ols.calculateBeta();
        
        // Estimate GLS parameters
        gls.newYSampleData(y);
        RealVector glsBeta = gls.calculateBeta();
        
        // Record deviations from "true" beta
        double dist = olsBeta.getDistance(b);
        olsBetaStats.addValue(dist * dist);
        dist = glsBeta.getDistance(b);
        glsBetaStats.addValue(dist * dist);
        
    }
    
    // Verify that GLS is on average more efficient, lower variance
    assert(olsBetaStats.getMean() > 1.5 * glsBetaStats.getMean());
    assert(olsBetaStats.getStandardDeviation() > glsBetaStats.getStandardDeviation());  
}
 

private double[] sumOverSegments(final JavaRDD<double[]> eZskBySeg) {
    final double[][] eZskBySeg2D = eZskBySeg.collect().stream().toArray(double[][]::new); // S x K
    final RealMatrix eZskBySegMatrix = MatrixUtils.createRealMatrix(eZskBySeg2D);
    return GATKProtectedMathUtils.columnSums(eZskBySegMatrix);
}
 

@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);

            }
        }
    }
 }
 

@Test
public void testSingularMatrix() {
   BiDiagonalTransformer transformer =
        new BiDiagonalTransformer(MatrixUtils.createRealMatrix(new double[][] {
            { 1.0, 2.0, 3.0 },
            { 2.0, 3.0, 4.0 },
            { 3.0, 5.0, 7.0 }
        }));
   final double s3  = FastMath.sqrt(3.0);
   final double s14 = FastMath.sqrt(14.0);
   final double s1553 = FastMath.sqrt(1553.0);
   RealMatrix uRef = MatrixUtils.createRealMatrix(new double[][] {
       {  -1.0 / s14,  5.0 / (s3 * s14),  1.0 / s3 },
       {  -2.0 / s14, -4.0 / (s3 * s14),  1.0 / s3 },
       {  -3.0 / s14,  1.0 / (s3 * s14), -1.0 / s3 }
   });
   RealMatrix bRef = MatrixUtils.createRealMatrix(new double[][] {
       { -s14, s1553 / s14,   0.0 },
       {  0.0, -87 * s3 / (s14 * s1553), -s3 * s14 / s1553 },
       {  0.0, 0.0, 0.0 }
   });
   RealMatrix vRef = MatrixUtils.createRealMatrix(new double[][] {
       { 1.0,   0.0,         0.0        },
       { 0.0,  -23 / s1553,  32 / s1553 },
       { 0.0,  -32 / s1553, -23 / s1553 }
   });

   // check values against known references
   RealMatrix u = transformer.getU();
   Assert.assertEquals(0, u.subtract(uRef).getNorm(), 1.0e-14);
   RealMatrix b = transformer.getB();
   Assert.assertEquals(0, b.subtract(bRef).getNorm(), 1.0e-14);
   RealMatrix v = transformer.getV();
   Assert.assertEquals(0, v.subtract(vRef).getNorm(), 1.0e-14);

   // check the same cached instance is returned the second time
   Assert.assertTrue(u == transformer.getU());
   Assert.assertTrue(b == transformer.getB());
   Assert.assertTrue(v == transformer.getV());

}
 

/**
 * {@inheritDoc}
 * @throws NumberIsTooSmallException if the number of observations
 * in a cell is < 2
 */
@Override
public RealMatrix getCovarianceMatrix() throws NumberIsTooSmallException {
    return MatrixUtils.createRealMatrix(getData());
}
 

/**
 * {@inheritDoc}
 * @throws NumberIsTooSmallException if the number of observations
 * in a cell is < 2
 */
@Override
public RealMatrix getCovarianceMatrix() throws NumberIsTooSmallException {
    return MatrixUtils.createRealMatrix(getData());
}