Example usage for org.apache.commons.math3.linear QRDecomposition getQ

Introduction

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

Prototype

public RealMatrix getQ() 

Source Link

Document

Returns the matrix Q of the decomposition.

Usage

From source file:edu.cudenver.bios.matrix.OrthogonalPolynomials.java

/**
 * Computes orthogonal polynomial contrasts for the specified data values.  Currently only
 * supports fitting (not prediction contrasts).  
 * /*from   ww w.  ja v a2 s  . c  o  m*/
 *    @param x the points at which the polynomials will be evaluated
 * @param maxDegree contrasts will be computed for degrees 1 to maxDegree
 * @return matrix containing 0th,1st, 2nd,...maxDegree-th degree contrasts in each column
 * @throws IllegalArgumentException
 */
public static RealMatrix orthogonalPolynomialCoefficients(double[] x, int maxDegree)
        throws IllegalArgumentException {
    if (x == null)
        throw new IllegalArgumentException("no data specified");
    if (maxDegree < 1)
        throw new IllegalArgumentException("max polynomial degree must be greater than 1");
    // count number of unique values
    HashSet<Double> s = new HashSet<Double>();
    for (double i : x)
        s.add(i);
    int uniqueCount = s.size();
    if (maxDegree >= uniqueCount)
        throw new IllegalArgumentException(
                "max polynomial degree must be less than the number of unique points");

    // center the data
    double xBar = StatUtils.mean(x);
    double[] xCentered = new double[x.length];
    for (int i = 0; i < x.length; i++)
        xCentered[i] = x[i] - xBar;
    // compute an "outer product" of the centered x vector and a vector 
    // containing the sequence 0 to maxDegree-1, but raise the x values
    // to the power in the sequence array
    double[][] xOuter = new double[x.length][maxDegree + 1];
    int row = 0;
    for (double xValue : xCentered) {
        for (int col = 0; col <= maxDegree; col++) {
            xOuter[row][col] = Math.pow(xValue, col);
        }
        row++;
    }
    // do some mysterious QR decomposition stuff.  See Emerson (1968)
    RealMatrix outerVector = new Array2DRowRealMatrix(xOuter);
    QRDecomposition qrDecomp = new QRDecomposition(outerVector);

    RealMatrix z = MatrixUtils.getDiagonalMatrix(qrDecomp.getR());
    RealMatrix raw = qrDecomp.getQ().multiply(z);

    // column sum of squared elements in raw
    double[] normalizingConstants = new double[raw.getColumnDimension()];
    for (int col = 0; col < raw.getColumnDimension(); col++) {
        normalizingConstants[col] = 0;
        for (row = 0; row < raw.getRowDimension(); row++) {
            double value = raw.getEntry(row, col);
            normalizingConstants[col] += value * value;
        }
    }

    // now normalize the raw values
    for (int col = 0; col < raw.getColumnDimension(); col++) {
        double normalConstantSqrt = Math.sqrt(normalizingConstants[col]);
        for (row = 0; row < raw.getRowDimension(); row++) {
            raw.setEntry(row, col, raw.getEntry(row, col) / normalConstantSqrt);
        }
    }

    return raw;
}

---

From source file:gamlss.utilities.MatrixFunctions.java

/**
* <p>Compute the "hat" matrix.//w ww .  ja  v  a  2  s. co  m
* </p>
* <p>The hat matrix is defined in terms of the design matrix X
*  by X(X<sup>T</sup>X)<sup>-1</sup>X<sup>T</sup>
* </p>
* <p>The implementation here uses the QR decomposition to compute the
* hat matrix as Q I<sub>p</sub>Q<sup>T</sup> where I<sub>p</sub> is the
* p-dimensional identity matrix augmented by 0's.  This computational
* formula is from "The Hat Matrix in Regression and ANOVA",
* David C. Hoaglin and Roy E. Welsch,
* <i>The American Statistician</i>, Vol. 32, No. 1 (Feb., 1978), pp. 17-22.
*@param m - matrix
* @return the hat matrix
*/
public static RealMatrix calculateHat(final BlockRealMatrix m) {
    QRDecomposition qr = new QRDecomposition(m);
    // Create augmented identity matrix
    RealMatrix qM = qr.getQ();
    final int p = qr.getR().getColumnDimension();
    final int n = qM.getColumnDimension();
    Array2DRowRealMatrix augI = new Array2DRowRealMatrix(n, n);
    double[][] augIData = augI.getDataRef();
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            if (i == j && i < p) {
                augIData[i][j] = 1d;
            } else {
                augIData[i][j] = 0d;
            }
        }
    }
    // Compute and return Hat matrix
    return qM.multiply(augI).multiply(qM.transpose());
}

---

From source file:com.joptimizer.util.UtilsTest.java

public void testQRdecomposition() throws Exception {
    log.debug("testQRdecomposition");
    double[][] A = new double[][] { { 1, 0, 0 }, { 2, 2, 0 }, { 3, 3, 3 } };
    QRDecomposition qrFact = new QRDecomposition(MatrixUtils.createRealMatrix(A));
    qrFact.getQ();
}

---

From source file:com.opengamma.strata.math.impl.linearalgebra.QRDecompositionCommonsResult.java

/**
 * @param qr The result of the QR decomposition, not null
 */// w  w  w  . j  a  v a2  s.  co m
public QRDecompositionCommonsResult(QRDecomposition qr) {
    ArgChecker.notNull(qr, "qr");
    _q = CommonsMathWrapper.unwrap(qr.getQ());
    _r = CommonsMathWrapper.unwrap(qr.getR());
    _qTranspose = _q.transpose();
    _solver = qr.getSolver();
}

---

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

public void testDecompose() throws Exception {
    log.debug("testDecompose");
    final double[][] A = new double[][] { { 1, 0, 0, 2 }, { 0, 0, 2, 0 }, { 2, 3, 0, 0 }, { 0, 0, 4, 4 } };

    double[][] EQ = { { 0.447214, -0.894427, 0., 0. }, { 0., 0., 0.447214, -0.894427 },
            { 0.894427, 0.447214, 0., 0. }, { 0., 0., 0.894427, 0.447214 } };
    double[][] ER = { { 2.23607, 2.68328, 0., 0.894427 }, { 0., 1.34164, 0., -1.78885 },
            { 0., 0., 4.47214, 3.57771 }, { 0., 0., 0., 1.78885 } };

    QRDecomposition dFact = new QRDecomposition(new Array2DRowRealMatrix(A));
    RealMatrix Q = dFact.getQ();
    RealMatrix R = dFact.getR();//w  ww . j a  va  2  s .c  o  m
    RealMatrix H = dFact.getH();
    log.debug("Q: " + ArrayUtils.toString(Q.getData()));
    log.debug("R: " + ArrayUtils.toString(R.getData()));
    //log.debug("H: " + ArrayUtils.toString(H.getData()));

    SparseDoubleMatrix2D S = new SparseDoubleMatrix2D(A);
    QRSparseFactorization qr = new QRSparseFactorization(S);
    qr.factorize();
    log.debug("R: " + ArrayUtils.toString(qr.getR().toArray()));
    for (int i = 0; i < R.getRowDimension(); i++) {
        for (int j = 0; j < R.getColumnDimension(); j++) {
            assertEquals(ER[i][j], qr.getR().getQuick(i, j), 1.e-5);
        }
    }
    assertTrue(qr.hasFullRank());
}

---

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"
 *//*from  w w  w  .jav a 2s  . c  o  m*/
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:com.itemanalysis.psychometrics.factoranalysis.GPArotation.java

private RealMatrix randomStart(int ncol) {
    NormalDistribution norm = new NormalDistribution(0.0, 1.0);
    RealMatrix T = new Array2DRowRealMatrix(ncol, ncol);
    for (int i = 0; i < ncol; i++) {
        for (int j = 0; j < ncol; j++) {
            T.setEntry(i, j, norm.sample());
        }//  w  ww.j a  v  a2  s .c o  m
    }
    QRDecomposition qr = new QRDecomposition(T);
    return qr.getQ();
}

---

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

/**
 * This test is with a rank-deficient matrix, and must reveal it.
 * It is a 381x694 matrix with rank=379 (as given by Mathematica).
 * Some doubt on the rank of this matrix.
 *//*w ww. j  ava2  s.c  o  m*/
public void testSimpleRankDeficient() throws Exception {
    log.debug("testSimpleRankDeficient");

    double[][] A = new double[][] { { 1, 2, 3 }, { 4, 5, 6 }, { 7, 8, 9 } };
    log.debug("A  : " + ArrayUtils.toString(A));

    //JAMA
    //      Matrix M1 = new Matrix(A);
    //      Jama.QRDecomposition qr1 = new Jama.QRDecomposition(M1);
    //      assertFalse(qr1.isFullRank());

    //COMMONS-MATH 3
    RealMatrix M2 = MatrixUtils.createRealMatrix(A);
    org.apache.commons.math3.linear.QRDecomposition qr2 = new org.apache.commons.math3.linear.QRDecomposition(
            M2, Utils.getDoubleMachineEpsilon());
    log.debug("Q: " + qr2.getQ());
    log.debug("R: " + qr2.getR());
    //assertFalse(qr2.getSolver().isNonSingular());//this fails

    //COLT
    //NOTE: COLT fails because it checks the diagonal elements of R againts 0, 
    //while JOptimizer tests against a threshold
    SparseDoubleMatrix2D M3 = new SparseDoubleMatrix2D(A);
    cern.colt.matrix.linalg.QRDecomposition qr3 = new cern.colt.matrix.linalg.QRDecomposition(M3);
    log.debug("Q: " + qr3.getQ());
    log.debug("R: " + qr3.getR());
    //assertFalse(qr3.hasFullRank());//this fails

    //PARALLEL-COLT
    //      SparseCCDoubleMatrix2D M4 = new SparseCCDoubleMatrix2D(A);
    //      cern.colt.matrix.tdouble.algo.decomposition.SparseDoubleQRDecomposition qr4 = new cern.colt.matrix.tdouble.algo.decomposition.SparseDoubleQRDecomposition(M4, 0);
    //      assertFalse(qr4.hasFullRank());

    //JOptimizer
    //NOTE: COLT fails because it checks the diagonal elements of R againts 0, 
    //while JOptimizer tests against a threshold
    SparseDoubleMatrix2D M5 = new SparseDoubleMatrix2D(A);
    QRSparseFactorization qr5 = new QRSparseFactorization(M5);
    qr5.factorize();
    assertFalse(qr5.hasFullRank());

}

---

From source file:gamlss.smoothing.PB.java

/**
  * GCV smoothing method./*from  w  w  w .j a v a2  s  . c o  m*/
  * @param lambda - smoothing parameter
  * @return fitted values
  */
private ArrayRealVector functionGCV(Double lambda) {

    //QR <-qr(sqrt(w)*X)
    //Rinv <- solve(qr.R(QR))
    QRDecomposition qr = new QRDecomposition(
            MatrixFunctions.multVectorMatrix(MatrixFunctions.sqrtVec(w), basisM));
    rM = (BlockRealMatrix) qr.getR();
    rM = rM.getSubMatrix(0, rM.getColumnDimension() - 1, 0, rM.getColumnDimension() - 1);
    rMinv = (BlockRealMatrix) new QRDecomposition(rM).getSolver().getInverse();

    //wy <- sqrt(w)*y
    ArrayRealVector wy = MatrixFunctions.sqrtVec(w).ebeMultiply(y);

    //y.y <- sum(wy^2)
    double y_y = MatrixFunctions.sumV(wy.ebeMultiply(wy));

    //S   <- t(D)%*%D
    sM = penaltyM.transpose().multiply(penaltyM);

    //UDU <- eigen(t(Rinv)%*%S%*%Rinv)
    uduM = new BlockRealMatrix(
            new EigenDecomposition(rMinv.transpose().multiply(sM).multiply(rMinv), Controls.SPLIT_TOLERANCE)
                    .getV().getData());

    //yy <- t(UDU$vectors)%*%t(qr.Q(QR))%*%wy
    BlockRealMatrix qM = (BlockRealMatrix) qr.getQ();
    //SingularValueDecomposition svd = new SingularValueDecomposition(
    //MatrixFunctions.multVectorMatrix(MatrixFunctions.sqrtVec(w), basisM));
    //BlockRealMatrix qM = new BlockRealMatrix(svd.getV().getData());

    //.... to finish !!!!!!!
    MatrixFunctions.matrixPrint(qM);
    return null;
}

---

From source file:org.meteoinfo.math.linalg.LinalgUtil.java

/**
 * Calculates the QR-decomposition of a matrix.
 * The QR-decomposition of a matrix A consists of two matrices Q and R that satisfy: A = QR, Q 
 * is orthogonal (QTQ = I), and R is upper triangular. If A is mn, Q is mm and R mn.
 * @param a Given matrix.//from  ww w  .  j  av a2  s .c o  m
 * @return Result Q/R arrays.
 */
public static Array[] qr(Array a) {
    int m = a.getShape()[0];
    int n = a.getShape()[1];
    Array Qa = Array.factory(DataType.DOUBLE, new int[] { m, m });
    Array Ra = Array.factory(DataType.DOUBLE, a.getShape());
    double[][] aa = (double[][]) ArrayUtil.copyToNDJavaArray(a);
    RealMatrix matrix = new Array2DRowRealMatrix(aa, false);
    QRDecomposition decomposition = new QRDecomposition(matrix);
    RealMatrix Q = decomposition.getQ();
    RealMatrix R = decomposition.getR();
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < m; j++) {
            Qa.setDouble(i * m + j, Q.getEntry(i, j));
        }
    }
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            Ra.setDouble(i * n + j, R.getEntry(i, j));
        }
    }

    return new Array[] { Qa, Ra };
}

---

• «
• 1
• 2
• »