[ 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.5.0, Dec 07 2006 ) */ 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 "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 upper 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 for(unsigned int k=0; k<acols; ++k) 00275 if(r(k,k) == NumericTraits<T>::zero()) 00276 return false; // a didn't have full rank. 00277 q.transpose(); 00278 reverseElimination(r, q * b, res); 00279 return true; 00280 } 00281 00282 //@} 00283 00284 } // namespace linalg 00285 00286 using linalg::inverse; 00287 using linalg::linearSolve; 00288 using linalg::qrDecomposition; 00289 using linalg::reverseElimination; 00290 00291 } // namespace vigra 00292 00293 00294 #endif // VIGRA_LINEAR_SOLVE_HXX
© Ullrich Köthe (koethe@informatik.uni-hamburg.de) |
html generated using doxygen and Python
|