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"
79 template <
class T,
class C1,
class C2,
class C3>
83 std::string method =
"QR")
121 template <
class T,
class C1,
class C2,
class C3,
class C4>
127 const unsigned int rows =
rowCount(A);
130 vigra_precondition(rows >= cols,
131 "weightedLeastSquares(): Input matrix A must be rectangular with rowCount >= columnCount.");
132 vigra_precondition(
rowCount(b) == rows,
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.");
146 for(
unsigned int l=0; l<cols; ++l)
147 wa(k,l) = w * A(k,l);
148 for(
unsigned int l=0; l<rhsCount; ++l)
149 wb(k,l) = w * b(k,l);
179 template <
class T,
class C1,
class C2,
class C3>
184 const unsigned int rows =
rowCount(A);
187 vigra_precondition(rows >= cols,
188 "ridgeRegression(): Input matrix A must be rectangular with rowCount >= columnCount.");
189 vigra_precondition(
rowCount(b) == rows,
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)
207 for(
unsigned int l=0; l<rhsCount; ++l)
208 t(k,l) *= s(k,0) / (
sq(s(k,0)) + lambda);
250 template <
class T,
class C1,
class C2,
class C3,
class C4>
256 const unsigned int rows =
rowCount(A);
259 vigra_precondition(rows >= cols,
260 "weightedRidgeRegression(): Input matrix A must be rectangular with rowCount >= columnCount.");
261 vigra_precondition(
rowCount(b) == rows,
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.");
277 for(
unsigned int l=0; l<cols; ++l)
278 wa(k,l) = w * A(k,l);
279 for(
unsigned int l=0; l<rhsCount; ++l)
280 wb(k,l) = w * b(k,l);
302 template <
class T,
class C1,
class C2,
class C3,
class Array>
307 const unsigned int rows =
rowCount(A);
309 const unsigned int lambdaCount = lambda.size();
310 vigra_precondition(rows >= cols,
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;
326 for(
unsigned int i=0; i<lambdaCount; ++i)
328 vigra_precondition(lambda[i] >= 0.0,
329 "ridgeRegressionSeries(): lambda >= 0.0 required.");
330 if(lambda[i] == 0.0 && rank < rows)
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;
445 template <
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;
482 template <
class T,
class C1,
class C2,
class Array1,
class Array2,
class Array3>
484 leastAngleRegressionMainLoop(LarsData<T, C1, C2> & d,
486 Array2 * lars_solutions, Array3 * lsq_solutions,
487 LeastAngleRegressionOptions
const & options)
489 using namespace vigra::functor;
492 typedef typename Matrix<T>::view_type Subarray;
493 typedef ArrayVector<MultiArrayIndex> Permutation;
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);
500 bool lasso_modification = (options.mode != LeastAngleRegressionOptions::LARS);
507 if(maxSolutionCount == 0)
508 maxSolutionCount = lasso_modification
512 bool needToRemoveColumn =
false;
515 while(currentSolutionCount < maxSolutionCount)
518 ColumnSet inactiveSet = d.columnPermutation.subarray(static_cast<unsigned int>(d.activeSetSize), static_cast<unsigned int>(cols));
521 Matrix<T> cLARS =
transpose(d.A) * (d.b - d.lars_prediction),
522 cLSQ =
transpose(d.A) * (d.b - d.next_lsq_prediction);
530 : argMax(
abs(cLARS));
531 T C =
abs(cLARS(cmaxIndex, 0));
533 Matrix<T> ac(cols - d.activeSetSize, 1);
536 T rho = cLSQ(inactiveSet[k], 0),
537 cc = C -
sign(rho)*cLARS(inactiveSet[k], 0);
542 ac(k,0) = cc / (cc + rho);
543 else if(enforce_positive)
546 ac(k,0) = cc / (cc - rho);
551 if(enforce_positive && needToRemoveColumn)
552 ac(columnToBeRemoved-d.activeSetSize,0) = 1.0;
557 columnToBeAdded =
argMin(ac);
560 T
gamma = (d.activeSetSize == maxRank)
562 : ac(columnToBeAdded, 0);
565 if(columnToBeAdded >= 0)
566 columnToBeAdded += d.activeSetSize;
569 needToRemoveColumn =
false;
570 if(lasso_modification)
573 Matrix<T> s(Shape(d.activeSetSize, 1), NumericTraits<T>::max());
576 if(( enforce_positive && d.next_lsq_solution(k,0) < 0.0) ||
577 (!enforce_positive &&
sign(d.lars_solution(k,0))*
sign(d.next_lsq_solution(k,0)) == -1.0))
578 s(k,0) = d.lars_solution(k,0) / (d.lars_solution(k,0) - d.next_lsq_solution(k,0));
581 columnToBeRemoved =
argMinIf(s, Arg1() <= Param(gamma));
582 if(columnToBeRemoved >= 0)
584 needToRemoveColumn =
true;
585 gamma = s(columnToBeRemoved, 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;
592 if(needToRemoveColumn)
593 d.lars_solution(columnToBeRemoved, 0) = 0.0;
596 ++currentSolutionCount;
597 activeSets.push_back(
typename Array1::value_type(d.columnPermutation.begin(), d.columnPermutation.begin()+d.activeSetSize));
599 if(lsq_solutions != 0)
603 ArrayVector<Matrix<T> > nnresults;
604 ArrayVector<ArrayVector<MultiArrayIndex> > nnactiveSets;
605 LarsData<T, C1, C2> nnd(d, d.activeSetSize);
607 leastAngleRegressionMainLoop(nnd, nnactiveSets, &nnresults, static_cast<Array3*>(0),
608 LeastAngleRegressionOptions().leastSquaresSolutions(
false).nnlasso());
610 typename Array2::value_type nnlsq_solution(Shape(d.activeSetSize, 1));
611 for(
unsigned int k=0; k<nnactiveSets.back().size(); ++k)
613 nnlsq_solution(nnactiveSets.back()[k],0) = nnresults.back()[k];
616 lsq_solutions->push_back(
typename Array3::value_type());
617 lsq_solutions->back() = nnlsq_solution;
622 lsq_solutions->push_back(
typename Array3::value_type());
623 lsq_solutions->back() = d.next_lsq_solution.subarray(Shape(0,0), Shape(d.activeSetSize, 1));
626 if(lars_solutions != 0)
629 lars_solutions->push_back(
typename Array2::value_type());
630 lars_solutions->back() = d.lars_solution.subarray(Shape(0,0), Shape(d.activeSetSize, 1));
637 if(needToRemoveColumn)
640 if(columnToBeRemoved != d.activeSetSize)
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));
649 columnToBeRemoved = d.activeSetSize;
651 d.lars_solution(d.activeSetSize,0) = 0.0;
652 d.next_lsq_solution(d.activeSetSize,0) = 0.0;
656 vigra_invariant(columnToBeAdded >= 0,
657 "leastAngleRegression(): internal error (columnToBeAdded < 0)");
659 if(d.activeSetSize != columnToBeAdded)
661 std::swap(d.columnPermutation[d.activeSetSize], d.columnPermutation[columnToBeAdded]);
663 columnToBeAdded = d.activeSetSize;
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));
677 Subarray qtbactive = d.qtb.subarray(Shape(0,0), Shape(d.activeSetSize, 1));
678 Subarray next_lsq_solution_view = d.next_lsq_solution.subarray(Shape(0,0), Shape(d.activeSetSize, 1));
682 d.next_lsq_prediction.init(0.0);
684 d.next_lsq_prediction += next_lsq_solution_view(k,0)*
columnVector(d.A, d.columnPermutation[k]);
687 return static_cast<unsigned int>(currentSolutionCount);
690 template <
class T,
class C1,
class C2,
class Array1,
class Array2>
692 leastAngleRegressionImpl(MultiArrayView<2, T, C1>
const & A, MultiArrayView<2, T, C2>
const &b,
693 Array1 & activeSets, Array2 * lasso_solutions, Array2 * lsq_solutions,
694 LeastAngleRegressionOptions
const & options)
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);
712 if(initialColumn == -1)
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]);
723 return leastAngleRegressionMainLoop(d, activeSets, lasso_solutions, lsq_solutions, options);
873 template <
class T,
class C1,
class C2,
class Array1,
class Array2>
876 Array1 & activeSets, Array2 & solutions,
877 LeastAngleRegressionOptions
const & options = LeastAngleRegressionOptions())
879 if(options.least_squares_solutions)
880 return detail::leastAngleRegressionImpl(A, b, activeSets, static_cast<Array2*>(0), &solutions, options);
882 return detail::leastAngleRegressionImpl(A, b, activeSets, &solutions, static_cast<Array2*>(0), options);
885 template <
class T,
class C1,
class C2,
class Array1,
class Array2>
888 Array1 & activeSets, Array2 & lasso_solutions, Array2 & lsq_solutions,
889 LeastAngleRegressionOptions
const & options = LeastAngleRegressionOptions())
891 return detail::leastAngleRegressionImpl(A, b, activeSets, &lasso_solutions, &lsq_solutions, options);
929 template <
class T,
class C1,
class C2,
class C3>
932 MultiArrayView<2, T, C2>
const &b, MultiArrayView<2, T, C3> &x)
935 "nonnegativeLeastSquares(): Matrix shape mismatch.");
937 "nonnegativeLeastSquares(): RHS and solution must be vectors (i.e. columnCount == 1).");
939 ArrayVector<ArrayVector<MultiArrayIndex> > activeSets;
940 ArrayVector<Matrix<T> > results;
943 LeastAngleRegressionOptions().leastSquaresSolutions(
false).nnlasso());
944 x.init(NumericTraits<T>::zero());
945 if(activeSets.size() > 0)
946 for(
unsigned int k=0; k<activeSets.back().size(); ++k)
947 x(activeSets.back()[k],0) = results.back()[k];
962 using linalg::LeastAngleRegressionOptions;
966 template <
class T,
class S>
973 template <
class T,
class S>
974 inline MultiArrayView<1, T>
977 return a.bindInner(i);
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;
1050 template <
unsigned int D,
class T,
class S1,
class S2,
1054 nonlinearLeastSquaresImpl(MultiArrayView<D, T, S1>
const & features,
1055 MultiArrayView<1, T, S2>
const & response,
1056 TinyVector<U, N> & p,
1058 NonlinearLSQOptions
const & options)
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
1068 ? std::max(epsilonT, epsilonU)
1071 linalg::Matrix<T> jj(N,N);
1072 TinyVector<U, N> jr, dp;
1075 bool didStep =
true;
1077 for(
int iter=0; iter<options.max_iter; ++iter)
1086 for(
int i=0; i<features.shape(0); ++i)
1088 autodiff::DualVector<U, N> res = model(detail::getRow(features, i), autodiff::dualMatrix(p));
1090 T r = response(i) - res.v;
1098 linalg::Matrix<T> djj(jj);
1099 djj.diagonal() *= 1.0 + l;
1102 TinyVector<U, N> p_new = p + dp;
1105 T residual_new = 0.0;
1106 for(
int i=0; i<features.shape(0); ++i)
1108 residual_new +=
sq(response(i) - model(detail::getRow(features, i), p_new));
1111 if(residual_new < residual)
1115 if(
std::abs((residual - residual_new) / residual) < epsilon)
1116 return residual_new;
1263 template <
class T,
class S1,
class S2,
1268 MultiArrayView<1, T, S2>
const & response,
1269 TinyVector<U, N> & p,
1271 NonlinearLSQOptions
const & options = NonlinearLSQOptions())
1273 return nonlinearLeastSquaresImpl(features, response, p, model, options);
1276 template <
class T,
class S1,
class S2,
1281 MultiArrayView<1, T, S2>
const & response,
1282 TinyVector<U, N> & p,
1284 NonlinearLSQOptions
const & options = NonlinearLSQOptions())
1286 return nonlinearLeastSquaresImpl(features, response, p, model, options);
1293 #endif // VIGRA_REGRESSION_HXX
LeastAngleRegressionOptions & maxSolutionCount(unsigned int n)
Definition: regression.hxx:365
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
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
MultiArrayView< 2, T, C > columnVector(MultiArrayView< 2, T, C > const &m, MultiArrayIndex d)
Definition: matrix.hxx:727
LeastAngleRegressionOptions & nnlasso()
Definition: regression.hxx:419
std::string tolower(std::string s)
Definition: utilities.hxx:96
Pass options to leastAngleRegression().
Definition: regression.hxx:344
NonlinearLSQOptions & maxIterations(int iter)
Set maximum number of iterations.
Definition: regression.hxx:1024
MultiArrayIndex rowCount(const MultiArrayView< 2, T, C > &x)
Definition: matrix.hxx:671
void transpose(const MultiArrayView< 2, T, C1 > &v, MultiArrayView< 2, T, C2 > &r)
Definition: matrix.hxx:965
int argMin(MultiArrayView< 2, T, C > const &a)
Find the index of the minimum element in a matrix.
Definition: matrix.hxx:1953
const difference_type & shape() const
Definition: multi_array.hxx:1648
NonlinearLSQOptions & dampingParamters(double lambda, double v)
Set damping parameters for Levenberg-Marquardt algorithm.
Definition: regression.hxx:1040
linalg::TemporaryMatrix< T > sign(MultiArrayView< 2, T, C > const &v)
LeastAngleRegressionOptions & leastSquaresSolutions(bool select=true)
Definition: regression.hxx:432
NonlinearLSQOptions()
Initialize options with default values.
Definition: regression.hxx:1000
int argMinIf(MultiArrayView< 2, T, C > const &a, UnaryFunctor condition)
Find the index of the minimum element in a matrix subject to a condition.
Definition: matrix.hxx:2021
LeastAngleRegressionOptions & lars()
Definition: regression.hxx:397
Pass options to nonlinearLeastSquares().
Definition: regression.hxx:991
int argMaxIf(MultiArrayView< 2, T, C > const &a, UnaryFunctor condition)
Find the index of the maximum element in a matrix subject to a condition.
Definition: matrix.hxx:2056
std::ptrdiff_t MultiArrayIndex
Definition: multi_fwd.hxx:60
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 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
unsigned int nonnegativeLeastSquares(...)
Definition: multi_fwd.hxx:63
LeastAngleRegressionOptions()
Definition: regression.hxx:351
linalg::TemporaryMatrix< T > sq(MultiArrayView< 2, T, C > const &v)
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
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
TinyVector< MultiArrayIndex, N > type
Definition: multi_shape.hxx:272
LeastAngleRegressionOptions & setMode(std::string mode)
Definition: regression.hxx:378
int argMax(MultiArrayView< 2, T, C > const &a)
Find the index of the maximum element in a matrix.
Definition: matrix.hxx:1986
double gamma(double x)
The gamma function.
Definition: mathutil.hxx:1587
void outer(const MultiArrayView< 2, T, C1 > &x, const MultiArrayView< 2, T, C2 > &y, MultiArrayView< 2, T, C3 > &r)
Definition: matrix.hxx:1459
MultiArrayIndex columnCount(const MultiArrayView< 2, T, C > &x)
Definition: matrix.hxx:684
LeastAngleRegressionOptions & lasso()
Definition: regression.hxx:408
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
bool linearSolveUpperTriangular(const MultiArrayView< 2, T, C1 > &r, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > x)
Definition: linear_solve.hxx:1068
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
void nonlinearLeastSquares(...)
Fit a non-linear model to given data by minimizing least squares loss.
linalg::TemporaryMatrix< T > abs(MultiArrayView< 2, T, C > const &v)
unsigned int leastAngleRegression(...)
NonlinearLSQOptions & tolerance(double eps)
Set minimum relative improvement in residual.
Definition: regression.hxx:1014
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616