001/* 002 * To change this template, choose Tools | Templates 003 * and open the template in the editor. 004 */ 005package org.jblas; 006 007import org.jblas.exceptions.LapackArgumentException; 008import org.jblas.exceptions.LapackPositivityException; 009import org.jblas.util.Permutations; 010import static org.jblas.util.Functions.min; 011 012/** 013 * Matrix which collects all kinds of decompositions. 014 */ 015public class Decompose { 016 /** 017 * Class to hold an LU decomposition result. 018 * 019 * Contains a lower matrix L, and upper matrix U, and a permutation matrix 020 * P such that P*L*U is the original matrix. 021 * @param <T> 022 */ 023 public static class LUDecomposition<T> { 024 025 public T l; 026 public T u; 027 public T p; 028 029 public LUDecomposition(T l, T u, T p) { 030 this.l = l; 031 this.u = u; 032 this.p = p; 033 } 034 035 @Override 036 public String toString() { 037 return String.format("<LUDecomposition L=%s U=%s P=%s>", l, u, p); 038 } 039 } 040 041 /** 042 * Compute LU Decomposition of a general matrix. 043 * 044 * Computes the LU decomposition using GETRF. Returns three matrices L, U, P, 045 * where L is lower diagonal, U is upper diagonal, and P is a permutation 046 * matrix such that A = P * L * U. 047 * 048 * @param A general matrix 049 * @return An LUDecomposition object. 050 */ 051 public static LUDecomposition<DoubleMatrix> lu(DoubleMatrix A) { 052 int[] ipiv = new int[min(A.rows, A.columns)]; 053 DoubleMatrix result = A.dup(); 054 NativeBlas.dgetrf(A.rows, A.columns, result.data, 0, A.rows, ipiv, 0); 055 056 // collect result 057 DoubleMatrix l = new DoubleMatrix(A.rows, min(A.rows, A.columns)); 058 DoubleMatrix u = new DoubleMatrix(min(A.columns, A.rows), A.columns); 059 decomposeLowerUpper(result, l, u); 060 DoubleMatrix p = Permutations.permutationDoubleMatrixFromPivotIndices(A.rows, ipiv); 061 return new LUDecomposition<DoubleMatrix>(l, u, p); 062 } 063 064 private static void decomposeLowerUpper(DoubleMatrix A, DoubleMatrix L, DoubleMatrix U) { 065 for (int i = 0; i < A.rows; i++) { 066 for (int j = 0; j < A.columns; j++) { 067 if (i < j) { 068 U.put(i, j, A.get(i, j)); 069 } else if (i == j) { 070 U.put(i, i, A.get(i, i)); 071 L.put(i, i, 1.0); 072 } else { 073 L.put(i, j, A.get(i, j)); 074 } 075 076 } 077 } 078 } 079 080 /**if (info ) 081 * Compute Cholesky decomposition of A 082 * 083 * @param A symmetric, positive definite matrix (only upper half is used) 084 * @return upper triangular matrix U such that A = U' * U 085 */ 086 public static FloatMatrix cholesky(FloatMatrix A) { 087 FloatMatrix result = A.dup(); 088 int info = NativeBlas.spotrf('U', A.rows, result.data, 0, A.rows); 089 if (info < 0) { 090 throw new LapackArgumentException("DPOTRF", -info); 091 } else if (info > 0) { 092 throw new LapackPositivityException("DPOTRF", "Minor " + info + " was negative. Matrix must be positive definite."); 093 } 094 clearLower(result); 095 return result; 096 } 097 098 private static void clearLower(FloatMatrix A) { 099 for (int j = 0; j < A.columns; j++) 100 for (int i = j + 1; i < A.rows; i++) 101 A.put(i, j, 0.0f); 102 } 103 104 /** 105 * Compute LU Decomposition of a general matrix. 106 * 107 * Computes the LU decomposition using GETRF. Returns three matrices L, U, P, 108 * where L is lower diagonal, U is upper diagonal, and P is a permutation 109 * matrix such that A = P * L * U. 110 * 111 * @param A general matrix 112 * @return An LUDecomposition object. 113 */ 114 public static LUDecomposition<FloatMatrix> lu(FloatMatrix A) { 115 int[] ipiv = new int[min(A.rows, A.columns)]; 116 FloatMatrix result = A.dup(); 117 NativeBlas.sgetrf(A.rows, A.columns, result.data, 0, A.rows, ipiv, 0); 118 119 // collect result 120 FloatMatrix l = new FloatMatrix(A.rows, min(A.rows, A.columns)); 121 FloatMatrix u = new FloatMatrix(min(A.columns, A.rows), A.columns); 122 decomposeLowerUpper(result, l, u); 123 FloatMatrix p = Permutations.permutationFloatMatrixFromPivotIndices(A.rows, ipiv); 124 return new LUDecomposition<FloatMatrix>(l, u, p); 125 } 126 127 private static void decomposeLowerUpper(FloatMatrix A, FloatMatrix L, FloatMatrix U) { 128 for (int i = 0; i < A.rows; i++) { 129 for (int j = 0; j < A.columns; j++) { 130 if (i < j) { 131 U.put(i, j, A.get(i, j)); 132 } else if (i == j) { 133 U.put(i, i, A.get(i, i)); 134 L.put(i, i, 1.0f); 135 } else { 136 L.put(i, j, A.get(i, j)); 137 } 138 139 } 140 } 141 } 142 143 /** 144 * Compute Cholesky decomposition of A 145 * 146 * @param A symmetric, positive definite matrix (only upper half is used) 147 * @return upper triangular matrix U such that A = U' * U 148 */ 149 public static DoubleMatrix cholesky(DoubleMatrix A) { 150 DoubleMatrix result = A.dup(); 151 int info = NativeBlas.dpotrf('U', A.rows, result.data, 0, A.rows); 152 if (info < 0) { 153 throw new LapackArgumentException("DPOTRF", -info); 154 } else if (info > 0) { 155 throw new LapackPositivityException("DPOTRF", "Minor " + info + " was negative. Matrix must be positive definite."); 156 } 157 clearLower(result); 158 return result; 159 } 160 161 private static void clearLower(DoubleMatrix A) { 162 for (int j = 0; j < A.columns; j++) 163 for (int i = j + 1; i < A.rows; i++) 164 A.put(i, j, 0.0); 165 } 166 167 /** 168 * Class to represent a QR decomposition. 169 * 170 * @param <T> 171 */ 172 public static class QRDecomposition<T> { 173 public T q; 174 public T r; 175 176 QRDecomposition(T q, T r) { 177 this.q = q; 178 this.r = r; 179 } 180 181 @Override 182 public String toString() { 183 return "<Q=" + q + " R=" + r + ">"; 184 } 185 } 186 187 /** 188 * QR decomposition. 189 * 190 * Decomposes (m,n) matrix A into a (m,m) matrix Q and an (m,n) matrix R such that 191 * Q is orthogonal, R is upper triangular and Q * R = A 192 * 193 * Note that if A has more rows than columns, then the lower rows of R will contain 194 * only zeros, such that the corresponding later columns of Q do not enter the computation 195 * at all. For some reason, LAPACK does not properly normalize those columns. 196 * 197 * @param A matrix 198 * @return QR decomposition 199 */ 200 public static QRDecomposition<DoubleMatrix> qr(DoubleMatrix A) { 201 int minmn = min(A.rows, A.columns); 202 DoubleMatrix result = A.dup(); 203 DoubleMatrix tau = new DoubleMatrix(minmn); 204 SimpleBlas.geqrf(result, tau); 205 DoubleMatrix R = new DoubleMatrix(A.rows, A.columns); 206 for (int i = 0; i < A.rows; i++) { 207 for (int j = i; j < A.columns; j++) { 208 R.put(i, j, result.get(i, j)); 209 } 210 } 211 DoubleMatrix Q = DoubleMatrix.eye(A.rows); 212 SimpleBlas.ormqr('L', 'N', result, tau, Q); 213 return new QRDecomposition<DoubleMatrix>(Q, R); 214 } 215 216 /** 217 * QR decomposition. 218 * 219 * Decomposes (m,n) matrix A into a (m,m) matrix Q and an (m,n) matrix R such that 220 * Q is orthogonal, R is upper triangular and Q * R = A 221 * 222 * @param A matrix 223 * @return QR decomposition 224 */ 225 public static QRDecomposition<FloatMatrix> qr(FloatMatrix A) { 226 int minmn = min(A.rows, A.columns); 227 FloatMatrix result = A.dup(); 228 FloatMatrix tau = new FloatMatrix(minmn); 229 SimpleBlas.geqrf(result, tau); 230 FloatMatrix R = new FloatMatrix(A.rows, A.columns); 231 for (int i = 0; i < A.rows; i++) { 232 for (int j = i; j < A.columns; j++) { 233 R.put(i, j, result.get(i, j)); 234 } 235 } 236 FloatMatrix Q = FloatMatrix.eye(A.rows); 237 SimpleBlas.ormqr('L', 'N', result, tau, Q); 238 return new QRDecomposition<FloatMatrix>(Q, R); 239 } 240}