Next: Singular Value Decomposition, Previous: QR Decomposition with Column Pivoting, Up: Linear Algebra [Index]

The complete orthogonal decomposition of a *M*-by-*N* matrix
*A* is a generalization of the QR decomposition with column pivoting, given by

A P = Q [ R11 0 ] Z [ 0 0 ]

where *P* is a *N*-by-*N* permutation matrix,
*Q* is *M*-by-*M* orthogonal, *R_{11}* is
*r*-by-*r* upper triangular, with *r = {\rm rank}(A)*,
and *Z* is *N*-by-*N* orthogonal. If *A*
has full rank, then *R_{11} = R*, *Z = I* and this reduces to the
QR decomposition with column pivoting. The advantage of using
the complete orthogonal decomposition for rank deficient matrices
is the ability to compute the minimum norm solution to the linear
least squares problem *Ax = b*, which is given by

x = P Z^T [ R11^-1 c1 ] [ 0 ]

and the vector *c_1* is the first *r* elements of *Q^T b*.

- Function:
*int***gsl_linalg_COD_decomp***(gsl_matrix **`A`, gsl_vector *`tau_Q`, gsl_vector *`tau_Z`, gsl_permutation *`p`, size_t *`rank`, gsl_vector *`work`) - Function:
*int***gsl_linalg_COD_decomp_e***(gsl_matrix **`A`, gsl_vector *`tau_Q`, gsl_vector *`tau_Z`, gsl_permutation *`p`, double`tol`, size_t *`rank`, gsl_vector *`work`) These functions factor the

*M*-by-*N*matrix`A`into the decomposition*A = Q R Z P^T*. The rank of`A`is computed as the number of diagonal elements of*R*greater than the tolerance`tol`and output in`rank`. If`tol`is not specified, a default value is used (see`gsl_linalg_QRPT_rank`

). On output, the permutation matrix*P*is stored in`p`. The matrix*R_{11}*is stored in the upper`rank`-by-`rank`block of`A`. The matrices*Q*and*Z*are encoded in packed storage in`A`on output. The vectors`tau_Q`and`tau_Z`contain the Householder scalars corresponding to the matrices*Q*and*Z*respectively and must be of length*k = \min(M,N)*. The vector`work`is additional workspace of length*N*.

- Function:
*int***gsl_linalg_COD_lssolve***(const gsl_matrix **`QRZ`, const gsl_vector *`tau_Q`, const gsl_vector *`tau_Z`, const gsl_permutation *`p`, const size_t`rank`, const gsl_vector *`b`, gsl_vector *`x`, gsl_vector *`residual`) This function finds the least squares solution to the overdetermined system

*A x = b*where the matrix`A`has more rows than columns. The least squares solution minimizes the Euclidean norm of the residual,*||b - A x||*. The routine requires as input the*QRZ*decomposition of*A*into (`QRZ`,`tau_Q`,`tau_Z`,`p`,`rank`) given by`gsl_linalg_COD_decomp`

. The solution is returned in`x`. The residual is computed as a by-product and stored in`residual`.

- Function:
*int***gsl_linalg_COD_unpack***(const gsl_matrix **`QRZ`, const gsl_vector *`tau_Q`, const gsl_vector *`tau_Z`, const size_t`rank`, gsl_matrix *`Q`, gsl_matrix *`R`, gsl_matrix *`Z`) This function unpacks the encoded

*QRZ*decomposition (`QRZ`,`tau_Q`,`tau_Z`,`rank`) into the matrices`Q`,`R`, and`Z`, where`Q`is*M*-by-*M*,`R`is*M*-by-*N*, and`Z`is*N*-by-*N*.

- Function:
*int***gsl_linalg_COD_matZ***(const gsl_matrix **`QRZ`, const gsl_vector *`tau_Z`, const size_t`rank`, gsl_matrix *`A`, gsl_vector *`work`) This function multiplies the input matrix

`A`on the right by`Z`,*A' = A Z*using the encoded*QRZ*decomposition (`QRZ`,`tau_Z`,`rank`).`A`must have*N*columns but may have any number of rows. Additional workspace of length*M*is provided in`work`.