001/* 002 * BioJava development code 003 * 004 * This code may be freely distributed and modified under the 005 * terms of the GNU Lesser General Public Licence. This should 006 * be distributed with the code. If you do not have a copy, 007 * see: 008 * 009 * http://www.gnu.org/copyleft/lesser.html 010 * 011 * Copyright for this code is held jointly by the individual 012 * authors. These should be listed in @author doc comments. 013 * 014 * For more information on the BioJava project and its aims, 015 * or to join the biojava-l mailing list, visit the home page 016 * at: 017 * 018 * http://www.biojava.org/ 019 * 020 */ 021package org.biojava.nbio.structure.jama; 022 023 024/** QR Decomposition. 025<P> 026 For an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n 027 orthogonal matrix Q and an n-by-n upper triangular matrix R so that 028 A = Q*R. 029<P> 030 The QR decompostion always exists, even if the matrix does not have 031 full rank, so the constructor will never fail. The primary use of the 032 QR decomposition is in the least squares solution of nonsquare systems 033 of simultaneous linear equations. This will fail if isFullRank() 034 returns false. 035 */ 036 037public class QRDecomposition implements java.io.Serializable { 038 039 static final long serialVersionUID = 10293720387423l; 040 041/* ------------------------ 042 Class variables 043 * ------------------------ */ 044 045 /** Array for internal storage of decomposition. 046 @serial internal array storage. 047 */ 048 private double[][] QR; 049 050 /** Row and column dimensions. 051 @serial column dimension. 052 @serial row dimension. 053 */ 054 private int m, n; 055 056 /** Array for internal storage of diagonal of R. 057 @serial diagonal of R. 058 */ 059 private double[] Rdiag; 060 061/* ------------------------ 062 Constructor 063 * ------------------------ */ 064 065 /** QR Decomposition, computed by Householder reflections. provides 066 a data structure to access R and the Householder vectors and 067 compute Q. 068 @param A Rectangular matrix 069 */ 070 071 public QRDecomposition (Matrix A) { 072 // Initialize. 073 QR = A.getArrayCopy(); 074 m = A.getRowDimension(); 075 n = A.getColumnDimension(); 076 Rdiag = new double[n]; 077 078 // Main loop. 079 for (int k = 0; k < n; k++) { 080 // Compute 2-norm of k-th column without under/overflow. 081 double nrm = 0; 082 for (int i = k; i < m; i++) { 083 nrm = Maths.hypot(nrm,QR[i][k]); 084 } 085 086 if (nrm != 0.0) { 087 // Form k-th Householder vector. 088 if (QR[k][k] < 0) { 089 nrm = -nrm; 090 } 091 for (int i = k; i < m; i++) { 092 QR[i][k] /= nrm; 093 } 094 QR[k][k] += 1.0; 095 096 // Apply transformation to remaining columns. 097 for (int j = k+1; j < n; j++) { 098 double s = 0.0; 099 for (int i = k; i < m; i++) { 100 s += QR[i][k]*QR[i][j]; 101 } 102 s = -s/QR[k][k]; 103 for (int i = k; i < m; i++) { 104 QR[i][j] += s*QR[i][k]; 105 } 106 } 107 } 108 Rdiag[k] = -nrm; 109 } 110 } 111 112/* ------------------------ 113 Public Methods 114 * ------------------------ */ 115 116 /** Is the matrix full rank? 117 @return true if R, and hence A, has full rank. 118 */ 119 120 public boolean isFullRank () { 121 for (int j = 0; j < n; j++) { 122 if (Rdiag[j] == 0) 123 return false; 124 } 125 return true; 126 } 127 128 /** Return the Householder vectors 129 @return Lower trapezoidal matrix whose columns define the reflections 130 */ 131 132 public Matrix getH () { 133 Matrix X = new Matrix(m,n); 134 double[][] H = X.getArray(); 135 for (int i = 0; i < m; i++) { 136 for (int j = 0; j < n; j++) { 137 if (i >= j) { 138 H[i][j] = QR[i][j]; 139 } else { 140 H[i][j] = 0.0; 141 } 142 } 143 } 144 return X; 145 } 146 147 /** Return the upper triangular factor 148 @return R 149 */ 150 151 public Matrix getR () { 152 Matrix X = new Matrix(n,n); 153 double[][] R = X.getArray(); 154 for (int i = 0; i < n; i++) { 155 for (int j = 0; j < n; j++) { 156 if (i < j) { 157 R[i][j] = QR[i][j]; 158 } else if (i == j) { 159 R[i][j] = Rdiag[i]; 160 } else { 161 R[i][j] = 0.0; 162 } 163 } 164 } 165 return X; 166 } 167 168 /** Generate and return the (economy-sized) orthogonal factor 169 @return Q 170 */ 171 172 public Matrix getQ () { 173 Matrix X = new Matrix(m,n); 174 double[][] Q = X.getArray(); 175 for (int k = n-1; k >= 0; k--) { 176 for (int i = 0; i < m; i++) { 177 Q[i][k] = 0.0; 178 } 179 Q[k][k] = 1.0; 180 for (int j = k; j < n; j++) { 181 if (QR[k][k] != 0) { 182 double s = 0.0; 183 for (int i = k; i < m; i++) { 184 s += QR[i][k]*Q[i][j]; 185 } 186 s = -s/QR[k][k]; 187 for (int i = k; i < m; i++) { 188 Q[i][j] += s*QR[i][k]; 189 } 190 } 191 } 192 } 193 return X; 194 } 195 196 /** Least squares solution of A*X = B 197 @param B A Matrix with as many rows as A and any number of columns. 198 @return X that minimizes the two norm of Q*R*X-B. 199 @exception IllegalArgumentException Matrix row dimensions must agree. 200 @exception RuntimeException Matrix is rank deficient. 201 */ 202 203 public Matrix solve (Matrix B) { 204 if (B.getRowDimension() != m) { 205 throw new IllegalArgumentException("Matrix row dimensions must agree."); 206 } 207 if (!this.isFullRank()) { 208 throw new RuntimeException("Matrix is rank deficient."); 209 } 210 211 // Copy right hand side 212 int nx = B.getColumnDimension(); 213 double[][] X = B.getArrayCopy(); 214 215 // Compute Y = transpose(Q)*B 216 for (int k = 0; k < n; k++) { 217 for (int j = 0; j < nx; j++) { 218 double s = 0.0; 219 for (int i = k; i < m; i++) { 220 s += QR[i][k]*X[i][j]; 221 } 222 s = -s/QR[k][k]; 223 for (int i = k; i < m; i++) { 224 X[i][j] += s*QR[i][k]; 225 } 226 } 227 } 228 // Solve R*X = Y; 229 for (int k = n-1; k >= 0; k--) { 230 for (int j = 0; j < nx; j++) { 231 X[k][j] /= Rdiag[k]; 232 } 233 for (int i = 0; i < k; i++) { 234 for (int j = 0; j < nx; j++) { 235 X[i][j] -= X[k][j]*QR[i][k]; 236 } 237 } 238 } 239 return (new Matrix(X,n,nx).getMatrix(0,n-1,0,nx-1)); 240 } 241}