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"
78 template <
class T,
class C1,
class C2,
class C3>
82 std::string method =
"QR")
120 template <
class T,
class C1,
class C2,
class C3,
class C4>
128 const unsigned int rows =
rowCount(A);
131 vigra_precondition(rows >= cols,
132 "weightedLeastSquares(): Input matrix A must be rectangular with rowCount >= columnCount.");
133 vigra_precondition(
rowCount(b) == rows,
134 "weightedLeastSquares(): Shape mismatch between matrices A and b.");
136 "weightedLeastSquares(): Weight matrix has wrong shape.");
138 "weightedLeastSquares(): Result matrix x has wrong shape.");
142 for(
unsigned int k=0; k<rows; ++k)
144 vigra_precondition(weights(k,0) >= 0,
145 "weightedLeastSquares(): Weights must be positive.");
147 for(
unsigned int l=0; l<cols; ++l)
148 wa(k,l) = w * A(k,l);
149 for(
unsigned int l=0; l<rhsCount; ++l)
150 wb(k,l) = w * b(k,l);
180 template <
class T,
class C1,
class C2,
class C3>
187 const unsigned int rows =
rowCount(A);
190 vigra_precondition(rows >= cols,
191 "ridgeRegression(): Input matrix A must be rectangular with rowCount >= columnCount.");
192 vigra_precondition(
rowCount(b) == rows,
193 "ridgeRegression(): Shape mismatch between matrices A and b.");
195 "ridgeRegression(): Result matrix x has wrong shape.");
196 vigra_precondition(lambda >= 0.0,
197 "ridgeRegression(): lambda >= 0.0 required.");
199 unsigned int m = rows;
200 unsigned int n = cols;
205 if(rank < n && lambda == 0.0)
209 for(
unsigned int k=0; k<cols; ++k)
210 for(
unsigned int l=0; l<rhsCount; ++l)
211 t(k,l) *= s(k,0) / (
sq(s(k,0)) + lambda);
253 template <
class T,
class C1,
class C2,
class C3,
class C4>
261 const unsigned int rows =
rowCount(A);
264 vigra_precondition(rows >= cols,
265 "weightedRidgeRegression(): Input matrix A must be rectangular with rowCount >= columnCount.");
266 vigra_precondition(
rowCount(b) == rows,
267 "weightedRidgeRegression(): Shape mismatch between matrices A and b.");
269 "weightedRidgeRegression(): Weight matrix has wrong shape.");
271 "weightedRidgeRegression(): Result matrix x has wrong shape.");
272 vigra_precondition(lambda >= 0.0,
273 "weightedRidgeRegression(): lambda >= 0.0 required.");
277 for(
unsigned int k=0; k<rows; ++k)
279 vigra_precondition(weights(k,0) >= 0,
280 "weightedRidgeRegression(): Weights must be positive.");
282 for(
unsigned int l=0; l<cols; ++l)
283 wa(k,l) = w * A(k,l);
284 for(
unsigned int l=0; l<rhsCount; ++l)
285 wb(k,l) = w * b(k,l);
307 template <
class T,
class C1,
class C2,
class C3,
class Array>
314 const unsigned int rows =
rowCount(A);
316 const unsigned int lambdaCount = lambda.size();
317 vigra_precondition(rows >= cols,
318 "ridgeRegressionSeries(): Input matrix A must be rectangular with rowCount >= columnCount.");
320 "ridgeRegressionSeries(): Shape mismatch between matrices A and b.");
322 "ridgeRegressionSeries(): Result matrix x has wrong shape.");
324 unsigned int m = rows;
325 unsigned int n = cols;
333 for(
unsigned int i=0; i<lambdaCount; ++i)
335 vigra_precondition(lambda[i] >= 0.0,
336 "ridgeRegressionSeries(): lambda >= 0.0 required.");
337 if(lambda[i] == 0.0 && rank < rows)
339 for(
unsigned int k=0; k<cols; ++k)
340 xt(k,0) = xl(k,0) * s(k,0) / (
sq(s(k,0)) + lambda[i]);
354 enum Mode { LARS, LASSO, NNLASSO };
359 : max_solution_count(0),
360 unconstrained_dimension_count(0),
362 least_squares_solutions(true)
374 max_solution_count = (int)n;
387 for(
unsigned int k=0; k<mode.size(); ++k)
388 mode[k] = (std::string::value_type)tolower(mode[k]);
391 else if(mode ==
"lasso")
393 else if(mode ==
"nnlasso")
396 vigra_fail(
"LeastAngleRegressionOptions.setMode(): Invalid mode.");
442 least_squares_solutions = select;
446 int max_solution_count, unconstrained_dimension_count;
448 bool least_squares_solutions;
453 template <
class T,
class C1,
class C2>
461 Matrix<T> R, qtb, lars_solution, lars_prediction, next_lsq_solution, next_lsq_prediction, searchVector;
467 A(Ai), b(bi), R(A), qtb(b),
468 lars_solution(A.shape(1), 1), lars_prediction(A.shape(0), 1),
469 next_lsq_solution(A.shape(1), 1), next_lsq_prediction(A.shape(0), 1), searchVector(A.shape(0), 1),
470 columnPermutation(A.shape(1))
472 for(
unsigned int k=0; k<columnPermutation.size(); ++k)
473 columnPermutation[k] = k;
477 LarsData(LarsData
const & d,
int asetSize)
478 : activeSetSize(asetSize),
479 A(d.R.subarray(Shape(0,0), Shape(d.A.shape(0), activeSetSize))), b(d.qtb), R(A), qtb(b),
480 lars_solution(d.lars_solution.subarray(Shape(0,0), Shape(activeSetSize, 1))), lars_prediction(d.lars_prediction),
481 next_lsq_solution(d.next_lsq_solution.subarray(Shape(0,0), Shape(activeSetSize, 1))),
482 next_lsq_prediction(d.next_lsq_prediction), searchVector(d.searchVector),
483 columnPermutation(A.shape(1))
485 for(
unsigned int k=0; k<columnPermutation.size(); ++k)
486 columnPermutation[k] = k;
490 template <
class T,
class C1,
class C2,
class Array1,
class Array2>
491 unsigned int leastAngleRegressionMainLoop(LarsData<T, C1, C2> & d,
492 Array1 & activeSets, Array2 * lars_solutions, Array2 * lsq_solutions,
493 LeastAngleRegressionOptions
const & options)
495 using namespace vigra::functor;
497 typedef typename MultiArrayShape<2>::type Shape;
498 typedef typename Matrix<T>::view_type Subarray;
499 typedef ArrayVector<MultiArrayIndex> Permutation;
500 typedef typename Permutation::view_type ColumnSet;
502 vigra_precondition(d.activeSetSize > 0,
503 "leastAngleRegressionMainLoop() must not be called with empty active set.");
505 bool enforce_positive = (options.mode == LeastAngleRegressionOptions::NNLASSO);
506 bool lasso_modification = (options.mode != LeastAngleRegressionOptions::LARS);
513 if(maxSolutionCount == 0)
514 maxSolutionCount = lasso_modification
518 bool needToRemoveColumn =
false;
521 while(currentSolutionCount < maxSolutionCount)
523 ColumnSet activeSet = d.columnPermutation.subarray(0, (
unsigned int)d.activeSetSize);
524 ColumnSet inactiveSet = d.columnPermutation.subarray((
unsigned int)d.activeSetSize, (
unsigned int)cols);
527 Matrix<T> cLARS =
transpose(d.A) * (d.b - d.lars_prediction),
528 cLSQ =
transpose(d.A) * (d.b - d.next_lsq_prediction);
536 : argMax(
abs(cLARS));
537 T C =
abs(cLARS(cmaxIndex, 0));
539 Matrix<T> ac(cols - d.activeSetSize, 1);
542 T rho = cLSQ(inactiveSet[k], 0),
543 cc = C -
sign(rho)*cLARS(inactiveSet[k], 0);
548 ac(k,0) = cc / (cc + rho);
549 else if(enforce_positive)
552 ac(k,0) = cc / (cc - rho);
557 if(enforce_positive && needToRemoveColumn)
558 ac(columnToBeRemoved-d.activeSetSize,0) = 1.0;
563 columnToBeAdded =
argMin(ac);
566 T gamma = (d.activeSetSize == maxRank)
568 : ac(columnToBeAdded, 0);
571 if(columnToBeAdded >= 0)
572 columnToBeAdded += d.activeSetSize;
575 needToRemoveColumn =
false;
576 if(lasso_modification)
579 Matrix<T> s(Shape(d.activeSetSize, 1), NumericTraits<T>::max());
582 if(( enforce_positive && d.next_lsq_solution(k,0) < 0.0) ||
583 (!enforce_positive &&
sign(d.lars_solution(k,0))*
sign(d.next_lsq_solution(k,0)) == -1.0))
584 s(k,0) = d.lars_solution(k,0) / (d.lars_solution(k,0) - d.next_lsq_solution(k,0));
587 columnToBeRemoved =
argMinIf(s, Arg1() <= Param(gamma));
588 if(columnToBeRemoved >= 0)
590 needToRemoveColumn =
true;
591 gamma = s(columnToBeRemoved, 0);
596 d.lars_prediction = gamma * d.next_lsq_prediction + (1.0 - gamma) * d.lars_prediction;
597 d.lars_solution = gamma * d.next_lsq_solution + (1.0 - gamma) * d.lars_solution;
598 if(needToRemoveColumn)
599 d.lars_solution(columnToBeRemoved, 0) = 0.0;
602 ++currentSolutionCount;
603 activeSets.push_back(
typename Array1::value_type(d.columnPermutation.begin(), d.columnPermutation.begin()+d.activeSetSize));
605 if(lsq_solutions != 0)
609 ArrayVector<Matrix<T> > nnresults;
610 ArrayVector<ArrayVector<MultiArrayIndex> > nnactiveSets;
611 LarsData<T, C1, C2> nnd(d, d.activeSetSize);
613 leastAngleRegressionMainLoop(nnd, nnactiveSets, &nnresults, (Array2*)0,
614 LeastAngleRegressionOptions().leastSquaresSolutions(
false).nnlasso());
615 Matrix<T> nnlsq_solution(d.activeSetSize, 1);
616 for(
unsigned int k=0; k<nnactiveSets.back().size(); ++k)
618 nnlsq_solution(nnactiveSets.back()[k],0) = nnresults.back()[k];
620 lsq_solutions->push_back(nnlsq_solution);
624 lsq_solutions->push_back(d.next_lsq_solution.subarray(Shape(0,0), Shape(d.activeSetSize, 1)));
627 if(lars_solutions != 0)
629 lars_solutions->push_back(d.lars_solution.subarray(Shape(0,0), Shape(d.activeSetSize, 1)));
636 if(needToRemoveColumn)
639 if(columnToBeRemoved != d.activeSetSize)
643 detail::upperTriangularSwapColumns(columnToBeRemoved, d.activeSetSize, d.R, d.qtb, d.columnPermutation);
646 std::swap(d.lars_solution(columnToBeRemoved, 0), d.lars_solution(d.activeSetSize,0));
647 std::swap(d.next_lsq_solution(columnToBeRemoved, 0), d.next_lsq_solution(d.activeSetSize,0));
648 columnToBeRemoved = d.activeSetSize;
650 d.lars_solution(d.activeSetSize,0) = 0.0;
651 d.next_lsq_solution(d.activeSetSize,0) = 0.0;
655 vigra_invariant(columnToBeAdded >= 0,
656 "leastAngleRegression(): internal error (columnToBeAdded < 0)");
658 if(d.activeSetSize != columnToBeAdded)
660 std::swap(d.columnPermutation[d.activeSetSize], d.columnPermutation[columnToBeAdded]);
662 columnToBeAdded = d.activeSetSize;
666 d.next_lsq_solution(d.activeSetSize,0) = 0.0;
667 d.lars_solution(d.activeSetSize,0) = 0.0;
670 detail::qrColumnHouseholderStep(d.activeSetSize, d.R, d.qtb);
675 Subarray Ractive = d.R.subarray(Shape(0,0), Shape(d.activeSetSize, d.activeSetSize));
676 Subarray qtbactive = d.qtb.subarray(Shape(0,0), Shape(d.activeSetSize, 1));
677 Subarray next_lsq_solution_view = d.next_lsq_solution.subarray(Shape(0,0), Shape(d.activeSetSize, 1));
681 d.next_lsq_prediction.init(0.0);
683 d.next_lsq_prediction += next_lsq_solution_view(k,0)*
columnVector(d.A, d.columnPermutation[k]);
686 return (
unsigned int)currentSolutionCount;
691 template <
class T,
class C1,
class C2,
class Array1,
class Array2>
692 unsigned int leastAngleRegressionMainLoop(LarsData<T, C1, C2> & d,
693 Array1 & activeSets, Array2 * lars_solutions, Array2 * lsq_solutions,
694 LeastAngleRegressionOptions
const & options)
696 using namespace vigra::functor;
698 typedef typename MultiArrayShape<2>::type Shape;
699 typedef typename Matrix<T>::view_type Subarray;
700 typedef ArrayVector<MultiArrayIndex> Permutation;
701 typedef typename Permutation::view_type ColumnSet;
703 vigra_precondition(d.activeSetSize > 0,
704 "leastAngleRegressionMainLoop() must not be called with empty active set.");
706 bool enforce_positive = (options.mode == LeastAngleRegressionOptions::NNLASSO);
707 bool lasso_modification = (options.mode != LeastAngleRegressionOptions::LARS);
714 if(maxSolutionCount == 0)
715 maxSolutionCount = lasso_modification
719 bool needToRemoveColumn =
false;
722 while(currentSolutionCount < maxSolutionCount)
724 ColumnSet activeSet = d.columnPermutation.subarray(0, (
unsigned int)d.activeSetSize);
725 ColumnSet inactiveSet = d.columnPermutation.subarray((
unsigned int)d.activeSetSize, (
unsigned int)cols);
728 Matrix<T> lars_residual = d.b - d.lars_prediction;
729 Matrix<T> lsq_residual = lars_residual - d.searchVector;
730 Matrix<T> c =
transpose(d.A) * lars_residual;
739 T C =
abs(c(cmaxIndex, 0));
741 Matrix<T> ac(cols - d.activeSetSize, 1);
745 a1 = (C - c(inactiveSet[k], 0)) / (C - a),
746 a2 = (C + c(inactiveSet[k], 0)) / (C + a);
750 else if(std::min(a1, a2) < 0.0)
751 ac(k, 0) = std::max(a1, a2);
753 ac(k, 0) = std::min(a1, a2);
758 if(enforce_positive && needToRemoveColumn)
759 ac(columnToBeRemoved-d.activeSetSize,0) = -1.0;
764 columnToBeAdded =
argMinIf(ac, Arg1() < Param(1.0) && Arg1() >= Param(NumericTraits<T>::zero()));
767 T gamma = (columnToBeAdded == -1 || d.activeSetSize == maxRank)
769 : ac(columnToBeAdded, 0);
772 if(columnToBeAdded >= 0)
773 columnToBeAdded += d.activeSetSize;
776 needToRemoveColumn =
false;
777 if(lasso_modification)
780 Matrix<T> s(Shape(d.activeSetSize, 1), NumericTraits<T>::max());
783 if(( enforce_positive && d.next_lsq_solution(k,0) < 0.0) ||
784 (!enforce_positive &&
sign(d.lars_solution(k,0))*
sign(d.next_lsq_solution(k,0)) == -1.0))
785 s(k,0) = d.lars_solution(k,0) / (d.lars_solution(k,0) - d.next_lsq_solution(k,0));
788 columnToBeRemoved =
argMinIf(s, Arg1() < Param(gamma) && Arg1() >= Param(NumericTraits<T>::zero()));
789 if(columnToBeRemoved >= 0)
791 needToRemoveColumn =
true;
792 gamma = s(columnToBeRemoved, 0);
797 d.lars_prediction += gamma * d.searchVector;
798 d.lars_solution = gamma * d.next_lsq_solution + (1.0 - gamma) * d.lars_solution;
799 if(needToRemoveColumn)
800 d.lars_solution(columnToBeRemoved, 0) = 0.0;
803 ++currentSolutionCount;
804 activeSets.push_back(
typename Array1::value_type(d.columnPermutation.begin(), d.columnPermutation.begin()+d.activeSetSize));
806 if(lsq_solutions != 0)
810 ArrayVector<Matrix<T> > nnresults;
811 ArrayVector<ArrayVector<MultiArrayIndex> > nnactiveSets;
812 LarsData<T, C1, C2> nnd(d, d.activeSetSize);
814 leastAngleRegressionMainLoop(nnd, nnactiveSets, &nnresults, (Array2*)0,
815 LeastAngleRegressionOptions().leastSquaresSolutions(
false).nnlasso());
816 Matrix<T> nnlsq_solution(d.activeSetSize, 1);
817 for(
unsigned int k=0; k<nnactiveSets.back().size(); ++k)
819 nnlsq_solution(nnactiveSets.back()[k],0) = nnresults.back()[k];
821 lsq_solutions->push_back(nnlsq_solution);
825 lsq_solutions->push_back(d.next_lsq_solution.subarray(Shape(0,0), Shape(d.activeSetSize, 1)));
828 if(lars_solutions != 0)
830 lars_solutions->push_back(d.lars_solution.subarray(Shape(0,0), Shape(d.activeSetSize, 1)));
837 if(needToRemoveColumn)
840 if(columnToBeRemoved != d.activeSetSize)
844 detail::upperTriangularSwapColumns(columnToBeRemoved, d.activeSetSize, d.R, d.qtb, d.columnPermutation);
847 std::swap(d.lars_solution(columnToBeRemoved, 0), d.lars_solution(d.activeSetSize,0));
848 std::swap(d.next_lsq_solution(columnToBeRemoved, 0), d.next_lsq_solution(d.activeSetSize,0));
849 columnToBeRemoved = d.activeSetSize;
851 d.lars_solution(d.activeSetSize,0) = 0.0;
852 d.next_lsq_solution(d.activeSetSize,0) = 0.0;
856 vigra_invariant(columnToBeAdded >= 0,
857 "leastAngleRegression(): internal error (columnToBeAdded < 0)");
859 if(d.activeSetSize != columnToBeAdded)
861 std::swap(d.columnPermutation[d.activeSetSize], d.columnPermutation[columnToBeAdded]);
863 columnToBeAdded = d.activeSetSize;
867 d.next_lsq_solution(d.activeSetSize,0) = 0.0;
868 d.lars_solution(d.activeSetSize,0) = 0.0;
871 detail::qrColumnHouseholderStep(d.activeSetSize, d.R, d.qtb);
876 Subarray Ractive = d.R.subarray(Shape(0,0), Shape(d.activeSetSize, d.activeSetSize));
877 Subarray qtbactive = d.qtb.subarray(Shape(0,0), Shape(d.activeSetSize, 1));
878 Subarray next_lsq_solution_view = d.next_lsq_solution.subarray(Shape(0,0), Shape(d.activeSetSize, 1));
882 Matrix<T> lsq_prediction(rows, 1);
884 lsq_prediction += next_lsq_solution_view(k,0)*
columnVector(d.A, d.columnPermutation[k]);
887 d.searchVector = lsq_prediction - d.lars_prediction;
890 return (
unsigned int)currentSolutionCount;
894 template <
class T,
class C1,
class C2,
class Array1,
class Array2>
896 leastAngleRegressionImpl(MultiArrayView<2, T, C1>
const & A, MultiArrayView<2, T, C2>
const &b,
897 Array1 & activeSets, Array2 * lasso_solutions, Array2 * lsq_solutions,
898 LeastAngleRegressionOptions
const & options)
900 using namespace vigra::functor;
905 "leastAngleRegression(): Shape mismatch between matrices A and b.");
907 bool enforce_positive = (options.mode == LeastAngleRegressionOptions::NNLASSO);
909 detail::LarsData<T, C1, C2> d(A, b);
916 if(initialColumn == -1)
920 std::swap(d.columnPermutation[0], d.columnPermutation[initialColumn]);
922 detail::qrColumnHouseholderStep(0, d.R, d.qtb);
923 d.next_lsq_solution(0,0) = d.qtb(0,0) / d.R(0,0);
924 d.next_lsq_prediction = d.next_lsq_solution(0,0) *
columnVector(A, d.columnPermutation[0]);
925 d.searchVector = d.next_lsq_solution(0,0) *
columnVector(A, d.columnPermutation[0]);
927 return leastAngleRegressionMainLoop(d, activeSets, lasso_solutions, lsq_solutions, options);
1077 template <
class T,
class C1,
class C2,
class Array1,
class Array2>
1080 Array1 & activeSets, Array2 & solutions,
1081 LeastAngleRegressionOptions
const & options = LeastAngleRegressionOptions())
1083 if(options.least_squares_solutions)
1084 return detail::leastAngleRegressionImpl(A, b, activeSets, (Array2*)0, &solutions, options);
1086 return detail::leastAngleRegressionImpl(A, b, activeSets, &solutions, (Array2*)0, options);
1089 template <
class T,
class C1,
class C2,
class Array1,
class Array2>
1092 Array1 & activeSets, Array2 & lasso_solutions, Array2 & lsq_solutions,
1093 LeastAngleRegressionOptions
const & options = LeastAngleRegressionOptions())
1095 return detail::leastAngleRegressionImpl(A, b, activeSets, &lasso_solutions, &lsq_solutions, options);
1116 template <
class T,
class C1,
class C2,
class C3>
1122 "nonnegativeLeastSquares(): Matrix shape mismatch.");
1124 "nonnegativeLeastSquares(): RHS and solution must be vectors (i.e. columnCount == 1).");
1131 x.
init(NumericTraits<T>::zero());
1132 if(activeSets.
size() > 0)
1133 for(
unsigned int k=0; k<activeSets.
back().size(); ++k)
1134 x(activeSets.
back()[k],0) = results.
back()[k];
1149 using linalg::LeastAngleRegressionOptions;
1153 #endif // VIGRA_REGRESSION_HXX