[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

vigra/regression.hxx

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)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
VIGRA 1.6.0 (5 Nov 2009)