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