Eigen::NNLS< MatrixType_ > Class Template Reference

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)

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_

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{\Vert Ax-b\Vert }_2^{2}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 ${\Vert A^N{x}^N-b\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
  };

Code

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()) {}

Code

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);
}

Code

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;
}

Code

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.

{ return info_; }

Code

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.

{ return iterations_; }

Code

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.

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

Code

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());
  }
}

Code

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());
}

Code

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;
  }

Code

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;
  }

Code

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);
    }
  }
}

Code

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());
}

Code

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.

{ return tolerance_; }

Code

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.

{ return x_; }

Code

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