[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 2008 by Ullrich Koethe */ 00004 /* */ 00005 /* This file is part of the VIGRA computer vision library. */ 00006 /* ( Version 1.6.0, Aug 13 2008 ) */ 00007 /* The VIGRA Website is */ 00008 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00009 /* Please direct questions, bug reports, and contributions to */ 00010 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00011 /* vigra@informatik.uni-hamburg.de */ 00012 /* */ 00013 /* Permission is hereby granted, free of charge, to any person */ 00014 /* obtaining a copy of this software and associated documentation */ 00015 /* files (the "Software"), to deal in the Software without */ 00016 /* restriction, including without limitation the rights to use, */ 00017 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00018 /* sell copies of the Software, and to permit persons to whom the */ 00019 /* Software is furnished to do so, subject to the following */ 00020 /* conditions: */ 00021 /* */ 00022 /* The above copyright notice and this permission notice shall be */ 00023 /* included in all copies or substantial portions of the */ 00024 /* Software. */ 00025 /* */ 00026 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00027 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00028 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00029 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00030 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00031 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00032 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00033 /* OTHER DEALINGS IN THE SOFTWARE. */ 00034 /* */ 00035 /************************************************************************/ 00036 00037 00038 #ifndef VIGRA_REGRESSION_HXX 00039 #define VIGRA_REGRESSION_HXX 00040 00041 #include "matrix.hxx" 00042 #include "linear_solve.hxx" 00043 #include "singular_value_decomposition.hxx" 00044 #include "numerictraits.hxx" 00045 #include "functorexpression.hxx" 00046 00047 00048 namespace vigra 00049 { 00050 00051 namespace linalg 00052 { 00053 00054 /** \addtogroup MatrixAlgebra 00055 */ 00056 //@{ 00057 /** Ordinary Least Squares Regression. 00058 00059 Given a matrix \a A with <tt>m</tt> rows and <tt>n</tt> columns (with <tt>m >= n</tt>), 00060 and a column vector \a b of length <tt>m</tt> rows, this function computes 00061 a column vector \a x of length <tt>n</tt> rows such that the residual 00062 00063 \f[ \left|\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right|^2 00064 \f] 00065 00066 is minimized. When \a b is a matrix with <tt>k</tt> columns, \a x must also have 00067 <tt>k</tt> columns, which will contain the solutions for the corresponding columns of 00068 \a b. Note that all matrices must already have the correct shape. 00069 00070 This function is just another name for \ref linearSolve(), perhaps 00071 leading to more readable code when \a A is a rectangular matrix. It returns 00072 <tt>false</tt> when the rank of \a A is less than <tt>n</tt>. 00073 See \ref linearSolve() for more documentation. 00074 00075 <b>\#include</b> <<a href="regression_8hxx-source.html">vigra/regression_8hxx.hxx</a>> or<br> 00076 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00077 Namespaces: vigra and vigra::linalg 00078 */ 00079 template <class T, class C1, class C2, class C3> 00080 inline bool 00081 leastSquares(MultiArrayView<2, T, C1> const & A, 00082 MultiArrayView<2, T, C2> const &b, MultiArrayView<2, T, C3> &x, 00083 std::string method = "QR") 00084 { 00085 return linearSolve(A, b, x, method); 00086 } 00087 00088 /** Weighted Least Squares Regression. 00089 00090 Given a matrix \a A with <tt>m</tt> rows and <tt>n</tt> columns (with <tt>m >= n</tt>), 00091 a vector \a b of length <tt>m</tt>, and a weight vector \a weights of length <tt>m</tt> 00092 with non-negative entries, this function computes a vector \a x of length <tt>n</tt> 00093 such that the weighted residual 00094 00095 \f[ \left(\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right)^T 00096 \textrm{diag}(\textrm{\bf weights}) 00097 \left(\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right) 00098 \f] 00099 00100 is minimized, where <tt>diag(weights)</tt> creates a diagonal matrix from \a weights. 00101 The algorithm calls \ref leastSquares() on the equivalent problem 00102 00103 \f[ \left|\textrm{diag}(\textrm{\bf weights})^{1/2}\textrm{\bf A} \textrm{\bf x} - 00104 \textrm{diag}(\textrm{\bf weights})^{1/2} \textrm{\bf b}\right|^2 00105 \f] 00106 00107 where the square root of \a weights is just taken element-wise. 00108 00109 When \a b is a matrix with <tt>k</tt> columns, \a x must also have 00110 <tt>k</tt> columns, which will contain the solutions for the corresponding columns of 00111 \a b. Note that all matrices must already have the correct shape. 00112 00113 The function returns 00114 <tt>false</tt> when the rank of the weighted matrix \a A is less than <tt>n</tt>. 00115 00116 <b>\#include</b> <<a href="regression_8hxx-source.html">vigra/regression.hxx</a>> or<br> 00117 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00118 Namespaces: vigra and vigra::linalg 00119 */ 00120 template <class T, class C1, class C2, class C3, class C4> 00121 bool 00122 weightedLeastSquares(MultiArrayView<2, T, C1> const & A, 00123 MultiArrayView<2, T, C2> const &b, MultiArrayView<2, T, C3> const &weights, 00124 MultiArrayView<2, T, C4> &x, std::string method = "QR") 00125 { 00126 typedef T Real; 00127 00128 const unsigned int rows = rowCount(A); 00129 const unsigned int cols = columnCount(A); 00130 const unsigned int rhsCount = columnCount(b); 00131 vigra_precondition(rows >= cols, 00132 "weightedLeastSquares(): Input matrix A must be rectangular with rowCount >= columnCount."); 00133 vigra_precondition(rowCount(b) == rows, 00134 "weightedLeastSquares(): Shape mismatch between matrices A and b."); 00135 vigra_precondition(rowCount(b) == rowCount(weights) && columnCount(weights) == 1, 00136 "weightedLeastSquares(): Weight matrix has wrong shape."); 00137 vigra_precondition(rowCount(x) == cols && columnCount(x) == rhsCount, 00138 "weightedLeastSquares(): Result matrix x has wrong shape."); 00139 00140 Matrix<T> wa(A.shape()), wb(b.shape()); 00141 00142 for(unsigned int k=0; k<rows; ++k) 00143 { 00144 vigra_precondition(weights(k,0) >= 0, 00145 "weightedLeastSquares(): Weights must be positive."); 00146 T w = std::sqrt(weights(k,0)); 00147 for(unsigned int l=0; l<cols; ++l) 00148 wa(k,l) = w * A(k,l); 00149 for(unsigned int l=0; l<rhsCount; ++l) 00150 wb(k,l) = w * b(k,l); 00151 } 00152 00153 return leastSquares(wa, wb, x, method); 00154 } 00155 00156 /** Ridge Regression. 00157 00158 Given a matrix \a A with <tt>m</tt> rows and <tt>n</tt> columns (with <tt>m >= n</tt>), 00159 a vector \a b of length <tt>m</tt>, and a regularization parameter <tt>lambda >= 0.0</tt>, 00160 this function computes a vector \a x of length <tt>n</tt> such that the residual 00161 00162 \f[ \left|\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right|^2 + 00163 \lambda \textrm{\bf x}^T\textrm{\bf x} 00164 \f] 00165 00166 is minimized. This is implemented by means of \ref singularValueDecomposition(). 00167 00168 When \a b is a matrix with <tt>k</tt> columns, \a x must also have 00169 <tt>k</tt> columns, which will contain the solutions for the corresponding columns of 00170 \a b. Note that all matrices must already have the correct shape. 00171 00172 The function returns <tt>false</tt> if the rank of \a A is less than <tt>n</tt> 00173 and <tt>lambda == 0.0</tt>. 00174 00175 <b>\#include</b> <<a href="regression_8hxx-source.html">vigra/regression.hxx</a>> or<br> 00176 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00177 Namespaces: vigra and vigra::linalg 00178 */ 00179 template <class T, class C1, class C2, class C3> 00180 bool 00181 ridgeRegression(MultiArrayView<2, T, C1> const & A, 00182 MultiArrayView<2, T, C2> const &b, MultiArrayView<2, T, C3> &x, double lambda) 00183 { 00184 typedef T Real; 00185 00186 const unsigned int rows = rowCount(A); 00187 const unsigned int cols = columnCount(A); 00188 const unsigned int rhsCount = columnCount(b); 00189 vigra_precondition(rows >= cols, 00190 "ridgeRegression(): Input matrix A must be rectangular with rowCount >= columnCount."); 00191 vigra_precondition(rowCount(b) == rows, 00192 "ridgeRegression(): Shape mismatch between matrices A and b."); 00193 vigra_precondition(rowCount(x) == cols && columnCount(x) == rhsCount, 00194 "ridgeRegression(): Result matrix x has wrong shape."); 00195 vigra_precondition(lambda >= 0.0, 00196 "ridgeRegression(): lambda >= 0.0 required."); 00197 00198 unsigned int m = rows; 00199 unsigned int n = cols; 00200 00201 Matrix<T> u(m, n), s(n, 1), v(n, n); 00202 00203 unsigned int rank = singularValueDecomposition(A, u, s, v); 00204 if(rank < n && lambda == 0.0) 00205 return false; 00206 00207 Matrix<T> t = transpose(u)*b; 00208 for(unsigned int k=0; k<cols; ++k) 00209 for(unsigned int l=0; l<rhsCount; ++l) 00210 t(k,l) *= s(k,0) / (sq(s(k,0)) + lambda); 00211 x = v*t; 00212 return true; 00213 } 00214 00215 /** Weighted ridge Regression. 00216 00217 Given a matrix \a A with <tt>m</tt> rows and <tt>n</tt> columns (with <tt>m >= n</tt>), 00218 a vector \a b of length <tt>m</tt>, a weight vector \a weights of length <tt>m</tt> 00219 with non-negative entries, and a regularization parameter <tt>lambda >= 0.0</tt> 00220 this function computes a vector \a x of length <tt>n</tt> such that the weighted residual 00221 00222 \f[ \left(\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right)^T 00223 \textrm{diag}(\textrm{\bf weights}) 00224 \left(\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right) + 00225 \lambda \textrm{\bf x}^T\textrm{\bf x} 00226 \f] 00227 00228 is minimized, where <tt>diag(weights)</tt> creates a diagonal matrix from \a weights. 00229 The algorithm calls \ref ridgeRegression() on the equivalent problem 00230 00231 \f[ \left|\textrm{diag}(\textrm{\bf weights})^{1/2}\textrm{\bf A} \textrm{\bf x} - 00232 \textrm{diag}(\textrm{\bf weights})^{1/2} \textrm{\bf b}\right|^2 + 00233 \lambda \textrm{\bf x}^T\textrm{\bf x} 00234 \f] 00235 00236 where the square root of \a weights is just taken element-wise. This solution is 00237 computed by means of \ref singularValueDecomposition(). 00238 00239 When \a b is a matrix with <tt>k</tt> columns, \a x must also have 00240 <tt>k</tt> columns, which will contain the solutions for the corresponding columns of 00241 \a b. Note that all matrices must already have the correct shape. 00242 00243 The function returns <tt>false</tt> if the rank of \a A is less than <tt>n</tt> 00244 and <tt>lambda == 0.0</tt>. 00245 00246 <b>\#include</b> <<a href="regression_8hxx-source.html">vigra/regression.hxx</a>> or<br> 00247 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00248 Namespaces: vigra and vigra::linalg 00249 */ 00250 template <class T, class C1, class C2, class C3, class C4> 00251 bool 00252 weightedRidgeRegression(MultiArrayView<2, T, C1> const & A, 00253 MultiArrayView<2, T, C2> const &b, MultiArrayView<2, T, C3> const &weights, 00254 MultiArrayView<2, T, C4> &x, double lambda) 00255 { 00256 typedef T Real; 00257 00258 const unsigned int rows = rowCount(A); 00259 const unsigned int cols = columnCount(A); 00260 const unsigned int rhsCount = columnCount(b); 00261 vigra_precondition(rows >= cols, 00262 "weightedRidgeRegression(): Input matrix A must be rectangular with rowCount >= columnCount."); 00263 vigra_precondition(rowCount(b) == rows, 00264 "weightedRidgeRegression(): Shape mismatch between matrices A and b."); 00265 vigra_precondition(rowCount(b) == rowCount(weights) && columnCount(weights) == 1, 00266 "weightedRidgeRegression(): Weight matrix has wrong shape."); 00267 vigra_precondition(rowCount(x) == cols && columnCount(x) == rhsCount, 00268 "weightedRidgeRegression(): Result matrix x has wrong shape."); 00269 vigra_precondition(lambda >= 0.0, 00270 "weightedRidgeRegression(): lambda >= 0.0 required."); 00271 00272 Matrix<T> wa(A.shape()), wb(b.shape()); 00273 00274 for(unsigned int k=0; k<rows; ++k) 00275 { 00276 vigra_precondition(weights(k,0) >= 0, 00277 "weightedRidgeRegression(): Weights must be positive."); 00278 T w = std::sqrt(weights(k,0)); 00279 for(unsigned int l=0; l<cols; ++l) 00280 wa(k,l) = w * A(k,l); 00281 for(unsigned int l=0; l<rhsCount; ++l) 00282 wb(k,l) = w * b(k,l); 00283 } 00284 00285 return ridgeRegression(wa, wb, x, lambda); 00286 } 00287 00288 /** Ridge Regression with many lambdas. 00289 00290 This executes \ref ridgeRegression() for a sequence of regularization parameters. This 00291 is implemented so that the \ref singularValueDecomposition() has to be executed only once. 00292 \a lambda must be an array conforming to the <tt>std::vector</tt> interface, i.e. must 00293 support <tt>lambda.size()</tt> and <tt>lambda[k]</tt>. The columns of the matrix \a x 00294 will contain the solutions for the corresponding lambda, so the number of columns of 00295 the matrix \a x must be equal to <tt>lambda.size()</tt>, and \a b must be a columns vector, 00296 i.e. cannot contain several right hand sides at once. 00297 00298 The function returns <tt>false</tt> when the matrix \a A is rank deficient. If this 00299 happens, and one of the lambdas is zero, the corresponding column of \a x will be skipped. 00300 00301 <b>\#include</b> <<a href="regression_8hxx-source.html">vigra/regression.hxx</a>> or<br> 00302 <b>\#include</b> <<a href="linear__algebra_8hxx-source.html">vigra/linear_algebra.hxx</a>><br> 00303 Namespaces: vigra and vigra::linalg 00304 */ 00305 template <class T, class C1, class C2, class C3, class Array> 00306 bool 00307 ridgeRegressionSeries(MultiArrayView<2, T, C1> const & A, 00308 MultiArrayView<2, T, C2> const &b, MultiArrayView<2, T, C3> &x, Array const & lambda) 00309 { 00310 typedef T Real; 00311 00312 const unsigned int rows = rowCount(A); 00313 const unsigned int cols = columnCount(A); 00314 const unsigned int lambdaCount = lambda.size(); 00315 vigra_precondition(rows >= cols, 00316 "ridgeRegressionSeries(): Input matrix A must be rectangular with rowCount >= columnCount."); 00317 vigra_precondition(rowCount(b) == rows && columnCount(b) == 1, 00318 "ridgeRegressionSeries(): Shape mismatch between matrices A and b."); 00319 vigra_precondition(rowCount(x) == cols && columnCount(x) == lambdaCount, 00320 "ridgeRegressionSeries(): Result matrix x has wrong shape."); 00321 00322 unsigned int m = rows; 00323 unsigned int n = cols; 00324 00325 Matrix<T> u(m, n), s(n, 1), v(n, n); 00326 00327 unsigned int rank = singularValueDecomposition(A, u, s, v); 00328 00329 Matrix<T> xl = transpose(u)*b; 00330 Matrix<T> xt(cols,1); 00331 for(unsigned int i=0; i<lambdaCount; ++i) 00332 { 00333 vigra_precondition(lambda[i] >= 0.0, 00334 "ridgeRegressionSeries(): lambda >= 0.0 required."); 00335 if(lambda == 0.0 && rank < rows) 00336 continue; 00337 for(unsigned int k=0; k<cols; ++k) 00338 xt(k,0) = xl(k,0) * s(k,0) / (sq(s(k,0)) + lambda[i]); 00339 columnVector(x, i) = v*xt; 00340 } 00341 return (rank < n); 00342 } 00343 00344 } // namespace linalg 00345 00346 using linalg::leastSquares; 00347 using linalg::weightedLeastSquares; 00348 using linalg::ridgeRegression; 00349 using linalg::weightedRidgeRegression; 00350 using linalg::ridgeRegressionSeries; 00351 00352 } // namespace vigra 00353 00354 #endif // VIGRA_REGRESSION_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|