IMSL_ITMAX, intitmax (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.
IMSL_DROP_MAX_POS_DUAL, intmaxdual (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 maxdual = 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: maxdual = 0
IMSL_DROP_TOLERANCE, float tol (Input) This is a rank-determination tolerance. A candidate column
has values eliminated below the first entry of . The resulting value must satisfy the relative condition
Otherwise the constraint remains satisfied because the column is linearly dependent on previously dropped columns. Default: tol=sqrt(imsl_f_machine(3));
IMSL_SUPPLY_WORK_ARRAYS , intlwork, floatwork[], intliwork, intiwork[] (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:
lworkmaxt*(m*(n+2) + n), and liworkmaxt*n.
Without the use of OpenMP and parallel threading, maxt=1.
IMSL_OPTIMIZED, int *flag (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.
flag
Description
0
the maximum number of iterations was reached.
1
the optimum residual norm was obtained.
IMSL_DUAL_SOLUTION, float**dual (Output) An array of length n containing the dual vector, . This may not be optimal (all components may not satisfy ), if the maximum number of iterations occurs first.
IMSL_DUAL_SOLUTION_USER, floatdual[] (Output) Storage for dual provided by the user. See IMSL_DUAL_SOLUTION.
IMSL_RESIDUAL_NORM, float*rnorm (Output) The value of the residual vector norm, ∥Ax-b∥2.
IMSL_RETURN_USER, floatx[] (Output) A user-allocated array of length n containing the approximate solution vector, .
Description
Function imsl_f_nonneg_least_squares computes the constrained least squares solution of , by minimizing ∥Ax-b∥2 subject to . 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.
Examples
Example 1
A model function of exponentials is
The exponential function argument parameters
are fixed. The coefficients
are estimated by sampling data values,
using non-negative least squares. The values used for the data are
with
#include <imsl.h>
#include <math.h>
#define M 21
#define N 3
int main() {
int i;
float a[M][N], b[M], *c;
for (i = 0; i < M; i++) {
/* 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 = imsl_f_nonneg_least_squares(M, N, &a[0][0], b, 0);
/* With noise level = 0, solution should be (1, 0.2, 0.3) */
imsl_f_write_matrix("Coefficients", 1, N, c, 0);
}
Output
Coefficients
1 2 3
1.0 0.2 0.3
Example 2
The model function of exponentials is
The values λ2,λ3 are the same as in Example 1. The function n (t) represents normally distributed random noise with a standard deviation . A simulation is done with ns = 10001 samples for n (t). The resulting problem is solved using OpenMP. To check that the OpenMP results are correct, a loop computes the solutions without OpenMP followed by the same loop using OpenMP. The residual norms agree, showing that the routine returns the same values using OpenMP as without using OpenMP.
#include <imsl.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>
#define M 21
#define N 3
#define NS 10001
int main() {
#define BS(i_,j_) bs[(i_)*M + (j_)]
#define X(i_,j_) x[(i_)*N + (j_)]
int thread_safe=1, seed=123457, i, *iwork, j, lwork, liwork, maxt;
