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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.7.0 (Thu Aug 25 2011)