nonnegLeastSquares¶

Compute the non-negative least squares (NNLS) solution of an m × n real linear least squares system, \(Ax \cong b\) , \(x \geq 0\) .

Synopsis¶

nonnegLeastSquares (a, b)

Required Arguments¶

float a[[]] (Input): An array of length m × n containing the matrix.
float b[] (Input): An array of length m containing the right-hand side vector.

Return Value¶

An array of length n containing the approximate solution vector, x ≥ 0.

Optional Arguments¶

itMax, int (Input)

The number of times a constraint is added or dropped should not exceed this maximum value. An approximate solution x ≥ 0 is returned when the maximum number is reached.

Default: itMax = 3 × n.

dropMaxPosDual, int (Input)

Indicates how a variable is moved from its constraint to a positive value, or dropped, when its current dual value is positive. By dropping the variable corresponding to the first computed positive dual value, instead of the maximum, better runtime efficiency usually results by avoiding work in the early stages of the algorithm.

If dropMaxPosDual = 0, the first encountered positive dual is used. Otherwise, the maximum positive dual, is used. The results for x ≥ 0 will usually vary slightly depending on the choice.

Default: dropMaxPosDual = 0

dropTolerance, float (Input)

This is a rank-determination tolerance. A candidate column

\[\begin{split}a = \begin{bmatrix} c \\ d \end{bmatrix}\end{split}\]

has values eliminated below the first entry of \(d\) . The resulting value must satisfy the relative condition

\[\|d\|_2 > \mathit{tol} \times \|c\|_2\]

Otherwise the constraint remains satisfied because the column \(a\) is linearly dependent on previously dropped columns.

Default: dropTolerance = sqrt(machine(3))

supplyWorkArrays , int lwork, float work[], int liwork, int iwork[] (Input/Output)

The use of this optional argument will increase efficiency and avoid memory fragmentation run-time failures for large problems by allowing the user to provide the sizes and locations of the working arrays work and iwork. With maxT as the maximum number of threads that will be active, it is required that:

lwork \(\geq\) maxT * (m * (n+2) + n), and liwork \(\geq\) maxT * n.

Without the use of OpenMP and parallel threading, maxT=1.

optimized (Output)

A 0-1 flag noting whether or not the optimum residual norm was obtained. A value of 1 indicates the optimum residual norm was obtained. A value of 0 occurs if the maximum number of iterations was reached.

optimized	Description
0	the maximum number of iterations was reached.
1	the optimum residual norm was obtained.

dualSolution, (Output): An array of length n containing the dual vector, \(w = A^T (Ax - b)\). This may not be optimal (all components may not satisfy \(w \leq 0\) ), if the maximum number of iterations occurs first.
residualNorm (Output): The value of the residual vector norm, \(\|Ax-b\|_2\).

Description¶

Function nonnegLeastSquares computes the constrained least squares solution of \(Ax \cong b\), by minimizing \(\|Ax-b\|_2\) subject to \(x \geq 0\) . It uses the algorithm NNLS found in Charles L. Lawson and Richard J. Hanson, Solving Least Squares Problems, SIAM Publications, Chap. 23, (1995). The functionality for multiple threads and the constraint dropping strategy are new features. The original NNLS algorithm was silent about multiple threads; all dual components were computed when only one was used. Using the first encountered eligible variable to make non-active usually improves performance. An optimum solution is obtained in either approach. There is no restriction on the relative sizes of m and n.

Example¶

A model function of exponentials is

\[f(t) = c_1 + c_2 \exp\left(-\lambda_2t\right) + c_3 \exp\left(-\lambda_3t\right), t \geq 0\]

The exponential function argument parameters

\[\lambda_2 = 1, \lambda_3 = 5\]

are fixed. The coefficients

\[c_j \geq 0, j = 1, 2, 3\]

are estimated by sampling data values,

\[f\left(t_i\right), i = 1, \ldots 21\]

using non-negative least squares. The values used for the data are

\[t_i = 0.25i, i = 0, \ldots 20\]

with

\[c_1 = 1, c_2 = 0.2, c_3 = 0.3\]

from numpy import *
from pyimsl.math.nonnegLeastSquares import nonnegLeastSquares
from pyimsl.math.writeMatrix import writeMatrix

m = 21
n = 3
a = zeros((m, n), dtype=double)
b = zeros((m), dtype=double)
for i in range(m):
    # Generate exponential values. This model is
    # y(t) = c_0 + c_1*exp(-t) + c_2*exp(-5*t)
    a[i][0] = 1.0
    a[i][1] = exp(-(i * 0.25))
    a[i][2] = exp(-(i * 0.25) * 5.0)
    # Compute sample values
    b[i] = a[i][0] + 0.2 * a[i][1] + 0.3 * a[i][2]

# Solve for coefficients, constraining values
# to be non-negative. */
c = nonnegLeastSquares(a, b)

# With noise level = 0, solution should be (1, 0.2, 0.3)
writeMatrix("Coefficients", c, False)

Output¶

 
            Coefficients
          1            2            3
        1.0          0.2          0.3

Warning Errors¶

IMSL_MAX_NNLS_ITER_REACHED The maximum number of iterations was reached. The best answer will be returned. “itMax” = # was used. A larger value may help the algorithm complete.