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 /** Cholesky Decomposition. 024 <P> 025 For a symmetric, positive definite matrix A, the Cholesky decomposition 026 is an lower triangular matrix L so that A = L*L'. 027 <P> 028 If the matrix is not symmetric or positive definite, the constructor 029 returns a partial decomposition and sets an internal flag that may 030 be queried by the isSPD() method. 031 */ 032 033public class CholeskyDecomposition implements java.io.Serializable { 034 035 static final long serialVersionUID = 224348942390823l; 036 037/* ------------------------ 038 Class variables 039 * ------------------------ */ 040 041 /** Array for internal storage of decomposition. 042 @serial internal array storage. 043 */ 044 private double[][] L; 045 046 /** Row and column dimension (square matrix). 047 @serial matrix dimension. 048 */ 049 private int n; 050 051 /** Symmetric and positive definite flag. 052 @serial is symmetric and positive definite flag. 053 */ 054 private boolean isspd; 055 056/* ------------------------ 057 Constructor 058 * ------------------------ */ 059 060 /** Cholesky algorithm for symmetric and positive definite matrix. 061 @param Arg Square, symmetric matrix. 062 */ 063 064 public CholeskyDecomposition (Matrix Arg) { 065 066 067 // Initialize. 068 double[][] A = Arg.getArray(); 069 n = Arg.getRowDimension(); 070 L = new double[n][n]; 071 isspd = (Arg.getColumnDimension() == n); 072 // Main loop. 073 for (int j = 0; j < n; j++) { 074 double[] Lrowj = L[j]; 075 double d = 0.0; 076 for (int k = 0; k < j; k++) { 077 double[] Lrowk = L[k]; 078 double s = 0.0; 079 for (int i = 0; i < k; i++) { 080 s += Lrowk[i]*Lrowj[i]; 081 } 082 Lrowj[k] = s = (A[j][k] - s)/L[k][k]; 083 d = d + s*s; 084 isspd = isspd && (A[k][j] == A[j][k]); 085 } 086 d = A[j][j] - d; 087 isspd = isspd && (d > 0.0); 088 L[j][j] = Math.sqrt(Math.max(d,0.0)); 089 for (int k = j+1; k < n; k++) { 090 L[j][k] = 0.0; 091 } 092 } 093 } 094 095/* ------------------------ 096 Temporary, experimental code. 097 * ------------------------ *\ 098 099 \** Right Triangular Cholesky Decomposition. 100 <P> 101 For a symmetric, positive definite matrix A, the Right Cholesky 102 decomposition is an upper triangular matrix R so that A = R'*R. 103 This constructor computes R with the Fortran inspired column oriented 104 algorithm used in LINPACK and MATLAB. In Java, we suspect a row oriented, 105 lower triangular decomposition is faster. We have temporarily included 106 this constructor here until timing experiments confirm this suspicion. 107 *\ 108 109 \** Array for internal storage of right triangular decomposition. **\ 110 private transient double[][] R; 111 112 \** Cholesky algorithm for symmetric and positive definite matrix. 113 @param A Square, symmetric matrix. 114 @param rightflag Actual value ignored. 115 @return Structure to access R and isspd flag. 116 *\ 117 118 public CholeskyDecomposition (Matrix Arg, int rightflag) { 119 // Initialize. 120 double[][] A = Arg.getArray(); 121 n = Arg.getColumnDimension(); 122 R = new double[n][n]; 123 isspd = (Arg.getColumnDimension() == n); 124 // Main loop. 125 for (int j = 0; j < n; j++) { 126 double d = 0.0; 127 for (int k = 0; k < j; k++) { 128 double s = A[k][j]; 129 for (int i = 0; i < k; i++) { 130 s = s - R[i][k]*R[i][j]; 131 } 132 R[k][j] = s = s/R[k][k]; 133 d = d + s*s; 134 isspd = isspd & (A[k][j] == A[j][k]); 135 } 136 d = A[j][j] - d; 137 isspd = isspd & (d > 0.0); 138 R[j][j] = Math.sqrt(Math.max(d,0.0)); 139 for (int k = j+1; k < n; k++) { 140 R[k][j] = 0.0; 141 } 142 } 143 } 144 145 \** Return upper triangular factor. 146 @return R 147 *\ 148 149 public Matrix getR () { 150 return new Matrix(R,n,n); 151 } 152 153\* ------------------------ 154 End of temporary code. 155 * ------------------------ */ 156 157/* ------------------------ 158 Public Methods 159 * ------------------------ */ 160 161 /** Is the matrix symmetric and positive definite? 162 @return true if A is symmetric and positive definite. 163 */ 164 165 public boolean isSPD () { 166 return isspd; 167 } 168 169 /** Return triangular factor. 170 @return L 171 */ 172 173 public Matrix getL () { 174 return new Matrix(L,n,n); 175 } 176 177 /** Solve A*X = B 178 @param B A Matrix with as many rows as A and any number of columns. 179 @return X so that L*L'*X = B 180 @exception IllegalArgumentException Matrix row dimensions must agree. 181 @exception RuntimeException Matrix is not symmetric positive definite. 182 */ 183 184 public Matrix solve (Matrix B) { 185 if (B.getRowDimension() != n) { 186 throw new IllegalArgumentException("Matrix row dimensions must agree."); 187 } 188 if (!isspd) { 189 throw new RuntimeException("Matrix is not symmetric positive definite."); 190 } 191 192 // Copy right hand side. 193 double[][] X = B.getArrayCopy(); 194 int nx = B.getColumnDimension(); 195 196 // Solve L*Y = B; 197 for (int k = 0; k < n; k++) { 198 for (int j = 0; j < nx; j++) { 199 for (int i = 0; i < k ; i++) { 200 X[k][j] -= X[i][j]*L[k][i]; 201 } 202 X[k][j] /= L[k][k]; 203 } 204 } 205 206 // Solve L'*X = Y; 207 for (int k = n-1; k >= 0; k--) { 208 for (int j = 0; j < nx; j++) { 209 for (int i = k+1; i < n ; i++) { 210 X[k][j] -= X[i][j]*L[i][k]; 211 } 212 X[k][j] /= L[k][k]; 213 } 214 } 215 216 217 return new Matrix(X,n,nx); 218 } 219} 220