List of usage examples for org.apache.commons.math3.linear QRDecomposition getR
public RealMatrix getR()
From source file:com.ibm.bi.dml.runtime.matrix.data.LibCommonsMath.java
/** * Function to perform QR decomposition on a given matrix. * /*from w w w . j ava 2 s .com*/ * @param in * @return * @throws DMLRuntimeException */ private static MatrixBlock[] computeQR(MatrixObject in) throws DMLRuntimeException { Array2DRowRealMatrix matrixInput = DataConverter.convertToArray2DRowRealMatrix(in); // Perform QR decomposition QRDecomposition qrdecompose = new QRDecomposition(matrixInput); RealMatrix H = qrdecompose.getH(); RealMatrix R = qrdecompose.getR(); // Read the results into native format MatrixBlock mbH = DataConverter.convertToMatrixBlock(H.getData()); MatrixBlock mbR = DataConverter.convertToMatrixBlock(R.getData()); return new MatrixBlock[] { mbH, mbR }; }
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). * //w w w .j av a2 s . co 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.//from ww w . j av a2 s . c o 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.opengamma.strata.math.impl.linearalgebra.QRDecompositionCommonsResult.java
/** * @param qr The result of the QR decomposition, not null *//* w ww.j a v a2 s. c o 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();//from ww w .j a v a 2 s. c o m RealMatrix R = dFact.getR(); 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 ww . j a v a2s . 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:edu.stanford.cfuller.imageanalysistools.fitting.BisquareLinearFit.java
/** * Calculates the leverages of data points for least squares fitting (assuming equal variances). * * @param indVarValues The values of the independent variable used for the fitting. * @return a RealVector containing a leverage value for each independent variable value. *///from w ww. ja v a2s . c o m protected RealVector calculateLeverages(RealVector indVarValues) { RealMatrix indVarMatrix = null; if (this.noIntercept) { indVarMatrix = new Array2DRowRealMatrix(indVarValues.getDimension(), 1); } else { indVarMatrix = new Array2DRowRealMatrix(indVarValues.getDimension(), 2); } indVarMatrix.setColumnVector(0, indVarValues); if (!this.noIntercept) { indVarMatrix.setColumnVector(1, indVarValues.mapMultiply(0).mapAdd(1)); } RealVector leverages = new ArrayRealVector(indVarValues.getDimension()); QRDecomposition xQR = new QRDecomposition(indVarMatrix); RealMatrix xR = xQR.getR(); int smallerDim = xR.getRowDimension() < xR.getColumnDimension() ? xR.getRowDimension() : xR.getColumnDimension(); RealMatrix xRSq = xR.getSubMatrix(0, smallerDim - 1, 0, smallerDim - 1); QRDecomposition xRQR = new QRDecomposition(xRSq); RealMatrix xRInv = xRQR.getSolver().getInverse(); RealMatrix xxRInv = indVarMatrix.multiply(xRInv); for (int i = 0; i < indVarValues.getDimension(); i++) { double sum = 0; for (int j = 0; j < xxRInv.getColumnDimension(); j++) { sum += Math.pow(xxRInv.getEntry(i, j), 2); } leverages.setEntry(i, sum); } return leverages; }
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 w w . j a v a 2s . 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:edu.cmu.tetrad.search.Lingam.java
/** * This is the method used in Patrik's code. *///from w w w. j a v a 2 s . c o m public TetradMatrix pruneEdgesByResampling(TetradMatrix data, int[] k) { if (k.length != data.columns()) { throw new IllegalArgumentException("Execting a permutation."); } Set<Integer> set = new LinkedHashSet<Integer>(); for (int i = 0; i < k.length; i++) { if (k[i] >= k.length) { throw new IllegalArgumentException("Expecting a permutation."); } if (set.contains(i)) { throw new IllegalArgumentException("Expecting a permutation."); } set.add(i); } TetradMatrix X = data.transpose(); int npieces = 10; int cols = X.columns(); int rows = X.rows(); int piecesize = (int) Math.floor(cols / npieces); List<TetradMatrix> bpieces = new ArrayList<TetradMatrix>(); // List<Vector> diststdpieces = new ArrayList<Vector>(); // List<Vector> cpieces = new ArrayList<Vector>(); for (int p = 0; p < npieces; p++) { // % Select subset of data, and permute the variables to the causal order // Xp = X(k,((p-1)*piecesize+1):(p*piecesize)); int p0 = (p) * piecesize; int p1 = (p + 1) * piecesize - 1; int[] range = range(p0, p1); TetradMatrix Xp = X.getSelection(k, range); // % Remember to subract out the mean // Xpm = mean(Xp,2); // Xp = Xp - Xpm*ones(1,size(Xp,2)); // // % Calculate covariance matrix // cov = (Xp*Xp')/size(Xp,2); double[] Xpm = new double[rows]; for (int i = 0; i < rows; i++) { double sum = 0.0; for (int j = 0; j < Xp.columns(); j++) { sum += Xp.get(i, j); } Xpm[i] = sum / Xp.columns(); } for (int i = 0; i < rows; i++) { for (int j = 0; j < Xp.columns(); j++) { Xp.set(i, j, Xp.get(i, j) - Xpm[i]); } } TetradMatrix cov = Xp.times(Xp.transpose()); // for (int i = 0; i < cov.rows(); i++) { // for (int j = 0; j < cov.columns(); j++) { // cov.set(i, j, cov.get(i, j) / Xp.columns()); // } // } // % Do QL decomposition on the inverse square root of cov // [Q,L] = tridecomp(cov^(-0.5),'ql'); boolean posDef = MatrixUtils.isPositiveDefinite(cov); // TetradLogger.getInstance().log("lingamDetails","Positive definite = " + posDef); if (!posDef) { System.out.println("Covariance matrix is not positive definite."); } TetradMatrix invSqrt = cov.sqrt().inverse(); QRDecomposition qr = new QRDecomposition(invSqrt.getRealMatrix()); RealMatrix r = qr.getR(); // % The estimated disturbance-stds are one over the abs of the diag of L // newestdisturbancestd = 1./diag(abs(L)); TetradVector newestdisturbancestd = new TetradVector(rows); for (int t = 0; t < rows; t++) { newestdisturbancestd.set(t, 1.0 / abs(r.getEntry(t, t))); } // % Normalize rows of L to unit diagonal // L = L./(diag(L)*ones(1,dims)); // for (int s = 0; s < rows; s++) { for (int t = 0; t < min(s, cols); t++) { r.setEntry(s, t, r.getEntry(s, t) / r.getEntry(s, s)); } } // % Calculate corresponding B // bnewest = eye(dims)-L; TetradMatrix bnewest = TetradMatrix.identity(rows); bnewest = bnewest.minus(new TetradMatrix(r)); // % Also calculate constants // cnewest = L*Xpm; // Vector cnewest = new DenseVector(rows); // cnewest = L.mult(new DenseVector(Xpm), cnewest); // % Permute back to original variable order // ik = iperm(k); // bnewest = bnewest(ik, ik); // newestdisturbancestd = newestdisturbancestd(ik); // cnewest = cnewest(ik); int[] ik = iperm(k); // System.out.println("ik = " + Arrays.toString(ik)); bnewest = bnewest.getSelection(ik, ik); // newestdisturbancestd = Matrices.getSubVector(newestdisturbancestd, ik); // cnewest = Matrices.getSubVector(cnewest, ik); // % Save results // Bpieces(:,:,p) = bnewest; // diststdpieces(:,p) = newestdisturbancestd; // cpieces(:,p) = cnewest; bpieces.add(bnewest); // diststdpieces.add(newestdisturbancestd); // cpieces.add(cnewest); // // end } TetradMatrix means = new TetradMatrix(rows, rows); TetradMatrix stds = new TetradMatrix(rows, rows); TetradMatrix BFinal = new TetradMatrix(rows, rows); for (int i = 0; i < rows; i++) { for (int j = 0; j < rows; j++) { double[] b = new double[npieces]; for (int y = 0; y < npieces; y++) { b[y] = bpieces.get(y).get(i, j); } double themean = StatUtils.mean(b); double thestd = StatUtils.sd(b); means.set(i, j, themean); stds.set(i, j, thestd); if (abs(themean) < getPruneFactor() * thestd) { BFinal.set(i, j, 0); } else { BFinal.set(i, j, themean); } } } return BFinal; }
From source file:gamlss.smoothing.PB.java
/** * GCV smoothing method.//from w w w . j av a 2s . 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; }