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