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