1dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/* 2dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Licensed to the Apache Software Foundation (ASF) under one or more 3dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * contributor license agreements. See the NOTICE file distributed with 4dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * this work for additional information regarding copyright ownership. 5dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The ASF licenses this file to You under the Apache License, Version 2.0 6dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * (the "License"); you may not use this file except in compliance with 7dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the License. You may obtain a copy of the License at 8dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 9dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * http://www.apache.org/licenses/LICENSE-2.0 10dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 11dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Unless required by applicable law or agreed to in writing, software 12dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * distributed under the License is distributed on an "AS IS" BASIS, 13dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * See the License for the specific language governing permissions and 15dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * limitations under the License. 16dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 17dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 18dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpackage org.apache.commons.math.linear; 19dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 20dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport java.util.Arrays; 21dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 22dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.MathRuntimeException; 23dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.exception.util.LocalizedFormats; 24dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondimport org.apache.commons.math.util.FastMath; 25dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 26dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 27dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond/** 28dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Calculates the QR-decomposition of a matrix. 29dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>The QR-decomposition of a matrix A consists of two matrices Q and R 30dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * that satisfy: A = QR, Q is orthogonal (Q<sup>T</sup>Q = I), and R is 31dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * upper triangular. If A is m×n, Q is m×m and R m×n.</p> 32dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>This class compute the decomposition using Householder reflectors.</p> 33dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>For efficiency purposes, the decomposition in packed form is transposed. 34dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * This allows inner loop to iterate inside rows, which is much more cache-efficient 35dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * in Java.</p> 36dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 37dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see <a href="http://mathworld.wolfram.com/QRDecomposition.html">MathWorld</a> 38dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @see <a href="http://en.wikipedia.org/wiki/QR_decomposition">Wikipedia</a> 39dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 40dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 août 2010) $ 41dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @since 1.2 42dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 43dee0849a9704d532af0b550146cbafbaa6ee1d19Raymondpublic class QRDecompositionImpl implements QRDecomposition { 44dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 45dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 46dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * A packed TRANSPOSED representation of the QR decomposition. 47dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>The elements BELOW the diagonal are the elements of the UPPER triangular 48dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors 49dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * from which an explicit form of Q can be recomputed if desired.</p> 50dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 51dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double[][] qrt; 52dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 53dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** The diagonal elements of R. */ 54dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private double[] rDiag; 55dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 56dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Cached value of Q. */ 57dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private RealMatrix cachedQ; 58dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 59dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Cached value of QT. */ 60dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private RealMatrix cachedQT; 61dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 62dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Cached value of R. */ 63dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private RealMatrix cachedR; 64dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 65dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Cached value of H. */ 66dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private RealMatrix cachedH; 67dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 68dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 69dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Calculates the QR-decomposition of the given matrix. 70dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param matrix The matrix to decompose. 71dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 72dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public QRDecompositionImpl(RealMatrix matrix) { 73dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 74dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int m = matrix.getRowDimension(); 75dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int n = matrix.getColumnDimension(); 76dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond qrt = matrix.transpose().getData(); 77dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond rDiag = new double[FastMath.min(m, n)]; 78dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedQ = null; 79dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedQT = null; 80dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedR = null; 81dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedH = null; 82dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 83dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /* 84dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The QR decomposition of a matrix A is calculated using Householder 85dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * reflectors by repeating the following operations to each minor 86dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * A(minor,minor) of A: 87dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 88dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int minor = 0; minor < FastMath.min(m, n); minor++) { 89dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 90dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] qrtMinor = qrt[minor]; 91dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 92dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /* 93dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Let x be the first column of the minor, and a^2 = |x|^2. 94dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * x will be in the positions qr[minor][minor] through qr[m][minor]. 95dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The first column of the transformed minor will be (a,0,0,..)' 96dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * The sign of a is chosen to be opposite to the sign of the first 97dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * component of x. Let's find a: 98dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 99dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double xNormSqr = 0; 100dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int row = minor; row < m; row++) { 101dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double c = qrtMinor[row]; 102dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond xNormSqr += c * c; 103dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 104dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double a = (qrtMinor[minor] > 0) ? -FastMath.sqrt(xNormSqr) : FastMath.sqrt(xNormSqr); 105dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond rDiag[minor] = a; 106dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 107dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (a != 0.0) { 108dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 109dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /* 110dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Calculate the normalized reflection vector v and transform 111dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * the first column. We know the norm of v beforehand: v = x-ae 112dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * so |v|^2 = <x-ae,x-ae> = <x,x>-2a<x,e>+a^2<e,e> = 113dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * a^2+a^2-2a<x,e> = 2a*(a - <x,e>). 114dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Here <x, e> is now qr[minor][minor]. 115dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * v = x-ae is stored in the column at qr: 116dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 117dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond qrtMinor[minor] -= a; // now |v|^2 = -2a*(qr[minor][minor]) 118dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 119dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /* 120dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Transform the rest of the columns of the minor: 121dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * They will be transformed by the matrix H = I-2vv'/|v|^2. 122dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * If x is a column vector of the minor, then 123dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Hx = (I-2vv'/|v|^2)x = x-2vv'x/|v|^2 = x - 2<x,v>/|v|^2 v. 124dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Therefore the transformation is easily calculated by 125dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * subtracting the column vector (2<x,v>/|v|^2)v from x. 126dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * 127dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Let 2<x,v>/|v|^2 = alpha. From above we have 128dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * |v|^2 = -2a*(qr[minor][minor]), so 129dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * alpha = -<x,v>/(a*qr[minor][minor]) 130dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 131dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int col = minor+1; col < n; col++) { 132dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] qrtCol = qrt[col]; 133dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double alpha = 0; 134dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int row = minor; row < m; row++) { 135dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond alpha -= qrtCol[row] * qrtMinor[row]; 136dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 137dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond alpha /= a * qrtMinor[minor]; 138dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 139dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // Subtract the column vector alpha*v from x. 140dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int row = minor; row < m; row++) { 141dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond qrtCol[row] -= alpha * qrtMinor[row]; 142dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 143dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 144dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 145dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 146dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 147dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 148dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 149dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public RealMatrix getR() { 150dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 151dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (cachedR == null) { 152dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 153dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // R is supposed to be m x n 154dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int n = qrt.length; 155dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int m = qrt[0].length; 156dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedR = MatrixUtils.createRealMatrix(m, n); 157dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 158dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // copy the diagonal from rDiag and the upper triangle of qr 159dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int row = FastMath.min(m, n) - 1; row >= 0; row--) { 160dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedR.setEntry(row, row, rDiag[row]); 161dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int col = row + 1; col < n; col++) { 162dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedR.setEntry(row, col, qrt[col][row]); 163dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 164dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 165dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 166dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 167dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 168dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // return the cached matrix 169dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return cachedR; 170dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 171dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 172dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 173dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 174dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public RealMatrix getQ() { 175dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (cachedQ == null) { 176dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedQ = getQT().transpose(); 177dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 178dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return cachedQ; 179dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 180dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 181dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 182dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public RealMatrix getQT() { 183dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 184dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (cachedQT == null) { 185dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 186dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // QT is supposed to be m x m 187dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int n = qrt.length; 188dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int m = qrt[0].length; 189dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedQT = MatrixUtils.createRealMatrix(m, m); 190dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 191dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /* 192dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Q = Q1 Q2 ... Q_m, so Q is formed by first constructing Q_m and then 193dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * applying the Householder transformations Q_(m-1),Q_(m-2),...,Q1 in 194dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * succession to the result 195dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 196dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int minor = m - 1; minor >= FastMath.min(m, n); minor--) { 197dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedQT.setEntry(minor, minor, 1.0); 198dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 199dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 200dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int minor = FastMath.min(m, n)-1; minor >= 0; minor--){ 201dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] qrtMinor = qrt[minor]; 202dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedQT.setEntry(minor, minor, 1.0); 203dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (qrtMinor[minor] != 0.0) { 204dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int col = minor; col < m; col++) { 205dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double alpha = 0; 206dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int row = minor; row < m; row++) { 207dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond alpha -= cachedQT.getEntry(col, row) * qrtMinor[row]; 208dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 209dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond alpha /= rDiag[minor] * qrtMinor[minor]; 210dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 211dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int row = minor; row < m; row++) { 212dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedQT.addToEntry(col, row, -alpha * qrtMinor[row]); 213dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 214dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 215dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 216dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 217dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 218dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 219dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 220dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // return the cached matrix 221dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return cachedQT; 222dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 223dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 224dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 225dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 226dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public RealMatrix getH() { 227dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 228dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (cachedH == null) { 229dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 230dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int n = qrt.length; 231dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int m = qrt[0].length; 232dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedH = MatrixUtils.createRealMatrix(m, n); 233dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < m; ++i) { 234dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = 0; j < FastMath.min(i + 1, n); ++j) { 235dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond cachedH.setEntry(i, j, qrt[j][i] / -rDiag[j]); 236dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 237dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 238dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 239dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 240dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 241dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // return the cached matrix 242dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return cachedH; 243dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 244dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 245dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 246dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 247dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public DecompositionSolver getSolver() { 248dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return new Solver(qrt, rDiag); 249dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 250dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 251dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Specialized solver. */ 252dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private static class Solver implements DecompositionSolver { 253dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 254dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 255dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * A packed TRANSPOSED representation of the QR decomposition. 256dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>The elements BELOW the diagonal are the elements of the UPPER triangular 257dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * matrix R, and the rows ABOVE the diagonal are the Householder reflector vectors 258dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * from which an explicit form of Q can be recomputed if desired.</p> 259dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 260dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double[][] qrt; 261dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 262dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** The diagonal elements of R. */ 263dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private final double[] rDiag; 264dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 265dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** 266dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * Build a solver from decomposed matrix. 267dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param qrt packed TRANSPOSED representation of the QR decomposition 268dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param rDiag diagonal elements of R 269dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 270dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond private Solver(final double[][] qrt, final double[] rDiag) { 271dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.qrt = qrt; 272dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond this.rDiag = rDiag; 273dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 274dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 275dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 276dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public boolean isNonSingular() { 277dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 278dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (double diag : rDiag) { 279dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (diag == 0) { 280dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return false; 281dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 282dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 283dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return true; 284dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 285dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 286dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 287dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 288dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public double[] solve(double[] b) 289dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws IllegalArgumentException, InvalidMatrixException { 290dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 291dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int n = qrt.length; 292dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int m = qrt[0].length; 293dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (b.length != m) { 294dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw MathRuntimeException.createIllegalArgumentException( 295dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond LocalizedFormats.VECTOR_LENGTH_MISMATCH, 296dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond b.length, m); 297dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 298dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (!isNonSingular()) { 299dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new SingularMatrixException(); 300dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 301dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 302dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] x = new double[n]; 303dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] y = b.clone(); 304dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 305dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // apply Householder transforms to solve Q.y = b 306dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int minor = 0; minor < FastMath.min(m, n); minor++) { 307dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 308dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] qrtMinor = qrt[minor]; 309dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond double dotProduct = 0; 310dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int row = minor; row < m; row++) { 311dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond dotProduct += y[row] * qrtMinor[row]; 312dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 313dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond dotProduct /= rDiag[minor] * qrtMinor[minor]; 314dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 315dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int row = minor; row < m; row++) { 316dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond y[row] += dotProduct * qrtMinor[row]; 317dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 318dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 319dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 320dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 321dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // solve triangular system R.x = y 322dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int row = rDiag.length - 1; row >= 0; --row) { 323dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond y[row] /= rDiag[row]; 324dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double yRow = y[row]; 325dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] qrtRow = qrt[row]; 326dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond x[row] = yRow; 327dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < row; i++) { 328dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond y[i] -= yRow * qrtRow[i]; 329dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 330dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 331dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 332dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return x; 333dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 334dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 335dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 336dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 337dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public RealVector solve(RealVector b) 338dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws IllegalArgumentException, InvalidMatrixException { 339dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond try { 340dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return solve((ArrayRealVector) b); 341dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } catch (ClassCastException cce) { 342dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return new ArrayRealVector(solve(b.getData()), false); 343dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 344dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 345dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 346dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** Solve the linear equation A × X = B. 347dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * <p>The A matrix is implicit here. It is </p> 348dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @param b right-hand side of the equation A × X = B 349dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @return a vector X that minimizes the two norm of A × X - B 350dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws IllegalArgumentException if matrices dimensions don't match 351dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond * @throws InvalidMatrixException if decomposed matrix is singular 352dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond */ 353dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public ArrayRealVector solve(ArrayRealVector b) 354dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws IllegalArgumentException, InvalidMatrixException { 355dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return new ArrayRealVector(solve(b.getDataRef()), false); 356dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 357dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 358dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 359dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public RealMatrix solve(RealMatrix b) 360dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws IllegalArgumentException, InvalidMatrixException { 361dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 362dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int n = qrt.length; 363dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int m = qrt[0].length; 364dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (b.getRowDimension() != m) { 365dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw MathRuntimeException.createIllegalArgumentException( 366dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond LocalizedFormats.DIMENSIONS_MISMATCH_2x2, 367dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond b.getRowDimension(), b.getColumnDimension(), m, "n"); 368dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 369dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond if (!isNonSingular()) { 370dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throw new SingularMatrixException(); 371dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 372dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 373dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int columns = b.getColumnDimension(); 374dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int blockSize = BlockRealMatrix.BLOCK_SIZE; 375dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int cBlocks = (columns + blockSize - 1) / blockSize; 376dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[][] xBlocks = BlockRealMatrix.createBlocksLayout(n, columns); 377dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[][] y = new double[b.getRowDimension()][blockSize]; 378dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] alpha = new double[blockSize]; 379dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 380dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int kBlock = 0; kBlock < cBlocks; ++kBlock) { 381dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int kStart = kBlock * blockSize; 382dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int kEnd = FastMath.min(kStart + blockSize, columns); 383dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int kWidth = kEnd - kStart; 384dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 385dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // get the right hand side vector 386dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond b.copySubMatrix(0, m - 1, kStart, kEnd - 1, y); 387dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 388dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // apply Householder transforms to solve Q.y = b 389dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int minor = 0; minor < FastMath.min(m, n); minor++) { 390dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] qrtMinor = qrt[minor]; 391dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double factor = 1.0 / (rDiag[minor] * qrtMinor[minor]); 392dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 393dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond Arrays.fill(alpha, 0, kWidth, 0.0); 394dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int row = minor; row < m; ++row) { 395dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double d = qrtMinor[row]; 396dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] yRow = y[row]; 397dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int k = 0; k < kWidth; ++k) { 398dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond alpha[k] += d * yRow[k]; 399dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 400dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 401dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int k = 0; k < kWidth; ++k) { 402dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond alpha[k] *= factor; 403dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 404dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 405dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int row = minor; row < m; ++row) { 406dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double d = qrtMinor[row]; 407dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] yRow = y[row]; 408dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int k = 0; k < kWidth; ++k) { 409dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yRow[k] += alpha[k] * d; 410dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 411dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 412dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 413dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 414dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 415dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond // solve triangular system R.x = y 416dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int j = rDiag.length - 1; j >= 0; --j) { 417dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int jBlock = j / blockSize; 418dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final int jStart = jBlock * blockSize; 419dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double factor = 1.0 / rDiag[j]; 420dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] yJ = y[j]; 421dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] xBlock = xBlocks[jBlock * cBlocks + kBlock]; 422dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond int index = (j - jStart) * kWidth; 423dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int k = 0; k < kWidth; ++k) { 424dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yJ[k] *= factor; 425dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond xBlock[index++] = yJ[k]; 426dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 427dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 428dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] qrtJ = qrt[j]; 429dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int i = 0; i < j; ++i) { 430dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double rIJ = qrtJ[i]; 431dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond final double[] yI = y[i]; 432dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond for (int k = 0; k < kWidth; ++k) { 433dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond yI[k] -= yJ[k] * rIJ; 434dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 435dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 436dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 437dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 438dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 439dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 440dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 441dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return new BlockRealMatrix(n, columns, xBlocks, false); 442dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 443dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 444dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 445dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond /** {@inheritDoc} */ 446dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond public RealMatrix getInverse() 447dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond throws InvalidMatrixException { 448dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond return solve(MatrixUtils.createRealIdentityMatrix(rDiag.length)); 449dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 450dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 451dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond } 452dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond 453dee0849a9704d532af0b550146cbafbaa6ee1d19Raymond} 454