NNLS (non-negative least squares) and required subroutines. More...

#include "libtpcmodel.h"

Functions
int	nnls (double *a, int m, int n, double b, double x, double rnorm, double wp, double zzp, int *indexp)

int	nnlsWght (int N, int M, double *A, double b, double *weight)

int	nnlsWghtSquared (int N, int M, double *A, double b, double *sweight)

Detailed Description

NNLS (non-negative least squares) and required subroutines.

Author: Vesa Oikonen, Kaisa Sederholm

This routine is based on the text and Fortran code in C.L. Lawson and R.J. Hanson, Solving Least Squares Problems, Prentice-Hall, Englewood Cliffs, New Jersey, 1974.

Definition in file nnls.c.

Function Documentation

◆ nnls()

int nnls	(	double **	a,
		int	m,
		int	n,
		double *	b,
		double *	x,
		double *	rnorm,
		double *	wp,
		double *	zzp,
		int *	indexp )

Algorithm NNLS (Non-negative least-squares).

Given an m by n matrix A, and an m-vector B, computes an n-vector X, that solves the least squares problem A * X = B , subject to X>=0

Instead of pointers for working space, NULL can be given to let this function to allocate and free the required memory.

See also: qrLSQ, nnlsWghtSquared

Returns: Function returns 0 if successful, 1, if iteration count exceeded 3*N, or 2 in case of invalid problem dimensions or memory allocation error.

Parameters

a	On entry, a[ 0... N ][ 0 ... M ] contains the M by N matrix A. On exit, a[][] contains the product matrix Q*A, where Q is an m by n orthogonal matrix generated implicitly by this function.
m	Matrix dimension m
n	Matrix dimension n
b	On entry, b[] must contain the m-vector B. On exit, b[] contains Q*B
x	On exit, x[] will contain the solution vector
rnorm	On exit, rnorm contains the squared Euclidean norm of the residual vector, R^2 (Sum-of-Squares). If NULL is given, no rnorm is calculated
wp	An n-array of working space, wp[]. On exit, wp[] will contain the dual solution vector. wp[i]=0.0 for all i in set p and wp[i]<=0.0 for all i in set z. Can be NULL, which causes this algorithm to allocate memory for it.
zzp	An m-array of working space, zz[]. Can be NULL, which causes this algorithm to allocate memory for it.
indexp	An n-array of working space, index[]. Can be NULL, which causes this algorithm to allocate memory for it.

Definition at line 37 of file nnls.c.

  {
  /* Check the parameters and data */
  if(m<=0 || n<=0 || a==NULL || b==NULL || x==NULL) return(2);
  if(rnorm!=NULL) *rnorm=nan("");
  /* Allocate memory for working space, if required */
  int *index=NULL, *indexl=NULL;
  double *w=NULL, *zz=NULL, *wl=NULL, *zzl=NULL;
  if(wp!=NULL) w=wp; else w=wl=(double*)calloc(n, sizeof(double));
  if(zzp!=NULL) zz=zzp; else zz=zzl=(double*)calloc(m, sizeof(double));
  if(indexp!=NULL) index=indexp; else index=indexl=(int*)calloc(n, sizeof(int));
  if(w==NULL || zz==NULL || index==NULL) return(2);
 
  /* Initialize the arrays INDEX[] and X[] */
  for(int ni=0; ni<n; ni++) {x[ni]=0.; index[ni]=ni;}
  int iz1=0;
  int iz2=n-1;
  int nsetp=0;
  int npp1=0;
 
  /* Main loop; quit if all coefficients are already in the solution or
     if M cols of A have been triangulated */
  double up=0.0;
  int itmax; if(n<3) itmax=n*3; else itmax=n*n;
  int iter=0; 
  int k, j=0, jj=0;
  while(iz1<=iz2 && nsetp<m) {
    /* Compute components of the dual (negative gradient) vector W[] */
    for(int iz=iz1; iz<=iz2; iz++) {
      int ni=index[iz];
      double sm=0.;
      for(int mi=npp1; mi<m; mi++) sm+=a[ni][mi]*b[mi];
      w[ni]=sm;
    }
 
    double wmax;
    int izmax=0;
    while(1) {
 
      /* Find largest positive W[j] */
      wmax=0.0;
      for(int iz=iz1; iz<=iz2; iz++) {
        int i=index[iz];
        if(w[i]>wmax) {wmax=w[i]; izmax=iz;}
      }
 
      /* Terminate if wmax<=0.; */
      /* it indicates satisfaction of the Kuhn-Tucker conditions */
      if(wmax<=0.0) break;
      j=index[izmax];
 
      /* The sign of W[j] is ok for j to be moved to set P.
         Begin the transformation and check new diagonal element to avoid
         near linear dependence. */
      double asave=a[j][npp1];
      up=0.0;
      _lss_h12(1, npp1, npp1+1, m, &a[j][0], 1, &up, NULL, 1, 1, 0);
      double unorm=0.0;
      if(nsetp!=0) for(int mi=0; mi<nsetp; mi++) unorm+=a[j][mi]*a[j][mi];
      unorm=sqrt(unorm);
      double d=unorm+fabs(a[j][npp1])*0.01;
      if((d-unorm)>0.0) {
        /* Col j is sufficiently independent. Copy B into ZZ, update ZZ
           and solve for ztest ( = proposed new value for X[j] ) */
        for(int mi=0; mi<m; mi++) zz[mi]=b[mi];
        _lss_h12(2, npp1, npp1+1, m, &a[j][0], 1, &up, zz, 1, 1, 1);
        double ztest=zz[npp1]/a[j][npp1];
        /* See if ztest is positive */
        if(ztest>0.) break;
      }
 
      /* Reject j as a candidate to be moved from set Z to set P. Restore
         A[npp1,j], set W[j]=0., and loop back to test dual coefficients again */
      a[j][npp1]=asave; w[j]=0.;
    } /* while(1) */
    if(wmax<=0.0) break;
 
    /* Index j=INDEX[izmax] has been selected to be moved from set Z to set P.
       Update B and indices, apply householder transformations to cols in
       new set Z, zero sub-diagonal elements in col j, set W[j]=0. */
    for(int mi=0; mi<m; mi++) b[mi]=zz[mi];
    index[izmax]=index[iz1]; index[iz1]=j; iz1++; nsetp=npp1+1; npp1++;
    if(iz1<=iz2)
      for(int jz=iz1; jz<=iz2; jz++) {
        jj=index[jz];
        _lss_h12(2, nsetp-1, npp1, m, &a[j][0], 1, &up, &a[jj][0], 1, m, 1);
      }
    if(nsetp!=m) for(int mi=npp1; mi<m; mi++) a[j][mi]=0.;
    w[j]=0.;
    
    /* Solve the triangular system; store the solution temporarily in Z[] */
    for(int mi=0; mi<nsetp; mi++) {
      int ip=nsetp-(mi+1);
      if(mi!=0) for(int ii=0; ii<=ip; ii++) zz[ii]-=a[jj][ii]*zz[ip+1];
      jj=index[ip]; zz[ip]/=a[jj][ip];
    }
 
    /* Secondary loop begins here */
    while(++iter<itmax) {
      /* See if all new constrained coefficients are feasible; if not, compute alpha */
      double alpha=2.0;
      for(int ip=0; ip<nsetp; ip++) {
        int ni=index[ip];
        if(zz[ip]<=0.) {
          double t=-x[ni]/(zz[ip]-x[ni]);
          if(alpha>t) {alpha=t; jj=ip-1;}
        }
      }
 
      /* If all new constrained coefficients are feasible then still alpha==2.
         If so, then exit from the secondary loop to main loop */
      if(alpha==2.0) break;
 
      /* Use alpha (0.<alpha<1.) to interpolate between old X and new ZZ */
      for(int ip=0; ip<nsetp; ip++) {
        int ni=index[ip]; x[ni]+=alpha*(zz[ip]-x[ni]);
      }
 
      /* Modify A and B and the INDEX arrays to move coefficient i from set P to set Z. */
      int pfeas=1;
      k=index[jj+1];
      do {
        x[k]=0.;
        if(jj!=(nsetp-1)) {
          jj++;
          for(int ni=jj+1; ni<nsetp; ni++) {
            int ii=index[ni]; index[ni-1]=ii;
            double ss, cc;
            _lss_g1(a[ii][ni-1], a[ii][ni], &cc, &ss, &a[ii][ni-1]);
            a[ii][ni]=0.0;
            for(int nj=0; nj<n; nj++) if(nj!=ii) {
              /* Apply procedure G2 (CC,SS,A(J-1,L),A(J,L)) */
              double temp=a[nj][ni-1];
              a[nj][ni-1]=cc*temp+ss*a[nj][ni];
              a[nj][ni]=-ss*temp+cc*a[nj][ni];
            }
            /* Apply procedure G2 (CC,SS,B(J-1),B(J)) */
            double temp=b[ni-1]; b[ni-1]=cc*temp+ss*b[ni]; b[ni]=-ss*temp+cc*b[ni];
          }
        }
        npp1=nsetp-1; nsetp--; iz1--; index[iz1]=k;
 
        /* See if the remaining coefficients in set P are feasible; they should be
           because of the way alpha was determined. If any are infeasible 
           it is due to round-off error. Any that are non-positive 
           will be set to zero and moved from set P to set Z. */
        for(jj=0, pfeas=1; jj<nsetp; jj++) {
          k=index[jj]; if(x[k]<=0.) {pfeas=0; break;}
        }
      } while(pfeas==0);
 
      /* Copy B[] into zz[], then solve again and loop back */
      for(int mi=0; mi<m; mi++) zz[mi]=b[mi];
      for(int mi=0; mi<nsetp; mi++) {
        int ip=nsetp-(mi+1);
        if(mi!=0) for(int ii=0; ii<=ip; ii++) zz[ii]-=a[jj][ii]*zz[ip+1];
        jj=index[ip]; zz[ip]/=a[jj][ip];
      }
    } /* end of secondary loop */
 
    if(iter>=itmax) break;
    for(int ip=0; ip<nsetp; ip++) {k=index[ip]; x[k]=zz[ip];}
  } /* end of main loop */
 
  if(wp != NULL) {
    if(npp1>=m) for(int ni=0; ni<n; ni++) w[ni]=0.;
  }
 
  /* Compute the norm of the final residual vector (sum-of-squares) */
  if(rnorm != NULL) {
    double ss=0.0;
    if(npp1<m) for(int mi=npp1; mi<m; mi++) ss+=(b[mi]*b[mi]);
    *rnorm=ss;
  }
 
  /* Free working space, if it was allocated here */
//w=NULL; zz=NULL; index=NULL;
  if(wl!=NULL) free(wl); 
  if(zzl!=NULL) free(zzl); 
  if(indexl!=NULL) free(indexl);
  if(iter>=itmax) return(1);
 
  return(0);
} /* nnls */

Referenced by fitExpDecayNNLS(), and img_k1_using_ki().

◆ nnlsWght()

int nnlsWght	(	int	N,
		int	M,
		double **	A,
		double *	b,
		double *	weight )

Algorithm for weighting the problem that is given to nnls-algorithm.

Square roots of weights are used because in nnls the difference w*A-w*b is squared.

Returns: Algorithm returns zero if successful, 1 if arguments are inappropriate.

See also: nnlsWghtSquared, nnls, qrLSQ, llsqWght

Parameters

N	NNLS dimension N (nr of parameters).
M	NNLS dimension M (nr of samples).
A	NNLS matrix A.
b	NNLS vector B.
weight	Weights for each sample (array of length M).

Definition at line 257 of file nnls.c.

  {
  int n, m;
  double *w;
 
  /* Check the arguments */
  if(N<1 || M<1 || A==NULL || b==NULL || weight==NULL) return(1);
 
  /* Allocate memory */
  w=(double*)malloc(M*sizeof(double)); if(w==NULL) return(2);
 
  /* Check that weights are not zero and get the square roots of them to w[] */
  for(m=0; m<M; m++) {
    if(weight[m]<=1.0e-20) w[m]=0.0;
    else w[m]=sqrt(weight[m]);
  }
 
  /* Multiply rows of matrix A and elements of vector b with weights*/
  for(m=0; m<M; m++) {
    for(n=0; n<N; n++) {
      A[n][m]*=w[m];
    }
    b[m]*=w[m];
  }
 
  free(w);
  return(0);
}

◆ nnlsWghtSquared()

int nnlsWghtSquared	(	int	N,
		int	M,
		double **	A,
		double *	b,
		double *	sweight )

Algorithm for weighting the problem that is given to nnls-algorithm.

Square roots of weights are used because in nnls the difference w*A-w*b is squared. Here user must give squared weights; this makes calculation faster, when this function needs to be called many times.

Returns: Algorithm returns zero if successful, 1 if arguments are inappropriate.

See also: nnlsWghtSquared, qrLSQ, nnls

Parameters

N	NNLS dimension N (nr of parameters).
M	NNLS dimension M (nr of samples).
A	NNLS matrix A.
b	NNLS vector B.
sweight	Squared weights for each sample (array of length M).

Definition at line 306 of file nnls.c.

  {
  int n, m;
 
  /* Check the arguments */
  if(N<1 || M<1 || A==NULL || b==NULL || sweight==NULL) return(1);
 
  /* Multiply rows of matrix A and elements of vector b with weights*/
  for(m=0; m<M; m++) {
    for(n=0; n<N; n++) {
      A[n][m]*=sweight[m];
    }
    b[m]*=sweight[m];
  }
 
  return(0);
}

Functions

Detailed Description

Function Documentation

◆ nnls()

◆ nnlsWght()

◆ nnlsWghtSquared()