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

| vigra/linear_solve.hxx |
--- | --- | ---

    00001 /************************************************************************/
    00002 /*                                                                      */
    00003 /*        Copyright 2004 by Gunnar Kedenburg and Ullrich Koethe         */
    00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
    00005 /*                                                                      */
    00006 /*    This file is part of the VIGRA computer vision library.           */
    00007 /*    ( Version 1.3.2, Jan 27 2005 )                                    */
    00008 /*    You may use, modify, and distribute this software according       */
    00009 /*    to the terms stated in the LICENSE file included in               */
    00010 /*    the VIGRA distribution.                                           */
    00011 /*                                                                      */
    00012 /*    The VIGRA Website is                                              */
    00013 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
    00014 /*    Please direct questions, bug reports, and contributions to        */
    00015 /*        koethe@informatik.uni-hamburg.de                              */
    00016 /*                                                                      */
    00017 /*  THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR          */
    00018 /*  IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED      */
    00019 /*  WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */
    00020 /*                                                                      */
    00021 /************************************************************************/
    00022
    00023
    00024 #ifndef VIGRA_LINEAR_SOLVE_HXX
    00025 #define VIGRA_LINEAR_SOLVE_HXX
    00026
    00027 #include "vigra/matrix.hxx"
    00028
    00029
    00030 namespace vigra
    00031 {
    00032
    00033 namespace linalg
    00034 {
    00035
    00036 /** \addtogroup LinearAlgebraFunctions Matrix functions
    00037  */
    00038 //@{
    00039     /** invert square matrix \a v.
    00040         The result is written into \a r which must have the same shape.
    00041         The inverse is calculated by means of QR decomposition. If \a v
    00042         is not invertible, <tt>vigra::PreconditionViolation</tt> exception is thrown.
    00043
    00044     <b>\#include</b> "<a href="linear__solve_8hxx-source.html">vigra/linear__solve.hxx</a>" or<br>
    00045     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
    00046         Namespaces: vigra and vigra::linalg
    00047      */
    00048 template <class T, class C1, class C2>
    00049 void inverse(const MultiArrayView<2, T, C1> &v, MultiArrayView<2, T, C2> &r)
    00050 {
    00051     const unsigned int n = rowCount(r);
    00052     vigra_precondition(n == columnCount(v) && n == rowCount(v) && n == columnCount(r),
    00053        "inverse(): matrices must be square.");
    00054     vigra_precondition(linearSolve(v, identityMatrix<T>(n), r),
    00055         "inverse(): matrix is not invertible.");
    00056 }
    00057
    00058     /** create the inverse of square matrix \a v.
    00059         The result is returned as a temporary matrix.
    00060         The inverse is calculated by means of QR decomposition. If \a v
    00061         is not invertible, <tt>vigra::PreconditionViolation</tt> exception is thrown.
    00062         Usage:
    00063
    00064         \code
    00065         vigra::Matrix<double> v(n, n);
    00066         v = ...;
    00067
    00068         vigra::Matrix<double> m = inverse(v);
    00069         \endcode
    00070
    00071     <b>\#include</b> "<a href="linear__solve_8hxx-source.html">vigra/linear__solve.hxx</a>" or<br>
    00072     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
    00073         Namespaces: vigra and vigra::linalg
    00074      */
    00075 template <class T, class C>
    00076 TemporaryMatrix<T> inverse(const MultiArrayView<2, T, C> &v)
    00077 {
    00078     const unsigned int n = rowCount(v);
    00079     vigra_precondition(n == columnCount(v),
    00080        "inverse(): matrix must be square.");
    00081     TemporaryMatrix<T> ret = identityMatrix<T>(n);
    00082     vigra_precondition(linearSolve(v, ret, ret),
    00083         "inverse(): matrix is not invertible.");
    00084     return ret;
    00085 }
    00086
    00087 template <class T>
    00088 TemporaryMatrix<T> inverse(const TemporaryMatrix<T> &v)
    00089 {
    00090     const unsigned int n = v.rowCount();
    00091     vigra_precondition(n == v.columnCount(),
    00092        "inverse(): matrix must be square.");
    00093     vigra_precondition(linearSolve(v, identityMatrix<T>(n), v),
    00094         "inverse(): matrix is not invertible.");
    00095     return v;
    00096 }
    00097
    00098     /** QR decomposition.
    00099
    00100         \a a contains the original matrix, results are returned in \a q and \a r, where
    00101         \a q is a orthogonal matrix, and \a r is an uppr triangular matrix, and
    00102         the following relation holds:
    00103
    00104         \code
    00105         assert(a == q * r);
    00106         \endcode
    00107
    00108         This implementation uses householder transformations. It can be applied in-place,
    00109         i.e. <tt>&a == &r</tt> is allowed.
    00110
    00111     <b>\#include</b> "<a href="linear__solve_8hxx-source.html">vigra/linear__solve.hxx</a>" or<br>
    00112     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
    00113         Namespaces: vigra and vigra::linalg
    00114      */
    00115 template <class T, class C1, class C2, class C3>
    00116 void qrDecomposition(MultiArrayView<2, T, C1> const & a,
    00117                      MultiArrayView<2, T, C2> &q, MultiArrayView<2, T, C3> &r)
    00118 {
    00119     typedef typename MultiArrayView<2, T, C2>::difference_type MatrixShape;
    00120     typedef typename MultiArray<1, T>::difference_type VectorShape;
    00121
    00122     // the orthogonal matrix q will have as many rows and columns as
    00123     // the original matrix has columns.
    00124     const unsigned int rows = rowCount(a);
    00125     const unsigned int cols = columnCount(a);
    00126     vigra_precondition(cols == columnCount(r) && cols == rowCount(r) &&
    00127                        cols == columnCount(q) && cols == rowCount(q),
    00128                        "qrDecomposition(): Matrix shape mismatch.");
    00129
    00130     identityMatrix(q);
    00131     r.copy(a);   // does nothing if &r == &a
    00132
    00133     for(unsigned int k = 0; (k < cols) && (k < rows - 1); ++k) {
    00134
    00135         const unsigned int rows_left = rows - k;
    00136         const unsigned int cols_left = cols - k;
    00137
    00138         // create a view on the remaining part of r
    00139         MatrixShape rul(k, k);
    00140         MultiArrayView<2, T, C2> rsub = r.subarray(rul, r.shape());
    00141
    00142         // decompose the first row
    00143         MultiArrayView <1, T, C2 > vec = rsub.bindOuter(0);
    00144
    00145         // defining householder vector
    00146         VectorShape ushape(rows_left);
    00147         MultiArray<1, T> u(ushape);
    00148         for(unsigned int i = 0; i < rows_left; ++i)
    00149             u(i) = vec(i);
    00150         u(0) += norm(vec);
    00151
    00152         const T divisor = squaredNorm(u);
    00153         const T scal = (divisor == 0) ? 0.0 : 2.0 / divisor;
    00154
    00155         // apply householder elimination on rsub
    00156         for(unsigned int i = 0; i < cols_left; ++i) {
    00157
    00158             // compute the inner product of the i'th column of rsub with u
    00159             T sum = dot(u, rsub.bindOuter(i));
    00160
    00161             // add rsub*(uu')/(u'u)
    00162             sum *= scal;
    00163             for(unsigned int j = 0; j < rows_left; ++j)
    00164                 rsub(j, i) -= sum * u(j);
    00165         }
    00166
    00167         MatrixShape qul(0, k);
    00168         MultiArrayView <2, T, C3 > qsub = q.subarray(qul, q.shape());
    00169
    00170         // apply the (self-inverse) householder matrix on q
    00171         for(unsigned int i = 0; i < cols; ++i) {
    00172
    00173             // compute the inner product of the i'th row of q with u
    00174             T sum = dot(qsub.bindInner(i), u);
    00175
    00176             // add q*(uu')/(u'u)
    00177             sum *= scal;
    00178             for(unsigned int j = 0; j < rows_left; ++j)
    00179                 qsub(i, j) -= sum * u(j);
    00180         }
    00181     }
    00182 }
    00183
    00184     /** Solve a linear system with right-triangular defining matrix.
    00185
    00186         The square matrix \a a must be a right-triangular coefficient matrix as can,
    00187         for example, be obtained by means of QR decomposition. The column vectors
    00188         in \a b are the right-hand sides of the equation (so, several equations
    00189         with the same coefficients can be solved in one go). The result is returned
    00190         int \a x, whose columns contain the solutions for the correspoinding
    00191         columns of \a b. The number of columns of \a a must equal the number of rows of
    00192         both \a b and \a x, and the number of columns of \a b and \a x must be
    00193         equal. This implementation can be applied in-place, i.e. <tt>&b == &x</tt> is allowed.
    00194
    00195     <b>\#include</b> "<a href="linear__solve_8hxx-source.html">vigra/linear__solve.hxx</a>" or<br>
    00196     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
    00197         Namespaces: vigra and vigra::linalg
    00198      */
    00199 template <class T, class C1, class C2, class C3>
    00200 void reverseElimination(const MultiArrayView<2, T, C1> &r, const MultiArrayView<2, T, C2> &b,
    00201                           MultiArrayView<2, T, C3> & x)
    00202 {
    00203     unsigned int m = columnCount(r);
    00204     unsigned int n = columnCount(b);
    00205     vigra_precondition(m == rowCount(r),
    00206         "reverseElimination(): square coefficient matrix required.");
    00207     vigra_precondition(m == rowCount(b) && m == rowCount(x) && n == columnCount(x),
    00208         "reverseElimination(): matrix shape mismatch.");
    00209
    00210     for(unsigned int k = 0; k < n; ++k)
    00211     {
    00212         x(m-1, k) = b(m-1, k) / r(m-1, m-1);
    00213         if(m >= 2)
    00214         {
    00215             for(int i = m-2; i >= 0; --i)
    00216             {
    00217                 // compute the i'th inner product, excluding the diagonal entry.
    00218                 T sum = NumericTraits<T>::zero();
    00219                 for(unsigned int j = i+1; j < m; ++j)
    00220                     sum += r(i, j) * x(j, k);
    00221                 if(r(i, i) != NumericTraits<T>::zero())
    00222                     x(i, k) = (b(i, k) - sum) / r(i, i);
    00223                 else
    00224                     x(i, k) = NumericTraits<T>::zero();
    00225             }
    00226         }
    00227     }
    00228 }
    00229
    00230     /** Solve a linear system.
    00231
    00232         The square matrix \a a is the coefficient matrix, and the column vectors
    00233         in \a b are the right-hand sides of the equation (so, several equations
    00234         with the same coefficients can be solved in one go). The result is returned
    00235         int \a res, whose columns contain the solutions for the correspoinding
    00236         columns of \a b. The number of columns of \a a must equal the number of rows of
    00237         both \a b and \a res, and the number of columns of \a b and \a res must be
    00238         equal. The algorithm uses QR decomposition of \a a. The algorithm returns
    00239         <tt>false</tt> if \a a doesn't have full rank. This implementation  can be
    00240         applied in-place, i.e. <tt>&b == &res</tt> or <tt>&a == &res</tt> are allowed.
    00241
    00242     <b>\#include</b> "<a href="linear__solve_8hxx-source.html">vigra/linear__solve.hxx</a>" or<br>
    00243     <b>\#include</b> "<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>"<br>
    00244         Namespaces: vigra and vigra::linalg
    00245      */
    00246 template <class T, class C1, class C2, class C3>
    00247 bool linearSolve(const MultiArrayView<2, T, C1> &a, const MultiArrayView<2, T, C2> &b,
    00248                  MultiArrayView<2, T, C3> & res)
    00249 {
    00250     unsigned int acols = columnCount(a);
    00251     unsigned int bcols = columnCount(b);
    00252     vigra_precondition(acols == rowCount(a),
    00253         "linearSolve(): square coefficient matrix required.");
    00254     vigra_precondition(acols == rowCount(b) && acols == rowCount(res) && bcols == columnCount(res),
    00255         "linearSolve(): matrix shape mismatch.");
    00256
    00257     Matrix<T> q(acols, acols), r(a);
    00258     qrDecomposition(r, q, r);
    00259     if(r(acols-1, acols-1) == NumericTraits<T>::zero())
    00260          return false; // a didn't have full rank.
    00261     q.transpose();
    00262     reverseElimination(r, q * b, res);
    00263     return true;
    00264 }
    00265
    00266 //@}
    00267
    00268 } // namespace linalg
    00269
    00270 using linalg::inverse;
    00271 using linalg::linearSolve;
    00272 using linalg::qrDecomposition;
    00273 using linalg::reverseElimination;
    00274
    00275 } // namespace vigra
    00276
    00277
    00278 #endif // VIGRA_LINEAR_SOLVE_HXX

© Ullrich Köthe (koethe@informatik.uni-hamburg.de) Cognitive Systems Group, University of Hamburg, Germany | html generated using doxygen and Python VIGRA 1.3.2 (27 Jan 2005)