14.3 QR Decomposition with Column Pivoting

The QR decomposition of an M-by-N matrix A can be extended to the rank deficient case by introducing a column permutation P,

A P = Q R

The first r columns of Q form an orthonormal basis for the range of A for a matrix with column rank r. This decomposition can also be used to convert the linear system A x = b into the triangular system R y = Q^T b, x = P y, which can be solved by back-substitution and permutation. We denote the QR decomposition with column pivoting by QRP^T since A = Q R P^T. When A is rank deficient with r = {\rm rank}(A), the matrix R can be partitioned as

R = [ R11 R12; 0 R22 ] =~ [ R11 R12; 0 0 ]

where R_{11} is r-by-r and nonsingular. In this case, a “basic” least squares solution for the overdetermined system A x = b can be obtained as

x = P [ R11^-1 c1 ; 0 ]

where c_1 consists of the first r elements of Q^T b. This basic solution is not guaranteed to be the minimum norm solution unless R_{12} = 0 (see Complete Orthogonal Decomposition).

Function: int gsl_linalg_QRPT_decomp (gsl_matrix * A, gsl_vector * tau, gsl_permutation * p, int * signum, gsl_vector * norm)

This function factorizes the M-by-N matrix A into the QRP^T decomposition A = Q R P^T. On output the diagonal and upper triangular part of the input matrix contain the matrix R. The permutation matrix P is stored in the permutation p. The sign of the permutation is given by signum. It has the value (-1)^n, where n is the number of interchanges in the permutation. The vector tau and the columns of the lower triangular part of the matrix A contain the Householder coefficients and vectors which encode the orthogonal matrix Q. The vector tau must be of length k=\min(M,N). The matrix Q is related to these components by, Q = Q_k ... Q_2 Q_1 where Q_i = I - \tau_i v_i v_i^T and v_i is the Householder vector v_i = (0,...,1,A(i+1,i),A(i+2,i),...,A(m,i)). This is the same storage scheme as used by LAPACK. The vector norm is a workspace of length N used for column pivoting.

The algorithm used to perform the decomposition is Householder QR with column pivoting (Golub & Van Loan, Matrix Computations, Algorithm 5.4.1).

Function: int gsl_linalg_QRPT_decomp2 (const gsl_matrix * A, gsl_matrix * q, gsl_matrix * r, gsl_vector * tau, gsl_permutation * p, int * signum, gsl_vector * norm)

This function factorizes the matrix A into the decomposition A = Q R P^T without modifying A itself and storing the output in the separate matrices q and r.

Function: int gsl_linalg_QRPT_solve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)

This function solves the square system A x = b using the QRP^T decomposition of A held in (QR, tau, p) which must have been computed previously by gsl_linalg_QRPT_decomp.

Function: int gsl_linalg_QRPT_svx (const gsl_matrix * QR, const gsl_vector * tau, const gsl_permutation * p, gsl_vector * x)

This function solves the square system A x = b in-place using the QRP^T decomposition of A held in (QR,tau,p). On input x should contain the right-hand side b, which is replaced by the solution on output.

Function: int gsl_linalg_QRPT_lssolve (const gsl_matrix * QR, const gsl_vector * tau, const gsl_permutation * p, 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 and is assumed to have full rank. The least squares solution minimizes the Euclidean norm of the residual, ||b - A x||. The routine requires as input the QR decomposition of A into (QR, tau, p) given by gsl_linalg_QRPT_decomp. The solution is returned in x. The residual is computed as a by-product and stored in residual. For rank deficient matrices, gsl_linalg_QRPT_lssolve2 should be used instead.

Function: int gsl_linalg_QRPT_lssolve2 (const gsl_matrix * QR, const gsl_vector * tau, const gsl_permutation * p, const gsl_vector * b, const size_t rank, 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 and has rank given by the input rank. If the user does not know the rank of A, the routine gsl_linalg_QRPT_rank can be called to estimate it. The least squares solution is the so-called “basic” solution discussed above and may not be the minimum norm solution. The routine requires as input the QR decomposition of A into (QR, tau, p) given by gsl_linalg_QRPT_decomp. The solution is returned in x. The residual is computed as a by-product and stored in residual.

Function: int gsl_linalg_QRPT_QRsolve (const gsl_matrix * Q, const gsl_matrix * R, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)

This function solves the square system R P^T x = Q^T b for x. It can be used when the QR decomposition of a matrix is available in unpacked form as (Q, R).

Function: int gsl_linalg_QRPT_update (gsl_matrix * Q, gsl_matrix * R, const gsl_permutation * p, gsl_vector * w, const gsl_vector * v)

This function performs a rank-1 update w v^T of the QRP^T decomposition (Q, R, p). The update is given by Q'R' = Q (R + w v^T P) where the output matrices Q' and R' are also orthogonal and right triangular. Note that w is destroyed by the update. The permutation p is not changed.

Function: int gsl_linalg_QRPT_Rsolve (const gsl_matrix * QR, const gsl_permutation * p, const gsl_vector * b, gsl_vector * x)

This function solves the triangular system R P^T x = b for the N-by-N matrix R contained in QR.

Function: int gsl_linalg_QRPT_Rsvx (const gsl_matrix * QR, const gsl_permutation * p, gsl_vector * x)

This function solves the triangular system R P^T x = b in-place for the N-by-N matrix R contained in QR. On input x should contain the right-hand side b, which is replaced by the solution on output.

Function: size_t gsl_linalg_QRPT_rank (const gsl_matrix * QR, const double tol)

This function estimates the rank of the triangular matrix R contained in QR. The algorithm simply counts the number of diagonal elements of R whose absolute value is greater than the specified tolerance tol. If the input tol is negative, a default value of 20 (M + N) eps(max(|diag(R)|)) is used.

Function: int gsl_linalg_QRPT_rcond (const gsl_matrix * QR, double * rcond, gsl_vector * work)

This function estimates the reciprocal condition number (using the 1-norm) of the R factor, stored in the upper triangle of QR. The reciprocal condition number estimate, defined as 1 / (||R||_1 \cdot ||R^{-1}||_1), is stored in rcond. Additional workspace of size 3 N is required in work.