nnls.c File Reference

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

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "tpcextensions.h"
#include "tpclinopt.h"

Go to the source code of this file.

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 Liukko

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 43 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;
  double *w=NULL, *zz=NULL;
  if(wp!=NULL) w=wp; else w=(double*)calloc(n, sizeof(double));
  if(zzp!=NULL) zz=zzp; else zz=(double*)calloc(m, sizeof(double));
  if(indexp!=NULL) index=indexp; else index=(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;
      nnls_lss_h12(1, npp1, npp1+1, m, a[j], &up, NULL);
      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];
        nnls_lss_h12(2, npp1, npp1+1, m, a[j], &up, zz);
        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];
        nnls_lss_h12(2, nsetp-1, npp1, m, &a[j][0], &up, a[jj]);
      }
    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;
            nnls_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 */
  if(wp==NULL) free(w); 
  if(zzp==NULL) free(zz); 
  if(indexp==NULL) free(index);
  if(iter>=itmax) return(1);
  return(0);
} /* nnls */

Referenced by spectralDExp(), spectralDMSurge(), and tacDelayCMFit().

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

Referenced by spectralDExp(), and spectralDMSurge().

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