qr.c File Reference

Routines needed in the use of QR decomposition when solving least squares problems. More...

#include "libtpcmodel.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	qr_weight (int N, int M, double **A, double *b, double *weight, double *ws)
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

Routines needed in the use of QR decomposition when solving least squares problems.

Author

Kaisa Sederholm, Vesa Oikonen

These routines are based on the code of Gerard Jungman and Brian Gough provided in the GSL library (https://sources.redhat.com/gsl/).

Definition in file qr.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.

Returns

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

See also

qr_decomp, qr_solve, qrLH

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 346 of file qr.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) 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.

Returns

Function returns 0 if ok.

See also

qr_solve

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 matrix of working space.
chain	m size array of working space.

Definition at line 427 of file qr.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++) 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) {
      for(m=i; m<M; m++)
        for(n=i+1; n<N; 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.

Returns

Function returns 0 if ok.

Precondition

qr_decomp

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 matrix of the working space.
chain	2m length array of the working space.

Definition at line 492 of file qr.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 MNmin; if(M<N) MNmin=M; else MNmin=N;
 
  int m, n;
  double **Rmatrix=cchain;
  double *h=chain;
  double *w=chain+M;
 
  /* 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 */
  /* Form the product Q^T residual from householder vectors saved in the lower triangle of 
     QR matrix and householder coefficients saved in vector tau. */
  for(int i=0; i<MNmin; i++) {
    for(m=i; m<M; m++) h[m-i]=QR[m][i]; 
    for(m=i; m<M; m++) w[m-i]=residual[m];
    if(householder_hv(tau[i], M-i, h, w)) return(2);
    for(m=i; m<M; m++) residual[m]=w[m-i];
  }
 
  /* Solve R x = b by computing x = R^-1 b */
  for(n=0; n<N; n++) x[n]=residual[n];
  /* back-substitution */
  x[N-1]=x[N-1]/Rmatrix[N-1][N-1];
  for(int i=N-2; i>=0; i--) {
    for(int j=i+1; j<N; j++) x[i]-=Rmatrix[i][j]*x[j];
    x[i]/=Rmatrix[i][i];
  }
 
  /* 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 */
  /* Form the product Q*residual from householder vectors saved in the lower triangle
     of QR matrix and householder coefficients saved in vector tau. */
  for(int i=MNmin-1; i>=0; i--) {
    for(int m=i; m<M; m++) h[m-i]=QR[m][i]; 
    for(int m=i; m<M; m++) w[m-i]=residual[m];
    if(householder_hv(tau[i], M-i, h, w)) return(2);
    for(int m=i; m<M; m++) residual[m]=w[m-i];
  }
 
  /* 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().

qr_weight()

int qr_weight	(	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.

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 576 of file qr.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);
}

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 633 of file qr.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 54 of file qr.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);
}