[ 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.1 (3 Dec 2010)