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

html generated using doxygen and Python
VIGRA 1.6.0 (13 Aug 2008)