[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

vigra/linear_solve.hxx

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*     Copyright 2003-2008 by Gunnar Kedenburg and Ullrich Koethe       */
00004 /*       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)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
VIGRA 1.6.0 (13 Aug 2008)