vigra - vigra/linear_solve.hxx Source File

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