qrlsq.c File Reference

QR decomposition for solving least squares problems. 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	qrLSQ (double **mat, double *rhs, double *sol, const unsigned int rows, const unsigned int cols, double *r2)
	QR least-squares solving routine.
int	qr (double **A, int m, int n, double *B, double *X, double *rnorm, double *tau, double *res, double **wws, double *ws)
int	qr_decomp (double **a, int M, int N, double *tau, double **cchain, double *chain)
int	qr_solve (double **QR, int M, int N, double *tau, double *b, double *x, double *residual, double *resNorm, double **cchain, double *chain)
int	qrWeight (int N, int M, double **A, double *b, double *weight, double *ws)
int	qrWeightRm (int N, int M, double **A, double *b, double *weight, double *ws)
int	qrSimpleLSQ (double **A, double *B, int M, int N, double *X, double *r2)
int	qrLH (const unsigned int m, const unsigned int n, double *a, double *b, double *x, double *r2)
	Solve over-determined least-squares problem A x ~ b using successive Householder rotations.

Detailed Description

QR decomposition for solving least squares problems.

Author

Kaisa Liukko, Vesa Oikonen

Definition in file qrlsq.c.

Function Documentation

qr()

int qr	(	double **	A,
		int	m,
		int	n,
		double *	B,
		double *	X,
		double *	rnorm,
		double *	tau,
		double *	res,
		double **	wws,
		double *	ws )

Algorithm QR.

Solves a matrix form least square problem min||A x - b|| => A x =~ b (A is m*n matrix, m>=n) using the QR decomposition for overdetermined systems. Based on GNU Scientific Library, edited by Kaisa Liukko.

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

See also

qrLSQ, nnls, qr_decomp

Returns

Function returns 0 if successful and 1 in case of invalid problem dimensions or memory allocation error.

Parameters

A	On entry, a[m][n] contains the m by n matrix A. On exit, a[][] contains the QR factorization.
m	Dimensions of matrix A are a[m][n].
n	Dimensions of matrix A are a[m][n].
B	B[] is an m-size vector containing the right-hand side vector b.
X	On exit, x[] will contain the solution vector x (size of n).
rnorm	On exit, rnorm (pointer to double) contains the squared Euclidean norm of the residual vector (R^2); enter NULL if not needed.
tau	On exit, tau[] will contain the householder coefficients (size of n); enter NULL, if not needed.
res	An m-size array of working space, res[]. On output contains residual b - Ax. Enter NULL to let qr() to handle it.
wws	m*n array of working space. Enter NULL to let qr() to handle it.
ws	2m-array of working space. Enter NULL to let qr() to handle it.

Definition at line 411 of file qrlsq.c.

  {
  int i;
  double *qrRes, *qrTau, **qrWws, *qrWs, *chain;
 
  /* Check the parameters and data */
  if(m<=0 || n<=0 || A==NULL || B==NULL || X==NULL) return(1);
  if(m<n) return(1);
 
  /* Allocate memory for working space, if required */
  if(tau!=NULL) qrTau=tau; else qrTau=(double*)calloc(n, sizeof(double));
  if(res!=NULL) qrRes=res; else qrRes=(double*)calloc(m, sizeof(double));
  if(wws!=NULL) {
    qrWws=wws; chain=(double*)NULL;
  } else {
    qrWws=(double**)malloc(m * sizeof(double*));
    chain=(double*)malloc(m*n * sizeof(double));
    for(i=0; i<m; i++) qrWws[i]=chain + i*n;
  }
  if(ws!=NULL) qrWs=ws; else qrWs=(double*)calloc(2*m, sizeof(double));
  if(qrTau==NULL || qrRes==NULL || qrWws==NULL || qrWs==NULL) return(1);
 
  /* Form the householder decomposition and solve the least square problem */
  if(qr_decomp(A, m, n, qrTau, qrWws, qrWs)) return(2);
  if(qr_solve(A, m, n, qrTau, B, X, qrRes, rnorm, qrWws, qrWs)) return(3);
 
  /* Free working space, if it was allocated here */
  //if(tau==NULL || res==NULL || wws==NULL || ws==NULL) printf("free in qr()\n");
  if(tau==NULL) free(qrTau);
  if(res==NULL) free(qrRes);
  if(wws==NULL) {free(qrWws); free(chain);}
  if(ws==NULL) free(qrWs);
  for(i=0; i<n; i++) if(isnan(X[i])) return(4);
  return(0);
} /* qr */

qr_decomp()

int qr_decomp	(	double **	a,
		int	M,
		int	N,
		double *	tau,
		double **	cchain,
		double *	chain )

Factorise a general M x N matrix A into A = Q R , where Q is orthogonal (M x M) and R is upper triangular (M x N).

Q is stored as a packed set of Householder vectors in the strict lower triangular part of the input matrix A and a set of coefficients in vector tau.

R is stored in the diagonal and upper triangle of the input matrix.

The full matrix for Q can be obtained as the product Q = Q_1 Q_2 .. Q_k and it's transform as the product Q^T = Q_k .. Q_2 Q_1 , where k = min(M,N) and Q_i = (I - tau_i * h_i * h_i^T) and where h_i is a Householder vector h_i = [1, A(i+1,i), A(i+2,i), ... , A(M,i)]. This storage scheme is the same as in LAPACK.

NOTICE! The calling program must take care that pointer tau is of size N or M, whichever is smaller.

See also

qr

Returns

Function returns 0 if ok.

Parameters

a	contains coefficient matrix A (m*n) as input and factorisation QR as output.
M	nr of rows in matrix A.
N	nr of columns in matrix A.
tau	Vector for householder coefficients, of length N or M, whichever is smaller.
cchain	m*n array of working space.
chain	m size array of working space.

Definition at line 493 of file qrlsq.c.

  {
  //printf("qr_decomp()\n");
  int i, m, n, MNmin;
  double *subvector, **submatrix;
 
  /* Local variables */
  if(M<N) MNmin=M; else MNmin=N;
  if(MNmin<1 || a==NULL || tau==NULL || cchain==NULL || chain==NULL) return(1);
  subvector=chain;
  submatrix=cchain;
 
  for(i=0; i<MNmin; i++) {
    //printf("i=%d (MNmin=%d)\n", i, MNmin);
    /* Compute the Householder transformation to reduce the j-th column of the matrix A to 
       a multiple of the j-th unit vector. Householder vector h_i is saved in the lower triangular 
       part of the column and Householder coefficient tau_i in the vector tau. */
    for(m=i; m<M; m++) {
      //printf("subvector[%d]=a[%d][%d]\n", m-1, m, i);
      subvector[m-i]=a[m][i];
    }
    tau[i] = householder_transform(subvector, M-i);
    for(m=i; m<M; m++) a[m][i]=subvector[m-i];
 
    /* Apply the transformation to the remaining columns to get upper triangular part of matrix R  */
    if(i+1 < N) {
      printf("     apply transformation\n");
      for(m=i; m<M; m++) for(n=i+1; n<N; n++) {
        if((m-i)<0 || (m-i)>=M || (n-i-1)<0 || (n-i-1)>=N) printf("OVERFLOW!\n"); 
        submatrix[m-i][n-i-1]=a[m][n];
      }
      if(householder_hm(tau[i], subvector, submatrix, M-i, N-i)) return(2);
      for(m=i; m<M; m++) for(n=i+1; n<N; n++)
        a[m][n]=submatrix[m-i][n-i-1];
    }
  }
 
  return(0);
}

Referenced by qr().

qr_solve()

int qr_solve	(	double **	QR,
		int	M,
		int	N,
		double *	tau,
		double *	b,
		double *	x,
		double *	residual,
		double *	resNorm,
		double **	cchain,
		double *	chain )

Find the least squares solution to the overdetermined system

A x = b

for m >= n using the QR factorisation A = Q R. qr_decomp() must be used prior to this function in order to form the QR factorisation of A. Solution is formed in the following order: QR x = b => R x = Q^T b => x = R^-1 (Q^T b)

NOTICE! The calling program must take care that pointers b, x and residual are of the right size.

See also

qr, qr_decomp

Returns

Function returns 0 if ok.

Parameters

QR	m*n matrix containing householder vectors of A.
M	nr of rows in matrix A.
N	nr of columns in matrix A.
tau	vector containing householder coefficients tau; length is M or N, whichever is smaller.
b	Contains m-size vector b of A x = b.
x	solution vector x of length n.
residual	residual vector of length m.
resNorm	norm^2 of the residual vector; enter NULL if not needed.
cchain	m*n array of the working space.
chain	2m length array of the working space.

Definition at line 563 of file qrlsq.c.

  {
  //printf("qr_solve()\n");
  if(QR==NULL || tau==NULL || b==NULL || x==NULL || residual==NULL) return(1);
  if(cchain==NULL || chain==NULL) return(1);
  if(M<1 || N<1) return(2);
 
  int m, n;
  double **Rmatrix;
 
  /* Local variable */
  Rmatrix=cchain;
 
  /* Get matrix R from the upper triangular part of QR
     First the rows N - M-1 are eliminated from R matrix*/
  for(m=0; m<N; m++) for(n=0; n<N; n++) Rmatrix[m][n]=QR[m][n];
  for(m=0; m<M; m++) residual[m]=b[m];
 
  /* Compute b = Q^T b */
  if(qr_QTvec(M, N, QR, tau, residual, chain)) return(3);
 
  /* Solve R x = b by computing x = R^-1 b */
  for(n=0; n<N; n++) x[n]=residual[n];
  if(qr_invRvec(N, Rmatrix, x)) return(4);
 
  /* Compute residual = b - A x = Q (Q^T b - R x) */
  for(n=0; n<N; n++) residual[n]=0.0;
  /* Compute residual= Q*residual */
  if(qr_Qvec(M, N, QR, tau, residual, chain)) return(4);
 
  /* Compute norm^2 of the residual vector, if needed */
  if(resNorm!=NULL)
    for(m=0, *resNorm=0.0; m<M; m++) *resNorm +=residual[m]*residual[m];
 
  return(0);
}

Referenced by qr().

qrLH()

int qrLH	(	const unsigned int	m,
		const unsigned int	n,
		double *	a,
		double *	b,
		double *	x,
		double *	r2 )

Solve over-determined least-squares problem A x ~ b using successive Householder rotations.

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, and Fortran code by R.L. Parker and P.B. Stark.

Returns

Returns 0 when successful, 1 if system is singular, and 2 in case of other errors, including that system is under-determined.

Parameters

m	Number of samples in matrix A and the length of vector b.
n	Number of parameters in matrix A and the length of vector x. The n must be smaller or equal to m.
a	Pointer to matrix A; matrix must be given as an n*m array, containing n consecutive m-length vectors. Contents of A are modified in this routine.
b	Pointer to vector b of length n. Contents of b are modified in this routine.
x	Pointer to the result vector x of length n.
r2	Pointer to a double value, in where the sum of squared residuals is written.

Definition at line 946 of file qrlsq.c.

  {
  /* Check the input */
  if(a==NULL || b==NULL || x==NULL || r2==NULL) return(2);
  if(n<1 || m<n) {*r2=nan(""); return(2);}
 
  /* Initiate output to zeroes, in case of exit because of singularity */
  for(unsigned int ni=0; ni<n; ni++) x[ni]=0.0;
  *r2=0.0;
 
  /* Rotates matrix A into upper triangular form */
  for(unsigned int ni=0; ni<n; ni++) {
    /* Find constants for rotation and diagonal entry */
    double sq=0.0;
    for(unsigned int mi=ni; mi<m; mi++) sq+=a[mi + ni*m]*a[mi + ni*m];
    if(sq==0.0) return(1);
    double qv1=-copysign(sqrt(sq), a[ni + ni*m]);
    double u1=a[ni + ni*m] - qv1;
    a[ni + ni*m]=qv1;
    unsigned int ni1=ni+1;
    /*  Rotate the remaining columns of sub-matrix. */
    for(unsigned int nj=ni1; nj<n; nj++) {
      double dot=u1*a[ni + nj*m];
      for(unsigned int mi=ni1; mi<m; mi++)
        dot+=a[mi + nj*m] * a[mi + ni*m];
      double c=dot/fabs(qv1*u1);
      for(unsigned int mi=ni1; mi<m; mi++)
        a[mi + nj*m]-=c*a[mi + ni*m];
      a[ni + nj*m]-=c*u1;
    }
    /* Rotate vector B */
    double dot=u1*b[ni];
    for(unsigned int mi=ni1; mi<m; mi++)
      dot+=b[mi]*a[mi + ni*m];
    double c=dot/fabs(qv1*u1);
    b[ni]-=c*u1;
    for(unsigned int mi=ni1; mi<m; mi++)
      b[mi]-=c*a[mi + ni*m];
  } // end of rotation loop
 
  /* Solve triangular system by back-substitution. */
  for(unsigned int ni=0; ni<n; ni++) {
    int k=n-ni-1;
    double s=b[k];
    for(unsigned int nj=k+1; nj<n; nj++) 
      s-=a[k + nj*m] * x[nj];
    if(a[k + k*m]==0.0) return(1);
    x[k]=s/a[k + k*m];
  }
 
  /* Calculate the sum of squared residuals. */
  *r2=0.0;
  for(unsigned int mi=n; mi<m; mi++)
    *r2 += b[mi]*b[mi];
 
  return(0);
}

Referenced by bvls().

qrLSQ()

int qrLSQ	(	double **	mat,
		double *	rhs,
		double *	sol,
		const unsigned int	rows,
		const unsigned int	cols,
		double *	r2 )

QR least-squares solving routine.

Solves a matrix form least square problem min||A x - b|| => A x =~ b (A is m*n matrix, m>=n) using the QR decomposition for overdetermined systems

Author

: Michael Mazack License: Public Domain. Redistribution and modification without restriction is granted. If you find this code helpful, please let the author know (https://mazack.org). Small editions by Vesa Oikonen when added to tpcclib.

See also

nnls, qr

Todo

Check that modification in one subroutine is ok, and add possibility to provide working space for the main and sub-functions.

Returns

Returns 0 if ok.

Parameters

mat	Pointer to the row-major matrix (2D array), A[rows][cols]; not modified.
rhs	Vector b[rows]; modified to contain the computed (fitted) b[], that is, matrix-vector product A*x.
sol	Solution vector x[cols]; solution x corresponds to the solution of both the modified and original systems A and B.
rows	Number of rows (samples) in matrix A and length of vector B.
cols	Number of cols (parameters) in matrix A and length of vector X.
r2	Pointer to value where R^2, the difference between original and computed right hand side (b); enter NULL if not needed.

Definition at line 72 of file qrlsq.c.

  {
  if(rows<1 || cols<1) return(1);
  if(mat==NULL || rhs==NULL || sol==NULL) return(2);
 
  unsigned int i, j, max_loc;
  unsigned int *p=NULL;
  double *buf=NULL, *ptr, *v=NULL, *orig_b=NULL, **A=NULL;
 
  /* Allocate memory for index vector and Householder transform vector */
  /* and original data matrix and vector. */
  p=malloc(sizeof(unsigned int)*cols);
  A=malloc(sizeof(double*)*rows);
  buf=malloc(sizeof(double)*(rows*cols+rows+rows));
  if(p==NULL || A==NULL || buf==NULL) {
    free(p); free(A); free(buf); 
    return(3);
  }
  ptr=buf; 
  for(i=0; i<rows; i++) {A[i]=ptr; ptr+=cols;}
  v=ptr; ptr+=rows; orig_b=ptr;
  /* copy the original data */
  for(i=0; i<rows; i++) {
    orig_b[i]=rhs[i];
    for(j=0; j<cols; j++) A[i][j]=mat[i][j];
  }
 
  /* Initial permutation vector. */
  for(i=0; i<cols; i++) p[i]=i;
  
  /* Apply rotators to make R and Q'*b */
  for(i=0; i<cols; i++) {
    if(qr_get_next_col(A, rows, cols, i, p, &max_loc)==0)
      qr_swap_cols(p, i, max_loc);
    qr_householder(A, rows, cols, i, p[i], v);
    qr_apply_householder(A, rhs, rows, cols, v, i, p);
  }
 
  /* Back solve Rx = Q'*b */
  qr_back_solve(A, rhs, rows, cols, sol, p);
 
  /* Compute fitted b[] (matrix-vector product A*x) using non-modified A */
  for(i=0; i<rows; i++) {
    rhs[i]=0.0;
    for(j=0; j<cols; j++) rhs[i]+=mat[i][j]*sol[j];
  }
 
  /* Compute R^2, if requested */
  if(r2!=NULL) {
    double ss=0, d;
    for(i=0; i<rows; i++) {
      d=orig_b[i]-rhs[i];
      ss+=d*d;
    }
    *r2=ss;
  }
 
  /* Collect garbage. */
  free(p); free(A); free(buf);
  return(0);
}

qrSimpleLSQ()

int qrSimpleLSQ	(	double **	A,
		double *	B,
		int	M,
		int	N,
		double *	X,
		double *	r2 )

Simplistic QR least-squares solving routine.

Reference: Máté A. Introduction to Numerical Analysis with C programs. Chapter 38. Overdetermined systems of linear equations. Brooklyn College of the City University of New York, 2014.

See also

qrLSQ

Returns

Function returns 0 when successful.

Parameters

A	Matrix A[M][N].
B	Vector B[M].
M	nr of rows in matrix A, must be >=N.
N	nr of columns in matrix A, must be <=M.
X	Vector X of length N for the results.
r2	Pointer to value where R^2, the difference between original and computed right hand side (b); enter NULL if not needed.

Definition at line 854 of file qrlsq.c.

  {
  //printf("qrSimpleLSQ()\n");
  if(A==NULL || B==NULL || X==NULL || N<1 || M<N) return(1);
 
  const double closeZero=1.0E-20;
 
  double cc[M], orig_B[M];
  double p[M][N], c[M][N];
 
  /* Householder transformation of matrix A into matrix P */
  for(int n=0; n<N; n++) {
    for(int nn=0; nn<n; nn++) p[nn][n]=0.0;
    p[n][n]=0.0;
    for(int m=n; m<M; m++) p[n][n]+=A[m][n]*A[m][n];
    p[n][n]=copysign(sqrt(p[n][n]), A[n][n]);
    p[n][n]+=A[n][n];
    if(fabs(p[n][n])<closeZero) return(2);
    for(int m=n+1; m<M; m++) p[m][n]=A[m][n];
    double norm=0.0;
    for(int m=n; m<M; m++) norm+=p[m][n]*p[m][n];
    norm=sqrt(norm);
    for(int m=n; m<M; m++) p[m][n]/=norm;
    for(int m=n; m<M; m++) {
      for(int nn=n; nn<N; nn++) {
        c[m][nn]=A[m][nn];
        for(int mm=n; mm<M; mm++) c[m][nn]-=2.0*p[m][n]*p[mm][n]*A[mm][nn];
      }
    }
    for(int m=n; m<M; m++)
      for(int nn=n; nn<N; nn++) A[m][nn]=c[m][nn];
  }
 
  /* Keep the original B[] for R^2 calculation */
  for(int m=0; m<M; m++) orig_B[m]=B[m];
 
  /* Compute P'b */
  for(int n=0; n<N; n++) {
    for(int m=n; m<M; m++) {
      cc[m]=B[m];
      for(int mm=n; mm<M; mm++) cc[m]-=2.0*p[m][n]*p[mm][n]*B[mm];
    }
    for(int m=n; m<M; m++) B[m]=cc[m];
  }
 
  /* Solve the linear system with backward substitution */
  X[N-1]=B[N-1]/A[N-1][N-1];
  for(int n=N-2; n>=0; n--) {
    X[n]=B[n];
    for(int nn=n+1; nn<N; nn++) X[n]-=A[n][nn]*X[nn];
    X[n]/=A[n][n];
  }
 
  /* Compute R^2, if requested */
  if(r2!=NULL) {
    double ss=0, d;
    for(int m=0; m<M; m++) {
      d=orig_B[m]-B[m];
      ss+=d*d;
    }
    *r2=ss;
  }
 
  return(0);
}

qrWeight()

int qrWeight	(	int	N,
		int	M,
		double **	A,
		double *	b,
		double *	weight,
		double *	ws )

Algorithm for weighting the problem that is given to QR algorithm. Square roots of weights are used because in QR the difference w*A-w*b is squared.

See also

qrWeightRm, qrLSQ, qr

Todo

Add tests.

Returns

Algorithm returns zero if successful, otherwise <>0.

Parameters

N	Dimensions of matrix A are a[m][n].
M	Dimensions of matrix A are a[m][n]; size of vector B is m.
A	Matrix a[m][n] for QR, contents will be weighted here.
b	B[] is an m-size vector for QR, contents will be weighted here.
weight	Pointer to array of size m, which contains sample weights, used here to weight matrix A and vector b.
ws	m-sized vector for working space; enter NULL to allocate locally.

Definition at line 745 of file qrlsq.c.

  {
  int n, m;
  double *w;
 
  /* Check the arguments */
  if(N<1 || M<1 || A==NULL || b==NULL || weight==NULL) return(1);
 
  /* Allocate memory, if necessary */
  if(ws==NULL) {
    w=(double*)malloc(M*sizeof(double)); if(w==NULL) return(2);
  } else {
    w=ws;
  }
 
  /* 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-100) w[m]=1.0e-50;
    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[m][n]*=w[m];
    b[m]*=w[m];
  }
 
  if(ws==NULL) free(w);
  return(0);
}

qrWeightRm()

int qrWeightRm	(	int	N,
		int	M,
		double **	A,
		double *	b,
		double *	weight,
		double *	ws )

Remove the weights from data provided to QR LSQ algorithms.

See also

qrWeight, qrLSQ, qr

Todo

Add tests.

Returns

Algorithm returns zero if successful, otherwise <>0.

Parameters

N	Dimensions of matrix A are a[m][n].
M	Dimensions of matrix A are a[m][n]; size of vector B is m.
A	Matrix a[m][n] for QR, weights will be removed here; enter NULL if not needed.
b	B[] is an m-size vector for QR, weights will be removed here; enter NULL if not needed here.
weight	Pointer to array of size m, which contains sample weights, used here to weight matrix A and vector b.
ws	m-sized vector for working space; enter NULL to allocate locally

Definition at line 796 of file qrlsq.c.

  {
  int n, m;
  double *w;
 
  /* Check the arguments */
  if(N<1 || M<1 || weight==NULL) return(1);
  if(A==NULL && b==NULL) return(0); // nothing to do
 
  /* Allocate memory, if necessary */
  if(ws==NULL) {
    w=(double*)malloc(M*sizeof(double)); if(w==NULL) return(2);
  } else {
    w=ws;
  }
 
  /* 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-100) w[m]=1.0e-50;
    else w[m]=sqrt(weight[m]);
  }
 
  /* Multiply rows of matrix A and elements of vector b with weights */
  for(m=0; m<M; m++) {
    if(A!=NULL) for(n=0; n<N; n++) A[m][n]/=w[m];
    if(b!=NULL) b[m]/=w[m];
  }
 
  if(ws==NULL) free(w);
  return(0);
}