[ 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.1 (3 Dec 2010)