[ 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 /*    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)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
VIGRA 1.6.0 (5 Nov 2009)