[ 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 method = tolower(method); 00823 00824 if(n == 1) 00825 return a(0,0); 00826 if(n == 2) 00827 return a(0,0)*a(1,1) - a(0,1)*a(1,0); 00828 if(method == "lu") 00829 { 00830 return detail::determinantByLUDecomposition(a); 00831 } 00832 else if(method == "cholesky") 00833 { 00834 Matrix<T> L(a.shape()); 00835 vigra_precondition(choleskyDecomposition(a, L), 00836 "determinant(): Cholesky method requires symmetric positive definite matrix."); 00837 T det = L(0,0); 00838 for(MultiArrayIndex k=1; k<n; ++k) 00839 det *= L(k,k); 00840 return sq(det); 00841 } 00842 else 00843 { 00844 vigra_precondition(false, "determinant(): Unknown solution method."); 00845 } 00846 return T(); 00847 } 00848 00849 /** Compute the logarithm of the determinant of a symmetric positive definite matrix. 00850 00851 This is useful to avoid multiplication of very large numbers in big matrices. 00852 It is implemented by means of Cholesky decomposition. 00853 00854 <b>\#include</b> <vigra/linear_solve.hxx> or<br> 00855 <b>\#include</b> <vigra/linear_algebra.hxx><br> 00856 Namespaces: vigra and vigra::linalg 00857 */ 00858 template <class T, class C1> 00859 T logDeterminant(MultiArrayView<2, T, C1> const & a) 00860 { 00861 MultiArrayIndex n = columnCount(a); 00862 vigra_precondition(rowCount(a) == n, 00863 "logDeterminant(): Square matrix required."); 00864 if(n == 1) 00865 { 00866 vigra_precondition(a(0,0) > 0.0, 00867 "logDeterminant(): Matrix not positive definite."); 00868 return std::log(a(0,0)); 00869 } 00870 if(n == 2) 00871 { 00872 T det = a(0,0)*a(1,1) - a(0,1)*a(1,0); 00873 vigra_precondition(det > 0.0, 00874 "logDeterminant(): Matrix not positive definite."); 00875 return std::log(det); 00876 } 00877 else 00878 { 00879 Matrix<T> L(a.shape()); 00880 vigra_precondition(choleskyDecomposition(a, L), 00881 "logDeterminant(): Matrix not positive definite."); 00882 T logdet = std::log(L(0,0)); 00883 for(MultiArrayIndex k=1; k<n; ++k) 00884 logdet += std::log(L(k,k)); // L(k,k) is guaranteed to be positive 00885 return 2.0*logdet; 00886 } 00887 } 00888 00889 /** Cholesky decomposition. 00890 00891 \a A must be a symmetric positive definite matrix, and \a L will be a lower 00892 triangular matrix, such that (up to round-off errors): 00893 00894 \code 00895 A == L * transpose(L); 00896 \endcode 00897 00898 This implementation cannot be applied in-place, i.e. <tt>&L == &A</tt> is an error. 00899 If \a A is not symmetric, a <tt>ContractViolation</tt> exception is thrown. If it 00900 is not positive definite, the function returns <tt>false</tt>. 00901 00902 <b>\#include</b> <vigra/linear_solve.hxx> or<br> 00903 <b>\#include</b> <vigra/linear_algebra.hxx><br> 00904 Namespaces: vigra and vigra::linalg 00905 */ 00906 template <class T, class C1, class C2> 00907 bool choleskyDecomposition(MultiArrayView<2, T, C1> const & A, 00908 MultiArrayView<2, T, C2> &L) 00909 { 00910 MultiArrayIndex n = columnCount(A); 00911 vigra_precondition(rowCount(A) == n, 00912 "choleskyDecomposition(): Input matrix must be square."); 00913 vigra_precondition(n == columnCount(L) && n == rowCount(L), 00914 "choleskyDecomposition(): Output matrix must have same shape as input matrix."); 00915 vigra_precondition(isSymmetric(A), 00916 "choleskyDecomposition(): Input matrix must be symmetric."); 00917 00918 for (MultiArrayIndex j = 0; j < n; ++j) 00919 { 00920 T d(0.0); 00921 for (MultiArrayIndex k = 0; k < j; ++k) 00922 { 00923 T s(0.0); 00924 for (MultiArrayIndex i = 0; i < k; ++i) 00925 { 00926 s += L(k, i)*L(j, i); 00927 } 00928 L(j, k) = s = (A(j, k) - s)/L(k, k); 00929 d = d + s*s; 00930 } 00931 d = A(j, j) - d; 00932 if(d <= 0.0) 00933 return false; // A is not positive definite 00934 L(j, j) = std::sqrt(d); 00935 for (MultiArrayIndex k = j+1; k < n; ++k) 00936 { 00937 L(j, k) = 0.0; 00938 } 00939 } 00940 return true; 00941 } 00942 00943 /** QR decomposition. 00944 00945 \a a contains the original matrix, results are returned in \a q and \a r, where 00946 \a q is a orthogonal matrix, and \a r is an upper triangular matrix, such that 00947 (up to round-off errors): 00948 00949 \code 00950 a == q * r; 00951 \endcode 00952 00953 If \a a doesn't have full rank, the function returns <tt>false</tt>. 00954 The decomposition is computed by householder transformations. It can be applied in-place, 00955 i.e. <tt>&a == &q</tt> or <tt>&a == &r</tt> are allowed. 00956 00957 <b>\#include</b> <vigra/linear_solve.hxx> or<br> 00958 <b>\#include</b> <vigra/linear_algebra.hxx><br> 00959 Namespaces: vigra and vigra::linalg 00960 */ 00961 template <class T, class C1, class C2, class C3> 00962 bool qrDecomposition(MultiArrayView<2, T, C1> const & a, 00963 MultiArrayView<2, T, C2> &q, MultiArrayView<2, T, C3> &r, 00964 double epsilon = 0.0) 00965 { 00966 const MultiArrayIndex m = rowCount(a); 00967 const MultiArrayIndex n = columnCount(a); 00968 vigra_precondition(n == columnCount(r) && m == rowCount(r) && 00969 m == columnCount(q) && m == rowCount(q), 00970 "qrDecomposition(): Matrix shape mismatch."); 00971 00972 q = identityMatrix<T>(m); 00973 MultiArrayView<2,T, StridedArrayTag> tq = transpose(q); 00974 r = a; 00975 ArrayVector<MultiArrayIndex> noPivoting; // intentionally empty 00976 return ((MultiArrayIndex)detail::qrTransformToUpperTriangular(r, tq, noPivoting, epsilon) == std::min(m,n)); 00977 } 00978 00979 /** Deprecated, use \ref linearSolveUpperTriangular(). 00980 */ 00981 template <class T, class C1, class C2, class C3> 00982 inline 00983 bool reverseElimination(const MultiArrayView<2, T, C1> &r, const MultiArrayView<2, T, C2> &b, 00984 MultiArrayView<2, T, C3> x) 00985 { 00986 return linearSolveUpperTriangular(r, b, x); 00987 } 00988 00989 /** Solve a linear system with upper-triangular coefficient matrix. 00990 00991 The square matrix \a r must be an upper-triangular coefficient matrix as can, 00992 for example, be obtained by means of QR decomposition. If \a r doesn't have full rank 00993 the function fails and returns <tt>false</tt>, otherwise it returns <tt>true</tt>. The 00994 lower triangular part of matrix \a r will not be touched, so it doesn't need to contain zeros. 00995 00996 The column vectors of matrix \a b are the right-hand sides of the equation (several equations 00997 with the same coefficients can thus be solved in one go). The result is returned 00998 int \a x, whose columns contain the solutions for the corresponding 00999 columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed. 01000 The following size requirements apply: 01001 01002 \code 01003 rowCount(r) == columnCount(r); 01004 rowCount(r) == rowCount(b); 01005 columnCount(r) == rowCount(x); 01006 columnCount(b) == columnCount(x); 01007 \endcode 01008 01009 <b>\#include</b> <vigra/linear_solve.hxx> or<br> 01010 <b>\#include</b> <vigra/linear_algebra.hxx><br> 01011 Namespaces: vigra and vigra::linalg 01012 */ 01013 template <class T, class C1, class C2, class C3> 01014 bool linearSolveUpperTriangular(const MultiArrayView<2, T, C1> &r, const MultiArrayView<2, T, C2> &b, 01015 MultiArrayView<2, T, C3> x) 01016 { 01017 typedef MultiArrayShape<2>::type Shape; 01018 MultiArrayIndex m = rowCount(r); 01019 MultiArrayIndex rhsCount = columnCount(b); 01020 vigra_precondition(m == columnCount(r), 01021 "linearSolveUpperTriangular(): square coefficient matrix required."); 01022 vigra_precondition(m == rowCount(b) && m == rowCount(x) && rhsCount == columnCount(x), 01023 "linearSolveUpperTriangular(): matrix shape mismatch."); 01024 01025 for(MultiArrayIndex k = 0; k < rhsCount; ++k) 01026 { 01027 for(int i=m-1; i>=0; --i) 01028 { 01029 if(r(i,i) == NumericTraits<T>::zero()) 01030 return false; // r doesn't have full rank 01031 T sum = b(i, k); 01032 for(MultiArrayIndex j=i+1; j<m; ++j) 01033 sum -= r(i, j) * x(j, k); 01034 x(i, k) = sum / r(i, i); 01035 } 01036 } 01037 return true; 01038 } 01039 01040 /** Solve a linear system with lower-triangular coefficient matrix. 01041 01042 The square matrix \a l must be a lower-triangular coefficient matrix. If \a l 01043 doesn't have full rank the function fails and returns <tt>false</tt>, 01044 otherwise it returns <tt>true</tt>. The upper triangular part of matrix \a l will not be touched, 01045 so it doesn't need to contain zeros. 01046 01047 The column vectors of matrix \a b are the right-hand sides of the equation (several equations 01048 with the same coefficients can thus be solved in one go). The result is returned 01049 in \a x, whose columns contain the solutions for the corresponding 01050 columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed. 01051 The following size requirements apply: 01052 01053 \code 01054 rowCount(l) == columnCount(l); 01055 rowCount(l) == rowCount(b); 01056 columnCount(l) == rowCount(x); 01057 columnCount(b) == columnCount(x); 01058 \endcode 01059 01060 <b>\#include</b> <vigra/linear_solve.hxx> or<br> 01061 <b>\#include</b> <vigra/linear_algebra.hxx><br> 01062 Namespaces: vigra and vigra::linalg 01063 */ 01064 template <class T, class C1, class C2, class C3> 01065 bool linearSolveLowerTriangular(const MultiArrayView<2, T, C1> &l, const MultiArrayView<2, T, C2> &b, 01066 MultiArrayView<2, T, C3> x) 01067 { 01068 MultiArrayIndex m = columnCount(l); 01069 MultiArrayIndex n = columnCount(b); 01070 vigra_precondition(m == rowCount(l), 01071 "linearSolveLowerTriangular(): square coefficient matrix required."); 01072 vigra_precondition(m == rowCount(b) && m == rowCount(x) && n == columnCount(x), 01073 "linearSolveLowerTriangular(): matrix shape mismatch."); 01074 01075 for(MultiArrayIndex k = 0; k < n; ++k) 01076 { 01077 for(MultiArrayIndex i=0; i<m; ++i) 01078 { 01079 if(l(i,i) == NumericTraits<T>::zero()) 01080 return false; // l doesn't have full rank 01081 T sum = b(i, k); 01082 for(MultiArrayIndex j=0; j<i; ++j) 01083 sum -= l(i, j) * x(j, k); 01084 x(i, k) = sum / l(i, i); 01085 } 01086 } 01087 return true; 01088 } 01089 01090 01091 /** Solve a linear system when the Cholesky decomposition of the left hand side is given. 01092 01093 The square matrix \a L must be a lower-triangular matrix resulting from Cholesky 01094 decomposition of some positive definite coefficient matrix. 01095 01096 The column vectors of matrix \a b are the right-hand sides of the equation (several equations 01097 with the same matrix \a L can thus be solved in one go). The result is returned 01098 in \a x, whose columns contain the solutions for the corresponding 01099 columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed. 01100 The following size requirements apply: 01101 01102 \code 01103 rowCount(L) == columnCount(L); 01104 rowCount(L) == rowCount(b); 01105 columnCount(L) == rowCount(x); 01106 columnCount(b) == columnCount(x); 01107 \endcode 01108 01109 <b>\#include</b> <vigra/linear_solve.hxx> or<br> 01110 <b>\#include</b> <vigra/linear_algebra.hxx><br> 01111 Namespaces: vigra and vigra::linalg 01112 */ 01113 template <class T, class C1, class C2, class C3> 01114 inline 01115 void choleskySolve(MultiArrayView<2, T, C1> & L, MultiArrayView<2, T, C2> const & b, MultiArrayView<2, T, C3> & x) 01116 { 01117 /* Solve L * y = b */ 01118 linearSolveLowerTriangular(L, b, x); 01119 /* Solve L^T * x = y */ 01120 linearSolveUpperTriangular(transpose(L), x, x); 01121 } 01122 01123 /** Solve a linear system. 01124 01125 \a A is the coefficient matrix, and the column vectors 01126 in \a b are the right-hand sides of the equation (so, several equations 01127 with the same coefficients can be solved in one go). The result is returned 01128 in \a res, whose columns contain the solutions for the corresponding 01129 columns of \a b. The number of columns of \a A must equal the number of rows of 01130 both \a b and \a res, and the number of columns of \a b and \a res must match. 01131 01132 \a method must be one of the following: 01133 <DL> 01134 <DT>"Cholesky"<DD> Compute the solution by means of Cholesky decomposition. The 01135 coefficient matrix \a A must by symmetric positive definite. If 01136 this is not the case, the function returns <tt>false</tt>. 01137 01138 <DT>"QR"<DD> (default) Compute the solution by means of QR decomposition. The 01139 coefficient matrix \a A can be square or rectangular. In the latter case, 01140 it must have more rows than columns, and the solution will be computed in the 01141 least squares sense. If \a A doesn't have full rank, the function 01142 returns <tt>false</tt>. 01143 01144 <DT>"SVD"<DD> Compute the solution by means of singular value decomposition. The 01145 coefficient matrix \a A can be square or rectangular. In the latter case, 01146 it must have more rows than columns, and the solution will be computed in the 01147 least squares sense. If \a A doesn't have full rank, the function 01148 returns <tt>false</tt>. 01149 01150 <DT>"NE"<DD> Compute the solution by means of the normal equations, i.e. by applying Cholesky 01151 decomposition to the equivalent problem <tt>A'*A*x = A'*b</tt>. This only makes sense 01152 when the equation is to be solved in the least squares sense, i.e. when \a A is a 01153 rectangular matrix with more rows than columns. If \a A doesn't have full column rank, 01154 the function returns <tt>false</tt>. 01155 </DL> 01156 01157 This function can be applied in-place, i.e. <tt>&b == &res</tt> or <tt>&A == &res</tt> are allowed 01158 (provided they have the required shapes). 01159 01160 The following size requirements apply: 01161 01162 \code 01163 rowCount(r) == rowCount(b); 01164 columnCount(r) == rowCount(x); 01165 columnCount(b) == columnCount(x); 01166 \endcode 01167 01168 <b>\#include</b> <vigra/linear_solve.hxx> or<br> 01169 <b>\#include</b> <vigra/linear_algebra.hxx><br> 01170 Namespaces: vigra and vigra::linalg 01171 */ 01172 template <class T, class C1, class C2, class C3> 01173 bool linearSolve(const MultiArrayView<2, T, C1> &A, const MultiArrayView<2, T, C2> &b, 01174 MultiArrayView<2, T, C3> & res, std::string method = "QR") 01175 { 01176 typedef typename Matrix<T>::difference_type Shape; 01177 typedef typename Matrix<T>::view_type SubMatrix; 01178 01179 const MultiArrayIndex n = columnCount(A); 01180 const MultiArrayIndex m = rowCount(A); 01181 01182 vigra_precondition(n <= m, 01183 "linearSolve(): Coefficient matrix A must have at least as many rows as columns."); 01184 vigra_precondition(n == rowCount(res) && 01185 m == rowCount(b) && columnCount(b) == columnCount(res), 01186 "linearSolve(): matrix shape mismatch."); 01187 01188 method = tolower(method); 01189 if(method == "cholesky") 01190 { 01191 vigra_precondition(columnCount(A) == rowCount(A), 01192 "linearSolve(): Cholesky method requires square coefficient matrix."); 01193 Matrix<T> L(A.shape()); 01194 if(!choleskyDecomposition(A, L)) 01195 return false; // false if A wasn't symmetric positive definite 01196 choleskySolve(L, b, res); 01197 } 01198 else if(method == "qr") 01199 { 01200 return (MultiArrayIndex)linearSolveQR(A, b, res) == n; 01201 } 01202 else if(method == "ne") 01203 { 01204 return linearSolve(transpose(A)*A, transpose(A)*b, res, "Cholesky"); 01205 } 01206 else if(method == "svd") 01207 { 01208 MultiArrayIndex rhsCount = columnCount(b); 01209 Matrix<T> u(A.shape()), s(n, 1), v(n, n); 01210 01211 MultiArrayIndex rank = (MultiArrayIndex)singularValueDecomposition(A, u, s, v); 01212 01213 Matrix<T> t = transpose(u)*b; 01214 for(MultiArrayIndex l=0; l<rhsCount; ++l) 01215 { 01216 for(MultiArrayIndex k=0; k<rank; ++k) 01217 t(k,l) /= s(k,0); 01218 for(MultiArrayIndex k=rank; k<n; ++k) 01219 t(k,l) = NumericTraits<T>::zero(); 01220 } 01221 res = v*t; 01222 01223 return (rank == n); 01224 } 01225 else 01226 { 01227 vigra_precondition(false, "linearSolve(): Unknown solution method."); 01228 } 01229 return true; 01230 } 01231 01232 //@} 01233 01234 } // namespace linalg 01235 01236 using linalg::inverse; 01237 using linalg::determinant; 01238 using linalg::logDeterminant; 01239 using linalg::linearSolve; 01240 using linalg::choleskySolve; 01241 using linalg::choleskyDecomposition; 01242 using linalg::qrDecomposition; 01243 using linalg::linearSolveUpperTriangular; 01244 using linalg::linearSolveLowerTriangular; 01245 01246 } // namespace vigra 01247 01248 01249 #endif // VIGRA_LINEAR_SOLVE_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|