[ 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) |
html generated using doxygen and Python
|