[ 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.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)