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

vigra/linear_solve.hxx VIGRA

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, \a res must have the transposed shape 
00708         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> <vigra/linear_solve.hxx> or<br>
00715     <b>\#include</b> <vigra/linear_algebra.hxx><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     typedef typename MultiArrayShape<2>::type Shape;
00722     
00723     const MultiArrayIndex n = columnCount(v);
00724     const MultiArrayIndex m = rowCount(v);
00725     vigra_precondition(n == rowCount(res) && m == columnCount(res),
00726        "inverse(): shape of output matrix must be the transpose of the input matrix' shape.");
00727     
00728     if(m < n)
00729     {
00730         MultiArrayView<2, T, StridedArrayTag> vt = transpose(v);
00731         Matrix<T> r(vt.shape()), q(n, n);
00732         if(!qrDecomposition(vt, q, r))
00733             return false; // a didn't have full rank
00734         linearSolveUpperTriangular(r.subarray(Shape(0,0), Shape(m,m)), 
00735                                    transpose(q).subarray(Shape(0,0), Shape(m,n)), 
00736                                    transpose(res)); 
00737     }
00738     else
00739     {
00740         Matrix<T> r(v.shape()), q(m, m);
00741         if(!qrDecomposition(v, q, r))
00742             return false; // a didn't have full rank
00743         linearSolveUpperTriangular(r.subarray(Shape(0,0), Shape(n,n)), 
00744                                    transpose(q).subarray(Shape(0,0), Shape(n,m)), 
00745                                    res); 
00746     }
00747     return true;
00748 }
00749 
00750     /** Create the inverse or pseudo-inverse of matrix \a v.
00751 
00752         The result is returned as a temporary matrix. If the matrix \a v is square, 
00753         the result will have the same shape and contains the inverse of \a v. 
00754         If \a v is rectangular, the result will have the transposed shape of \a v. 
00755         The inverse is then computed in the least-squares 
00756         sense, i.e. \a res will be the pseudo-inverse (Moore-Penrose inverse).
00757         The inverse is computed by means of QR decomposition. If \a v
00758         is not invertible, <tt>vigra::PreconditionViolation</tt> exception is thrown.
00759         Usage:
00760 
00761         \code
00762         vigra::Matrix<double> v(n, n);
00763         v = ...;
00764 
00765         vigra::Matrix<double> m = inverse(v);
00766         \endcode
00767 
00768     <b>\#include</b> <vigra/linear_solve.hxx> or<br>
00769     <b>\#include</b> <vigra/linear_algebra.hxx><br>
00770         Namespaces: vigra and vigra::linalg
00771      */
00772 template <class T, class C>
00773 TemporaryMatrix<T> inverse(const MultiArrayView<2, T, C> &v)
00774 {
00775     TemporaryMatrix<T> ret(columnCount(v), rowCount(v));  // transpose shape
00776     vigra_precondition(inverse(v, ret),
00777         "inverse(): matrix is not invertible.");
00778     return ret;
00779 }
00780 
00781 template <class T>
00782 TemporaryMatrix<T> inverse(const TemporaryMatrix<T> &v)
00783 {
00784     if(columnCount(v) == rowCount(v))
00785     {
00786         vigra_precondition(inverse(v, const_cast<TemporaryMatrix<T> &>(v)),
00787             "inverse(): matrix is not invertible.");
00788         return v;      
00789     }
00790     else
00791     {
00792         TemporaryMatrix<T> ret(columnCount(v), rowCount(v));  // transpose shape
00793         vigra_precondition(inverse(v, ret),
00794             "inverse(): matrix is not invertible.");
00795         return ret;
00796     }
00797 }
00798 
00799     /** Compute the determinant of a square matrix.
00800 
00801         \a method must be one of the following:
00802         <DL>
00803         <DT>"Cholesky"<DD> Compute the solution by means of Cholesky decomposition. This
00804                            method is faster than "LU", but requires the matrix \a a 
00805                            to be symmetric positive definite. If this is 
00806                            not the case, a <tt>ContractViolation</tt> exception is thrown.
00807                            
00808         <DT>"LU"<DD> (default) Compute the solution by means of LU decomposition.
00809         </DL>
00810 
00811     <b>\#include</b> <vigra/linear_solve.hxx> or<br>
00812     <b>\#include</b> <vigra/linear_algebra.hxx><br>
00813         Namespaces: vigra and vigra::linalg
00814      */
00815 template <class T, class C1>
00816 T determinant(MultiArrayView<2, T, C1> const & a, std::string method = "LU")
00817 {
00818     MultiArrayIndex n = columnCount(a);
00819     vigra_precondition(rowCount(a) == n,
00820                "determinant(): Square matrix required.");    
00821 
00822     for(unsigned int k=0; k<method.size(); ++k)
00823         method[k] = tolower(method[k]);
00824     
00825     if(n == 1)
00826         return a(0,0);
00827     if(n == 2)
00828         return a(0,0)*a(1,1) - a(0,1)*a(1,0);
00829     if(method == "lu")
00830     {
00831         return detail::determinantByLUDecomposition(a);
00832     }
00833     else if(method == "cholesky")
00834     {
00835         Matrix<T> L(a.shape());
00836         vigra_precondition(choleskyDecomposition(a, L),
00837            "determinant(): Cholesky method requires symmetric positive definite matrix.");
00838         T det = L(0,0);
00839         for(MultiArrayIndex k=1; k<n; ++k)
00840             det *= L(k,k);
00841         return sq(det);
00842     }
00843     else
00844     {
00845         vigra_precondition(false, "determinant(): Unknown solution method.");
00846     }
00847     return T();
00848 }
00849 
00850     /** Compute the logarithm of the determinant of a symmetric positive definite matrix.
00851 
00852         This is useful to avoid multiplication of very large numbers in big matrices.
00853         It is implemented by means of Cholesky decomposition.
00854 
00855     <b>\#include</b> <vigra/linear_solve.hxx> or<br>
00856     <b>\#include</b> <vigra/linear_algebra.hxx><br>
00857         Namespaces: vigra and vigra::linalg
00858      */
00859 template <class T, class C1>
00860 T logDeterminant(MultiArrayView<2, T, C1> const & a)
00861 {
00862     MultiArrayIndex n = columnCount(a);
00863     vigra_precondition(rowCount(a) == n,
00864                "logDeterminant(): Square matrix required.");    
00865     if(n == 1)
00866     {
00867         vigra_precondition(a(0,0) > 0.0,
00868                    "logDeterminant(): Matrix not positive definite.");    
00869         return std::log(a(0,0));
00870     }
00871     if(n == 2)
00872     {
00873         T det = a(0,0)*a(1,1) - a(0,1)*a(1,0);
00874         vigra_precondition(det > 0.0,
00875                    "logDeterminant(): Matrix not positive definite.");    
00876         return std::log(det);
00877     }
00878     else
00879     {
00880         Matrix<T> L(a.shape());
00881         vigra_precondition(choleskyDecomposition(a, L),
00882                 "logDeterminant(): Matrix not positive definite.");  
00883         T logdet = std::log(L(0,0)); 
00884         for(MultiArrayIndex k=1; k<n; ++k)
00885             logdet += std::log(L(k,k));  // L(k,k) is guaranteed to be positive
00886         return 2.0*logdet;
00887     }
00888 }
00889 
00890     /** Cholesky decomposition.
00891 
00892         \a A must be a symmetric positive definite matrix, and \a L will be a lower
00893         triangular matrix, such that (up to round-off errors):
00894 
00895         \code
00896         A == L * transpose(L);
00897         \endcode
00898 
00899         This implementation cannot be applied in-place, i.e. <tt>&L == &A</tt> is an error.
00900         If \a A is not symmetric, a <tt>ContractViolation</tt> exception is thrown. If it
00901         is not positive definite, the function returns <tt>false</tt>.
00902 
00903     <b>\#include</b> <vigra/linear_solve.hxx> or<br>
00904     <b>\#include</b> <vigra/linear_algebra.hxx><br>
00905         Namespaces: vigra and vigra::linalg
00906      */
00907 template <class T, class C1, class C2>
00908 bool choleskyDecomposition(MultiArrayView<2, T, C1> const & A,
00909                            MultiArrayView<2, T, C2> &L)
00910 {
00911     MultiArrayIndex n = columnCount(A); 
00912     vigra_precondition(rowCount(A) == n,
00913                        "choleskyDecomposition(): Input matrix must be square.");
00914     vigra_precondition(n == columnCount(L) && n == rowCount(L),
00915                        "choleskyDecomposition(): Output matrix must have same shape as input matrix.");
00916     vigra_precondition(isSymmetric(A),
00917                        "choleskyDecomposition(): Input matrix must be symmetric.");
00918 
00919     for (MultiArrayIndex j = 0; j < n; ++j) 
00920     {
00921         T d(0.0);
00922         for (MultiArrayIndex k = 0; k < j; ++k) 
00923         {
00924             T s(0.0);
00925             for (MultiArrayIndex i = 0; i < k; ++i) 
00926             {
00927                s += L(k, i)*L(j, i);
00928             }
00929             L(j, k) = s = (A(j, k) - s)/L(k, k);
00930             d = d + s*s;
00931         }
00932         d = A(j, j) - d;
00933         if(d <= 0.0)
00934             return false;  // A is not positive definite
00935         L(j, j) = std::sqrt(d);
00936         for (MultiArrayIndex k = j+1; k < n; ++k) 
00937         {
00938            L(j, k) = 0.0;
00939         }
00940     }
00941     return true;
00942 }
00943 
00944     /** QR decomposition.
00945 
00946         \a a contains the original matrix, results are returned in \a q and \a r, where
00947         \a q is a orthogonal matrix, and \a r is an upper triangular matrix, such that 
00948         (up to round-off errors):
00949 
00950         \code
00951         a == q * r;
00952         \endcode
00953 
00954         If \a a doesn't have full rank, the function returns <tt>false</tt>. 
00955         The decomposition is computed by householder transformations. It can be applied in-place,
00956         i.e. <tt>&a == &q</tt> or <tt>&a == &r</tt> are allowed.
00957 
00958     <b>\#include</b> <vigra/linear_solve.hxx> or<br>
00959     <b>\#include</b> <vigra/linear_algebra.hxx><br>
00960         Namespaces: vigra and vigra::linalg
00961      */
00962 template <class T, class C1, class C2, class C3>
00963 bool qrDecomposition(MultiArrayView<2, T, C1> const & a,
00964                      MultiArrayView<2, T, C2> &q, MultiArrayView<2, T, C3> &r,
00965                      double epsilon = 0.0)
00966 {
00967     const MultiArrayIndex m = rowCount(a);
00968     const MultiArrayIndex n = columnCount(a);
00969     vigra_precondition(n == columnCount(r) && m == rowCount(r) &&
00970                        m == columnCount(q) && m == rowCount(q),
00971                        "qrDecomposition(): Matrix shape mismatch.");
00972 
00973     q = identityMatrix<T>(m);
00974     MultiArrayView<2,T, StridedArrayTag> tq = transpose(q);
00975     r = a;
00976     ArrayVector<MultiArrayIndex> noPivoting; // intentionally empty
00977     return ((MultiArrayIndex)detail::qrTransformToUpperTriangular(r, tq, noPivoting, epsilon) == std::min(m,n));
00978 }
00979 
00980     /** Deprecated, use \ref linearSolveUpperTriangular().
00981      */
00982 template <class T, class C1, class C2, class C3>
00983 inline 
00984 bool reverseElimination(const MultiArrayView<2, T, C1> &r, const MultiArrayView<2, T, C2> &b,
00985                         MultiArrayView<2, T, C3> x)
00986 {
00987     return linearSolveUpperTriangular(r, b, x);
00988 }
00989 
00990     /** Solve a linear system with upper-triangular coefficient matrix.
00991 
00992         The square matrix \a r must be an upper-triangular coefficient matrix as can,
00993         for example, be obtained by means of QR decomposition. If \a r doesn't have full rank
00994         the function fails and returns <tt>false</tt>, otherwise it returns <tt>true</tt>. The 
00995         lower triangular part of matrix \a r will not be touched, so it doesn't need to contain zeros.
00996         
00997         The column vectors of matrix \a b are the right-hand sides of the equation (several equations
00998         with the same coefficients can thus be solved in one go). The result is returned
00999         int \a x, whose columns contain the solutions for the corresponding
01000         columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed.
01001         The following size requirements apply:
01002         
01003         \code
01004         rowCount(r) == columnCount(r);
01005         rowCount(r) == rowCount(b);
01006         columnCount(r) == rowCount(x);
01007         columnCount(b) == columnCount(x);
01008         \endcode
01009 
01010     <b>\#include</b> <vigra/linear_solve.hxx> or<br>
01011     <b>\#include</b> <vigra/linear_algebra.hxx><br>
01012         Namespaces: vigra and vigra::linalg
01013      */
01014 template <class T, class C1, class C2, class C3>
01015 bool linearSolveUpperTriangular(const MultiArrayView<2, T, C1> &r, const MultiArrayView<2, T, C2> &b,
01016                                 MultiArrayView<2, T, C3> x)
01017 {
01018     typedef MultiArrayShape<2>::type Shape;
01019     MultiArrayIndex m = rowCount(r);
01020     MultiArrayIndex rhsCount = columnCount(b);
01021     vigra_precondition(m == columnCount(r),
01022         "linearSolveUpperTriangular(): square coefficient matrix required.");
01023     vigra_precondition(m == rowCount(b) && m == rowCount(x) && rhsCount == columnCount(x),
01024         "linearSolveUpperTriangular(): matrix shape mismatch.");
01025 
01026     for(MultiArrayIndex k = 0; k < rhsCount; ++k)
01027     {
01028         for(int i=m-1; i>=0; --i)
01029         {
01030             if(r(i,i) == NumericTraits<T>::zero())
01031                 return false;  // r doesn't have full rank
01032             T sum = b(i, k);
01033             for(MultiArrayIndex j=i+1; j<m; ++j)
01034                  sum -= r(i, j) * x(j, k);
01035             x(i, k) = sum / r(i, i);
01036         }
01037     }
01038     return true;
01039 }
01040 
01041     /** Solve a linear system with lower-triangular coefficient matrix.
01042 
01043         The square matrix \a l must be a lower-triangular coefficient matrix. If \a l 
01044         doesn't have full rank the function fails and returns <tt>false</tt>, 
01045         otherwise it returns <tt>true</tt>. The upper triangular part of matrix \a l will not be touched, 
01046         so it doesn't need to contain zeros.
01047         
01048         The column vectors of matrix \a b are the right-hand sides of the equation (several equations
01049         with the same coefficients can thus be solved in one go). The result is returned
01050         in \a x, whose columns contain the solutions for the corresponding
01051         columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed.
01052         The following size requirements apply:
01053         
01054         \code
01055         rowCount(l) == columnCount(l);
01056         rowCount(l) == rowCount(b);
01057         columnCount(l) == rowCount(x);
01058         columnCount(b) == columnCount(x);
01059         \endcode
01060 
01061     <b>\#include</b> <vigra/linear_solve.hxx> or<br>
01062     <b>\#include</b> <vigra/linear_algebra.hxx><br>
01063         Namespaces: vigra and vigra::linalg
01064      */
01065 template <class T, class C1, class C2, class C3>
01066 bool linearSolveLowerTriangular(const MultiArrayView<2, T, C1> &l, const MultiArrayView<2, T, C2> &b,
01067                             MultiArrayView<2, T, C3> x)
01068 {
01069     MultiArrayIndex m = columnCount(l);
01070     MultiArrayIndex n = columnCount(b);
01071     vigra_precondition(m == rowCount(l),
01072         "linearSolveLowerTriangular(): square coefficient matrix required.");
01073     vigra_precondition(m == rowCount(b) && m == rowCount(x) && n == columnCount(x),
01074         "linearSolveLowerTriangular(): matrix shape mismatch.");
01075 
01076     for(MultiArrayIndex k = 0; k < n; ++k)
01077     {
01078         for(MultiArrayIndex i=0; i<m; ++i)
01079         {
01080             if(l(i,i) == NumericTraits<T>::zero())
01081                 return false;  // l doesn't have full rank
01082             T sum = b(i, k);
01083             for(MultiArrayIndex j=0; j<i; ++j)
01084                  sum -= l(i, j) * x(j, k);
01085             x(i, k) = sum / l(i, i);
01086         }
01087     }
01088     return true;
01089 }
01090 
01091 
01092     /** Solve a linear system when the Cholesky decomposition of the left hand side is given.
01093 
01094         The square matrix \a L must be a lower-triangular matrix resulting from Cholesky
01095         decomposition of some positive definite coefficient matrix.
01096         
01097         The column vectors of matrix \a b are the right-hand sides of the equation (several equations
01098         with the same matrix \a L can thus be solved in one go). The result is returned
01099         in \a x, whose columns contain the solutions for the corresponding
01100         columns of \a b. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed.
01101         The following size requirements apply:
01102         
01103         \code
01104         rowCount(L) == columnCount(L);
01105         rowCount(L) == rowCount(b);
01106         columnCount(L) == rowCount(x);
01107         columnCount(b) == columnCount(x);
01108         \endcode
01109 
01110     <b>\#include</b> <vigra/linear_solve.hxx> or<br>
01111     <b>\#include</b> <vigra/linear_algebra.hxx><br>
01112         Namespaces: vigra and vigra::linalg
01113      */
01114 template <class T, class C1, class C2, class C3>
01115 inline 
01116 void choleskySolve(MultiArrayView<2, T, C1> & L, MultiArrayView<2, T, C2> const & b, MultiArrayView<2, T, C3> & x)
01117 {
01118     /* Solve L * y = b */
01119     linearSolveLowerTriangular(L, b, x);
01120     /* Solve L^T * x = y */
01121     linearSolveUpperTriangular(transpose(L), x, x);
01122 }
01123 
01124     /** Solve a linear system.
01125 
01126         \a A is the coefficient matrix, and the column vectors
01127         in \a b are the right-hand sides of the equation (so, several equations
01128         with the same coefficients can be solved in one go). The result is returned 
01129         in \a res, whose columns contain the solutions for the corresponding
01130         columns of \a b. The number of columns of \a A must equal the number of rows of
01131         both \a b and \a res, and the number of columns of \a b and \a res must match. 
01132         
01133         \a method must be one of the following:
01134         <DL>
01135         <DT>"Cholesky"<DD> Compute the solution by means of Cholesky decomposition. The 
01136                            coefficient matrix \a A must by symmetric positive definite. If
01137                            this is not the case, the function returns <tt>false</tt>.
01138                            
01139         <DT>"QR"<DD> (default) Compute the solution by means of QR decomposition.  The 
01140                            coefficient matrix \a A can be square or rectangular. In the latter case,
01141                            it must have more rows than columns, and the solution will be computed in the 
01142                            least squares sense. If \a A doesn't have full rank, the function 
01143                            returns <tt>false</tt>.
01144 
01145         <DT>"SVD"<DD> Compute the solution by means of singular value decomposition.  The 
01146                            coefficient matrix \a A can be square or rectangular. In the latter case,
01147                            it must have more rows than columns, and the solution will be computed in the 
01148                            least squares sense. If \a A doesn't have full rank, the function 
01149                            returns <tt>false</tt>.
01150 
01151         <DT>"NE"<DD> Compute the solution by means of the normal equations, i.e. by applying Cholesky
01152                            decomposition to the equivalent problem <tt>A'*A*x = A'*b</tt>. This only makes sense
01153                            when the equation is to be solved in the least squares sense, i.e. when \a A is a 
01154                            rectangular matrix with more rows than columns. If \a A doesn't have full column rank, 
01155                            the function returns <tt>false</tt>.
01156         </DL>
01157         
01158         This function can be applied in-place, i.e. <tt>&b == &res</tt> or <tt>&A == &res</tt> are allowed
01159         (provided they have the required shapes).
01160 
01161         The following size requirements apply:
01162         
01163         \code
01164         rowCount(r) == rowCount(b);
01165         columnCount(r) == rowCount(x);
01166         columnCount(b) == columnCount(x);
01167         \endcode
01168 
01169     <b>\#include</b> <vigra/linear_solve.hxx> or<br>
01170     <b>\#include</b> <vigra/linear_algebra.hxx><br>
01171         Namespaces: vigra and vigra::linalg
01172      */
01173 template <class T, class C1, class C2, class C3>
01174 bool linearSolve(const MultiArrayView<2, T, C1> &A, const MultiArrayView<2, T, C2> &b,
01175                  MultiArrayView<2, T, C3> & res, std::string method = "QR")
01176 {
01177     typedef typename Matrix<T>::difference_type Shape;
01178     typedef typename Matrix<T>::view_type SubMatrix;
01179     
01180     const MultiArrayIndex n = columnCount(A);
01181     const MultiArrayIndex m = rowCount(A);
01182 
01183     vigra_precondition(n <= m,
01184         "linearSolve(): Coefficient matrix A must have at least as many rows as columns.");
01185     vigra_precondition(n == rowCount(res) && 
01186                        m == rowCount(b) && columnCount(b) == columnCount(res),
01187         "linearSolve(): matrix shape mismatch.");
01188 
01189     for(unsigned int k=0; k<method.size(); ++k)
01190         method[k] = (std::string::value_type)tolower(method[k]);
01191     
01192     if(method == "cholesky")
01193     {
01194         vigra_precondition(columnCount(A) == rowCount(A),
01195             "linearSolve(): Cholesky method requires square coefficient matrix.");
01196         Matrix<T> L(A.shape());
01197         if(!choleskyDecomposition(A, L))
01198             return false; // false if A wasn't symmetric positive definite
01199         choleskySolve(L, b, res);
01200     }
01201     else if(method == "qr")
01202     {
01203         return (MultiArrayIndex)linearSolveQR(A, b, res) == n;
01204     }
01205     else if(method == "ne")
01206     {
01207         return linearSolve(transpose(A)*A, transpose(A)*b, res, "Cholesky");
01208     }
01209     else if(method == "svd")
01210     {
01211         MultiArrayIndex rhsCount = columnCount(b);
01212         Matrix<T> u(A.shape()), s(n, 1), v(n, n);
01213 
01214         MultiArrayIndex rank = (MultiArrayIndex)singularValueDecomposition(A, u, s, v);
01215 
01216         Matrix<T> t = transpose(u)*b;
01217         for(MultiArrayIndex l=0; l<rhsCount; ++l)
01218         {
01219             for(MultiArrayIndex k=0; k<rank; ++k)
01220                 t(k,l) /= s(k,0);
01221             for(MultiArrayIndex k=rank; k<n; ++k)
01222                 t(k,l) = NumericTraits<T>::zero();
01223         }
01224         res = v*t;
01225 
01226         return (rank == n);
01227     }
01228     else
01229     {
01230         vigra_precondition(false, "linearSolve(): Unknown solution method.");
01231     }
01232     return true;
01233 }
01234 
01235 //@}
01236 
01237 } // namespace linalg
01238 
01239 using linalg::inverse;
01240 using linalg::determinant;
01241 using linalg::logDeterminant;
01242 using linalg::linearSolve;
01243 using linalg::choleskySolve;
01244 using linalg::choleskyDecomposition;
01245 using linalg::qrDecomposition;
01246 using linalg::linearSolveUpperTriangular;
01247 using linalg::linearSolveLowerTriangular;
01248 
01249 } // namespace vigra
01250 
01251 
01252 #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.8.0 (20 Sep 2011)