[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
vigra/linear_solve.hxx | ![]() |
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 2003-2008 by Gunnar Kedenburg and Ullrich Koethe */ 00004 /* */ 00005 /* This file is part of the VIGRA computer vision library. */ 00006 /* The VIGRA Website is */ 00007 /* http://hci.iwr.uni-heidelberg.de/vigra/ */ 00008 /* Please direct questions, bug reports, and contributions to */ 00009 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00010 /* vigra@informatik.uni-hamburg.de */ 00011 /* */ 00012 /* Permission is hereby granted, free of charge, to any person */ 00013 /* obtaining a copy of this software and associated documentation */ 00014 /* files (the "Software"), to deal in the Software without */ 00015 /* restriction, including without limitation the rights to use, */ 00016 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00017 /* sell copies of the Software, and to permit persons to whom the */ 00018 /* Software is furnished to do so, subject to the following */ 00019 /* conditions: */ 00020 /* */ 00021 /* The above copyright notice and this permission notice shall be */ 00022 /* included in all copies or substantial portions of the */ 00023 /* Software. */ 00024 /* */ 00025 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00026 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00027 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00028 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00029 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00030 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00031 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00032 /* OTHER DEALINGS IN THE SOFTWARE. */ 00033 /* */ 00034 /************************************************************************/ 00035 00036 00037 #ifndef VIGRA_LINEAR_SOLVE_HXX 00038 #define VIGRA_LINEAR_SOLVE_HXX 00039 00040 #include <ctype.h> 00041 #include <string> 00042 #include "mathutil.hxx" 00043 #include "matrix.hxx" 00044 #include "singular_value_decomposition.hxx" 00045 00046 00047 namespace vigra 00048 { 00049 00050 namespace linalg 00051 { 00052 00053 namespace detail { 00054 00055 template <class T, class C1> 00056 T determinantByLUDecomposition(MultiArrayView<2, T, C1> const & a) 00057 { 00058 typedef MultiArrayShape<2>::type Shape; 00059 00060 MultiArrayIndex m = rowCount(a), n = columnCount(a); 00061 vigra_precondition(n == m, 00062 "determinant(): square matrix required."); 00063 00064 Matrix<T> LU(a); 00065 T det = 1.0; 00066 00067 for (MultiArrayIndex j = 0; j < n; ++j) 00068 { 00069 // Apply previous transformations. 00070 for (MultiArrayIndex i = 0; i < m; ++i) 00071 { 00072 MultiArrayIndex end = std::min(i, j); 00073 T s = dot(rowVector(LU, Shape(i,0), end), columnVector(LU, Shape(0,j), end)); 00074 LU(i,j) = LU(i,j) -= s; 00075 } 00076 00077 // Find pivot and exchange if necessary. 00078 MultiArrayIndex p = j + argMax(abs(columnVector(LU, Shape(j,j), m))); 00079 if (p != j) 00080 { 00081 rowVector(LU, p).swapData(rowVector(LU, j)); 00082 det = -det; 00083 } 00084 00085 det *= LU(j,j); 00086 00087 // Compute multipliers. 00088 if (LU(j,j) != 0.0) 00089 columnVector(LU, Shape(j+1,j), m) /= LU(j,j); 00090 else 00091 break; // det is zero 00092 } 00093 return det; 00094 } 00095 00096 // returns the new value of 'a' (when this Givens rotation is applied to 'a' and 'b') 00097 // the new value of 'b' is zero, of course 00098 template <class T> 00099 T givensCoefficients(T a, T b, T & c, T & s) 00100 { 00101 if(abs(a) < abs(b)) 00102 { 00103 T t = a/b, 00104 r = std::sqrt(1.0 + t*t); 00105 s = 1.0 / r; 00106 c = t*s; 00107 return r*b; 00108 } 00109 else if(a != 0.0) 00110 { 00111 T t = b/a, 00112 r = std::sqrt(1.0 + t*t); 00113 c = 1.0 / r; 00114 s = t*c; 00115 return r*a; 00116 } 00117 else // a == b == 0.0 00118 { 00119 c = 1.0; 00120 s = 0.0; 00121 return 0.0; 00122 } 00123 } 00124 00125 // see Golub, van Loan: Algorithm 5.1.3 (p. 216) 00126 template <class T> 00127 bool givensRotationMatrix(T a, T b, Matrix<T> & gTranspose) 00128 { 00129 if(b == 0.0) 00130 return false; // no rotation needed 00131 givensCoefficients(a, b, gTranspose(0,0), gTranspose(0,1)); 00132 gTranspose(1,1) = gTranspose(0,0); 00133 gTranspose(1,0) = -gTranspose(0,1); 00134 return true; 00135 } 00136 00137 // reflections are symmetric matrices and can thus be applied to rows 00138 // and columns in the same way => code simplification relative to rotations 00139 template <class T> 00140 inline bool 00141 givensReflectionMatrix(T a, T b, Matrix<T> & g) 00142 { 00143 if(b == 0.0) 00144 return false; // no reflection needed 00145 givensCoefficients(a, b, g(0,0), g(0,1)); 00146 g(1,1) = -g(0,0); 00147 g(1,0) = g(0,1); 00148 return true; 00149 } 00150 00151 // see Golub, van Loan: Algorithm 5.2.2 (p. 227) and Section 12.5.2 (p. 608) 00152 template <class T, class C1, class C2> 00153 bool 00154 qrGivensStepImpl(MultiArrayIndex i, MultiArrayView<2, T, C1> r, MultiArrayView<2, T, C2> rhs) 00155 { 00156 typedef typename Matrix<T>::difference_type Shape; 00157 00158 const MultiArrayIndex m = rowCount(r); 00159 const MultiArrayIndex n = columnCount(r); 00160 const MultiArrayIndex rhsCount = columnCount(rhs); 00161 vigra_precondition(m == rowCount(rhs), 00162 "qrGivensStepImpl(): Matrix shape mismatch."); 00163 00164 Matrix<T> givens(2,2); 00165 for(int k=m-1; k>(int)i; --k) 00166 { 00167 if(!givensReflectionMatrix(r(k-1,i), r(k,i), givens)) 00168 continue; // r(k,i) was already zero 00169 00170 r(k-1,i) = givens(0,0)*r(k-1,i) + givens(0,1)*r(k,i); 00171 r(k,i) = 0.0; 00172 00173 r.subarray(Shape(k-1,i+1), Shape(k+1,n)) = givens*r.subarray(Shape(k-1,i+1), Shape(k+1,n)); 00174 rhs.subarray(Shape(k-1,0), Shape(k+1,rhsCount)) = givens*rhs.subarray(Shape(k-1,0), Shape(k+1,rhsCount)); 00175 } 00176 return r(i,i) != 0.0; 00177 } 00178 00179 // see Golub, van Loan: Section 12.5.2 (p. 608) 00180 template <class T, class C1, class C2, class Permutation> 00181 void 00182 upperTriangularCyclicShiftColumns(MultiArrayIndex i, MultiArrayIndex j, 00183 MultiArrayView<2, T, C1> &r, MultiArrayView<2, T, C2> &rhs, Permutation & permutation) 00184 { 00185 typedef typename Matrix<T>::difference_type Shape; 00186 00187 const MultiArrayIndex m = rowCount(r); 00188 const MultiArrayIndex n = columnCount(r); 00189 const MultiArrayIndex rhsCount = columnCount(rhs); 00190 vigra_precondition(i < n && j < n, 00191 "upperTriangularCyclicShiftColumns(): Shift indices out of range."); 00192 vigra_precondition(m == rowCount(rhs), 00193 "upperTriangularCyclicShiftColumns(): Matrix shape mismatch."); 00194 00195 if(j == i) 00196 return; 00197 if(j < i) 00198 std::swap(j,i); 00199 00200 Matrix<T> t = columnVector(r, i); 00201 MultiArrayIndex ti = permutation[i]; 00202 for(MultiArrayIndex k=i; k<j;++k) 00203 { 00204 columnVector(r, k) = columnVector(r, k+1); 00205 permutation[k] = permutation[k+1]; 00206 } 00207 columnVector(r, j) = t; 00208 permutation[j] = ti; 00209 00210 Matrix<T> givens(2,2); 00211 for(MultiArrayIndex k=i; k<j; ++k) 00212 { 00213 if(!givensReflectionMatrix(r(k,k), r(k+1,k), givens)) 00214 continue; // r(k+1,k) was already zero 00215 00216 r(k,k) = givens(0,0)*r(k,k) + givens(0,1)*r(k+1,k); 00217 r(k+1,k) = 0.0; 00218 00219 r.subarray(Shape(k,k+1), Shape(k+2,n)) = givens*r.subarray(Shape(k,k+1), Shape(k+2,n)); 00220 rhs.subarray(Shape(k,0), Shape(k+2,rhsCount)) = givens*rhs.subarray(Shape(k,0), Shape(k+2,rhsCount)); 00221 } 00222 } 00223 00224 // see Golub, van Loan: Section 12.5.2 (p. 608) 00225 template <class T, class C1, class C2, class Permutation> 00226 void 00227 upperTriangularSwapColumns(MultiArrayIndex i, MultiArrayIndex j, 00228 MultiArrayView<2, T, C1> &r, MultiArrayView<2, T, C2> &rhs, Permutation & permutation) 00229 { 00230 typedef typename Matrix<T>::difference_type Shape; 00231 00232 const MultiArrayIndex m = rowCount(r); 00233 const MultiArrayIndex n = columnCount(r); 00234 const MultiArrayIndex rhsCount = columnCount(rhs); 00235 vigra_precondition(i < n && j < n, 00236 "upperTriangularSwapColumns(): Swap indices out of range."); 00237 vigra_precondition(m == rowCount(rhs), 00238 "upperTriangularSwapColumns(): Matrix shape mismatch."); 00239 00240 if(j == i) 00241 return; 00242 if(j < i) 00243 std::swap(j,i); 00244 00245 columnVector(r, i).swapData(columnVector(r, j)); 00246 std::swap(permutation[i], permutation[j]); 00247 00248 Matrix<T> givens(2,2); 00249 for(int k=m-1; k>(int)i; --k) 00250 { 00251 if(!givensReflectionMatrix(r(k-1,i), r(k,i), givens)) 00252 continue; // r(k,i) was already zero 00253 00254 r(k-1,i) = givens(0,0)*r(k-1,i) + givens(0,1)*r(k,i); 00255 r(k,i) = 0.0; 00256 00257 r.subarray(Shape(k-1,i+1), Shape(k+1,n)) = givens*r.subarray(Shape(k-1,i+1), Shape(k+1,n)); 00258 rhs.subarray(Shape(k-1,0), Shape(k+1,rhsCount)) = givens*rhs.subarray(Shape(k-1,0), Shape(k+1,rhsCount)); 00259 } 00260 MultiArrayIndex end = std::min(j, m-1); 00261 for(MultiArrayIndex k=i+1; k<end; ++k) 00262 { 00263 if(!givensReflectionMatrix(r(k,k), r(k+1,k), givens)) 00264 continue; // r(k+1,k) was already zero 00265 00266 r(k,k) = givens(0,0)*r(k,k) + givens(0,1)*r(k+1,k); 00267 r(k+1,k) = 0.0; 00268 00269 r.subarray(Shape(k,k+1), Shape(k+2,n)) = givens*r.subarray(Shape(k,k+1), Shape(k+2,n)); 00270 rhs.subarray(Shape(k,0), Shape(k+2,rhsCount)) = givens*rhs.subarray(Shape(k,0), Shape(k+2,rhsCount)); 00271 } 00272 } 00273 00274 // see Lawson & Hanson: Algorithm H1 (p. 57) 00275 template <class T, class C1, class C2, class U> 00276 bool householderVector(MultiArrayView<2, T, C1> const & v, MultiArrayView<2, T, C2> & u, U & vnorm) 00277 { 00278 vnorm = (v(0,0) > 0.0) 00279 ? -norm(v) 00280 : norm(v); 00281 U f = std::sqrt(vnorm*(vnorm - v(0,0))); 00282 00283 if(f == NumericTraits<U>::zero()) 00284 { 00285 u.init(NumericTraits<T>::zero()); 00286 return false; 00287 } 00288 else 00289 { 00290 u(0,0) = (v(0,0) - vnorm) / f; 00291 for(MultiArrayIndex k=1; k<rowCount(u); ++k) 00292 u(k,0) = v(k,0) / f; 00293 return true; 00294 } 00295 } 00296 00297 // see Lawson & Hanson: Algorithm H1 (p. 57) 00298 template <class T, class C1, class C2, class C3> 00299 bool 00300 qrHouseholderStepImpl(MultiArrayIndex i, MultiArrayView<2, T, C1> & r, 00301 MultiArrayView<2, T, C2> & rhs, MultiArrayView<2, T, C3> & householderMatrix) 00302 { 00303 typedef typename Matrix<T>::difference_type Shape; 00304 00305 const MultiArrayIndex m = rowCount(r); 00306 const MultiArrayIndex n = columnCount(r); 00307 const MultiArrayIndex rhsCount = columnCount(rhs); 00308 00309 vigra_precondition(i < n && i < m, 00310 "qrHouseholderStepImpl(): Index i out of range."); 00311 00312 Matrix<T> u(m-i,1); 00313 T vnorm; 00314 bool nontrivial = householderVector(columnVector(r, Shape(i,i), m), u, vnorm); 00315 00316 r(i,i) = vnorm; 00317 columnVector(r, Shape(i+1,i), m).init(NumericTraits<T>::zero()); 00318 00319 if(columnCount(householderMatrix) == n) 00320 columnVector(householderMatrix, Shape(i,i), m) = u; 00321 00322 if(nontrivial) 00323 { 00324 for(MultiArrayIndex k=i+1; k<n; ++k) 00325 columnVector(r, Shape(i,k), m) -= dot(columnVector(r, Shape(i,k), m), u) * u; 00326 for(MultiArrayIndex k=0; k<rhsCount; ++k) 00327 columnVector(rhs, Shape(i,k), m) -= dot(columnVector(rhs, Shape(i,k), m), u) * u; 00328 } 00329 return r(i,i) != 0.0; 00330 } 00331 00332 template <class T, class C1, class C2> 00333 bool 00334 qrColumnHouseholderStep(MultiArrayIndex i, MultiArrayView<2, T, C1> &r, MultiArrayView<2, T, C2> &rhs) 00335 { 00336 Matrix<T> dontStoreHouseholderVectors; // intentionally empty 00337 return qrHouseholderStepImpl(i, r, rhs, dontStoreHouseholderVectors); 00338 } 00339 00340 template <class T, class C1, class C2> 00341 bool 00342 qrRowHouseholderStep(MultiArrayIndex i, MultiArrayView<2, T, C1> &r, MultiArrayView<2, T, C2> & householderMatrix) 00343 { 00344 Matrix<T> dontTransformRHS; // intentionally empty 00345 MultiArrayView<2, T, StridedArrayTag> rt = transpose(r), 00346 ht = transpose(householderMatrix); 00347 return qrHouseholderStepImpl(i, rt, dontTransformRHS, ht); 00348 } 00349 00350 // O(n) algorithm due to Bischof: Incremental Condition Estimation, 1990 00351 template <class T, class C1, class C2, class SNType> 00352 void 00353 incrementalMaxSingularValueApproximation(MultiArrayView<2, T, C1> const & newColumn, 00354 MultiArrayView<2, T, C2> & z, SNType & v) 00355 { 00356 typedef typename Matrix<T>::difference_type Shape; 00357 MultiArrayIndex n = rowCount(newColumn) - 1; 00358 00359 SNType vneu = squaredNorm(newColumn); 00360 T yv = dot(columnVector(newColumn, Shape(0,0),n), columnVector(z, Shape(0,0),n)); 00361 // use atan2 as it is robust against overflow/underflow 00362 T t = 0.5*std::atan2(T(2.0*yv), T(sq(v)-vneu)), 00363 s = std::sin(t), 00364 c = std::cos(t); 00365 v = std::sqrt(sq(c*v) + sq(s)*vneu + 2.0*s*c*yv); 00366 columnVector(z, Shape(0,0),n) = c*columnVector(z, Shape(0,0),n) + s*columnVector(newColumn, Shape(0,0),n); 00367 z(n,0) = s*newColumn(n,0); 00368 } 00369 00370 // O(n) algorithm due to Bischof: Incremental Condition Estimation, 1990 00371 template <class T, class C1, class C2, class SNType> 00372 void 00373 incrementalMinSingularValueApproximation(MultiArrayView<2, T, C1> const & newColumn, 00374 MultiArrayView<2, T, C2> & z, SNType & v, double tolerance) 00375 { 00376 typedef typename Matrix<T>::difference_type Shape; 00377 00378 if(v <= tolerance) 00379 { 00380 v = 0.0; 00381 return; 00382 } 00383 00384 MultiArrayIndex n = rowCount(newColumn) - 1; 00385 00386 T gamma = newColumn(n,0); 00387 if(gamma == 0.0) 00388 { 00389 v = 0.0; 00390 return; 00391 } 00392 00393 T yv = dot(columnVector(newColumn, Shape(0,0),n), columnVector(z, Shape(0,0),n)); 00394 // use atan2 as it is robust against overflow/underflow 00395 T t = 0.5*std::atan2(T(-2.0*yv), T(squaredNorm(gamma / v) + squaredNorm(yv) - 1.0)), 00396 s = std::sin(t), 00397 c = std::cos(t); 00398 columnVector(z, Shape(0,0),n) *= c; 00399 z(n,0) = (s - c*yv) / gamma; 00400 v *= norm(gamma) / hypot(c*gamma, v*(s - c*yv)); 00401 } 00402 00403 // QR algorithm with optional column pivoting 00404 template <class T, class C1, class C2, class C3> 00405 unsigned int 00406 qrTransformToTriangularImpl(MultiArrayView<2, T, C1> & r, MultiArrayView<2, T, C2> & rhs, MultiArrayView<2, T, C3> & householder, 00407 ArrayVector<MultiArrayIndex> & permutation, double epsilon) 00408 { 00409 typedef typename Matrix<T>::difference_type Shape; 00410 typedef typename NormTraits<MultiArrayView<2, T, C1> >::NormType NormType; 00411 typedef typename NormTraits<MultiArrayView<2, T, C1> >::SquaredNormType SNType; 00412 00413 const MultiArrayIndex m = rowCount(r); 00414 const MultiArrayIndex n = columnCount(r); 00415 const MultiArrayIndex maxRank = std::min(m, n); 00416 00417 vigra_precondition(m >= n, 00418 "qrTransformToTriangularImpl(): Coefficient matrix with at least as many rows as columns required."); 00419 00420 const MultiArrayIndex rhsCount = columnCount(rhs); 00421 bool transformRHS = rhsCount > 0; 00422 vigra_precondition(!transformRHS || m == rowCount(rhs), 00423 "qrTransformToTriangularImpl(): RHS matrix shape mismatch."); 00424 00425 bool storeHouseholderSteps = columnCount(householder) > 0; 00426 vigra_precondition(!storeHouseholderSteps || r.shape() == householder.shape(), 00427 "qrTransformToTriangularImpl(): Householder matrix shape mismatch."); 00428 00429 bool pivoting = permutation.size() > 0; 00430 vigra_precondition(!pivoting || n == (MultiArrayIndex)permutation.size(), 00431 "qrTransformToTriangularImpl(): Permutation array size mismatch."); 00432 00433 if(n == 0) 00434 return 0; // trivial solution 00435 00436 Matrix<SNType> columnSquaredNorms; 00437 if(pivoting) 00438 { 00439 columnSquaredNorms.reshape(Shape(1,n)); 00440 for(MultiArrayIndex k=0; k<n; ++k) 00441 columnSquaredNorms[k] = squaredNorm(columnVector(r, k)); 00442 00443 int pivot = argMax(columnSquaredNorms); 00444 if(pivot != 0) 00445 { 00446 columnVector(r, 0).swapData(columnVector(r, pivot)); 00447 std::swap(columnSquaredNorms[0], columnSquaredNorms[pivot]); 00448 std::swap(permutation[0], permutation[pivot]); 00449 } 00450 } 00451 00452 qrHouseholderStepImpl(0, r, rhs, householder); 00453 00454 MultiArrayIndex rank = 1; 00455 NormType maxApproxSingularValue = norm(r(0,0)), 00456 minApproxSingularValue = maxApproxSingularValue; 00457 00458 double tolerance = (epsilon == 0.0) 00459 ? m*maxApproxSingularValue*NumericTraits<T>::epsilon() 00460 : epsilon; 00461 00462 bool simpleSingularValueApproximation = (n < 4); 00463 Matrix<T> zmax, zmin; 00464 if(minApproxSingularValue <= tolerance) 00465 { 00466 rank = 0; 00467 pivoting = false; 00468 simpleSingularValueApproximation = true; 00469 } 00470 if(!simpleSingularValueApproximation) 00471 { 00472 zmax.reshape(Shape(m,1)); 00473 zmin.reshape(Shape(m,1)); 00474 zmax(0,0) = r(0,0); 00475 zmin(0,0) = 1.0 / r(0,0); 00476 } 00477 00478 for(MultiArrayIndex k=1; k<maxRank; ++k) 00479 { 00480 if(pivoting) 00481 { 00482 for(MultiArrayIndex l=k; l<n; ++l) 00483 columnSquaredNorms[l] -= squaredNorm(r(k, l)); 00484 int pivot = k + argMax(rowVector(columnSquaredNorms, Shape(0,k), n)); 00485 if(pivot != (int)k) 00486 { 00487 columnVector(r, k).swapData(columnVector(r, pivot)); 00488 std::swap(columnSquaredNorms[k], columnSquaredNorms[pivot]); 00489 std::swap(permutation[k], permutation[pivot]); 00490 } 00491 } 00492 00493 qrHouseholderStepImpl(k, r, rhs, householder); 00494 00495 if(simpleSingularValueApproximation) 00496 { 00497 NormType nv = norm(r(k,k)); 00498 maxApproxSingularValue = std::max(nv, maxApproxSingularValue); 00499 minApproxSingularValue = std::min(nv, minApproxSingularValue); 00500 } 00501 else 00502 { 00503 incrementalMaxSingularValueApproximation(columnVector(r, Shape(0,k),k+1), zmax, maxApproxSingularValue); 00504 incrementalMinSingularValueApproximation(columnVector(r, Shape(0,k),k+1), zmin, minApproxSingularValue, tolerance); 00505 } 00506 00507 #if 0 00508 Matrix<T> u(k+1,k+1), s(k+1, 1), v(k+1,k+1); 00509 singularValueDecomposition(r.subarray(Shape(0,0), Shape(k+1,k+1)), u, s, v); 00510 std::cerr << "estimate, svd " << k << ": " << minApproxSingularValue << " " << s(k,0) << "\n"; 00511 #endif 00512 00513 if(epsilon == 0.0) 00514 tolerance = m*maxApproxSingularValue*NumericTraits<T>::epsilon(); 00515 00516 if(minApproxSingularValue > tolerance) 00517 ++rank; 00518 else 00519 pivoting = false; // matrix doesn't have full rank, triangulize the rest without pivoting 00520 } 00521 return (unsigned int)rank; 00522 } 00523 00524 template <class T, class C1, class C2> 00525 unsigned int 00526 qrTransformToUpperTriangular(MultiArrayView<2, T, C1> & r, MultiArrayView<2, T, C2> & rhs, 00527 ArrayVector<MultiArrayIndex> & permutation, double epsilon = 0.0) 00528 { 00529 Matrix<T> dontStoreHouseholderVectors; // intentionally empty 00530 return qrTransformToTriangularImpl(r, rhs, dontStoreHouseholderVectors, permutation, epsilon); 00531 } 00532 00533 // QR algorithm with optional row pivoting 00534 template <class T, class C1, class C2, class C3> 00535 unsigned int 00536 qrTransformToLowerTriangular(MultiArrayView<2, T, C1> & r, MultiArrayView<2, T, C2> & rhs, MultiArrayView<2, T, C3> & householderMatrix, 00537 double epsilon = 0.0) 00538 { 00539 ArrayVector<MultiArrayIndex> permutation((unsigned int)rowCount(rhs)); 00540 for(MultiArrayIndex k=0; k<(MultiArrayIndex)permutation.size(); ++k) 00541 permutation[k] = k; 00542 Matrix<T> dontTransformRHS; // intentionally empty 00543 MultiArrayView<2, T, StridedArrayTag> rt = transpose(r), 00544 ht = transpose(householderMatrix); 00545 unsigned int rank = qrTransformToTriangularImpl(rt, dontTransformRHS, ht, permutation, epsilon); 00546 00547 // apply row permutation to RHS 00548 Matrix<T> tempRHS(rhs); 00549 for(MultiArrayIndex k=0; k<(MultiArrayIndex)permutation.size(); ++k) 00550 rowVector(rhs, k) = rowVector(tempRHS, permutation[k]); 00551 return rank; 00552 } 00553 00554 // QR algorithm without column pivoting 00555 template <class T, class C1, class C2> 00556 inline bool 00557 qrTransformToUpperTriangular(MultiArrayView<2, T, C1> & r, MultiArrayView<2, T, C2> & rhs, 00558 double epsilon = 0.0) 00559 { 00560 ArrayVector<MultiArrayIndex> noPivoting; // intentionally empty 00561 00562 return (qrTransformToUpperTriangular(r, rhs, noPivoting, epsilon) == 00563 (unsigned int)columnCount(r)); 00564 } 00565 00566 // QR algorithm without row pivoting 00567 template <class T, class C1, class C2> 00568 inline bool 00569 qrTransformToLowerTriangular(MultiArrayView<2, T, C1> & r, MultiArrayView<2, T, C2> & householder, 00570 double epsilon = 0.0) 00571 { 00572 Matrix<T> noPivoting; // intentionally empty 00573 00574 return (qrTransformToLowerTriangular(r, noPivoting, householder, epsilon) == 00575 (unsigned int)rowCount(r)); 00576 } 00577 00578 // restore ordering of result vector elements after QR solution with column pivoting 00579 template <class T, class C1, class C2, class Permutation> 00580 void inverseRowPermutation(MultiArrayView<2, T, C1> &permuted, MultiArrayView<2, T, C2> &res, 00581 Permutation const & permutation) 00582 { 00583 for(MultiArrayIndex k=0; k<columnCount(permuted); ++k) 00584 for(MultiArrayIndex l=0; l<rowCount(permuted); ++l) 00585 res(permutation[l], k) = permuted(l,k); 00586 } 00587 00588 template <class T, class C1, class C2> 00589 void applyHouseholderColumnReflections(MultiArrayView<2, T, C1> const &householder, MultiArrayView<2, T, C2> &res) 00590 { 00591 typedef typename Matrix<T>::difference_type Shape; 00592 MultiArrayIndex n = rowCount(householder); 00593 MultiArrayIndex m = columnCount(householder); 00594 MultiArrayIndex rhsCount = columnCount(res); 00595 00596 for(int k = m-1; k >= 0; --k) 00597 { 00598 MultiArrayView<2, T, C1> u = columnVector(householder, Shape(k,k), n); 00599 for(MultiArrayIndex l=0; l<rhsCount; ++l) 00600 columnVector(res, Shape(k,l), n) -= dot(columnVector(res, Shape(k,l), n), u) * u; 00601 } 00602 } 00603 00604 } // namespace detail 00605 00606 template <class T, class C1, class C2, class C3> 00607 unsigned int 00608 linearSolveQRReplace(MultiArrayView<2, T, C1> &A, MultiArrayView<2, T, C2> &b, 00609 MultiArrayView<2, T, C3> & res, 00610 double epsilon = 0.0) 00611 { 00612 typedef typename Matrix<T>::difference_type Shape; 00613 00614 MultiArrayIndex n = columnCount(A); 00615 MultiArrayIndex m = rowCount(A); 00616 MultiArrayIndex rhsCount = columnCount(res); 00617 MultiArrayIndex rank = std::min(m,n); 00618 Shape ul(MultiArrayIndex(0), MultiArrayIndex(0)); 00619 00620 00621 vigra_precondition(rhsCount == columnCount(b), 00622 "linearSolveQR(): RHS and solution must have the same number of columns."); 00623 vigra_precondition(m == rowCount(b), 00624 "linearSolveQR(): Coefficient matrix and RHS must have the same number of rows."); 00625 vigra_precondition(n == rowCount(res), 00626 "linearSolveQR(): Mismatch between column count of coefficient matrix and row count of solution."); 00627 vigra_precondition(epsilon >= 0.0, 00628 "linearSolveQR(): 'epsilon' must be non-negative."); 00629 00630 if(m < n) 00631 { 00632 // minimum norm solution of underdetermined system 00633 Matrix<T> householderMatrix(n, m); 00634 MultiArrayView<2, T, StridedArrayTag> ht = transpose(householderMatrix); 00635 rank = (MultiArrayIndex)detail::qrTransformToLowerTriangular(A, b, ht, epsilon); 00636 res.subarray(Shape(rank,0), Shape(n, rhsCount)).init(NumericTraits<T>::zero()); 00637 if(rank < m) 00638 { 00639 // system is also rank-deficient => compute minimum norm least squares solution 00640 MultiArrayView<2, T, C1> Asub = A.subarray(ul, Shape(m,rank)); 00641 detail::qrTransformToUpperTriangular(Asub, b, epsilon); 00642 linearSolveUpperTriangular(A.subarray(ul, Shape(rank,rank)), 00643 b.subarray(ul, Shape(rank,rhsCount)), 00644 res.subarray(ul, Shape(rank, rhsCount))); 00645 } 00646 else 00647 { 00648 // system has full rank => compute minimum norm solution 00649 linearSolveLowerTriangular(A.subarray(ul, Shape(rank,rank)), 00650 b.subarray(ul, Shape(rank, rhsCount)), 00651 res.subarray(ul, Shape(rank, rhsCount))); 00652 } 00653 detail::applyHouseholderColumnReflections(householderMatrix.subarray(ul, Shape(n, rank)), res); 00654 } 00655 else 00656 { 00657 // solution of well-determined or overdetermined system 00658 ArrayVector<MultiArrayIndex> permutation((unsigned int)n); 00659 for(MultiArrayIndex k=0; k<n; ++k) 00660 permutation[k] = k; 00661 00662 rank = (MultiArrayIndex)detail::qrTransformToUpperTriangular(A, b, permutation, epsilon); 00663 00664 Matrix<T> permutedSolution(n, rhsCount); 00665 if(rank < n) 00666 { 00667 // system is rank-deficient => compute minimum norm solution 00668 Matrix<T> householderMatrix(n, rank); 00669 MultiArrayView<2, T, StridedArrayTag> ht = transpose(householderMatrix); 00670 MultiArrayView<2, T, C1> Asub = A.subarray(ul, Shape(rank,n)); 00671 detail::qrTransformToLowerTriangular(Asub, ht, epsilon); 00672 linearSolveLowerTriangular(A.subarray(ul, Shape(rank,rank)), 00673 b.subarray(ul, Shape(rank, rhsCount)), 00674 permutedSolution.subarray(ul, Shape(rank, rhsCount))); 00675 detail::applyHouseholderColumnReflections(householderMatrix, permutedSolution); 00676 } 00677 else 00678 { 00679 // system has full rank => compute exact or least squares solution 00680 linearSolveUpperTriangular(A.subarray(ul, Shape(rank,rank)), 00681 b.subarray(ul, Shape(rank,rhsCount)), 00682 permutedSolution); 00683 } 00684 detail::inverseRowPermutation(permutedSolution, res, permutation); 00685 } 00686 return (unsigned int)rank; 00687 } 00688 00689 template <class T, class C1, class C2, class C3> 00690 unsigned int linearSolveQR(MultiArrayView<2, T, C1> const & A, MultiArrayView<2, T, C2> const & b, 00691 MultiArrayView<2, T, C3> & res) 00692 { 00693 Matrix<T> r(A), rhs(b); 00694 return linearSolveQRReplace(r, rhs, res); 00695 } 00696 00697 /** \defgroup MatrixAlgebra Advanced Matrix Algebra 00698 00699 \brief Solution of linear systems, eigen systems, linear least squares etc. 00700 00701 \ingroup LinearAlgebraModule 00702 */ 00703 //@{ 00704 /** Create the inverse or pseudo-inverse of matrix \a v. 00705 00706 If the matrix \a v is square, \a res must have the same shape and will contain the 00707 inverse of \a v. If \a v is rectangular, \a res must have the transposed shape 00708 of \a v. The inverse is then computed in the least-squares 00709 sense, i.e. \a res will be the pseudo-inverse (Moore-Penrose inverse). 00710 The function returns <tt>true</tt> upon success, and <tt>false</tt> if \a v 00711 is not invertible (has not full rank). The inverse is computed by means of QR 00712 decomposition. This function can be applied in-place. 00713 00714 <b>\#include</b> <vigra/linear_solve.hxx> or<br> 00715 <b>\#include</b> <vigra/linear_algebra.hxx><br> 00716 Namespaces: vigra and vigra::linalg 00717 */ 00718 template <class T, class C1, class C2> 00719 bool inverse(const MultiArrayView<2, T, C1> &v, MultiArrayView<2, T, C2> &res) 00720 { 00721 typedef typename MultiArrayShape<2>::type Shape; 00722 00723 const MultiArrayIndex n = columnCount(v); 00724 const MultiArrayIndex m = rowCount(v); 00725 vigra_precondition(n == rowCount(res) && m == columnCount(res), 00726 "inverse(): shape of output matrix must be the transpose of the input matrix' shape."); 00727 00728 if(m < n) 00729 { 00730 MultiArrayView<2, T, StridedArrayTag> vt = transpose(v); 00731 Matrix<T> r(vt.shape()), q(n, n); 00732 if(!qrDecomposition(vt, q, r)) 00733 return false; // a didn't have full rank 00734 linearSolveUpperTriangular(r.subarray(Shape(0,0), Shape(m,m)), 00735 transpose(q).subarray(Shape(0,0), Shape(m,n)), 00736 transpose(res)); 00737 } 00738 else 00739 { 00740 Matrix<T> r(v.shape()), q(m, m); 00741 if(!qrDecomposition(v, q, r)) 00742 return false; // a didn't have full rank 00743 linearSolveUpperTriangular(r.subarray(Shape(0,0), Shape(n,n)), 00744 transpose(q).subarray(Shape(0,0), Shape(n,m)), 00745 res); 00746 } 00747 return true; 00748 } 00749 00750 /** Create the inverse or pseudo-inverse of matrix \a v. 00751 00752 The result is returned as a temporary matrix. If the matrix \a v is square, 00753 the result will have the same shape and contains the inverse of \a v. 00754 If \a v is rectangular, the result will have the transposed shape of \a v. 00755 The inverse is then computed in the least-squares 00756 sense, i.e. \a res will be the pseudo-inverse (Moore-Penrose inverse). 00757 The inverse is computed by means of QR decomposition. If \a v 00758 is not invertible, <tt>vigra::PreconditionViolation</tt> exception is thrown. 00759 Usage: 00760 00761 \code 00762 vigra::Matrix<double> v(n, n); 00763 v = ...; 00764 00765 vigra::Matrix<double> m = inverse(v); 00766 \endcode 00767 00768 <b>\#include</b> <vigra/linear_solve.hxx> or<br> 00769 <b>\#include</b> <vigra/linear_algebra.hxx><br> 00770 Namespaces: vigra and vigra::linalg 00771 */ 00772 template <class T, class C> 00773 TemporaryMatrix<T> inverse(const MultiArrayView<2, T, C> &v) 00774 { 00775 TemporaryMatrix<T> ret(columnCount(v), rowCount(v)); // transpose shape 00776 vigra_precondition(inverse(v, ret), 00777 "inverse(): matrix is not invertible."); 00778 return ret; 00779 } 00780 00781 template <class T> 00782 TemporaryMatrix<T> inverse(const TemporaryMatrix<T> &v) 00783 { 00784 if(columnCount(v) == rowCount(v)) 00785 { 00786 vigra_precondition(inverse(v, const_cast<TemporaryMatrix<T> &>(v)), 00787 "inverse(): matrix is not invertible."); 00788 return v; 00789 } 00790 else 00791 { 00792 TemporaryMatrix<T> ret(columnCount(v), rowCount(v)); // transpose shape 00793 vigra_precondition(inverse(v, ret), 00794 "inverse(): matrix is not invertible."); 00795 return ret; 00796 } 00797 } 00798 00799 /** Compute the determinant of a square matrix. 00800 00801 \a method must be one of the following: 00802 <DL> 00803 <DT>"Cholesky"<DD> Compute the solution by means of Cholesky decomposition. This 00804 method is faster than "LU", but requires the matrix \a a 00805 to be symmetric positive definite. If this is 00806 not the case, a <tt>ContractViolation</tt> exception is thrown. 00807 00808 <DT>"LU"<DD> (default) Compute the solution by means of LU decomposition. 00809 </DL> 00810 00811 <b>\#include</b> <vigra/linear_solve.hxx> or<br> 00812 <b>\#include</b> <vigra/linear_algebra.hxx><br> 00813 Namespaces: vigra and vigra::linalg 00814 */ 00815 template <class T, class C1> 00816 T determinant(MultiArrayView<2, T, C1> const & a, std::string method = "LU") 00817 { 00818 MultiArrayIndex n = columnCount(a); 00819 vigra_precondition(rowCount(a) == n, 00820 "determinant(): Square matrix required."); 00821 00822 for(unsigned int k=0; k<method.size(); ++k) 00823 method[k] = tolower(method[k]); 00824 00825 if(n == 1) 00826 return a(0,0); 00827 if(n == 2) 00828 return a(0,0)*a(1,1) - a(0,1)*a(1,0); 00829 if(method == "lu") 00830 { 00831 return detail::determinantByLUDecomposition(a); 00832 } 00833 else if(method == "cholesky") 00834 { 00835 Matrix<T> L(a.shape()); 00836 vigra_precondition(choleskyDecomposition(a, L), 00837 "determinant(): Cholesky method requires symmetric positive definite matrix."); 00838 T det = L(0,0); 00839 for(MultiArrayIndex k=1; k<n; ++k) 00840 det *= L(k,k); 00841 return sq(det); 00842 } 00843 else 00844 { 00845 vigra_precondition(false, "determinant(): Unknown solution method."); 00846 } 00847 return T(); 00848 } 00849 00850 /** Compute the logarithm of the determinant of a symmetric positive definite matrix. 00851 00852 This is useful to avoid multiplication of very large numbers in big matrices. 00853 It is implemented by means of Cholesky decomposition. 00854 00855 <b>\#include</b> <vigra/linear_solve.hxx> or<br> 00856 <b>\#include</b> <vigra/linear_algebra.hxx><br> 00857 Namespaces: vigra and vigra::linalg 00858 */ 00859 template <class T, class C1> 00860 T logDeterminant(MultiArrayView<2, T, C1> const & a) 00861 { 00862 MultiArrayIndex n = columnCount(a); 00863 vigra_precondition(rowCount(a) == n, 00864 "logDeterminant(): Square matrix required."); 00865 if(n == 1) 00866 { 00867 vigra_precondition(a(0,0) > 0.0, 00868 "logDeterminant(): Matrix not positive definite."); 00869 return std::log(a(0,0)); 00870 } 00871 if(n == 2) 00872 { 00873 T det = a(0,0)*a(1,1) - a(0,1)*a(1,0); 00874 vigra_precondition(det > 0.0, 00875 "logDeterminant(): Matrix not positive definite."); 00876 return std::log(det); 00877 } 00878 else 00879 { 00880 Matrix<T> L(a.shape()); 00881 vigra_precondition(choleskyDecomposition(a, L), 00882 "logDeterminant(): Matrix not positive definite."); 00883 T logdet = std::log(L(0,0)); 00884 for(MultiArrayIndex k=1; k<n; ++k) 00885 logdet += std::log(L(k,k)); // L(k,k) is guaranteed to be positive 00886 return 2.0*logdet; 00887 } 00888 } 00889 00890 /** Cholesky decomposition. 00891 00892 \a A must be a symmetric positive definite matrix, and \a L will be a lower 00893 triangular matrix, such that (up to round-off errors): 00894 00895 \code 00896 A == L * transpose(L); 00897 \endcode 00898 00899 This implementation cannot be applied in-place, i.e. <tt>&L == &A</tt> is an error. 00900 If \a A is not symmetric, a <tt>ContractViolation</tt> exception is thrown. If it 00901 is not positive definite, the function returns <tt>false</tt>. 00902 00903 <b>\#include</b> <vigra/linear_solve.hxx> or<br> 00904 <b>\#include</b> <vigra/linear_algebra.hxx><br> 00905 Namespaces: vigra and vigra::linalg 00906 */ 00907 template <class T, class C1, class C2> 00908 bool choleskyDecomposition(MultiArrayView<2, T, C1> const & A, 00909 MultiArrayView<2, T, C2> &L) 00910 { 00911 MultiArrayIndex n = columnCount(A); 00912 vigra_precondition(rowCount(A) == n, 00913 "choleskyDecomposition(): Input matrix must be square."); 00914 vigra_precondition(n == columnCount(L) && n == rowCount(L), 00915 "choleskyDecomposition(): Output matrix must have same shape as input matrix."); 00916 vigra_precondition(isSymmetric(A), 00917 "choleskyDecomposition(): Input matrix must be symmetric."); 00918 00919 for (MultiArrayIndex j = 0; j < n; ++j) 00920 { 00921 T d(0.0); 00922 for (MultiArrayIndex k = 0; k < j; ++k) 00923 { 00924 T s(0.0); 00925 for (MultiArrayIndex i = 0; i < k; ++i) 00926 { 00927 s += L(k, i)*L(j, i); 00928 } 00929 L(j, k) = s = (A(j, k) - s)/L(k, k); 00930 d = d + s*s; 00931 } 00932 d = A(j, j) - d; 00933 if(d <= 0.0) 00934 return false; // A is not positive definite 00935 L(j, j) = std::sqrt(d); 00936 for (MultiArrayIndex k = j+1; k < n; ++k) 00937 { 00938 L(j, k) = 0.0; 00939 } 00940 } 00941 return true; 00942 } 00943 00944 /** QR decomposition. 00945 00946 \a a contains the original matrix, results are returned in \a q and \a r, where 00947 \a q is a orthogonal matrix, and \a r is an upper triangular matrix, such that 00948 (up to round-off errors): 00949 00950 \code 00951 a == q * r; 00952 \endcode 00953 00954 If \a a doesn't have full rank, the function returns <tt>false</tt>. 00955 The decomposition is computed by householder transformations. It can be applied in-place, 00956 i.e. <tt>&a == &q</tt> or <tt>&a == &r</tt> are allowed. 00957 00958 <b>\#include</b> <vigra/linear_solve.hxx> or<br> 00959 <b>\#include</b> <vigra/linear_algebra.hxx><br> 00960 Namespaces: vigra and vigra::linalg 00961 */ 00962 template <class T, class C1, class C2, class C3> 00963 bool qrDecomposition(MultiArrayView<2, T, C1> const & a, 00964 MultiArrayView<2, T, C2> &q, MultiArrayView<2, T, C3> &r, 00965 double epsilon = 0.0) 00966 { 00967 const MultiArrayIndex m = rowCount(a); 00968 const MultiArrayIndex n = columnCount(a); 00969 vigra_precondition(n == columnCount(r) && m == rowCount(r) && 00970 m == columnCount(q) && m == rowCount(q), 00971 "qrDecomposition(): Matrix shape mismatch."); 00972 00973 q = identityMatrix<T>(m); 00974 MultiArrayView<2,T, StridedArrayTag> tq = transpose(q); 00975 r = a; 00976 ArrayVector<MultiArrayIndex> noPivoting; // intentionally empty 00977 return ((MultiArrayIndex)detail::qrTransformToUpperTriangular(r, tq, noPivoting, epsilon) == std::min(m,n)); 00978 } 00979 00980 /** Deprecated, use \ref linearSolveUpperTriangular(). 00981 */ 00982 template <class T, class C1, class C2, class C3> 00983 inline 00984 bool reverseElimination(const MultiArrayView<2, T, C1> &r, const MultiArrayView<2, T, C2> &b, 00985 MultiArrayView<2, T, C3> x) 00986 { 00987 return linearSolveUpperTriangular(r, b, x); 00988 } 00989 00990 /** Solve a linear system with upper-triangular coefficient matrix. 00991 00992 The square matrix \a r must be an upper-triangular coefficient matrix as can, 00993 for example, be obtained by means of QR decomposition. If \a r doesn't have full rank 00994 the function fails and returns <tt>false</tt>, otherwise it returns <tt>true</tt>. The 00995 lower triangular part of matrix \a r will not be touched, so it doesn't need to contain zeros. 00996 00997 The column vectors of matrix \a b are the right-hand sides of the equation (several equations 00998 with the same coefficients can thus be solved in one go). The result is returned 00999 int \a x, whose columns contain the solutions for the corresponding 01000 columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed. 01001 The following size requirements apply: 01002 01003 \code 01004 rowCount(r) == columnCount(r); 01005 rowCount(r) == rowCount(b); 01006 columnCount(r) == rowCount(x); 01007 columnCount(b) == columnCount(x); 01008 \endcode 01009 01010 <b>\#include</b> <vigra/linear_solve.hxx> or<br> 01011 <b>\#include</b> <vigra/linear_algebra.hxx><br> 01012 Namespaces: vigra and vigra::linalg 01013 */ 01014 template <class T, class C1, class C2, class C3> 01015 bool linearSolveUpperTriangular(const MultiArrayView<2, T, C1> &r, const MultiArrayView<2, T, C2> &b, 01016 MultiArrayView<2, T, C3> x) 01017 { 01018 typedef MultiArrayShape<2>::type Shape; 01019 MultiArrayIndex m = rowCount(r); 01020 MultiArrayIndex rhsCount = columnCount(b); 01021 vigra_precondition(m == columnCount(r), 01022 "linearSolveUpperTriangular(): square coefficient matrix required."); 01023 vigra_precondition(m == rowCount(b) && m == rowCount(x) && rhsCount == columnCount(x), 01024 "linearSolveUpperTriangular(): matrix shape mismatch."); 01025 01026 for(MultiArrayIndex k = 0; k < rhsCount; ++k) 01027 { 01028 for(int i=m-1; i>=0; --i) 01029 { 01030 if(r(i,i) == NumericTraits<T>::zero()) 01031 return false; // r doesn't have full rank 01032 T sum = b(i, k); 01033 for(MultiArrayIndex j=i+1; j<m; ++j) 01034 sum -= r(i, j) * x(j, k); 01035 x(i, k) = sum / r(i, i); 01036 } 01037 } 01038 return true; 01039 } 01040 01041 /** Solve a linear system with lower-triangular coefficient matrix. 01042 01043 The square matrix \a l must be a lower-triangular coefficient matrix. If \a l 01044 doesn't have full rank the function fails and returns <tt>false</tt>, 01045 otherwise it returns <tt>true</tt>. The upper triangular part of matrix \a l will not be touched, 01046 so it doesn't need to contain zeros. 01047 01048 The column vectors of matrix \a b are the right-hand sides of the equation (several equations 01049 with the same coefficients can thus be solved in one go). The result is returned 01050 in \a x, whose columns contain the solutions for the corresponding 01051 columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed. 01052 The following size requirements apply: 01053 01054 \code 01055 rowCount(l) == columnCount(l); 01056 rowCount(l) == rowCount(b); 01057 columnCount(l) == rowCount(x); 01058 columnCount(b) == columnCount(x); 01059 \endcode 01060 01061 <b>\#include</b> <vigra/linear_solve.hxx> or<br> 01062 <b>\#include</b> <vigra/linear_algebra.hxx><br> 01063 Namespaces: vigra and vigra::linalg 01064 */ 01065 template <class T, class C1, class C2, class C3> 01066 bool linearSolveLowerTriangular(const MultiArrayView<2, T, C1> &l, const MultiArrayView<2, T, C2> &b, 01067 MultiArrayView<2, T, C3> x) 01068 { 01069 MultiArrayIndex m = columnCount(l); 01070 MultiArrayIndex n = columnCount(b); 01071 vigra_precondition(m == rowCount(l), 01072 "linearSolveLowerTriangular(): square coefficient matrix required."); 01073 vigra_precondition(m == rowCount(b) && m == rowCount(x) && n == columnCount(x), 01074 "linearSolveLowerTriangular(): matrix shape mismatch."); 01075 01076 for(MultiArrayIndex k = 0; k < n; ++k) 01077 { 01078 for(MultiArrayIndex i=0; i<m; ++i) 01079 { 01080 if(l(i,i) == NumericTraits<T>::zero()) 01081 return false; // l doesn't have full rank 01082 T sum = b(i, k); 01083 for(MultiArrayIndex j=0; j<i; ++j) 01084 sum -= l(i, j) * x(j, k); 01085 x(i, k) = sum / l(i, i); 01086 } 01087 } 01088 return true; 01089 } 01090 01091 01092 /** Solve a linear system when the Cholesky decomposition of the left hand side is given. 01093 01094 The square matrix \a L must be a lower-triangular matrix resulting from Cholesky 01095 decomposition of some positive definite coefficient matrix. 01096 01097 The column vectors of matrix \a b are the right-hand sides of the equation (several equations 01098 with the same matrix \a L can thus be solved in one go). The result is returned 01099 in \a x, whose columns contain the solutions for the corresponding 01100 columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed. 01101 The following size requirements apply: 01102 01103 \code 01104 rowCount(L) == columnCount(L); 01105 rowCount(L) == rowCount(b); 01106 columnCount(L) == rowCount(x); 01107 columnCount(b) == columnCount(x); 01108 \endcode 01109 01110 <b>\#include</b> <vigra/linear_solve.hxx> or<br> 01111 <b>\#include</b> <vigra/linear_algebra.hxx><br> 01112 Namespaces: vigra and vigra::linalg 01113 */ 01114 template <class T, class C1, class C2, class C3> 01115 inline 01116 void choleskySolve(MultiArrayView<2, T, C1> & L, MultiArrayView<2, T, C2> const & b, MultiArrayView<2, T, C3> & x) 01117 { 01118 /* Solve L * y = b */ 01119 linearSolveLowerTriangular(L, b, x); 01120 /* Solve L^T * x = y */ 01121 linearSolveUpperTriangular(transpose(L), x, x); 01122 } 01123 01124 /** Solve a linear system. 01125 01126 \a A is the coefficient matrix, and the column vectors 01127 in \a b are the right-hand sides of the equation (so, several equations 01128 with the same coefficients can be solved in one go). The result is returned 01129 in \a res, whose columns contain the solutions for the corresponding 01130 columns of \a b. The number of columns of \a A must equal the number of rows of 01131 both \a b and \a res, and the number of columns of \a b and \a res must match. 01132 01133 \a method must be one of the following: 01134 <DL> 01135 <DT>"Cholesky"<DD> Compute the solution by means of Cholesky decomposition. The 01136 coefficient matrix \a A must by symmetric positive definite. If 01137 this is not the case, the function returns <tt>false</tt>. 01138 01139 <DT>"QR"<DD> (default) Compute the solution by means of QR decomposition. The 01140 coefficient matrix \a A can be square or rectangular. In the latter case, 01141 it must have more rows than columns, and the solution will be computed in the 01142 least squares sense. If \a A doesn't have full rank, the function 01143 returns <tt>false</tt>. 01144 01145 <DT>"SVD"<DD> Compute the solution by means of singular value decomposition. The 01146 coefficient matrix \a A can be square or rectangular. In the latter case, 01147 it must have more rows than columns, and the solution will be computed in the 01148 least squares sense. If \a A doesn't have full rank, the function 01149 returns <tt>false</tt>. 01150 01151 <DT>"NE"<DD> Compute the solution by means of the normal equations, i.e. by applying Cholesky 01152 decomposition to the equivalent problem <tt>A'*A*x = A'*b</tt>. This only makes sense 01153 when the equation is to be solved in the least squares sense, i.e. when \a A is a 01154 rectangular matrix with more rows than columns. If \a A doesn't have full column rank, 01155 the function returns <tt>false</tt>. 01156 </DL> 01157 01158 This function can be applied in-place, i.e. <tt>&b == &res</tt> or <tt>&A == &res</tt> are allowed 01159 (provided they have the required shapes). 01160 01161 The following size requirements apply: 01162 01163 \code 01164 rowCount(r) == rowCount(b); 01165 columnCount(r) == rowCount(x); 01166 columnCount(b) == columnCount(x); 01167 \endcode 01168 01169 <b>\#include</b> <vigra/linear_solve.hxx> or<br> 01170 <b>\#include</b> <vigra/linear_algebra.hxx><br> 01171 Namespaces: vigra and vigra::linalg 01172 */ 01173 template <class T, class C1, class C2, class C3> 01174 bool linearSolve(const MultiArrayView<2, T, C1> &A, const MultiArrayView<2, T, C2> &b, 01175 MultiArrayView<2, T, C3> & res, std::string method = "QR") 01176 { 01177 typedef typename Matrix<T>::difference_type Shape; 01178 typedef typename Matrix<T>::view_type SubMatrix; 01179 01180 const MultiArrayIndex n = columnCount(A); 01181 const MultiArrayIndex m = rowCount(A); 01182 01183 vigra_precondition(n <= m, 01184 "linearSolve(): Coefficient matrix A must have at least as many rows as columns."); 01185 vigra_precondition(n == rowCount(res) && 01186 m == rowCount(b) && columnCount(b) == columnCount(res), 01187 "linearSolve(): matrix shape mismatch."); 01188 01189 for(unsigned int k=0; k<method.size(); ++k) 01190 method[k] = (std::string::value_type)tolower(method[k]); 01191 01192 if(method == "cholesky") 01193 { 01194 vigra_precondition(columnCount(A) == rowCount(A), 01195 "linearSolve(): Cholesky method requires square coefficient matrix."); 01196 Matrix<T> L(A.shape()); 01197 if(!choleskyDecomposition(A, L)) 01198 return false; // false if A wasn't symmetric positive definite 01199 choleskySolve(L, b, res); 01200 } 01201 else if(method == "qr") 01202 { 01203 return (MultiArrayIndex)linearSolveQR(A, b, res) == n; 01204 } 01205 else if(method == "ne") 01206 { 01207 return linearSolve(transpose(A)*A, transpose(A)*b, res, "Cholesky"); 01208 } 01209 else if(method == "svd") 01210 { 01211 MultiArrayIndex rhsCount = columnCount(b); 01212 Matrix<T> u(A.shape()), s(n, 1), v(n, n); 01213 01214 MultiArrayIndex rank = (MultiArrayIndex)singularValueDecomposition(A, u, s, v); 01215 01216 Matrix<T> t = transpose(u)*b; 01217 for(MultiArrayIndex l=0; l<rhsCount; ++l) 01218 { 01219 for(MultiArrayIndex k=0; k<rank; ++k) 01220 t(k,l) /= s(k,0); 01221 for(MultiArrayIndex k=rank; k<n; ++k) 01222 t(k,l) = NumericTraits<T>::zero(); 01223 } 01224 res = v*t; 01225 01226 return (rank == n); 01227 } 01228 else 01229 { 01230 vigra_precondition(false, "linearSolve(): Unknown solution method."); 01231 } 01232 return true; 01233 } 01234 01235 //@} 01236 01237 } // namespace linalg 01238 01239 using linalg::inverse; 01240 using linalg::determinant; 01241 using linalg::logDeterminant; 01242 using linalg::linearSolve; 01243 using linalg::choleskySolve; 01244 using linalg::choleskyDecomposition; 01245 using linalg::qrDecomposition; 01246 using linalg::linearSolveUpperTriangular; 01247 using linalg::linearSolveLowerTriangular; 01248 01249 } // namespace vigra 01250 01251 01252 #endif // VIGRA_LINEAR_SOLVE_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|