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)
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;
689 template <
class T,
class C1,
class C2,
class Array1,
class Array2>
691 leastAngleRegressionImpl(MultiArrayView<2, T, C1>
const & A, MultiArrayView<2, T, C2>
const &b,
692 Array1 & activeSets, Array2 * lasso_solutions, Array2 * lsq_solutions,
693 LeastAngleRegressionOptions
const & options)
695 using namespace vigra::functor;
700 "leastAngleRegression(): Shape mismatch between matrices A and b.");
702 bool enforce_positive = (options.mode == LeastAngleRegressionOptions::NNLASSO);
704 detail::LarsData<T, C1, C2> d(A, b);
711 if(initialColumn == -1)
715 std::swap(d.columnPermutation[0], d.columnPermutation[initialColumn]);
717 detail::qrColumnHouseholderStep(0, d.R, d.qtb);
718 d.next_lsq_solution(0,0) = d.qtb(0,0) / d.R(0,0);
719 d.next_lsq_prediction = d.next_lsq_solution(0,0) *
columnVector(A, d.columnPermutation[0]);
720 d.searchVector = d.next_lsq_solution(0,0) *
columnVector(A, d.columnPermutation[0]);
722 return leastAngleRegressionMainLoop(d, activeSets, lasso_solutions, lsq_solutions, options);
872 template <
class T,
class C1,
class C2,
class Array1,
class Array2>
875 Array1 & activeSets, Array2 & solutions,
876 LeastAngleRegressionOptions
const & options = LeastAngleRegressionOptions())
878 if(options.least_squares_solutions)
879 return detail::leastAngleRegressionImpl(A, b, activeSets, (Array2*)0, &solutions, options);
881 return detail::leastAngleRegressionImpl(A, b, activeSets, &solutions, (Array2*)0, options);
884 template <
class T,
class C1,
class C2,
class Array1,
class Array2>
887 Array1 & activeSets, Array2 & lasso_solutions, Array2 & lsq_solutions,
888 LeastAngleRegressionOptions
const & options = LeastAngleRegressionOptions())
890 return detail::leastAngleRegressionImpl(A, b, activeSets, &lasso_solutions, &lsq_solutions, options);
911 template <
class T,
class C1,
class C2,
class C3>
917 "nonnegativeLeastSquares(): Matrix shape mismatch.");
919 "nonnegativeLeastSquares(): RHS and solution must be vectors (i.e. columnCount == 1).");
926 x.
init(NumericTraits<T>::zero());
927 if(activeSets.
size() > 0)
928 for(
unsigned int k=0; k<activeSets.
back().size(); ++k)
929 x(activeSets.
back()[k],0) = results.
back()[k];
944 using linalg::LeastAngleRegressionOptions;
948 #endif // VIGRA_REGRESSION_HXX