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

Classes | Functions
Optimization and Regression

Classes

class  LeastAngleRegressionOptions
 Pass options to leastAngleRegression(). More...

Functions

template<... >
unsigned int leastAngleRegression (...)
template<class T , class C1 , class C2 , class C3 >
bool leastSquares (MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > const &b, MultiArrayView< 2, T, C3 > &x, std::string method="QR")
template<class T , class C1 , class C2 , class C3 >
void nonnegativeLeastSquares (MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > const &b, MultiArrayView< 2, T, C3 > &x)
template<class T , class C1 , class C2 , class C3 , class C4 , class C5 , class C6 , class C7 >
quadraticProgramming (MultiArrayView< 2, T, C1 > const &G, MultiArrayView< 2, T, C2 > const &g, MultiArrayView< 2, T, C3 > const &CE, MultiArrayView< 2, T, C4 > const &ce, MultiArrayView< 2, T, C5 > const &CI, MultiArrayView< 2, T, C6 > const &ci, MultiArrayView< 2, T, C7 > &x)
template<class T , class C1 , class C2 , class C3 >
bool ridgeRegression (MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > const &b, MultiArrayView< 2, T, C3 > &x, double lambda)
template<class T , class C1 , class C2 , class C3 , class Array >
bool ridgeRegressionSeries (MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > const &b, MultiArrayView< 2, T, C3 > &x, Array const &lambda)
template<class T , class C1 , class C2 , class C3 , class C4 >
bool weightedLeastSquares (MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > const &b, MultiArrayView< 2, T, C3 > const &weights, MultiArrayView< 2, T, C4 > &x, std::string method="QR")
template<class T , class C1 , class C2 , class C3 , class C4 >
bool weightedRidgeRegression (MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > const &b, MultiArrayView< 2, T, C3 > const &weights, MultiArrayView< 2, T, C4 > &x, double lambda)

Function Documentation

T vigra::quadraticProgramming ( MultiArrayView< 2, T, C1 > const &  G,
MultiArrayView< 2, T, C2 > const &  g,
MultiArrayView< 2, T, C3 > const &  CE,
MultiArrayView< 2, T, C4 > const &  ce,
MultiArrayView< 2, T, C5 > const &  CI,
MultiArrayView< 2, T, C6 > const &  ci,
MultiArrayView< 2, T, C7 > &  x 
)

Solve Quadratic Programming Problem.

The quadraticProgramming() function implements the algorithm described in

D. Goldfarb, A. Idnani: "A numerically stable dual method for solving strictly convex quadratic programs", Mathematical Programming 27:1-33, 1983.

for the solution of (convex) quadratic programming problems by means of a primal-dual method.

#include <vigra/quadprog.hxx> Namespaces: vigra

Declaration:

     namespace vigra { 
         template <class T, class C1, class C2, class C3, class C4, class C5, class C6, class C7>
         T 
         quadraticProgramming(MultiArrayView<2, T, C1> const & GG, MultiArrayView<2, T, C2> const & g,  
                              MultiArrayView<2, T, C3> const & CE, MultiArrayView<2, T, C4> const & ce,  
                              MultiArrayView<2, T, C5> const & CI, MultiArrayView<2, T, C6> const & ci, 
                              MultiArrayView<2, T, C7> & x);
     }

The problem must be specified in the form:

\begin{eqnarray*} \mbox{minimize } &\,& \frac{1}{2} \mbox{\bf x}'\,\mbox{\bf G}\, \mbox{\bf x} + \mbox{\bf g}'\,\mbox{\bf x} \\ \mbox{subject to} &\,& \mbox{\bf C}_E\, \mbox{\bf x} = \mbox{\bf c}_e \\ &\,& \mbox{\bf C}_I\,\mbox{\bf x} \ge \mbox{\bf c}_i \end{eqnarray*}

Matrix G G must be symmetric positive definite, and matrix CE must have full row rank. Matrix and vector dimensions must be as follows:

  • G: [n * n], g: [n * 1]
  • CE: [me * n], ce: [me * 1]
  • CI: [mi * n], ci: [mi * 1]
  • x: [n * 1]

The function writes the optimal solution into the vector x and returns the cost of this solution. If the problem is infeasible, std::numeric_limits::infinity() is returned. In this case the value of vector x is undefined.

Usage:

Minimize f = 0.5 * x'*G*x + g'*x subject to -1 <= x <= 1. The solution is x' = [1.0, 0.5, -1.0] with f = -22.625.

      double Gdata[] = {13.0, 12.0, -2.0,
                        12.0, 17.0,  6.0,
                        -2.0,  6.0, 12.0};

      double gdata[] = {-22.0, -14.5, 13.0};

      double CIdata[] = { 1.0,  0.0,  0.0,
                          0.0,  1.0,  0.0,
                          0.0,  0.0,  1.0,
                         -1.0,  0.0,  0.0,
                          0.0, -1.0,  0.0,
                          0.0,  0.0, -1.0};
                        
      double cidata[] = {-1.0, -1.0, -1.0, -1.0, -1.0, -1.0};

      Matrix<double> G(3,3, Gdata), 
                     g(3,1, gdata), 
                     CE,             // empty since there are no equality constraints
                     ce,             // likewise
                     CI(7,3, CIdata), 
                     ci(7,1, cidata), 
                     x(3,1);
                   
      double f = quadraticProgramming(G, g, CE, ce, CI, ci, x);
bool vigra::linalg::leastSquares ( MultiArrayView< 2, T, C1 > const &  A,
MultiArrayView< 2, T, C2 > const &  b,
MultiArrayView< 2, T, C3 > &  x,
std::string  method = "QR" 
)

Ordinary Least Squares Regression.

Given a matrix A with m rows and n columns (with m >= n), and a column vector b of length m rows, this function computes the column vector x of length n rows that solves the optimization problem

\[ \tilde \textrm{\bf x} = \textrm{argmin} \left|\left|\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right|\right|_2^2 \]

When b is a matrix with k columns, x must also have k columns, which will contain the solutions for the corresponding columns of b. Note that all matrices must already have the correct shape.

This function is just another name for linearSolve(), perhaps leading to more readable code when A is a rectangular matrix. It returns false when the rank of A is less than n. See linearSolve() for more documentation.

#include <vigra/regression.hxx> Namespaces: vigra and vigra::linalg

bool vigra::linalg::weightedLeastSquares ( MultiArrayView< 2, T, C1 > const &  A,
MultiArrayView< 2, T, C2 > const &  b,
MultiArrayView< 2, T, C3 > const &  weights,
MultiArrayView< 2, T, C4 > &  x,
std::string  method = "QR" 
)

Weighted Least Squares Regression.

Given a matrix A with m rows and n columns (with m >= n), a vector b of length m, and a weight vector weights of length m with non-negative entries, this function computes the vector x of length n that solves the optimization problem

\[ \tilde \textrm{\bf x} = \textrm{argmin} \left(\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right)^T \textrm{diag}(\textrm{\bf weights}) \left(\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right) \]

where diag(weights) creates a diagonal matrix from weights. The algorithm calls leastSquares() on the equivalent problem

\[ \tilde \textrm{\bf x} = \textrm{argmin} \left|\left|\textrm{diag}(\textrm{\bf weights})^{1/2}\textrm{\bf A} \textrm{\bf x} - \textrm{diag}(\textrm{\bf weights})^{1/2} \textrm{\bf b}\right|\right|_2^2 \]

where the square root of weights is just taken element-wise.

When b is a matrix with k columns, x must also have k columns, which will contain the solutions for the corresponding columns of b. Note that all matrices must already have the correct shape.

The function returns false when the rank of the weighted matrix A is less than n.

#include <vigra/regression.hxx> Namespaces: vigra and vigra::linalg

bool vigra::linalg::ridgeRegression ( MultiArrayView< 2, T, C1 > const &  A,
MultiArrayView< 2, T, C2 > const &  b,
MultiArrayView< 2, T, C3 > &  x,
double  lambda 
)

Ridge Regression.

Given a matrix A with m rows and n columns (with m >= n), a vector b of length m, and a regularization parameter lambda >= 0.0, this function computes the vector x of length n that solves the optimization problem

\[ \tilde \textrm{\bf x} = \textrm{argmin} \left|\left|\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right|\right|_2^2 + \lambda \textrm{\bf x}^T\textrm{\bf x} \]

This is implemented by means of singularValueDecomposition().

When b is a matrix with k columns, x must also have k columns, which will contain the solutions for the corresponding columns of b. Note that all matrices must already have the correct shape.

The function returns false if the rank of A is less than n and lambda == 0.0.

#include <vigra/regression.hxx> Namespaces: vigra and vigra::linalg

bool vigra::linalg::weightedRidgeRegression ( MultiArrayView< 2, T, C1 > const &  A,
MultiArrayView< 2, T, C2 > const &  b,
MultiArrayView< 2, T, C3 > const &  weights,
MultiArrayView< 2, T, C4 > &  x,
double  lambda 
)

Weighted ridge Regression.

Given a matrix A with m rows and n columns (with m >= n), a vector b of length m, a weight vector weights of length m with non-negative entries, and a regularization parameter lambda >= 0.0 this function computes the vector x of length n that solves the optimization problem

\[ \tilde \textrm{\bf x} = \textrm{argmin} \left(\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right)^T \textrm{diag}(\textrm{\bf weights}) \left(\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right) + \lambda \textrm{\bf x}^T\textrm{\bf x} \]

where diag(weights) creates a diagonal matrix from weights. The algorithm calls ridgeRegression() on the equivalent problem

\[ \tilde \textrm{\bf x} = \textrm{argmin} \left|\left|\textrm{diag}(\textrm{\bf weights})^{1/2}\textrm{\bf A} \textrm{\bf x} - \textrm{diag}(\textrm{\bf weights})^{1/2} \textrm{\bf b}\right|\right|_2^2 + \lambda \textrm{\bf x}^T\textrm{\bf x} \]

where the square root of weights is just taken element-wise. This solution is computed by means of singularValueDecomposition().

When b is a matrix with k columns, x must also have k columns, which will contain the solutions for the corresponding columns of b. Note that all matrices must already have the correct shape.

The function returns false if the rank of A is less than n and lambda == 0.0.

#include <vigra/regression.hxx> Namespaces: vigra and vigra::linalg

bool vigra::linalg::ridgeRegressionSeries ( MultiArrayView< 2, T, C1 > const &  A,
MultiArrayView< 2, T, C2 > const &  b,
MultiArrayView< 2, T, C3 > &  x,
Array const &  lambda 
)

Ridge Regression with many lambdas.

This executes ridgeRegression() for a sequence of regularization parameters. This is implemented so that the singularValueDecomposition() has to be executed only once. lambda must be an array conforming to the std::vector interface, i.e. must support lambda.size() and lambda[k]. The columns of the matrix x will contain the solutions for the corresponding lambda, so the number of columns of the matrix x must be equal to lambda.size(), and b must be a columns vector, i.e. cannot contain several right hand sides at once.

The function returns false when the matrix A is rank deficient. If this happens, and one of the lambdas is zero, the corresponding column of x will be skipped.

#include <vigra/regression.hxx> Namespaces: vigra and vigra::linalg

unsigned int vigra::linalg::leastAngleRegression (   ...)

Least Angle Regression.

#include <vigra/regression.hxx> Namespaces: vigra and vigra::linalg

Declarations:

    namespace vigra {
      namespace linalg {
        // compute either LASSO or least sqaures solutions
        template <class T, class C1, class C2, class Array1, class Array2>
        unsigned int
        leastAngleRegression(MultiArrayView<2, T, C1> const & A, MultiArrayView<2, T, C2> const &b,
                             Array1 & activeSets, Array2 & solutions,
                             LeastAngleRegressionOptions const & options = LeastAngleRegressionOptions());

        // compute LASSO and least sqaures solutions
        template <class T, class C1, class C2, class Array1, class Array2>
        unsigned int
        leastAngleRegression(MultiArrayView<2, T, C1> const & A, MultiArrayView<2, T, C2> const &b,
                             Array1 & activeSets, Array2 & lasso_solutions, Array2 & lsq_solutions,
                             LeastAngleRegressionOptions const & options = LeastAngleRegressionOptions());
      }
      using linalg::leastAngleRegression;
    }

This function implements Least Angle Regression (LARS) as described in

        B.Efron, T.Hastie, I.Johnstone, and R.Tibshirani: "Least Angle Regression", Annals of Statistics 32(2):407-499, 2004.

It is an efficient algorithm to solve the L1-regularized least squares (LASSO) problem

\[ \tilde \textrm{\bf x} = \textrm{argmin} \left|\left|\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right|\right|_2^2 \textrm{ subject to } \left|\left|\textrm{\bf x}\right|\right|_1\le s \]

and the L1-regularized non-negative least squares (NN-LASSO) problem

\[ \tilde \textrm{\bf x} = \textrm{argmin} \left|\left|\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right|\right|_2^2 \textrm{ subject to } \left|\left|\textrm{\bf x}\right|\right|_1\le s \textrm{ and } \textrm{\bf x}\ge \textrm{\bf 0} \]

where A is a matrix with m rows and n columns (often with m < n), b a vector of length m, and a regularization parameter s >= 0.0. L1-regularization has the desirable effect that it causes the solution x to be sparse, i.e. only the most important variables (called the active set) have non-zero values. The key insight of the LARS algorithm is the following: When the solution vector is considered as a function of the regularization parameter s, then x(s) is a piecewise linear function, i.e. a polyline in n-dimensional space. The knots of the polyline occur precisely at those values of s where one variable enters or leaves the active set, and can be efficiently computed.

Therefore, leastAngleRegression() returns the entire solution path as a sequence of knot points, starting at $\textrm{\bf x}(s=0)$ (where the only feasible solution is obviously x = 0) and ending at $\textrm{\bf x}(s=\infty)$ (where the solution becomes the ordinary least squares solution). Actually, the initial null solution is not explicitly returned, i.e. the sequence starts at the first non-zero solution with one variable in the active set. The function leastAngleRegression() returns the number of solutions( i.e. knot points) computed.

The sequences of active sets and corresponding variable weights are returned in activeSets and solutions respectively. That is, activeSets[i] is an ArrayVector<int> containing the indices of the variables that are active at the i-th knot, and solutions is a Matrix<T> containing the weights of those variables, in the same order (see example below). Variables not contained in activeSets[i] are zero at this solution.

The behavior of the algorithm can be adapted by LeastAngleRegressionOptions:

options.lasso() (active by default)
Compute the LASSO solution as described above.
options.nnlasso() (inactive by default)
Compute non-negative LASSO solutions, i.e. use the additional constraint that x >= 0 in all solutions.
options.lars() (inactive by default)
Compute a solution path according to the plain LARS rule, i.e. never remove a variable from the active set once it entered.
options.leastSquaresSolutions(bool) (default: true)
Use the algorithm mode selected above to determine the sequence of active sets, but then compute and return an ordinary (unconstrained) least squares solution for every active set.
Note: The second form of leastAngleRegression() ignores this option and does always compute both constrained and unconstrained solutions (returned in lasso_solutions and lsq_solutions respectively).
maxSolutionCount(unsigned int n) (default: n = 0, i.e. compute all solutions)
Compute at most n solutions.

Usage:

        int m = ..., n = ...;
        Matrix<double> A(m, n), b(m, 1);
        ... // fill A and b

        // normalize the input
        Matrix<double> offset(1,n), scaling(1,n);
        prepareColumns(A, A, offset, scaling, DataPreparationGoals(ZeroMean|UnitVariance));
        prepareColumns(b, b, DataPreparationGoals(ZeroMean));

        // arrays to hold the output
        ArrayVector<ArrayVector<int> > activeSets;
        ArrayVector<Matrix<double> > solutions;

        // run leastAngleRegression() in non-negative LASSO mode
        int numSolutions = leastAngleRegression(A, b, activeSets, solutions,
                                    LeastAngleRegressionOptions().nnlasso());

        // print results
        Matrix<double> denseSolution(1, n);
        for (MultiArrayIndex k = 0; k < numSolutions; ++k)
        {
            // transform the sparse solution into a dense vector
            denseSolution.init(0.0); // ensure that inactive variables are zero
            for (unsigned int i = 0; i < activeSets[k].size(); ++i)
            {
                // set the values of the active variables;
                // activeSets[k][i] is the true index of the i-th variable in the active set
                denseSolution(0, activeSets[k][i]) = solutions[k](i,0);
            }

            // invert the input normalization
            denseSolution = denseSolution * pointWise(scaling);

            // output the solution
            std::cout << "solution " << k << ":\n" << denseSolution << std::endl;
        }

Required Interface:

  • T must be numeric type (compatible to double)
  • Array1 a1;
    a1.push_back(ArrayVector<int>());
  • Array2 a2;
    a2.push_back(Matrix<T>());
void vigra::linalg::nonnegativeLeastSquares ( MultiArrayView< 2, T, C1 > const &  A,
MultiArrayView< 2, T, C2 > const &  b,
MultiArrayView< 2, T, C3 > &  x 
)

Non-negative Least Squares Regression.

Given a matrix A with m rows and n columns (with m >= n), and a column vector b of length m rows, this function computes a column vector x of length n with non-negative entries that solves the optimization problem

\[ \tilde \textrm{\bf x} = \textrm{argmin} \left|\left|\textrm{\bf A} \textrm{\bf x} - \textrm{\bf b}\right|\right|_2^2 \textrm{ subject to } \textrm{\bf x} \ge \textrm{\bf 0} \]

Both b and x must be column vectors (i.e. matrices with 1 column). Note that all matrices must already have the correct shape. The solution is computed by means of leastAngleRegression() with non-negativity constraint.

#include <vigra/regression.hxx> Namespaces: vigra and vigra::linalg

© 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)