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

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

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

Introduction

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

Prototype

RealMatrix transpose();

Source Link

Document

Returns the transpose of this matrix.

Usage

From source file:com.joptimizer.algebra.CholeskyComparisonTest.java

/**
 * Test with a diagonal matrix./*www.j  a v a2s.c  om*/
 */
public void testCompareFactorizations2() throws Exception {
    log.debug("testCompareFactorizations2");

    int iterations = 4;
    int m = 1000;
    int dim = m * m;

    DoubleMatrix2D QMatrix = DoubleFactory2D.sparse.diagonal(DoubleFactory1D.sparse.make(m, 1.));
    double[][] QMatrixData = QMatrix.toArray();
    //log.debug("QMatrix: " + Utils.toString(QMatrix.toArray()));
    RealMatrix A = MatrixUtils.createRealMatrix(QMatrix.toArray());
    log.debug("cardinality: " + QMatrix.cardinality());
    int nz = dim - QMatrix.cardinality();
    log.debug("sparsity index: " + 100 * new Double(nz) / dim + " %");

    //try Cholesky 1
    long t1F = System.currentTimeMillis();
    CholeskyFactorization myc1 = null;
    for (int i = 0; i < iterations; i++) {
        //factorization
        myc1 = new CholeskyFactorization(QMatrix);
        myc1.factorize();
    }
    log.debug("Cholesky standard factorization time: " + (System.currentTimeMillis() - t1F));
    RealMatrix L1 = MatrixUtils.createRealMatrix(myc1.getL().toArray());
    double norm1 = A.subtract(L1.multiply(L1.transpose())).getNorm();
    log.debug("norm1                               : " + norm1);
    assertEquals(0., norm1, 1.E-12 * Math.sqrt(dim));
    long t1I = System.currentTimeMillis();
    for (int i = 0; i < iterations; i++) {
        //inversion
        myc1.getInverse();
    }
    log.debug("Cholesky standard inversion time    : " + (System.currentTimeMillis() - t1I));
    RealMatrix AInv1 = MatrixUtils.createRealMatrix(myc1.getInverse().toArray());

    //try Cholesky 5
    long t5 = System.currentTimeMillis();
    CholeskySparseFactorization myc5 = null;
    for (int i = 0; i < iterations; i++) {
        myc5 = new CholeskySparseFactorization(new SparseDoubleMatrix2D(QMatrix.toArray()));
        myc5.factorize();
    }
    log.debug("Cholesky sparse factorization time  : " + (System.currentTimeMillis() - t5));
    RealMatrix L5 = MatrixUtils.createRealMatrix(myc5.getL().toArray());
    double norm5 = A.subtract(L5.multiply(L5.transpose())).getNorm();
    log.debug("norm5                               : " + norm5);
    assertEquals(0., norm5, 1.E-12 * Math.sqrt(dim));
    assertEquals(0., L1.subtract(L5).getNorm(), 1.E-10 * Math.sqrt(dim));
}

From source file:com.joptimizer.optimizers.OptimizationRequestHandler.java

/**
 * Find a solution of the linear (equalities) system A.x = b.
 * @see "S.Boyd and L.Vandenberghe, Convex Optimization, p. 682"
 *//*  w w w  . j a  v  a  2 s .c  om*/
protected double[] findEqFeasiblePoint(double[][] A, double[] b) throws Exception {
    RealMatrix AT = new Array2DRowRealMatrix(A).transpose();
    int p = A.length;

    SingularValueDecomposition dFact1 = new SingularValueDecomposition(AT);
    int rangoAT = dFact1.getRank();
    if (rangoAT != p) {
        throw new RuntimeException("Equalities matrix A must have full rank");
    }

    QRDecomposition dFact = new QRDecomposition(AT);
    //A = QR
    RealMatrix Q1Q2 = dFact.getQ();
    RealMatrix R0 = dFact.getR();
    RealMatrix Q1 = Q1Q2.getSubMatrix(0, AT.getRowDimension() - 1, 0, p - 1);
    RealMatrix R = R0.getSubMatrix(0, p - 1, 0, p - 1);
    double[][] rData = R.copy().getData();
    //inversion
    double[][] rInvData = Utils.upperTriangularMatrixUnverse(rData);
    RealMatrix RInv = new Array2DRowRealMatrix(rInvData);

    //w = Q1 *   Inv([R]T) . b
    double[] w = Q1.operate(RInv.transpose().operate(new ArrayRealVector(b))).toArray();
    return w;
}

From source file:eagle.security.userprofile.model.eigen.UserProfileEigenModeler.java

private RealMatrix computeMaxDistanceOnPCs(int index) {

    RealMatrix pc = new Array2DRowRealMatrix(principalComponents[index].toArray());
    RealMatrix transposePC1 = pc.transpose();
    RealMatrix finalDataTranspose = finalMatrixWithoutLowVariantCmds.transpose();

    RealMatrix trainingData = pc.multiply(transposePC1).multiply(finalDataTranspose);
    RealMatrix trainingDataTranspose = trainingData.transpose();

    double maxDistance = 0.0;
    double minDistance = 0.0;
    RealVector p1 = null, p2 = null;//from  www  . j  av a2  s. c o m
    if (LOG.isDebugEnabled())
        LOG.debug("Training data transpose size: " + trainingDataTranspose.getRowDimension() + "x"
                + trainingDataTranspose.getColumnDimension());
    for (int i = 0; i < trainingDataTranspose.getRowDimension(); i++) {
        RealVector iRowVector = new ArrayRealVector(trainingDataTranspose.getRow(i));
        RealVector transposePC1Vect = transposePC1.getRowVector(0);
        double distance = iRowVector.getDistance(transposePC1Vect);
        if (distance > maxDistance) {
            maxDistance = distance;
            p1 = iRowVector;
            p2 = transposePC1Vect;
        }
        if (distance < minDistance)
            minDistance = distance;
        //}
    }
    maximumL2Norm.setEntry(index, maxDistance);
    minimumL2Norm.setEntry(index, minDistance);
    minVector = p1;
    maxVector = p2;
    return trainingDataTranspose;
}

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

@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);// ww  w . j av a 2 s . com
    }

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

From source file:com.joptimizer.algebra.CholeskyFactorizationTest.java

public void testInvert1() throws Exception {
    log.debug("testInvert1");
    double[][] QData = new double[][] { { 1, .12, .13, .14, .15 }, { .12, 2, .23, .24, .25 },
            { .13, .23, 3, 0, 0 }, { .14, .24, 0, 4, 0 }, { .15, .25, 0, 0, 5 } };
    RealMatrix Q = MatrixUtils.createRealMatrix(QData);

    CholeskyFactorization myc = new CholeskyFactorization(DoubleFactory2D.dense.make(QData));
    myc.factorize();//  w  ww .j av  a2 s  .  c  o m
    RealMatrix L = new Array2DRowRealMatrix(myc.getL().toArray());
    RealMatrix LT = new Array2DRowRealMatrix(myc.getLT().toArray());
    log.debug("L: " + ArrayUtils.toString(L.getData()));
    log.debug("LT: " + ArrayUtils.toString(LT.getData()));
    log.debug("L.LT: " + ArrayUtils.toString(L.multiply(LT).getData()));
    log.debug("LT.L: " + ArrayUtils.toString(LT.multiply(L).getData()));

    // check Q = L.LT
    double norm = L.multiply(LT).subtract(Q).getNorm();
    log.debug("norm: " + norm);
    assertTrue(norm < 1.E-15);

    RealMatrix LInv = new SingularValueDecomposition(L).getSolver().getInverse();
    log.debug("LInv: " + ArrayUtils.toString(LInv.getData()));
    RealMatrix LInvT = LInv.transpose();
    log.debug("LInvT: " + ArrayUtils.toString(LInvT.getData()));
    RealMatrix LTInv = new SingularValueDecomposition(LT).getSolver().getInverse();
    log.debug("LTInv: " + ArrayUtils.toString(LTInv.getData()));
    RealMatrix LTInvT = LTInv.transpose();
    log.debug("LTInvT: " + ArrayUtils.toString(LTInvT.getData()));
    log.debug("LInv.LInvT: " + ArrayUtils.toString(LInv.multiply(LInvT).getData()));
    log.debug("LTInv.LTInvT: " + ArrayUtils.toString(LTInv.multiply(LTInvT).getData()));

    RealMatrix Id = MatrixUtils.createRealIdentityMatrix(Q.getRowDimension());
    //check Q.(LTInv * LInv) = 1
    norm = Q.multiply(LTInv.multiply(LInv)).subtract(Id).getNorm();
    log.debug("norm: " + norm);
    assertTrue(norm < 5.E-15);

    // check Q.QInv = 1
    RealMatrix QInv = MatrixUtils.createRealMatrix(myc.getInverse().toArray());
    norm = Q.multiply(QInv).subtract(Id).getNorm();
    log.debug("norm: " + norm);
    assertTrue(norm < 1.E-15);
}

From source file:com.joptimizer.algebra.CholeskyFactorizationTest.java

/**
 * The same as before, but with rescaling.
 *//*from  ww  w . java  2s.com*/
public void testInvert2() throws Exception {
    log.debug("testInvert2");
    double[][] QData = new double[][] { { 1, .12, .13, .14, .15 }, { .12, 2, .23, .24, .25 },
            { .13, .23, 3, 0, 0 }, { .14, .24, 0, 4, 0 }, { .15, .25, 0, 0, 5 } };
    RealMatrix Q = MatrixUtils.createRealMatrix(QData);

    CholeskyFactorization myc = new CholeskyFactorization(DoubleFactory2D.dense.make(QData),
            new Matrix1NormRescaler());
    myc.factorize();
    RealMatrix L = new Array2DRowRealMatrix(myc.getL().toArray());
    RealMatrix LT = new Array2DRowRealMatrix(myc.getLT().toArray());
    log.debug("L: " + ArrayUtils.toString(L.getData()));
    log.debug("LT: " + ArrayUtils.toString(LT.getData()));
    log.debug("L.LT: " + ArrayUtils.toString(L.multiply(LT).getData()));
    log.debug("LT.L: " + ArrayUtils.toString(LT.multiply(L).getData()));

    // check Q = L.LT
    double norm = L.multiply(LT).subtract(Q).getNorm();
    log.debug("norm: " + norm);
    assertTrue(norm < 1.E-15);

    RealMatrix LInv = new SingularValueDecomposition(L).getSolver().getInverse();
    log.debug("LInv: " + ArrayUtils.toString(LInv.getData()));
    RealMatrix LInvT = LInv.transpose();
    log.debug("LInvT: " + ArrayUtils.toString(LInvT.getData()));
    RealMatrix LTInv = new SingularValueDecomposition(LT).getSolver().getInverse();
    log.debug("LTInv: " + ArrayUtils.toString(LTInv.getData()));
    RealMatrix LTInvT = LTInv.transpose();
    log.debug("LTInvT: " + ArrayUtils.toString(LTInvT.getData()));
    log.debug("LInv.LInvT: " + ArrayUtils.toString(LInv.multiply(LInvT).getData()));
    log.debug("LTInv.LTInvT: " + ArrayUtils.toString(LTInv.multiply(LTInvT).getData()));

    RealMatrix Id = MatrixUtils.createRealIdentityMatrix(Q.getRowDimension());
    //check Q.(LTInv * LInv) = 1
    norm = Q.multiply(LTInv.multiply(LInv)).subtract(Id).getNorm();
    log.debug("norm: " + norm);
    assertTrue(norm < 5.E-15);

    // check Q.QInv = 1
    RealMatrix QInv = MatrixUtils.createRealMatrix(myc.getInverse().toArray());
    norm = Q.multiply(QInv).subtract(Id).getNorm();
    log.debug("norm: " + norm);
    assertTrue(norm < 1.E-15);
}

From source file:edu.cudenver.bios.power.glmm.GLMMTest.java

/**
 * Create a statistical test for to compute power
 * @param params GLMM input parameters/*from  w  w w .j a  v a  2 s .com*/
 */
public GLMMTest(FApproximation fMethod, RealMatrix Xessence, RealMatrix XtXInverse, int perGroupN, int rank,
        FixedRandomMatrix C, RealMatrix U, RealMatrix thetaNull, RealMatrix beta, RealMatrix sigmaError) {
    this.fMethod = fMethod;
    this.Xessence = Xessence;
    this.XtXInverse = XtXInverse;
    int xrd = Xessence.getRowDimension();
    if (xrd != 0 && perGroupN > Integer.MAX_VALUE / xrd) {
        throw new GLMMTestException(OVERFLOW_ERROR_MESSAGE);
    }
    this.totalN = (double) xrd * perGroupN;
    this.rank = rank;
    this.CFixed = C.getFixedMatrix();
    this.C = C.getCombinedMatrix();
    this.U = U;
    this.thetaNull = thetaNull;
    this.beta = beta;
    this.sigmaError = sigmaError;

    // check if uni/multivariate design
    multivariate = (beta.getColumnDimension() > 1);

    // cache the value of M
    RealMatrix CFixed = C.getFixedMatrix();
    RealMatrix cxxcEssence = CFixed.multiply((XtXInverse).multiply(CFixed.transpose()));
    RealMatrix cxxcEssenceInverse = new LUDecomposition(cxxcEssence).getSolver().getInverse();
    this.M = cxxcEssenceInverse.scalarMultiply(perGroupN);
}

From source file:com.joptimizer.algebra.CholeskyComparisonTest.java

/**
 * Test with a generic sparse matrix./*from  w ww.j a  v a 2  s  .c  om*/
 */
public void testCompareFactorizations1() throws Exception {
    log.debug("testCompareFactorizations1");

    int iterations = 4;
    int m = 1000;
    int dim = m * m;

    DoubleMatrix2D sMatrix = Utils.randomValuesSparseMatrix(m, m, -10, 10, 0.97, 12345L);
    //log.debug("sMatrix: " + Utils.toString(sMatrix.toArray()));
    DoubleMatrix2D QMatrix = Algebra.DEFAULT.mult(sMatrix, Algebra.DEFAULT.transpose(sMatrix));//positive and symmetric
    double[][] QMatrixData = QMatrix.toArray();
    //log.debug("QMatrix: " + Utils.toString(QMatrix.toArray()));
    RealMatrix A = MatrixUtils.createRealMatrix(QMatrix.toArray());
    log.debug("cardinality: " + QMatrix.cardinality());
    int nz = dim - QMatrix.cardinality();
    log.debug("sparsity index: " + 100 * new Double(nz) / dim + " %");

    //try Cholesky 1
    long t1F = System.currentTimeMillis();
    CholeskyFactorization myc1 = null;
    for (int i = 0; i < iterations; i++) {
        //factorization
        myc1 = new CholeskyFactorization(QMatrix);
        myc1.factorize();
    }
    log.debug("Cholesky standard factorization time: " + (System.currentTimeMillis() - t1F));
    RealMatrix L1 = MatrixUtils.createRealMatrix(myc1.getL().toArray());
    double norm1 = A.subtract(L1.multiply(L1.transpose())).getNorm();
    log.debug("norm1                               : " + norm1);
    assertEquals(0., norm1, 1.E-12 * Math.sqrt(dim));
    long t1I = System.currentTimeMillis();
    for (int i = 0; i < iterations; i++) {
        //inversion
        myc1.getInverse();
    }
    log.debug("Cholesky standard inversion time    : " + (System.currentTimeMillis() - t1I));
    RealMatrix AInv1 = MatrixUtils.createRealMatrix(myc1.getInverse().toArray());

    //try Cholesky 2
    long t2F = System.currentTimeMillis();
    CholeskyRCTFactorization myc2 = null;
    for (int i = 0; i < iterations; i++) {
        //factorization
        myc2 = new CholeskyRCTFactorization(QMatrix);
        myc2.factorize();
    }
    log.debug("Cholesky RCT factorization time     : " + (System.currentTimeMillis() - t2F));
    RealMatrix L2 = MatrixUtils.createRealMatrix(myc2.getL().toArray());
    double norm2 = A.subtract(L2.multiply(L2.transpose())).getNorm();
    log.debug("norm2                               : " + norm2);
    assertEquals(0., norm2, 1.E-12 * Math.sqrt(dim));
    assertEquals(0., L1.subtract(L2).getNorm(), 1.E-10 * Math.sqrt(dim));
    long t2I = System.currentTimeMillis();
    for (int i = 0; i < iterations; i++) {
        //inversion
        myc2.getInverse();
    }
    log.debug("Cholesky RCT inversion time         : " + (System.currentTimeMillis() - t2I));
    RealMatrix AInv2 = MatrixUtils.createRealMatrix(myc2.getInverse().toArray());
    assertEquals(0., AInv1.subtract(AInv2).getNorm(), 1.E-10 * Math.sqrt(dim));

    //try Cholesky 3
    long t3F = System.currentTimeMillis();
    CholeskyRCFactorization myc3 = null;
    for (int i = 0; i < iterations; i++) {
        //factorization
        myc3 = new CholeskyRCFactorization(QMatrix);
        myc3.factorize();
    }
    log.debug("Cholesky RC factorization time      : " + (System.currentTimeMillis() - t3F));
    RealMatrix L3 = MatrixUtils.createRealMatrix(myc3.getL().toArray());
    double norm3 = A.subtract(L3.multiply(L3.transpose())).getNorm();
    log.debug("norm3                               : " + norm3);
    assertEquals(0., norm3, 1.E-12 * Math.sqrt(dim));
    assertEquals(0., L1.subtract(L3).getNorm(), 1.E-10 * Math.sqrt(dim));
    long t3I = System.currentTimeMillis();
    for (int i = 0; i < iterations; i++) {
        //inversion
        myc3.getInverse();
    }
    log.debug("Cholesky RC inversion time          : " + (System.currentTimeMillis() - t3I));
    RealMatrix AInv3 = MatrixUtils.createRealMatrix(myc3.getInverse().toArray());
    assertEquals(0., AInv1.subtract(AInv3).getNorm(), 1.E-10 * Math.sqrt(dim));

    //try Cholesky 4
    long t4 = System.currentTimeMillis();
    LDLTFactorization myc4 = null;
    for (int i = 0; i < iterations; i++) {
        myc4 = new LDLTFactorization(new DenseDoubleMatrix2D(QMatrixData));
        myc4.factorize();
    }
    log.debug("Cholesky LDLT factorization time    : " + (System.currentTimeMillis() - t4));
    RealMatrix L4 = MatrixUtils.createRealMatrix(myc4.getL().toArray());
    RealMatrix D4 = MatrixUtils.createRealMatrix(myc4.getD().toArray());
    double norm4 = A.subtract(L4.multiply(D4.multiply(L4.transpose()))).getNorm();
    log.debug("norm4                               : " + norm4);
    assertEquals(0., norm4, 1.E-12 * Math.sqrt(dim));

    //try Cholesky 5
    long t5 = System.currentTimeMillis();
    CholeskySparseFactorization myc5 = null;
    for (int i = 0; i < iterations; i++) {
        myc5 = new CholeskySparseFactorization(new SparseDoubleMatrix2D(QMatrix.toArray()));
        myc5.factorize();
    }
    log.debug("Cholesky sparse factorization time  : " + (System.currentTimeMillis() - t5));
    RealMatrix L5 = MatrixUtils.createRealMatrix(myc5.getL().toArray());
    double norm5 = A.subtract(L5.multiply(L5.transpose())).getNorm();
    log.debug("norm5                               : " + norm5);
    assertEquals(0., norm5, 1.E-12 * Math.sqrt(dim));
    assertEquals(0., L1.subtract(L5).getNorm(), 1.E-10 * Math.sqrt(dim));
}

From source file:com.joptimizer.algebra.CholeskySparseFactorizationTest.java

public void testSimple1_norescaling() throws Exception {
    log.debug("testSimple1_norescaling");
    double[][] A = new double[][] { { 4, 0, 2, 2 }, { 0, 4, 2, -2 }, { 2, 2, 6, 0 }, { 2, -2, 0, 6 } };
    //expected L//from w  w w .  j  av  a2s.  c  o  m
    double[][] EL = new double[][] { { 2, 0, 0, 0 }, { 0, 2, 0, 0 }, { 1, 1, 2, 0 }, { 1, -1, 0, 2 } };

    SparseDoubleMatrix2D ASparse = new SparseDoubleMatrix2D(A);
    CholeskySparseFactorization cs = new CholeskySparseFactorization(ASparse);
    cs.factorize();
    DoubleMatrix2D L = cs.getL();
    DoubleMatrix2D LT = cs.getLT();
    log.debug("L : " + ArrayUtils.toString(L.toArray()));
    log.debug("LT: " + ArrayUtils.toString(LT.toArray()));

    RealMatrix ELMatrix = MatrixUtils.createRealMatrix(EL);
    RealMatrix LMatrix = MatrixUtils.createRealMatrix(L.toArray());
    RealMatrix LTMatrix = MatrixUtils.createRealMatrix(LT.toArray());
    assertEquals((ELMatrix.subtract(LMatrix).getNorm()), 0.);
    assertEquals((ELMatrix.subtract(LTMatrix.transpose()).getNorm()), 0.);

    //A.x = b
    double[] b = new double[] { 1, 2, 3, 4 };
    double[] x = cs.solve(F1.make(b)).toArray();

    //check the norm ||A.x-b||
    double norm = new Array2DRowRealMatrix(A).operate(new ArrayRealVector(x)).subtract(new ArrayRealVector(b))
            .getNorm();
    log.debug("norm: " + norm);
    assertEquals(0., norm, Utils.getDoubleMachineEpsilon());

    //check the scaled residual
    double residual = Utils.calculateScaledResidual(A, x, b);
    log.debug("residual: " + residual);
    assertEquals(0., residual, Utils.getDoubleMachineEpsilon());
}

From source file:edu.ucdenver.bios.power.ConditionalOrthogonalPolynomial2FactorTest.java

/**
 * /* www  .j  ava  2  s  .  c  om*/
 * @param checker
 * @param params
 */
private void checkPower(PowerChecker checker, String title, GLMMPowerParameters params) {
    /* 
     * get orthogonal contrasts for within subject factors
     * Log base 2 spacing Clip (2,4,16) and Region(2,8,32) 
     */
    ArrayList<Factor> factorList = new ArrayList<Factor>();
    factorList.add(factorA);
    factorList.add(factorB);
    OrthogonalPolynomialContrastCollection collection = OrthogonalPolynomials.withinSubjectContrast(factorList);
    params.setWithinSubjectContrast(collection.getInteractionContrast(factorList).getContrastMatrix());

    // theta critical matrix used to back-tranform beta from U
    //    THETA = {.25}#{.5 1 -1 .5}; * =Theta(cr) from 1st sentence *after* 
    //  equation 7, Coffey and Muller (2003); 
    double[][] thetaData = { { 0.5, 1, -1, 0.5 } };
    RealMatrix thetaCr = (new Array2DRowRealMatrix(thetaData)).scalarMultiply(0.25);

    // loop over the various sigma matrices
    //Test[] testList = Test.values();
    GLMMTestFactory.Test[] testList = { GLMMTestFactory.Test.UNIREP, GLMMTestFactory.Test.UNIREP_BOX,
            GLMMTestFactory.Test.UNIREP_GEISSER_GREENHOUSE, GLMMTestFactory.Test.UNIREP_HUYNH_FELDT };
    for (GLMMTestFactory.Test test : testList) {
        params.clearTestList();
        params.addTest(test);

        for (RealMatrix sigStar : sigmaStars) {
            RealMatrix U = params.getWithinSubjectContrast();
            // 1st paragraph in section 2.4, Coffey and Muller 2003 *;
            RealMatrix sigmaTemp = U.multiply(sigStar).multiply(U.transpose());
            int dimension = sigmaTemp.getRowDimension();
            RealMatrix Uother = MatrixUtils.getRealMatrixWithFilledValue(dimension, 1,
                    1 / Math.sqrt(U.getColumnDimension()));
            Uother = MatrixUtils.getHorizontalAppend(Uother,
                    collection.getMainEffectContrast(factorA).getContrastMatrix());
            Uother = MatrixUtils.getHorizontalAppend(Uother,
                    collection.getMainEffectContrast(factorB).getContrastMatrix());
            double varianceMean = (double) sigStar.getTrace() / (double) sigStar.getColumnDimension();
            RealMatrix sigmaError = sigmaTemp
                    .add(Uother.multiply(Uother.transpose()).scalarMultiply(varianceMean));

            if (verbose)
                printMatrix("Sigma Error", sigmaError);
            // 1st paragraph in section 2.4, Coffey and Muller 2003 *; 
            RealMatrix beta = thetaCr.multiply(U.transpose());
            params.setSigmaError(sigmaError);
            params.setBeta(new FixedRandomMatrix(beta.getData(), null, false));

            checker.checkPower(params);
        }

    }

    // output the results
    ValidationReportBuilder reportBuilder = new ValidationReportBuilder();
    reportBuilder.createValidationReportAsStdout(checker, title, false);

    assertTrue("results outside tolerance: " + TOLERANCE, checker.isSASDeviationBelowTolerance(TOLERANCE));
}