[ 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://hci.iwr.uni-heidelberg.de/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 Optimization Optimization and Regression
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        the column vector \a x of length <tt>n</tt> rows that solves the optimization problem
00061 
00062         \f[ \tilde \textrm{\bf x} = \textrm{argmin}
00063             \left|\left|\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right|\right|_2^2
00064         \f]
00065 
00066        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.hxx</a>>
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 the vector \a x of length <tt>n</tt>
00092        that solves the optimization problem
00093 
00094         \f[ \tilde \textrm{\bf x} = \textrm{argmin}
00095             \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        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[ \tilde \textrm{\bf x} = \textrm{argmin}
00104              \left|\left|\textrm{diag}(\textrm{\bf weights})^{1/2}\textrm{\bf A} \textrm{\bf x} -
00105                   \textrm{diag}(\textrm{\bf weights})^{1/2} \textrm{\bf b}\right|\right|_2^2
00106         \f]
00107 
00108        where the square root of \a weights is just taken element-wise.
00109 
00110        When \a b is a matrix with <tt>k</tt> columns, \a x must also have
00111        <tt>k</tt> columns, which will contain the solutions for the corresponding columns of
00112        \a b. Note that all matrices must already have the correct shape.
00113 
00114        The function returns
00115        <tt>false</tt> when the rank of the weighted matrix \a A is less than <tt>n</tt>.
00116 
00117     <b>\#include</b> <<a href="regression_8hxx-source.html">vigra/regression.hxx</a>>
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 the vector \a x of length <tt>n</tt>
00161        that solves the optimization problem
00162 
00163         \f[ \tilde \textrm{\bf x} = \textrm{argmin}
00164             \left|\left|\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right|\right|_2^2 +
00165             \lambda \textrm{\bf x}^T\textrm{\bf x}
00166         \f]
00167 
00168        This is implemented by means of \ref singularValueDecomposition().
00169 
00170        When \a b is a matrix with <tt>k</tt> columns, \a x must also have
00171        <tt>k</tt> columns, which will contain the solutions for the corresponding columns of
00172        \a b. Note that all matrices must already have the correct shape.
00173 
00174        The function returns <tt>false</tt> if the rank of \a A is less than <tt>n</tt>
00175        and <tt>lambda == 0.0</tt>.
00176 
00177     <b>\#include</b> <<a href="regression_8hxx-source.html">vigra/regression.hxx</a>>
00178         Namespaces: vigra and vigra::linalg
00179    */
00180 template <class T, class C1, class C2, class C3>
00181 bool
00182 ridgeRegression(MultiArrayView<2, T, C1> const & A,
00183                 MultiArrayView<2, T, C2> const &b, MultiArrayView<2, T, C3> &x, double lambda)
00184 {
00185     typedef T Real;
00186 
00187     const unsigned int rows = rowCount(A);
00188     const unsigned int cols = columnCount(A);
00189     const unsigned int rhsCount = columnCount(b);
00190     vigra_precondition(rows >= cols,
00191        "ridgeRegression(): Input matrix A must be rectangular with rowCount >= columnCount.");
00192     vigra_precondition(rowCount(b) == rows,
00193        "ridgeRegression(): Shape mismatch between matrices A and b.");
00194     vigra_precondition(rowCount(x) == cols && columnCount(x) == rhsCount,
00195        "ridgeRegression(): Result matrix x has wrong shape.");
00196     vigra_precondition(lambda >= 0.0,
00197        "ridgeRegression(): lambda >= 0.0 required.");
00198 
00199     unsigned int m = rows;
00200     unsigned int n = cols;
00201 
00202     Matrix<T> u(m, n), s(n, 1), v(n, n);
00203 
00204     unsigned int rank = singularValueDecomposition(A, u, s, v);
00205     if(rank < n && lambda == 0.0)
00206         return false;
00207 
00208     Matrix<T> t = transpose(u)*b;
00209     for(unsigned int k=0; k<cols; ++k)
00210         for(unsigned int l=0; l<rhsCount; ++l)
00211             t(k,l) *= s(k,0) / (sq(s(k,0)) + lambda);
00212     x = v*t;
00213     return true;
00214 }
00215 
00216    /** Weighted ridge Regression.
00217 
00218        Given a matrix \a A with <tt>m</tt> rows and <tt>n</tt> columns (with <tt>m >= n</tt>),
00219        a vector \a b of length <tt>m</tt>, a weight vector \a weights of length <tt>m</tt>
00220        with non-negative entries, and a regularization parameter <tt>lambda >= 0.0</tt>
00221        this function computes the vector \a x of length <tt>n</tt>
00222        that solves the optimization problem
00223 
00224         \f[ \tilde \textrm{\bf x} = \textrm{argmin}
00225             \left(\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right)^T
00226              \textrm{diag}(\textrm{\bf weights})
00227              \left(\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right) +
00228              \lambda \textrm{\bf x}^T\textrm{\bf x}
00229         \f]
00230 
00231        where <tt>diag(weights)</tt> creates a diagonal matrix from \a weights.
00232        The algorithm calls \ref ridgeRegression() on the equivalent problem
00233 
00234         \f[ \tilde \textrm{\bf x} = \textrm{argmin}
00235             \left|\left|\textrm{diag}(\textrm{\bf weights})^{1/2}\textrm{\bf A} \textrm{\bf x} -
00236                   \textrm{diag}(\textrm{\bf weights})^{1/2} \textrm{\bf b}\right|\right|_2^2 +
00237              \lambda \textrm{\bf x}^T\textrm{\bf x}
00238         \f]
00239 
00240        where the square root of \a weights is just taken element-wise.  This solution is
00241        computed by means of \ref singularValueDecomposition().
00242 
00243        When \a b is a matrix with <tt>k</tt> columns, \a x must also have
00244        <tt>k</tt> columns, which will contain the solutions for the corresponding columns of
00245        \a b. Note that all matrices must already have the correct shape.
00246 
00247        The function returns <tt>false</tt> if the rank of \a A is less than <tt>n</tt>
00248        and <tt>lambda == 0.0</tt>.
00249 
00250     <b>\#include</b> <<a href="regression_8hxx-source.html">vigra/regression.hxx</a>>
00251         Namespaces: vigra and vigra::linalg
00252    */
00253 template <class T, class C1, class C2, class C3, class C4>
00254 bool
00255 weightedRidgeRegression(MultiArrayView<2, T, C1> const & A,
00256              MultiArrayView<2, T, C2> const &b, MultiArrayView<2, T, C3> const &weights,
00257              MultiArrayView<2, T, C4> &x, double lambda)
00258 {
00259     typedef T Real;
00260 
00261     const unsigned int rows = rowCount(A);
00262     const unsigned int cols = columnCount(A);
00263     const unsigned int rhsCount = columnCount(b);
00264     vigra_precondition(rows >= cols,
00265        "weightedRidgeRegression(): Input matrix A must be rectangular with rowCount >= columnCount.");
00266     vigra_precondition(rowCount(b) == rows,
00267        "weightedRidgeRegression(): Shape mismatch between matrices A and b.");
00268     vigra_precondition(rowCount(b) == rowCount(weights) && columnCount(weights) == 1,
00269        "weightedRidgeRegression(): Weight matrix has wrong shape.");
00270     vigra_precondition(rowCount(x) == cols && columnCount(x) == rhsCount,
00271        "weightedRidgeRegression(): Result matrix x has wrong shape.");
00272     vigra_precondition(lambda >= 0.0,
00273        "weightedRidgeRegression(): lambda >= 0.0 required.");
00274 
00275     Matrix<T> wa(A.shape()), wb(b.shape());
00276 
00277     for(unsigned int k=0; k<rows; ++k)
00278     {
00279         vigra_precondition(weights(k,0) >= 0,
00280            "weightedRidgeRegression(): Weights must be positive.");
00281         T w = std::sqrt(weights(k,0));
00282         for(unsigned int l=0; l<cols; ++l)
00283             wa(k,l) = w * A(k,l);
00284         for(unsigned int l=0; l<rhsCount; ++l)
00285             wb(k,l) = w * b(k,l);
00286     }
00287 
00288     return ridgeRegression(wa, wb, x, lambda);
00289 }
00290 
00291    /** Ridge Regression with many lambdas.
00292 
00293        This executes \ref ridgeRegression() for a sequence of regularization parameters. This
00294        is implemented so that the \ref singularValueDecomposition() has to be executed only once.
00295        \a lambda must be an array conforming to the <tt>std::vector</tt> interface, i.e. must
00296        support <tt>lambda.size()</tt> and <tt>lambda[k]</tt>. The columns of the matrix \a x
00297        will contain the solutions for the corresponding lambda, so the  number of columns of
00298        the matrix \a x must be equal to <tt>lambda.size()</tt>, and \a b must be a columns vector,
00299        i.e. cannot contain several right hand sides at once.
00300 
00301        The function returns <tt>false</tt> when the matrix \a A is rank deficient. If this
00302        happens, and one of the lambdas is zero, the corresponding column of \a x will be skipped.
00303 
00304     <b>\#include</b> <<a href="regression_8hxx-source.html">vigra/regression.hxx</a>>
00305         Namespaces: vigra and vigra::linalg
00306    */
00307 template <class T, class C1, class C2, class C3, class Array>
00308 bool
00309 ridgeRegressionSeries(MultiArrayView<2, T, C1> const & A,
00310           MultiArrayView<2, T, C2> const &b, MultiArrayView<2, T, C3> &x, Array const & lambda)
00311 {
00312     typedef T Real;
00313 
00314     const unsigned int rows = rowCount(A);
00315     const unsigned int cols = columnCount(A);
00316     const unsigned int lambdaCount = lambda.size();
00317     vigra_precondition(rows >= cols,
00318        "ridgeRegressionSeries(): Input matrix A must be rectangular with rowCount >= columnCount.");
00319     vigra_precondition(rowCount(b) == rows && columnCount(b) == 1,
00320        "ridgeRegressionSeries(): Shape mismatch between matrices A and b.");
00321     vigra_precondition(rowCount(x) == cols && columnCount(x) == lambdaCount,
00322        "ridgeRegressionSeries(): Result matrix x has wrong shape.");
00323 
00324     unsigned int m = rows;
00325     unsigned int n = cols;
00326 
00327     Matrix<T> u(m, n), s(n, 1), v(n, n);
00328 
00329     unsigned int rank = singularValueDecomposition(A, u, s, v);
00330 
00331     Matrix<T> xl = transpose(u)*b;
00332     Matrix<T> xt(cols,1);
00333     for(unsigned int i=0; i<lambdaCount; ++i)
00334     {
00335         vigra_precondition(lambda[i] >= 0.0,
00336            "ridgeRegressionSeries(): lambda >= 0.0 required.");
00337         if(lambda[i] == 0.0 && rank < rows)
00338             continue;
00339         for(unsigned int k=0; k<cols; ++k)
00340             xt(k,0) = xl(k,0) * s(k,0) / (sq(s(k,0)) + lambda[i]);
00341         columnVector(x, i) = v*xt;
00342     }
00343     return (rank == n);
00344 }
00345 
00346 /** \brief Pass options to leastAngleRegression().
00347 
00348     <b>\#include</b> <<a href="regression_8hxx-source.html">vigra/regression.hxx</a>>
00349         Namespaces: vigra and vigra::linalg
00350 */
00351 class LeastAngleRegressionOptions
00352 {
00353   public:
00354     enum Mode { LARS, LASSO, NNLASSO };
00355 
00356         /** Initialize all options with default values.
00357         */
00358     LeastAngleRegressionOptions()
00359     : max_solution_count(0),
00360       unconstrained_dimension_count(0),
00361       mode(LASSO),
00362       least_squares_solutions(true)
00363     {}
00364 
00365         /** Maximum number of solutions to be computed.
00366 
00367             If \a n is 0 (the default), the number of solutions is determined by the length
00368             of the solution array. Otherwise, the minimum of maxSolutionCount() and that
00369             length is taken.<br>
00370             Default: 0 (use length of solution array)
00371         */
00372     LeastAngleRegressionOptions & maxSolutionCount(unsigned int n)
00373     {
00374         max_solution_count = (int)n;
00375         return *this;
00376     }
00377 
00378         /** Set the mode of the algorithm.
00379 
00380             Mode must be one of "lars", "lasso", "nnlasso". The function just calls
00381             the member function of the corresponding name to set the mode.
00382 
00383             Default: "lasso"
00384         */
00385     LeastAngleRegressionOptions & setMode(std::string mode)
00386     {
00387         for(unsigned int k=0; k<mode.size(); ++k)
00388             mode[k] = (std::string::value_type)tolower(mode[k]);
00389         if(mode == "lars")
00390             this->lars();
00391         else if(mode == "lasso")
00392             this->lasso();
00393         else if(mode == "nnlasso")
00394             this->nnlasso();
00395         else
00396             vigra_fail("LeastAngleRegressionOptions.setMode(): Invalid mode.");
00397         return *this;
00398     }
00399 
00400 
00401         /** Use the plain LARS algorithm.
00402 
00403             Default: inactive
00404         */
00405     LeastAngleRegressionOptions & lars()
00406     {
00407         mode = LARS;
00408         return *this;
00409     }
00410 
00411         /** Use the LASSO modification of the LARS algorithm.
00412 
00413             This allows features to be removed from the active set under certain conditions.<br>
00414             Default: active
00415         */
00416     LeastAngleRegressionOptions & lasso()
00417     {
00418         mode = LASSO;
00419         return *this;
00420     }
00421 
00422         /** Use the non-negative LASSO modification of the LARS algorithm.
00423 
00424             This enforces all non-zero entries in the solution to be positive.<br>
00425             Default: inactive
00426         */
00427     LeastAngleRegressionOptions & nnlasso()
00428     {
00429         mode = NNLASSO;
00430         return *this;
00431     }
00432 
00433         /** Compute least squares solutions.
00434 
00435             Use least angle regression to determine active sets, but
00436             return least squares solutions for the features in each active set,
00437             instead of constrained solutions.<br>
00438             Default: <tt>true</tt>
00439         */
00440     LeastAngleRegressionOptions & leastSquaresSolutions(bool select = true)
00441     {
00442         least_squares_solutions = select;
00443         return *this;
00444     }
00445 
00446     int max_solution_count, unconstrained_dimension_count;
00447     Mode mode;
00448     bool least_squares_solutions;
00449 };
00450 
00451 namespace detail {
00452 
00453 template <class T, class C1, class C2>
00454 struct LarsData
00455 {
00456     typedef typename MultiArrayShape<2>::type Shape;
00457 
00458     int activeSetSize;
00459     MultiArrayView<2, T, C1> A;
00460     MultiArrayView<2, T, C2> b;
00461     Matrix<T> R, qtb, lars_solution, lars_prediction, next_lsq_solution, next_lsq_prediction, searchVector;
00462     ArrayVector<MultiArrayIndex> columnPermutation;
00463 
00464     // init data for a new run
00465     LarsData(MultiArrayView<2, T, C1> const & Ai, MultiArrayView<2, T, C2> const & bi)
00466     : activeSetSize(1),
00467       A(Ai), b(bi), R(A), qtb(b),
00468       lars_solution(A.shape(1), 1), lars_prediction(A.shape(0), 1),
00469       next_lsq_solution(A.shape(1), 1), next_lsq_prediction(A.shape(0), 1), searchVector(A.shape(0), 1),
00470       columnPermutation(A.shape(1))
00471     {
00472         for(unsigned int k=0; k<columnPermutation.size(); ++k)
00473             columnPermutation[k] = k;
00474     }
00475 
00476     // copy data for the recursive call in nnlassolsq
00477     LarsData(LarsData const & d, int asetSize)
00478     : activeSetSize(asetSize),
00479       A(d.R.subarray(Shape(0,0), Shape(d.A.shape(0), activeSetSize))), b(d.qtb), R(A), qtb(b),
00480       lars_solution(d.lars_solution.subarray(Shape(0,0), Shape(activeSetSize, 1))), lars_prediction(d.lars_prediction),
00481       next_lsq_solution(d.next_lsq_solution.subarray(Shape(0,0), Shape(activeSetSize, 1))), 
00482       next_lsq_prediction(d.next_lsq_prediction), searchVector(d.searchVector),
00483       columnPermutation(A.shape(1))
00484     {
00485         for(unsigned int k=0; k<columnPermutation.size(); ++k)
00486             columnPermutation[k] = k;
00487     }
00488 };
00489 
00490 template <class T, class C1, class C2, class Array1, class Array2>
00491 unsigned int leastAngleRegressionMainLoop(LarsData<T, C1, C2> & d,
00492                   Array1 & activeSets, Array2 * lars_solutions, Array2 * lsq_solutions,
00493                   LeastAngleRegressionOptions const & options)
00494 {
00495     using namespace vigra::functor;
00496 
00497     typedef typename MultiArrayShape<2>::type Shape;
00498     typedef typename Matrix<T>::view_type Subarray;
00499     typedef ArrayVector<MultiArrayIndex> Permutation;
00500     typedef typename Permutation::view_type ColumnSet;
00501 
00502     vigra_precondition(d.activeSetSize > 0,
00503        "leastAngleRegressionMainLoop() must not be called with empty active set.");
00504 
00505     bool enforce_positive = (options.mode == LeastAngleRegressionOptions::NNLASSO);
00506     bool lasso_modification = (options.mode != LeastAngleRegressionOptions::LARS);
00507 
00508     const MultiArrayIndex rows = rowCount(d.R);
00509     const MultiArrayIndex cols = columnCount(d.R);
00510     const MultiArrayIndex maxRank = std::min(rows, cols);
00511 
00512     MultiArrayIndex maxSolutionCount = options.max_solution_count;
00513     if(maxSolutionCount == 0)
00514         maxSolutionCount = lasso_modification
00515                                 ? 10*maxRank
00516                                 : maxRank;
00517 
00518     bool needToRemoveColumn = false;
00519     MultiArrayIndex columnToBeAdded, columnToBeRemoved;
00520     MultiArrayIndex currentSolutionCount = 0;
00521     while(currentSolutionCount < maxSolutionCount)
00522     {
00523         ColumnSet activeSet = d.columnPermutation.subarray(0, (unsigned int)d.activeSetSize);
00524         ColumnSet inactiveSet = d.columnPermutation.subarray((unsigned int)d.activeSetSize, (unsigned int)cols);
00525 
00526         // find next dimension to be activated
00527         Matrix<T> cLARS = transpose(d.A) * (d.b - d.lars_prediction),      // correlation with LARS residual
00528                   cLSQ  = transpose(d.A) * (d.b - d.next_lsq_prediction);  // correlation with LSQ residual
00529 
00530         // In theory, all vectors in the active set should have the same correlation C, and
00531         // the correlation of all others should not exceed this. In practice, we may find the
00532         // maximum correlation in any variable due to tiny numerical inaccuracies. Therefore, we
00533         // determine C from the entire set of variables.
00534         MultiArrayIndex cmaxIndex = enforce_positive
00535                                       ? argMax(cLARS)
00536                                       : argMax(abs(cLARS));
00537         T C = abs(cLARS(cmaxIndex, 0));
00538 
00539         Matrix<T> ac(cols - d.activeSetSize, 1);
00540         for(MultiArrayIndex k = 0; k<cols-d.activeSetSize; ++k)
00541         {
00542             T rho = cLSQ(inactiveSet[k], 0), 
00543               cc  = C - sign(rho)*cLARS(inactiveSet[k], 0);
00544 
00545             if(rho == 0.0)  // make sure that 0/0 cannot happen in the other cases
00546                 ac(k,0) = 1.0; // variable k is linearly dependent on the active set
00547             else if(rho > 0.0)
00548                 ac(k,0) = cc / (cc + rho); // variable k would enter the active set with positive sign
00549             else if(enforce_positive)
00550                 ac(k,0) = 1.0; // variable k cannot enter the active set because it would be negative
00551             else
00552                 ac(k,0) = cc / (cc - rho); // variable k would enter the active set with negative sign
00553         }
00554 
00555         // in the non-negative case: make sure that a column just removed cannot re-enter right away
00556         // (in standard LASSO, this is allowed, because the variable may re-enter with opposite sign)
00557         if(enforce_positive && needToRemoveColumn)
00558                 ac(columnToBeRemoved-d.activeSetSize,0) = 1.0;
00559 
00560         // find candidate
00561         // Note: R uses Arg1() > epsilon, but this is only possible because it allows several variables to
00562         //       join the active set simultaneously, so that gamma = 0 cannot occur.
00563         columnToBeAdded = argMin(ac);
00564 
00565         // if no new column can be added, we do a full step gamma = 1.0 and then stop, unless a column is removed below
00566         T gamma = (d.activeSetSize == maxRank)
00567                      ? 1.0
00568                      : ac(columnToBeAdded, 0);
00569 
00570         // adjust columnToBeAdded: we skipped the active set
00571         if(columnToBeAdded >= 0)
00572             columnToBeAdded += d.activeSetSize;
00573 
00574         // check whether we have to remove a column from the active set
00575         needToRemoveColumn = false;
00576         if(lasso_modification)
00577         {
00578             // find dimensions whose weight changes sign below gamma*searchDirection
00579             Matrix<T> s(Shape(d.activeSetSize, 1), NumericTraits<T>::max());
00580             for(MultiArrayIndex k=0; k<d.activeSetSize; ++k)
00581             {
00582                 if(( enforce_positive && d.next_lsq_solution(k,0) < 0.0) ||
00583                    (!enforce_positive && sign(d.lars_solution(k,0))*sign(d.next_lsq_solution(k,0)) == -1.0))
00584                         s(k,0) = d.lars_solution(k,0) / (d.lars_solution(k,0) - d.next_lsq_solution(k,0));
00585             }
00586 
00587             columnToBeRemoved = argMinIf(s, Arg1() <= Param(gamma));
00588             if(columnToBeRemoved >= 0)
00589             {
00590                 needToRemoveColumn = true; // remove takes precedence over add
00591                 gamma = s(columnToBeRemoved, 0);
00592             }
00593         }
00594 
00595         // compute the current solutions
00596         d.lars_prediction  = gamma * d.next_lsq_prediction + (1.0 - gamma) * d.lars_prediction;
00597         d.lars_solution    = gamma * d.next_lsq_solution   + (1.0 - gamma) * d.lars_solution;
00598         if(needToRemoveColumn)
00599             d.lars_solution(columnToBeRemoved, 0) = 0.0;  // turn possible epsilon into an exact zero
00600 
00601         // write the current solution
00602         ++currentSolutionCount;
00603         activeSets.push_back(typename Array1::value_type(d.columnPermutation.begin(), d.columnPermutation.begin()+d.activeSetSize));
00604 
00605         if(lsq_solutions != 0)
00606         {
00607             if(enforce_positive)
00608             {
00609                 ArrayVector<Matrix<T> > nnresults;
00610                 ArrayVector<ArrayVector<MultiArrayIndex> > nnactiveSets;
00611                 LarsData<T, C1, C2> nnd(d, d.activeSetSize);
00612 
00613                 leastAngleRegressionMainLoop(nnd, nnactiveSets, &nnresults, (Array2*)0,
00614                                              LeastAngleRegressionOptions().leastSquaresSolutions(false).nnlasso());
00615                 Matrix<T> nnlsq_solution(d.activeSetSize, 1);
00616                 for(unsigned int k=0; k<nnactiveSets.back().size(); ++k)
00617                 {
00618                     nnlsq_solution(nnactiveSets.back()[k],0) = nnresults.back()[k];
00619                 }
00620                 lsq_solutions->push_back(nnlsq_solution);
00621             }
00622             else
00623             {
00624                 lsq_solutions->push_back(d.next_lsq_solution.subarray(Shape(0,0), Shape(d.activeSetSize, 1)));
00625             }
00626         }
00627         if(lars_solutions != 0)
00628         {
00629             lars_solutions->push_back(d.lars_solution.subarray(Shape(0,0), Shape(d.activeSetSize, 1)));
00630         }
00631 
00632         // no further solutions possible
00633         if(gamma == 1.0)
00634             break;
00635 
00636         if(needToRemoveColumn)
00637         {
00638             --d.activeSetSize;
00639             if(columnToBeRemoved != d.activeSetSize)
00640             {
00641                 // remove column 'columnToBeRemoved' and restore triangular form of R
00642                 // note: columnPermutation is automatically swapped here
00643                 detail::upperTriangularSwapColumns(columnToBeRemoved, d.activeSetSize, d.R, d.qtb, d.columnPermutation);
00644 
00645                 // swap solution entries
00646                 std::swap(d.lars_solution(columnToBeRemoved, 0), d.lars_solution(d.activeSetSize,0));
00647                 std::swap(d.next_lsq_solution(columnToBeRemoved, 0), d.next_lsq_solution(d.activeSetSize,0));
00648                 columnToBeRemoved = d.activeSetSize; // keep track of removed column
00649             }
00650             d.lars_solution(d.activeSetSize,0) = 0.0;
00651             d.next_lsq_solution(d.activeSetSize,0) = 0.0;
00652         }
00653         else
00654         {
00655             vigra_invariant(columnToBeAdded >= 0,
00656                 "leastAngleRegression(): internal error (columnToBeAdded < 0)");
00657             // add column 'columnToBeAdded'
00658             if(d.activeSetSize != columnToBeAdded)
00659             {
00660                 std::swap(d.columnPermutation[d.activeSetSize], d.columnPermutation[columnToBeAdded]);
00661                 columnVector(d.R, d.activeSetSize).swapData(columnVector(d.R, columnToBeAdded));
00662                 columnToBeAdded = d.activeSetSize; // keep track of added column
00663             }
00664 
00665             // zero the corresponding entries of the solutions
00666             d.next_lsq_solution(d.activeSetSize,0) = 0.0;
00667             d.lars_solution(d.activeSetSize,0) = 0.0;
00668 
00669             // reduce R (i.e. its newly added column) to triangular form
00670             detail::qrColumnHouseholderStep(d.activeSetSize, d.R, d.qtb);
00671             ++d.activeSetSize;
00672         }
00673 
00674         // compute the LSQ solution of the new active set
00675         Subarray Ractive = d.R.subarray(Shape(0,0), Shape(d.activeSetSize, d.activeSetSize));
00676         Subarray qtbactive = d.qtb.subarray(Shape(0,0), Shape(d.activeSetSize, 1));
00677         Subarray next_lsq_solution_view = d.next_lsq_solution.subarray(Shape(0,0), Shape(d.activeSetSize, 1));
00678         linearSolveUpperTriangular(Ractive, qtbactive, next_lsq_solution_view);
00679 
00680         // compute the LSQ prediction of the new active set
00681         d.next_lsq_prediction.init(0.0);
00682         for(MultiArrayIndex k=0; k<d.activeSetSize; ++k)
00683             d.next_lsq_prediction += next_lsq_solution_view(k,0)*columnVector(d.A, d.columnPermutation[k]);
00684     }
00685 
00686     return (unsigned int)currentSolutionCount;
00687 }
00688 
00689 #if 0
00690 // old version, keep it until we are sure that the new version works
00691 template <class T, class C1, class C2, class Array1, class Array2>
00692 unsigned int leastAngleRegressionMainLoop(LarsData<T, C1, C2> & d,
00693                   Array1 & activeSets, Array2 * lars_solutions, Array2 * lsq_solutions,
00694                   LeastAngleRegressionOptions const & options)
00695 {
00696     using namespace vigra::functor;
00697 
00698     typedef typename MultiArrayShape<2>::type Shape;
00699     typedef typename Matrix<T>::view_type Subarray;
00700     typedef ArrayVector<MultiArrayIndex> Permutation;
00701     typedef typename Permutation::view_type ColumnSet;
00702 
00703     vigra_precondition(d.activeSetSize > 0,
00704        "leastAngleRegressionMainLoop() must not be called with empty active set.");
00705 
00706     bool enforce_positive = (options.mode == LeastAngleRegressionOptions::NNLASSO);
00707     bool lasso_modification = (options.mode != LeastAngleRegressionOptions::LARS);
00708 
00709     const MultiArrayIndex rows = rowCount(d.R);
00710     const MultiArrayIndex cols = columnCount(d.R);
00711     const MultiArrayIndex maxRank = std::min(rows, cols);
00712 
00713     MultiArrayIndex maxSolutionCount = options.max_solution_count;
00714     if(maxSolutionCount == 0)
00715         maxSolutionCount = lasso_modification
00716                                 ? 10*maxRank
00717                                 : maxRank;
00718 
00719     bool needToRemoveColumn = false;
00720     MultiArrayIndex columnToBeAdded, columnToBeRemoved;
00721     MultiArrayIndex currentSolutionCount = 0;
00722     while(currentSolutionCount < maxSolutionCount)
00723     {
00724         ColumnSet activeSet = d.columnPermutation.subarray(0, (unsigned int)d.activeSetSize);
00725         ColumnSet inactiveSet = d.columnPermutation.subarray((unsigned int)d.activeSetSize, (unsigned int)cols);
00726 
00727         // find next dimension to be activated
00728         Matrix<T> lars_residual = d.b - d.lars_prediction;
00729         Matrix<T> lsq_residual  = lars_residual - d.searchVector;
00730         Matrix<T> c = transpose(d.A) * lars_residual;
00731 
00732         // In theory, all vecors in the active set should have the same correlation C, and
00733         // the correlation of all others should not exceed this. In practice, we may find the
00734         // maximum correlation in any variable due to tiny numerical inaccuracies. Therefore, we
00735         // determine C from the entire set of variables.
00736         MultiArrayIndex cmaxIndex = enforce_positive
00737                                       ? argMax(c)
00738                                       : argMax(abs(c));
00739         T C = abs(c(cmaxIndex, 0));
00740 
00741         Matrix<T> ac(cols - d.activeSetSize, 1);
00742         for(MultiArrayIndex k = 0; k<cols-d.activeSetSize; ++k)
00743         {
00744             T a   = dot(columnVector(d.A, inactiveSet[k]), d.searchVector),
00745               a1  = (C - c(inactiveSet[k], 0)) / (C - a),
00746               a2  = (C + c(inactiveSet[k], 0)) / (C + a);
00747 
00748             if(enforce_positive)
00749                 ac(k, 0) = a1;
00750             else if(std::min(a1, a2) < 0.0)
00751                 ac(k, 0) = std::max(a1, a2);
00752             else
00753                 ac(k, 0) = std::min(a1, a2);
00754         }
00755 
00756         // in the non-negative case: make sure that a column just removed cannot re-enter right away
00757         // (in standard LASSO, this is allowed, because the variable may re-enter with opposite sign)
00758         if(enforce_positive && needToRemoveColumn)
00759                 ac(columnToBeRemoved-d.activeSetSize,0) = -1.0;
00760 
00761         // find candidate
00762         // Note: R uses Arg1() > epsilon, but this is only possible because it allows several variables to
00763         //       join the active set simultaneously, so that gamma = 0 cannot occur.
00764         columnToBeAdded = argMinIf(ac, Arg1() < Param(1.0) && Arg1() >= Param(NumericTraits<T>::zero()));
00765 
00766         // if no new column can be added, we do a full step gamma = 1.0 and then stop, unless a column is removed below
00767         T gamma = (columnToBeAdded == -1 || d.activeSetSize == maxRank)
00768                      ? 1.0
00769                      : ac(columnToBeAdded, 0);
00770 
00771         // adjust columnToBeAdded: we skipped the active set
00772         if(columnToBeAdded >= 0)
00773             columnToBeAdded += d.activeSetSize;
00774 
00775         // check whether we have to remove a column from the active set
00776         needToRemoveColumn = false;
00777         if(lasso_modification)
00778         {
00779             // find dimensions whose weight changes sign below gamma*searchDirection
00780             Matrix<T> s(Shape(d.activeSetSize, 1), NumericTraits<T>::max());
00781             for(MultiArrayIndex k=0; k<d.activeSetSize; ++k)
00782             {
00783                 if(( enforce_positive && d.next_lsq_solution(k,0) < 0.0) ||
00784                    (!enforce_positive && sign(d.lars_solution(k,0))*sign(d.next_lsq_solution(k,0)) == -1.0))
00785                         s(k,0) = d.lars_solution(k,0) / (d.lars_solution(k,0) - d.next_lsq_solution(k,0));
00786             }
00787 
00788             columnToBeRemoved = argMinIf(s, Arg1() < Param(gamma) && Arg1() >= Param(NumericTraits<T>::zero()));
00789             if(columnToBeRemoved >= 0)
00790             {
00791                 needToRemoveColumn = true; // remove takes precedence over add
00792                 gamma = s(columnToBeRemoved, 0);
00793             }
00794         }
00795 
00796         // compute the current solutions
00797         d.lars_prediction += gamma * d.searchVector;
00798         d.lars_solution    = gamma * d.next_lsq_solution + (1.0 - gamma) * d.lars_solution;
00799         if(needToRemoveColumn)
00800             d.lars_solution(columnToBeRemoved, 0) = 0.0;  // turn possible epsilon into an exact zero
00801 
00802         // write the current solution
00803         ++currentSolutionCount;
00804         activeSets.push_back(typename Array1::value_type(d.columnPermutation.begin(), d.columnPermutation.begin()+d.activeSetSize));
00805 
00806         if(lsq_solutions != 0)
00807         {
00808             if(enforce_positive)
00809             {
00810                 ArrayVector<Matrix<T> > nnresults;
00811                 ArrayVector<ArrayVector<MultiArrayIndex> > nnactiveSets;
00812                 LarsData<T, C1, C2> nnd(d, d.activeSetSize);
00813 
00814                 leastAngleRegressionMainLoop(nnd, nnactiveSets, &nnresults, (Array2*)0,
00815                                              LeastAngleRegressionOptions().leastSquaresSolutions(false).nnlasso());
00816                 Matrix<T> nnlsq_solution(d.activeSetSize, 1);
00817                 for(unsigned int k=0; k<nnactiveSets.back().size(); ++k)
00818                 {
00819                     nnlsq_solution(nnactiveSets.back()[k],0) = nnresults.back()[k];
00820                 }
00821                 lsq_solutions->push_back(nnlsq_solution);
00822             }
00823             else
00824             {
00825                 lsq_solutions->push_back(d.next_lsq_solution.subarray(Shape(0,0), Shape(d.activeSetSize, 1)));
00826             }
00827         }
00828         if(lars_solutions != 0)
00829         {
00830             lars_solutions->push_back(d.lars_solution.subarray(Shape(0,0), Shape(d.activeSetSize, 1)));
00831         }
00832 
00833         // no further solutions possible
00834         if(gamma == 1.0)
00835             break;
00836 
00837         if(needToRemoveColumn)
00838         {
00839             --d.activeSetSize;
00840             if(columnToBeRemoved != d.activeSetSize)
00841             {
00842                 // remove column 'columnToBeRemoved' and restore triangular form of R
00843                 // note: columnPermutation is automatically swapped here
00844                 detail::upperTriangularSwapColumns(columnToBeRemoved, d.activeSetSize, d.R, d.qtb, d.columnPermutation);
00845 
00846                 // swap solution entries
00847                 std::swap(d.lars_solution(columnToBeRemoved, 0), d.lars_solution(d.activeSetSize,0));
00848                 std::swap(d.next_lsq_solution(columnToBeRemoved, 0), d.next_lsq_solution(d.activeSetSize,0));
00849                 columnToBeRemoved = d.activeSetSize; // keep track of removed column
00850             }
00851             d.lars_solution(d.activeSetSize,0) = 0.0;
00852             d.next_lsq_solution(d.activeSetSize,0) = 0.0;
00853         }
00854         else
00855         {
00856             vigra_invariant(columnToBeAdded >= 0,
00857                 "leastAngleRegression(): internal error (columnToBeAdded < 0)");
00858             // add column 'columnToBeAdded'
00859             if(d.activeSetSize != columnToBeAdded)
00860             {
00861                 std::swap(d.columnPermutation[d.activeSetSize], d.columnPermutation[columnToBeAdded]);
00862                 columnVector(d.R, d.activeSetSize).swapData(columnVector(d.R, columnToBeAdded));
00863                 columnToBeAdded = d.activeSetSize; // keep track of added column
00864             }
00865 
00866             // zero the corresponding entries of the solutions
00867             d.next_lsq_solution(d.activeSetSize,0) = 0.0;
00868             d.lars_solution(d.activeSetSize,0) = 0.0;
00869 
00870             // reduce R (i.e. its newly added column) to triangular form
00871             detail::qrColumnHouseholderStep(d.activeSetSize, d.R, d.qtb);
00872             ++d.activeSetSize;
00873         }
00874 
00875         // compute the LSQ solution of the new active set
00876         Subarray Ractive = d.R.subarray(Shape(0,0), Shape(d.activeSetSize, d.activeSetSize));
00877         Subarray qtbactive = d.qtb.subarray(Shape(0,0), Shape(d.activeSetSize, 1));
00878         Subarray next_lsq_solution_view = d.next_lsq_solution.subarray(Shape(0,0), Shape(d.activeSetSize, 1));
00879         linearSolveUpperTriangular(Ractive, qtbactive, next_lsq_solution_view);
00880 
00881         // compute the LSQ prediction of the new active set
00882         Matrix<T> lsq_prediction(rows, 1);
00883         for(MultiArrayIndex k=0; k<d.activeSetSize; ++k)
00884             lsq_prediction += next_lsq_solution_view(k,0)*columnVector(d.A, d.columnPermutation[k]);
00885 
00886         // finally, the new search direction
00887         d.searchVector = lsq_prediction - d.lars_prediction;
00888     }
00889 
00890     return (unsigned int)currentSolutionCount;
00891 }
00892 #endif
00893 
00894 template <class T, class C1, class C2, class Array1, class Array2>
00895 unsigned int
00896 leastAngleRegressionImpl(MultiArrayView<2, T, C1> const & A, MultiArrayView<2, T, C2> const &b,
00897                          Array1 & activeSets, Array2 * lasso_solutions, Array2 * lsq_solutions,
00898                          LeastAngleRegressionOptions const & options)
00899 {
00900     using namespace vigra::functor;
00901 
00902     const MultiArrayIndex rows = rowCount(A);
00903 
00904     vigra_precondition(rowCount(b) == rows && columnCount(b) == 1,
00905        "leastAngleRegression(): Shape mismatch between matrices A and b.");
00906 
00907     bool enforce_positive = (options.mode == LeastAngleRegressionOptions::NNLASSO);
00908 
00909     detail::LarsData<T, C1, C2> d(A, b);
00910 
00911     // find dimension with largest correlation
00912     Matrix<T> c = transpose(A)*b;
00913     MultiArrayIndex initialColumn = enforce_positive
00914                                        ? argMaxIf(c, Arg1() > Param(0.0))
00915                                        : argMax(abs(c));
00916     if(initialColumn == -1)
00917         return 0; // no solution found
00918 
00919     // prepare initial active set and search direction etc.
00920     std::swap(d.columnPermutation[0], d.columnPermutation[initialColumn]);
00921     columnVector(d.R, 0).swapData(columnVector(d.R, initialColumn));
00922     detail::qrColumnHouseholderStep(0, d.R, d.qtb);
00923     d.next_lsq_solution(0,0) = d.qtb(0,0) / d.R(0,0);
00924     d.next_lsq_prediction = d.next_lsq_solution(0,0) * columnVector(A, d.columnPermutation[0]);
00925     d.searchVector = d.next_lsq_solution(0,0) * columnVector(A, d.columnPermutation[0]);
00926 
00927     return leastAngleRegressionMainLoop(d, activeSets, lasso_solutions, lsq_solutions, options);
00928 }
00929 
00930 } // namespace detail
00931 
00932    /** Least Angle Regression.
00933 
00934     <b>\#include</b> <<a href="regression_8hxx-source.html">vigra/regression.hxx</a>>
00935         Namespaces: vigra and vigra::linalg
00936 
00937    <b> Declarations:</b>
00938 
00939     \code
00940     namespace vigra {
00941       namespace linalg {
00942         // compute either LASSO or least sqaures solutions
00943         template <class T, class C1, class C2, class Array1, class Array2>
00944         unsigned int
00945         leastAngleRegression(MultiArrayView<2, T, C1> const & A, MultiArrayView<2, T, C2> const &b,
00946                              Array1 & activeSets, Array2 & solutions,
00947                              LeastAngleRegressionOptions const & options = LeastAngleRegressionOptions());
00948 
00949         // compute LASSO and least sqaures solutions
00950         template <class T, class C1, class C2, class Array1, class Array2>
00951         unsigned int
00952         leastAngleRegression(MultiArrayView<2, T, C1> const & A, MultiArrayView<2, T, C2> const &b,
00953                              Array1 & activeSets, Array2 & lasso_solutions, Array2 & lsq_solutions,
00954                              LeastAngleRegressionOptions const & options = LeastAngleRegressionOptions());
00955       }
00956       using linalg::leastAngleRegression;
00957     }
00958     \endcode
00959 
00960        This function implements Least Angle Regression (LARS) as described in
00961 
00962        &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;
00963        B.Efron, T.Hastie, I.Johnstone, and R.Tibshirani: <em>"Least Angle Regression"</em>,
00964        Annals of Statistics 32(2):407-499, 2004.
00965 
00966        It is an efficient algorithm to solve the L1-regularized least squares (LASSO) problem
00967 
00968         \f[ \tilde \textrm{\bf x} = \textrm{argmin}
00969              \left|\left|\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right|\right|_2^2
00970            \textrm{ subject to } \left|\left|\textrm{\bf x}\right|\right|_1\le s
00971         \f]
00972 
00973        and the L1-regularized non-negative least squares (NN-LASSO) problem
00974 
00975         \f[ \tilde \textrm{\bf x} = \textrm{argmin} \left|\left|\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right|\right|_2^2
00976            \textrm{ subject to } \left|\left|\textrm{\bf x}\right|\right|_1\le s \textrm{ and } \textrm{\bf x}\ge \textrm{\bf 0}
00977         \f]
00978 
00979        where \a A is a matrix with <tt>m</tt> rows and <tt>n</tt> columns (often with <tt>m < n</tt>),
00980        \a b a vector of length <tt>m</tt>, and a regularization parameter s >= 0.0.
00981        L1-regularization has the desirable effect that it causes the solution \a x to be sparse, i.e. only
00982        the most important variables (called the <em>active set</em>) have non-zero values. The
00983        key insight of the LARS algorithm is the following: When the solution vector is considered
00984        as a function of the regularization parameter s, then <b>x</b>(s) is a piecewise
00985        linear function, i.e. a polyline in n-dimensional space. The knots of the polyline
00986        occur precisely at those values of s where one variable enters or leaves the active set,
00987        and can be efficiently computed.
00988 
00989        Therefore, leastAngleRegression() returns the entire solution path as a sequence of knot points, starting
00990        at \f$\textrm{\bf x}(s=0)\f$ (where the only feasible solution is obviously <b>x</b> = 0) and ending at
00991        \f$\textrm{\bf x}(s=\infty)\f$ (where the solution becomes the ordinary least squares solution). Actually,
00992        the initial null solution is not explicitly returned, i.e. the sequence starts at the first non-zero
00993        solution with one variable in the active set. The function leastAngleRegression() returns the number
00994        of solutions( i.e. knot points) computed.
00995 
00996        The sequences of active sets and corresponding variable weights are returned in \a activeSets and
00997        \a solutions respectively. That is, <tt>activeSets[i]</tt> is an \ref vigra::ArrayVector "ArrayVector<int>"
00998        containing the indices of the variables that are active at the i-th knot, and <tt>solutions</tt> is a
00999        \ref vigra::linalg::Matrix "Matrix<T>" containing the weights of those variables, in the same order (see
01000        example below). Variables not contained in <tt>activeSets[i]</tt> are zero at this solution.
01001 
01002        The behavior of the algorithm can be adapted by \ref vigra::linalg::LeastAngleRegressionOptions
01003        "LeastAngleRegressionOptions":
01004         <DL>
01005         <DT><b>options.lasso()</b> (active by default)
01006                           <DD> Compute the LASSO solution as described above.
01007         <DT><b>options.nnlasso()</b> (inactive by default)
01008                           <DD> Compute non-negative LASSO solutions, i.e. use the additional constraint that
01009                                <b>x</b> >= 0 in all solutions.
01010         <DT><b>options.lars()</b> (inactive by default)
01011                           <DD> Compute a solution path according to the plain LARS rule, i.e. never remove
01012                                a variable from the active set once it entered.
01013         <DT><b>options.leastSquaresSolutions(bool)</b> (default: true)
01014                           <DD> Use the algorithm mode selected above
01015                                to determine the sequence of active sets, but then compute and return an
01016                                ordinary (unconstrained) least squares solution for every active set.<br>
01017                                <b>Note:</b> The second form of leastAngleRegression() ignores this option and
01018                                does always compute both constrained and unconstrained solutions (returned in
01019                                \a lasso_solutions and \a lsq_solutions respectively).
01020         <DT><b>maxSolutionCount(unsigned int n)</b> (default: n = 0, i.e. compute all solutions)
01021                           <DD> Compute at most <tt>n</tt> solutions.
01022         </DL>
01023 
01024         <b>Usage:</b>
01025 
01026         \code
01027         int m = ..., n = ...;
01028         Matrix<double> A(m, n), b(m, 1);
01029         ... // fill A and b
01030 
01031         // normalize the input
01032         Matrix<double> offset(1,n), scaling(1,n);
01033         prepareColumns(A, A, offset, scaling, DataPreparationGoals(ZeroMean|UnitVariance));
01034         prepareColumns(b, b, DataPreparationGoals(ZeroMean));
01035 
01036         // arrays to hold the output
01037         ArrayVector<ArrayVector<int> > activeSets;
01038         ArrayVector<Matrix<double> > solutions;
01039 
01040         // run leastAngleRegression() in non-negative LASSO mode
01041         int numSolutions = leastAngleRegression(A, b, activeSets, solutions,
01042                                     LeastAngleRegressionOptions().nnlasso());
01043 
01044         // print results
01045         Matrix<double> denseSolution(1, n);
01046         for (MultiArrayIndex k = 0; k < numSolutions; ++k)
01047         {
01048             // transform the sparse solution into a dense vector
01049             denseSolution.init(0.0); // ensure that inactive variables are zero
01050             for (unsigned int i = 0; i < activeSets[k].size(); ++i)
01051             {
01052                 // set the values of the active variables;
01053                 // activeSets[k][i] is the true index of the i-th variable in the active set
01054                 denseSolution(0, activeSets[k][i]) = solutions[k](i,0);
01055             }
01056 
01057             // invert the input normalization
01058             denseSolution = denseSolution * pointWise(scaling);
01059 
01060             // output the solution
01061             std::cout << "solution " << k << ":\n" << denseSolution << std::endl;
01062         }
01063         \endcode
01064 
01065         <b>Required Interface:</b>
01066 
01067         <ul>
01068         <li> <tt>T</tt> must be numeric type (compatible to double)
01069         <li> <tt>Array1 a1;</tt><br>
01070              <tt>a1.push_back(ArrayVector<int>());</tt>
01071         <li> <tt>Array2 a2;</tt><br>
01072              <tt>a2.push_back(Matrix<T>());</tt>
01073         </ul>
01074    */
01075 doxygen_overloaded_function(template <...> unsigned int leastAngleRegression)
01076 
01077 template <class T, class C1, class C2, class Array1, class Array2>
01078 inline unsigned int
01079 leastAngleRegression(MultiArrayView<2, T, C1> const & A, MultiArrayView<2, T, C2> const &b,
01080                      Array1 & activeSets, Array2 & solutions,
01081                      LeastAngleRegressionOptions const & options = LeastAngleRegressionOptions())
01082 {
01083     if(options.least_squares_solutions)
01084         return detail::leastAngleRegressionImpl(A, b, activeSets, (Array2*)0, &solutions, options);
01085     else
01086         return detail::leastAngleRegressionImpl(A, b, activeSets, &solutions, (Array2*)0, options);
01087 }
01088 
01089 template <class T, class C1, class C2, class Array1, class Array2>
01090 inline unsigned int
01091 leastAngleRegression(MultiArrayView<2, T, C1> const & A, MultiArrayView<2, T, C2> const &b,
01092                      Array1 & activeSets, Array2 & lasso_solutions, Array2 & lsq_solutions,
01093                      LeastAngleRegressionOptions const & options = LeastAngleRegressionOptions())
01094 {
01095     return detail::leastAngleRegressionImpl(A, b, activeSets, &lasso_solutions, &lsq_solutions, options);
01096 }
01097 
01098    /** Non-negative Least Squares Regression.
01099 
01100        Given a matrix \a A with <tt>m</tt> rows and <tt>n</tt> columns (with <tt>m >= n</tt>),
01101        and a column vector \a b of length <tt>m</tt> rows, this function computes
01102        a column vector \a x of length <tt>n</tt> with <b>non-negative entries</b> that solves the optimization problem
01103 
01104         \f[ \tilde \textrm{\bf x} = \textrm{argmin}
01105             \left|\left|\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right|\right|_2^2
01106             \textrm{ subject to } \textrm{\bf x} \ge \textrm{\bf 0}
01107         \f]
01108 
01109        Both \a b and \a x must be column vectors (i.e. matrices with <tt>1</tt> column).
01110        Note that all matrices must already have the correct shape. The solution is computed by means
01111        of \ref leastAngleRegression() with non-negativity constraint.
01112 
01113     <b>\#include</b> <<a href="regression_8hxx-source.html">vigra/regression.hxx</a>>
01114         Namespaces: vigra and vigra::linalg
01115    */
01116 template <class T, class C1, class C2, class C3>
01117 inline void
01118 nonnegativeLeastSquares(MultiArrayView<2, T, C1> const & A,
01119                         MultiArrayView<2, T, C2> const &b, MultiArrayView<2, T, C3> &x)
01120 {
01121     vigra_precondition(columnCount(A) == rowCount(x) && rowCount(A) == rowCount(b),
01122         "nonnegativeLeastSquares(): Matrix shape mismatch.");
01123     vigra_precondition(columnCount(b) == 1 && columnCount(x) == 1,
01124         "nonnegativeLeastSquares(): RHS and solution must be vectors (i.e. columnCount == 1).");
01125 
01126     ArrayVector<ArrayVector<MultiArrayIndex> > activeSets;
01127     ArrayVector<Matrix<T> > results;
01128 
01129     leastAngleRegression(A, b, activeSets, results,
01130                          LeastAngleRegressionOptions().leastSquaresSolutions(false).nnlasso());
01131     x.init(NumericTraits<T>::zero());
01132     if(activeSets.size() > 0)
01133         for(unsigned int k=0; k<activeSets.back().size(); ++k)
01134             x(activeSets.back()[k],0) = results.back()[k];
01135 }
01136 
01137 
01138 //@}
01139 
01140 } // namespace linalg
01141 
01142 using linalg::leastSquares;
01143 using linalg::weightedLeastSquares;
01144 using linalg::ridgeRegression;
01145 using linalg::weightedRidgeRegression;
01146 using linalg::ridgeRegressionSeries;
01147 using linalg::nonnegativeLeastSquares;
01148 using linalg::leastAngleRegression;
01149 using linalg::LeastAngleRegressionOptions;
01150 
01151 } // namespace vigra
01152 
01153 #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.7.0 (Thu Aug 25 2011)