37#ifndef VIGRA_REGRESSION_HXX
38#define VIGRA_REGRESSION_HXX
41#include "linear_solve.hxx"
42#include "singular_value_decomposition.hxx"
43#include "numerictraits.hxx"
44#include "functorexpression.hxx"
45#include "autodiff.hxx"
79template <
class T,
class C1,
class C2,
class C3>
83 std::string method =
"QR")
121template <
class T,
class C1,
class C2,
class C3,
class C4>
131 "weightedLeastSquares(): Input matrix A must be rectangular with rowCount >= columnCount.");
133 "weightedLeastSquares(): Shape mismatch between matrices A and b.");
135 "weightedLeastSquares(): Weight matrix has wrong shape.");
137 "weightedLeastSquares(): Result matrix x has wrong shape.");
141 for(
unsigned int k=0;
k<
rows; ++
k)
143 vigra_precondition(weights(
k,0) >= 0,
144 "weightedLeastSquares(): Weights must be positive.");
145 T w = std::sqrt(weights(
k,0));
146 for(
unsigned int l=0;
l<
cols; ++
l)
179template <
class T,
class C1,
class C2,
class C3>
188 "ridgeRegression(): Input matrix A must be rectangular with rowCount >= columnCount.");
190 "ridgeRegression(): Shape mismatch between matrices A and b.");
192 "ridgeRegression(): Result matrix x has wrong shape.");
193 vigra_precondition(lambda >= 0.0,
194 "ridgeRegression(): lambda >= 0.0 required.");
196 unsigned int m =
rows;
197 unsigned int n =
cols;
202 if(
rank < n && lambda == 0.0)
206 for(
unsigned int k=0;
k<
cols; ++
k)
208 t(
k,
l) *= s(
k,0) / (sq(s(
k,0)) + lambda);
250template <
class T,
class C1,
class C2,
class C3,
class C4>
260 "weightedRidgeRegression(): Input matrix A must be rectangular with rowCount >= columnCount.");
262 "weightedRidgeRegression(): Shape mismatch between matrices A and b.");
264 "weightedRidgeRegression(): Weight matrix has wrong shape.");
266 "weightedRidgeRegression(): Result matrix x has wrong shape.");
267 vigra_precondition(lambda >= 0.0,
268 "weightedRidgeRegression(): lambda >= 0.0 required.");
272 for(
unsigned int k=0;
k<
rows; ++
k)
274 vigra_precondition(weights(
k,0) >= 0,
275 "weightedRidgeRegression(): Weights must be positive.");
276 T w = std::sqrt(weights(
k,0));
277 for(
unsigned int l=0;
l<
cols; ++
l)
302template <
class T,
class C1,
class C2,
class C3,
class Array>
311 "ridgeRegressionSeries(): Input matrix A must be rectangular with rowCount >= columnCount.");
313 "ridgeRegressionSeries(): Shape mismatch between matrices A and b.");
315 "ridgeRegressionSeries(): Result matrix x has wrong shape.");
317 unsigned int m =
rows;
318 unsigned int n =
cols;
328 vigra_precondition(lambda[
i] >= 0.0,
329 "ridgeRegressionSeries(): lambda >= 0.0 required.");
332 for(
unsigned int k=0;
k<
cols; ++
k)
333 xt(
k,0) =
xl(
k,0) * s(
k,0) / (sq(s(
k,0)) + lambda[
i]);
347 enum Mode { LARS, LASSO, NNLASSO };
352 : max_solution_count(0),
353 unconstrained_dimension_count(0),
355 least_squares_solutions(
true)
367 max_solution_count =
static_cast<int>(n);
383 else if(mode ==
"lasso")
385 else if(mode ==
"nnlasso")
388 vigra_fail(
"LeastAngleRegressionOptions.setMode(): Invalid mode.");
434 least_squares_solutions =
select;
438 int max_solution_count, unconstrained_dimension_count;
440 bool least_squares_solutions;
445template <
class T,
class C1,
class C2>
453 Matrix<T> R, qtb, lars_solution, lars_prediction, next_lsq_solution, next_lsq_prediction, searchVector;
459 A(
Ai), b(
bi), R(A), qtb(b),
460 lars_solution(A.shape(1), 1), lars_prediction(A.shape(0), 1),
461 next_lsq_solution(A.shape(1), 1), next_lsq_prediction(A.shape(0), 1), searchVector(A.shape(0), 1),
462 columnPermutation(A.shape(1))
464 for(
unsigned int k=0;
k<columnPermutation.
size(); ++
k)
465 columnPermutation[
k] =
k;
469 LarsData(LarsData
const & d,
int asetSize)
470 : activeSetSize(asetSize),
471 A(d.R.subarray(Shape(0,0), Shape(d.A.shape(0), activeSetSize))), b(d.qtb), R(A), qtb(b),
472 lars_solution(d.lars_solution.subarray(Shape(0,0), Shape(activeSetSize, 1))), lars_prediction(d.lars_prediction),
473 next_lsq_solution(d.next_lsq_solution.subarray(Shape(0,0), Shape(activeSetSize, 1))),
474 next_lsq_prediction(d.next_lsq_prediction), searchVector(d.searchVector),
475 columnPermutation(A.shape(1))
477 for(
unsigned int k=0; k<columnPermutation.
size(); ++k)
478 columnPermutation[k] = k;
482template <
class T,
class C1,
class C2,
class Array1,
class Array2,
class Array3>
484leastAngleRegressionMainLoop(LarsData<T, C1, C2> & d,
486 Array2 * lars_solutions, Array3 * lsq_solutions,
487 LeastAngleRegressionOptions
const & options)
489 using namespace vigra::functor;
494 typedef typename Permutation::view_type
ColumnSet;
496 vigra_precondition(d.activeSetSize > 0,
497 "leastAngleRegressionMainLoop() must not be called with empty active set.");
499 bool enforce_positive = (options.mode == LeastAngleRegressionOptions::NNLASSO);
507 if(maxSolutionCount == 0)
518 ColumnSet inactiveSet = d.columnPermutation.subarray(
static_cast<unsigned int>(d.activeSetSize),
static_cast<unsigned int>(
cols));
573 Matrix<T> s(Shape(d.activeSetSize, 1), NumericTraits<T>::max());
578 s(
k,0) = d.lars_solution(
k,0) / (d.lars_solution(
k,0) - d.next_lsq_solution(
k,0));
590 d.lars_prediction =
gamma * d.next_lsq_prediction + (1.0 -
gamma) * d.lars_prediction;
591 d.lars_solution =
gamma * d.next_lsq_solution + (1.0 -
gamma) * d.lars_solution;
597 activeSets.push_back(
typename Array1::value_type(d.columnPermutation.begin(), d.columnPermutation.begin()+d.activeSetSize));
610 typename Array2::value_type
nnlsq_solution(Shape(d.activeSetSize, 1));
623 lsq_solutions->back() = d.next_lsq_solution.subarray(Shape(0,0), Shape(d.activeSetSize, 1));
630 lars_solutions->back() = d.lars_solution.subarray(Shape(0,0), Shape(d.activeSetSize, 1));
644 detail::upperTriangularSwapColumns(
columnToBeRemoved, d.activeSetSize, d.R, d.qtb, d.columnPermutation);
647 std::swap(d.lars_solution(
columnToBeRemoved, 0), d.lars_solution(d.activeSetSize,0));
648 std::swap(d.next_lsq_solution(
columnToBeRemoved, 0), d.next_lsq_solution(d.activeSetSize,0));
651 d.lars_solution(d.activeSetSize,0) = 0.0;
652 d.next_lsq_solution(d.activeSetSize,0) = 0.0;
657 "leastAngleRegression(): internal error (columnToBeAdded < 0)");
661 std::swap(d.columnPermutation[d.activeSetSize], d.columnPermutation[
columnToBeAdded]);
667 d.next_lsq_solution(d.activeSetSize,0) = 0.0;
668 d.lars_solution(d.activeSetSize,0) = 0.0;
671 detail::qrColumnHouseholderStep(d.activeSetSize, d.R, d.qtb);
676 Subarray Ractive = d.R.subarray(Shape(0,0), Shape(d.activeSetSize, d.activeSetSize));
682 d.next_lsq_prediction.init(0.0);
690template <
class T,
class C1,
class C2,
class Array1,
class Array2>
696 using namespace vigra::functor;
701 "leastAngleRegression(): Shape mismatch between matrices A and b.");
703 bool enforce_positive = (options.mode == LeastAngleRegressionOptions::NNLASSO);
705 detail::LarsData<T, C1, C2> d(A, b);
716 std::swap(d.columnPermutation[0], d.columnPermutation[
initialColumn]);
718 detail::qrColumnHouseholderStep(0, d.R, d.qtb);
719 d.next_lsq_solution(0,0) = d.qtb(0,0) / d.R(0,0);
720 d.next_lsq_prediction = d.next_lsq_solution(0,0) *
columnVector(A, d.columnPermutation[0]);
721 d.searchVector = d.next_lsq_solution(0,0) *
columnVector(A, d.columnPermutation[0]);
873template <
class T,
class C1,
class C2,
class Array1,
class Array2>
879 if(options.least_squares_solutions)
885template <
class T,
class C1,
class C2,
class Array1,
class Array2>
929template <
class T,
class C1,
class C2,
class C3>
935 "nonnegativeLeastSquares(): Matrix shape mismatch.");
937 "nonnegativeLeastSquares(): RHS and solution must be vectors (i.e. columnCount == 1).");
944 x.
init(NumericTraits<T>::zero());
966template <
class T,
class S>
973template <
class T,
class S>
995 double epsilon, lambda, tau;
1042 vigra_precondition(lambda > 0.0 && v > 0.0,
1043 "NonlinearLSQOptions::dampingParamters(): parameters must be positive.");
1044 this->lambda = lambda;
1050template <
unsigned int D,
class T,
class S1,
class S2,
1060 vigra_precondition(features.
shape(0) == response.
shape(0),
1061 "nonlinearLeastSquares(): shape mismatch between features and response.");
1063 double t = options.tau,
l = options.lambda;
1065 double epsilonT = NumericTraits<T>::epsilon()*10.0,
1066 epsilonU = NumericTraits<U>::epsilon()*10.0,
1067 epsilon = options.epsilon <= 0.0
1077 for(
int iter=0; iter<options.max_iter; ++iter)
1086 for(
int i=0;
i<features.
shape(0); ++
i)
1090 T r = response(
i) - res.v;
1099 djj.diagonal() *= 1.0 +
l;
1106 for(
int i=0;
i<features.
shape(0); ++
i)
1263template <
class T,
class S1,
class S2,
1273 return nonlinearLeastSquaresImpl(features, response, p, model, options);
1276template <
class T,
class S1,
class S2,
1286 return nonlinearLeastSquaresImpl(features, response, p, model, options);
Base class for, and view to, vigra::MultiArray.
Definition multi_array.hxx:705
const difference_type & shape() const
Definition multi_array.hxx:1648
MultiArrayView< N-M, T, StridedArrayTag > bindInner(const TinyVector< Index, M > &d) const
Pass options to nonlinearLeastSquares().
Definition regression.hxx:992
NonlinearLSQOptions & dampingParamters(double lambda, double v)
Set damping parameters for Levenberg-Marquardt algorithm.
Definition regression.hxx:1040
NonlinearLSQOptions & maxIterations(int iter)
Set maximum number of iterations.
Definition regression.hxx:1024
NonlinearLSQOptions()
Initialize options with default values.
Definition regression.hxx:1000
NonlinearLSQOptions & tolerance(double eps)
Set minimum relative improvement in residual.
Definition regression.hxx:1014
Class for a single RGB value.
Definition rgbvalue.hxx:128
void init(Iterator i, Iterator end)
Definition tinyvector.hxx:708
size_type size() const
Definition tinyvector.hxx:913
Class for fixed size vectors.
Definition tinyvector.hxx:1008
Definition autodiff.hxx:88
Pass options to leastAngleRegression().
Definition regression.hxx:345
LeastAngleRegressionOptions & nnlasso()
Definition regression.hxx:419
LeastAngleRegressionOptions & lars()
Definition regression.hxx:397
LeastAngleRegressionOptions & lasso()
Definition regression.hxx:408
LeastAngleRegressionOptions()
Definition regression.hxx:351
LeastAngleRegressionOptions & maxSolutionCount(unsigned int n)
Definition regression.hxx:365
LeastAngleRegressionOptions & leastSquaresSolutions(bool select=true)
Definition regression.hxx:432
LeastAngleRegressionOptions & setMode(std::string mode)
Definition regression.hxx:378
Definition matrix.hxx:125
unsigned int singularValueDecomposition(MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > &U, MultiArrayView< 2, T, C3 > &S, MultiArrayView< 2, T, C4 > &V)
Definition singular_value_decomposition.hxx:75
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")
Definition regression.hxx:123
void transpose(const MultiArrayView< 2, T, C1 > &v, MultiArrayView< 2, T, C2 > &r)
Definition matrix.hxx:965
bool leastSquares(MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > const &b, MultiArrayView< 2, T, C3 > &x, std::string method="QR")
Definition regression.hxx:81
bool linearSolveUpperTriangular(const MultiArrayView< 2, T, C1 > &r, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > x)
Definition linear_solve.hxx:1068
MultiArrayIndex columnCount(const MultiArrayView< 2, T, C > &x)
Definition matrix.hxx:684
bool ridgeRegressionSeries(MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > const &b, MultiArrayView< 2, T, C3 > &x, Array const &lambda)
Definition regression.hxx:304
bool ridgeRegression(MultiArrayView< 2, T, C1 > const &A, MultiArrayView< 2, T, C2 > const &b, MultiArrayView< 2, T, C3 > &x, double lambda)
Definition regression.hxx:181
MultiArrayIndex rowCount(const MultiArrayView< 2, T, C > &x)
Definition matrix.hxx:671
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)
Definition regression.hxx:252
MultiArrayView< 2, T, C > columnVector(MultiArrayView< 2, T, C > const &m, MultiArrayIndex d)
Definition matrix.hxx:727
unsigned int leastAngleRegression(...)
unsigned int nonnegativeLeastSquares(...)
linalg::TemporaryMatrix< T > sign(MultiArrayView< 2, T, C > const &v)
Iterator argMaxIf(Iterator first, Iterator last, UnaryFunctor condition)
Find the maximum element in a sequence conforming to a condition.
Definition algorithm.hxx:165
std::string tolower(std::string s)
Definition utilities.hxx:96
Iterator argMin(Iterator first, Iterator last)
Find the minimum element in a sequence.
Definition algorithm.hxx:68
double gamma(double x)
The gamma function.
Definition mathutil.hxx:1587
Iterator argMax(Iterator first, Iterator last)
Find the maximum element in a sequence.
Definition algorithm.hxx:96
void nonlinearLeastSquares(...)
Fit a non-linear model to given data by minimizing least squares loss.
std::ptrdiff_t MultiArrayIndex
Definition multi_fwd.hxx:60
Iterator argMinIf(Iterator first, Iterator last, UnaryFunctor condition)
Find the minimum element in a sequence conforming to a condition.
Definition algorithm.hxx:129