14.4 Complete Orthogonal Decomposition

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.