Implementation of the Non-Negative Least Squares (NNLS) algorithm. More...

Public Types
enum	{ RowsAtCompileTime , ColsAtCompileTime , Options , MaxRowsAtCompileTime , MaxColsAtCompileTime }

typedef MatrixType::Index	Index

typedef Matrix< Index, ColsAtCompileTime, 1 >	IndicesType

typedef MatrixType_	MatrixType

typedef MatrixType::RealScalar	RealScalar

typedef Matrix< Scalar, RowsAtCompileTime, 1 >	RhsVectorType

typedef MatrixType::Scalar	Scalar

typedef Matrix< Scalar, ColsAtCompileTime, 1 >	SolutionVectorType

Public Member Functions
template<typename MatrixDerived >
NNLS< MatrixType > &	compute (const EigenBase< MatrixDerived > &A)

ComputationInfo	info () const

Index	iterations () const

Index	maxIterations () const

	NNLS ()

	NNLS (const MatrixType &A, Index max_iter=-1, Scalar tol=NumTraits< Scalar >::dummy_precision())
	Constructs a NNLS sovler and initializes it with the given system matrix `A`. More...

NNLS< MatrixType > &	setMaxIterations (Index maxIters)

NNLS< MatrixType > &	setTolerance (const Scalar &tolerance)

const SolutionVectorType &	solve (const RhsVectorType &b)
	Solves the NNLS problem. More...

Scalar	tolerance () const

const SolutionVectorType &	x () const
	Returns the solution if a problem was solved. If not, an uninitialized vector may be returned. More...

Private Types
typedef Matrix< Scalar, ColsAtCompileTime, ColsAtCompileTime >	MatrixAtAType

Private Member Functions
void	moveToActiveSet_ (Index idx)

void	moveToInactiveSet_ (Index idx)

void	solveInactiveSet_ (const RhsVectorType &b)

Detailed Description

template<class MatrixType_>
class Eigen::NNLS< MatrixType_ >

Implementation of the Non-Negative Least Squares (NNLS) algorithm.

Template Parameters

MatrixType The type of the system matrix \(A\).

This class implements the NNLS algorithm as described in "SOLVING LEAST SQUARES PROBLEMS", Charles L. Lawson and Richard J. Hanson, Prentice-Hall, 1974. This algorithm solves a least squares problem iteratively and ensures that the solution is non-negative. I.e.

\[ \min \left\Vert Ax-b\right\Vert_2^2\quad s.t.\, x\ge 0 \]

The algorithm solves the constrained least-squares problem above by iteratively improving an estimate of which constraints are active (elements of \(x\) equal to zero) and which constraints are inactive (elements of \(x\) greater than zero). Each iteration, an unconstrained linear least-squares problem solves for the components of \(x\) in the (estimated) inactive set and the sets are updated. The unconstrained problem minimizes \(\left\Vert A^Nx^N-b\right\Vert_2^2\), where \(A^N\) is a matrix formed by selecting all columns of A which are in the inactive set \(N\).

See the wikipedia page on non-negative least squares for more background information.

Note: Please note that it is possible to construct an NNLS problem for which the algorithm does not converge. In practice these cases are extremely rare.

Definition at line 66 of file NNLS.

Member Typedef Documentation

◆ Index

template<class MatrixType_ >

typedef MatrixType::Index Eigen::NNLS< MatrixType_ >::Index

Definition at line 80 of file NNLS.

◆ IndicesType

template<class MatrixType_ >

typedef Matrix<Index, ColsAtCompileTime, 1> Eigen::NNLS< MatrixType_ >::IndicesType

Definition at line 86 of file NNLS.

◆ MatrixAtAType

template<class MatrixType_ >

typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime> Eigen::NNLS< MatrixType_ >::MatrixAtAType

private

Definition at line 167 of file NNLS.

◆ MatrixType

template<class MatrixType_ >

typedef MatrixType_ Eigen::NNLS< MatrixType_ >::MatrixType

Definition at line 68 of file NNLS.

◆ RealScalar

template<class MatrixType_ >

typedef MatrixType::RealScalar Eigen::NNLS< MatrixType_ >::RealScalar

Definition at line 79 of file NNLS.

◆ RhsVectorType

template<class MatrixType_ >

typedef Matrix<Scalar, RowsAtCompileTime, 1> Eigen::NNLS< MatrixType_ >::RhsVectorType

Type of a column vector of the system matrix \(A\).

Definition at line 85 of file NNLS.

◆ Scalar

template<class MatrixType_ >

typedef MatrixType::Scalar Eigen::NNLS< MatrixType_ >::Scalar

Definition at line 78 of file NNLS.

◆ SolutionVectorType

template<class MatrixType_ >

typedef Matrix<Scalar, ColsAtCompileTime, 1> Eigen::NNLS< MatrixType_ >::SolutionVectorType

Type of a row vector of the system matrix \(A\).

Definition at line 83 of file NNLS.

Member Enumeration Documentation

◆ anonymous enum

template<class MatrixType_ >

anonymous enum

Enumerator
RowsAtCompileTime
ColsAtCompileTime
Options
MaxRowsAtCompileTime
MaxColsAtCompileTime

Definition at line 70 of file NNLS.

        {
     RowsAtCompileTime = MatrixType::RowsAtCompileTime,
     ColsAtCompileTime = MatrixType::ColsAtCompileTime,
     Options = MatrixType::Options,
     MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
     MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
   };

Constructor & Destructor Documentation

◆ NNLS() [1/2]

template<typename MatrixType >

Eigen::NNLS< MatrixType >::NNLS

Definition at line 209 of file NNLS.

     : max_iter_(-1),
       iterations_(0),
       info_(ComputationInfo::InvalidInput),
       numInactive_(0),
       tolerance_(NumTraits<Scalar>::dummy_precision()) {}

◆ NNLS() [2/2]

template<typename MatrixType >

Eigen::NNLS< MatrixType >::NNLS	(	const MatrixType &	A,
		Index	max_iter = `-1`,
		Scalar	tol = `NumTraits<Scalar>::dummy_precision()`
	)

Constructs a NNLS sovler and initializes it with the given system matrix A.

Parameters

A	Specifies the system matrix.
max_iter	Specifies the maximum number of iterations to solve the system.
tol	Specifies the precision of the optimum. This is an absolute tolerance on the gradient of the Lagrangian, \(A^T(Ax-b)-\lambda\) (with Lagrange multipliers \(\lambda\)).

Definition at line 217 of file NNLS.

                                                                       : max_iter_(max_iter), tolerance_(tol) {
   compute(A);
 }

Member Function Documentation

◆ compute()

template<typename MatrixType >

template<typename MatrixDerived >

NNLS< MatrixType > & Eigen::NNLS< MatrixType >::compute ( const EigenBase< MatrixDerived > & A )

Initializes the solver with the matrix A for further solving NNLS problems.

This function mostly initializes/computes the preconditioner. In the future we might, for instance, implement column reordering for faster matrix vector products.

Definition at line 223 of file NNLS.

                                                                              {
   // Ensure Scalar type is real. The non-negativity constraint doesn't obviously extend to complex numbers.
   EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
  
   // max_iter_: unchanged
   iterations_ = 0;
   info_ = ComputationInfo::Success;
   numInactive_ = 0;
   // tolerance: unchanged
   A_ = A.derived();
   AtA_.noalias() = A_.transpose() * A_;
   x_.resize(A_.cols());
   gradient_.resize(A_.cols());
   y_.resize(A_.cols());
   Atb_.resize(A_.cols());
   index_sets_.resize(A_.cols());
   QR_.resize(A_.rows(), A_.cols());
   qrCoeffs_.resize(A_.cols());
   tempSolutionVector_.resize(A_.cols());
   tempRhsVector_.resize(A_.rows());
  
   return *this;
 }

◆ info()

template<class MatrixType_ >

ComputationInfo Eigen::NNLS< MatrixType_ >::info ( ) const

inline

Returns: Success if the iterations converged, and an error values otherwise.

Definition at line 154 of file NNLS.

154 { return info_; }

◆ iterations()

template<class MatrixType_ >

Index Eigen::NNLS< MatrixType_ >::iterations ( ) const

inline

Returns: the number of iterations (least-squares solves) performed during the last solve

Definition at line 151 of file NNLS.

151 { return iterations_; }

◆ maxIterations()

template<class MatrixType_ >

Index Eigen::NNLS< MatrixType_ >::maxIterations ( ) const

inline

Returns: the max number of iterations. It is either the value set by setMaxIterations or, by default, twice the number of columns of the matrix.

Definition at line 139 of file NNLS.

139 { return max_iter_ < 0 ? 2 * A_.cols() : max_iter_; }

◆ moveToActiveSet_()

template<typename MatrixType >

void Eigen::NNLS< MatrixType >::moveToActiveSet_ ( Index idx )

private

Definition at line 350 of file NNLS.

                                                  {
   // swap index with last inactive one & reduce number of inactive columns
   std::swap(index_sets_(idx), index_sets_(numInactive_ - 1));
   numInactive_--;
   // Update QR decomposition starting from the removed index up to the end [idx, ..., numInactive_]
   for (Index i = idx; i < numInactive_; i++) {
     Index col = index_sets_(i);
     internal::householder_qr_inplace_update(QR_, qrCoeffs_, A_.col(col), i, tempSolutionVector_.data());
   }
 }

◆ moveToInactiveSet_()

template<typename MatrixType >

void Eigen::NNLS< MatrixType >::moveToInactiveSet_ ( Index idx )

private

Definition at line 339 of file NNLS.

                                                    {
   // Update permutation matrix:
   std::swap(index_sets_(idx), index_sets_(numInactive_));
   numInactive_++;
  
   // Perform rank-1 update of the QR decomposition stored in QR_ & qrCoeff_
   internal::householder_qr_inplace_update(QR_, qrCoeffs_, A_.col(index_sets_(numInactive_ - 1)), numInactive_ - 1,
                                           tempSolutionVector_.data());
 }

◆ setMaxIterations()

template<class MatrixType_ >

NNLS<MatrixType>& Eigen::NNLS< MatrixType_ >::setMaxIterations ( Index maxIters )

inline

Sets the max number of iterations. Default is twice the number of columns of the matrix. The algorithm requires at least k iterations to produce a solution vector with k non-zero entries.

Definition at line 145 of file NNLS.

                                                      {
     max_iter_ = maxIters;
     return *this;
   }

◆ setTolerance()

template<class MatrixType_ >

NNLS<MatrixType>& Eigen::NNLS< MatrixType_ >::setTolerance ( const Scalar & tolerance )

inline

Sets the tolerance threshold used by the stopping criteria.

This is an absolute tolerance on the gradient of the Lagrangian, \(A^T(Ax-b)-\lambda\) (with Lagrange multipliers \(\lambda\)).

Definition at line 131 of file NNLS.

                                                           {
     tolerance_ = tolerance;
     return *this;
   }

◆ solve()

template<typename MatrixType >

const NNLS< MatrixType >::SolutionVectorType & Eigen::NNLS< MatrixType >::solve ( const RhsVectorType & b )

Solves the NNLS problem.

The dimension of b must be equal to the number of rows of A, given to the constructor.

Returns: The approximate solution vector \( x \). Use info() to determine if the solve was a success or not.

See also: info()

Definition at line 248 of file NNLS.

                                                                                                  {
   // Initialize solver
   iterations_ = 0;
   info_ = ComputationInfo::NumericalIssue;
   x_.setZero();
  
   index_sets_ = IndicesType::LinSpaced(A_.cols(), 0, A_.cols() - 1);  // Identity permutation.
   numInactive_ = 0;
  
   // Precompute A^T*b
   Atb_.noalias() = A_.transpose() * b;
  
   const Index maxIterations = this->maxIterations();
  
   // OUTER LOOP
   while (true) {
     // Early exit if all variables are inactive, which breaks 'maxCoeff' below.
     if (A_.cols() == numInactive_) {
       info_ = ComputationInfo::Success;
       return x_;
     }
  
     // Find the maximum element of the gradient in the active set.
     // If it is small or negative, then we have converged.
     // Else, we move that variable to the inactive set.
     gradient_.noalias() = Atb_ - AtA_ * x_;
  
     const Index numActive = A_.cols() - numInactive_;
     Index argmaxGradient = -1;
     const Scalar maxGradient = gradient_(index_sets_.tail(numActive)).maxCoeff(&argmaxGradient);
     argmaxGradient += numInactive_;  // beacause tail() skipped the first numInactive_ elements
  
     if (maxGradient < tolerance_) {
       info_ = ComputationInfo::Success;
       return x_;
     }
  
     moveToInactiveSet_(argmaxGradient);
  
     // INNER LOOP
     while (true) {
       // Check if max. number of iterations is reached
       if (iterations_ >= maxIterations) {
         info_ = ComputationInfo::NoConvergence;
         return x_;
       }
  
       // Solve least-squares problem in inactive set only,
       // this step is rather trivial as moveToInactiveSet_ & moveToActiveSet_
       // updates the QR decomposition of inactive columns A^N.
       // solveInactiveSet_ puts the solution in y_
       solveInactiveSet_(b);
       ++iterations_;  // The solve is expensive, so that is what we count as an iteration.
  
       // Check feasability...
       bool feasible = true;
       Scalar alpha = NumTraits<Scalar>::highest();
       Index infeasibleIdx = -1;  // Which variable became infeasible first.
       for (Index i = 0; i < numInactive_; i++) {
         Index idx = index_sets_[i];
         if (y_(idx) < 0) {
           // t should always be in [0,1].
           Scalar t = -x_(idx) / (y_(idx) - x_(idx));
           if (alpha > t) {
             alpha = t;
             infeasibleIdx = i;
             feasible = false;
           }
         }
       }
       eigen_assert(feasible || 0 <= infeasibleIdx);
  
       // If solution is feasible, exit to outer loop
       if (feasible) {
         x_ = y_;
         break;
       }
  
       // Infeasible solution -> interpolate to feasible one
       for (Index i = 0; i < numInactive_; i++) {
         Index idx = index_sets_[i];
         x_(idx) += alpha * (y_(idx) - x_(idx));
       }
  
       // Remove these indices from the inactive set and update QR decomposition
       moveToActiveSet_(infeasibleIdx);
     }
   }
 }

◆ solveInactiveSet_()

template<typename MatrixType >

void Eigen::NNLS< MatrixType >::solveInactiveSet_ ( const RhsVectorType & b )

private

Definition at line 362 of file NNLS.

                                                                {
   eigen_assert(numInactive_ > 0);
  
   tempRhsVector_ = b;
  
   // tmpRHS(0:numInactive_-1) := Q'*b
   // tmpRHS(numInactive_:end) := useless stuff we would rather not compute at all.
   tempRhsVector_.applyOnTheLeft(
       householderSequence(QR_.leftCols(numInactive_), qrCoeffs_.head(numInactive_)).transpose());
  
   // tempSol(0:numInactive_-1) := inv(R) * Q' * b
   //  = the least-squares solution for the inactive variables.
   tempSolutionVector_.head(numInactive_) =            //
       QR_.topLeftCorner(numInactive_, numInactive_)   //
           .template triangularView<Upper>()           //
           .solve(tempRhsVector_.head(numInactive_));  //
  
   // tempSol(numInactive_:end) := 0 = the value for the constrained variables.
   tempSolutionVector_.tail(y_.size() - numInactive_).setZero();
  
   // Back permute into original column order of A
   y_.noalias() = index_sets_.asPermutation() * tempSolutionVector_.head(y_.size());
 }

◆ tolerance()

template<class MatrixType_ >

Scalar Eigen::NNLS< MatrixType_ >::tolerance ( ) const

inline

Returns: the tolerance threshold used by the stopping criteria.

See also: setTolerance()

Definition at line 124 of file NNLS.

124 { return tolerance_; }

◆ x()

template<class MatrixType_ >

const SolutionVectorType& Eigen::NNLS< MatrixType_ >::x ( ) const

inline

Returns the solution if a problem was solved. If not, an uninitialized vector may be returned.

Definition at line 119 of file NNLS.

119 { return x_; }

Member Data Documentation

◆ A_

template<class MatrixType_ >

MatrixType Eigen::NNLS< MatrixType_ >::A_

private

Definition at line 181 of file NNLS.

◆ AtA_

template<class MatrixType_ >

MatrixAtAType Eigen::NNLS< MatrixType_ >::AtA_

private

Definition at line 183 of file NNLS.

◆ Atb_

template<class MatrixType_ >

SolutionVectorType Eigen::NNLS< MatrixType_ >::Atb_

private

Definition at line 191 of file NNLS.

◆ gradient_

template<class MatrixType_ >

SolutionVectorType Eigen::NNLS< MatrixType_ >::gradient_

private

Definition at line 187 of file NNLS.

◆ index_sets_

template<class MatrixType_ >

IndicesType Eigen::NNLS< MatrixType_ >::index_sets_

private

Definition at line 194 of file NNLS.

◆ info_

template<class MatrixType_ >

ComputationInfo Eigen::NNLS< MatrixType_ >::info_

private

Definition at line 175 of file NNLS.

◆ iterations_

template<class MatrixType_ >

Index Eigen::NNLS< MatrixType_ >::iterations_

private

Definition at line 173 of file NNLS.

◆ max_iter_

template<class MatrixType_ >

Index Eigen::NNLS< MatrixType_ >::max_iter_

private

Definition at line 171 of file NNLS.

◆ numInactive_

template<class MatrixType_ >

Index Eigen::NNLS< MatrixType_ >::numInactive_

private

Definition at line 177 of file NNLS.

◆ QR_

template<class MatrixType_ >

MatrixType Eigen::NNLS< MatrixType_ >::QR_

private

Definition at line 196 of file NNLS.

◆ qrCoeffs_

template<class MatrixType_ >

SolutionVectorType Eigen::NNLS< MatrixType_ >::qrCoeffs_

private

Definition at line 198 of file NNLS.

◆ tempRhsVector_

template<class MatrixType_ >

RhsVectorType Eigen::NNLS< MatrixType_ >::tempRhsVector_

private

Definition at line 201 of file NNLS.

◆ tempSolutionVector_

template<class MatrixType_ >

SolutionVectorType Eigen::NNLS< MatrixType_ >::tempSolutionVector_

private

Definition at line 200 of file NNLS.

◆ tolerance_

template<class MatrixType_ >

Scalar Eigen::NNLS< MatrixType_ >::tolerance_

private

Definition at line 179 of file NNLS.

◆ x_

template<class MatrixType_ >

SolutionVectorType Eigen::NNLS< MatrixType_ >::x_

private

Definition at line 185 of file NNLS.

◆ y_

template<class MatrixType_ >

SolutionVectorType Eigen::NNLS< MatrixType_ >::y_

private

Definition at line 189 of file NNLS.

The documentation for this class was generated from the following file:

NNLS

Private Attributes
MatrixType	A_

MatrixAtAType	AtA_

SolutionVectorType	Atb_

SolutionVectorType	gradient_

IndicesType	index_sets_

ComputationInfo	info_

Index	iterations_

Index	max_iter_

Index	numInactive_

MatrixType	QR_

SolutionVectorType	qrCoeffs_

RhsVectorType	tempRhsVector_

SolutionVectorType	tempSolutionVector_

Scalar	tolerance_

SolutionVectorType	x_

SolutionVectorType	y_

Public Types

Public Member Functions

Private Types

Private Member Functions

Private Attributes

Detailed Description

template<class MatrixType_> class Eigen::NNLS< MatrixType_ >

Member Typedef Documentation

◆ Index

◆ IndicesType

◆ MatrixAtAType

◆ MatrixType

◆ RealScalar

◆ RhsVectorType

◆ Scalar

◆ SolutionVectorType

Member Enumeration Documentation

◆ anonymous enum

Constructor & Destructor Documentation

◆ NNLS() [1/2]

◆ NNLS() [2/2]

Member Function Documentation

◆ compute()

◆ info()

◆ iterations()

◆ maxIterations()

◆ moveToActiveSet_()

◆ moveToInactiveSet_()

◆ setMaxIterations()

◆ setTolerance()

◆ solve()

◆ solveInactiveSet_()

◆ tolerance()

◆ x()

Member Data Documentation

◆ A_

◆ AtA_

◆ Atb_

◆ gradient_

◆ index_sets_

◆ info_

◆ iterations_

◆ max_iter_

◆ numInactive_

◆ QR_

◆ qrCoeffs_

◆ tempRhsVector_

◆ tempSolutionVector_

◆ tolerance_

◆ x_

◆ y_

template<class MatrixType_>
class Eigen::NNLS< MatrixType_ >