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

details vigra/linear_solve.hxx VIGRA

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