COMPORT Call
provides complete orthogonal decomposition
by Householder transformations
- CALL COMPORT( q, r, p, piv, lindep, a<, b><,
sing>);
The COMPORT subroutine returns the following values:
- q
- If b is not specified, q is the m ×m
orthogonal matrix Q that is the product of the
min(m,n) separate Householder transformations.
If b is specified, q is the m ×p
matrix Q' B that has the transposed
Householder transformations Q' applied
on the p columns of the argument matrix B.
- r
- is the n ×n upper triangular matrix R that
contains the r ×r nonsingular upper triangular
matrix L' of the complete orthogonal
decomposition, where is the rank of A.
The full m ×n upper triangular matrix R of the
orthogonal decomposition of matrix A can be obtained
by vertical concatenation (IML operator //) of
the (m-n) ×n zero matrix to the result r.
- piv
- is an n ×1 vector of permutations of the columns of A.
That is, the QR decomposition is computed, not of A, but of
the matrix with columns [Apiv[1] ... Apiv[n]].
The vector piv corresponds to an n ×n
permutation matrix, , of the pivoted QR
decomposition in the first step of orthogonal decomposition.
- lindep
- specifies the number of linearly dependent columns in the
matrix A detected by applying the r Householder
transformation in the order specified by the argument piv.
That is, lindep=n-r.
The inputs to the COMPORT subroutine are as follows:
- a
- specifies the m ×n matrix A, with ,
which is to be decomposed into the product of the m ×m
orthogonal matrix Q, the n ×n upper triangular
matrix R, and the n ×n orthogonal matrix P,
- b
- specifies an optional m ×p matrix
B that is to be left multiplied by the
transposed m ×m matrix Q'.
- sing
- is an optional scalar specifying a singularity criterion.
The complete orthogonal decomposition of the
singular matrix A can be used to compute the
Moore-Penrose inverse A-, r = rank(A) < n,
or to compute the minimum 2-norm solution of the (rank
deficient) least-squares problem |Ax-b|22.
- Use the QR decomposition of A with column pivoting,
where
is upper trapezoidal,
is upper triangular and invertible,
,
is orthogonal,
,
,
and permutes the columns of A.
- Use the transpose L12 of the upper trapezoidal
matrix ,
with rank(L12) = rank(L1) = r,
lower triangular,
. The lower trapezoidal
matrix is premultiplied
with r Householder transformations P1, ... ,Pr:
each zeroing out one of the r columns of
L2 and producing the nonsingular lower
triangular matrix .
Therefore, you obtain
with and upper triangular L'.
This second step is described in Golub
and Van Loan (1989, p. 220 and p. 236).
- Compute the Moore-Penrose Inverse A- explicitly.
- (a)
- Obtain Y in
explicitly by applying the r Householder
transformations obtained in the first step to
.
- (b)
- Solve the r ×r lower triangular system
(L')-1Y' with t right
hand sides using backward substitution, which
yields an r ×t intermediate matrix.
- (c)
- Left-apply the r Householder transformations in
P on the r ×t intermediate matrix
, which results in the
symmetric matrix .
The GINV function computes the Moore-Penrose inverse
A- using the singular value decomposition of A.
Using complete orthogonal decomposition to compute
A- usually needs far fewer floating point operations.
However, it may be slightly more sensitive to rounding
errors, which can disturb the detection of the true
rank of A, than singular value decomposition.
Copyright © 1999 by SAS Institute Inc., Cary, NC, USA. All rights reserved.