Next: Complete Orthogonal Decomposition, Previous: QR Decomposition, Up: Linear Algebra [Index]

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`.

Next: Complete Orthogonal Decomposition, Previous: QR Decomposition, Up: Linear Algebra [Index]