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