Chapter Contents
Chapter Contents
Previous
Previous
Next
Next
Language Reference

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 r\leq n 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, {{\Pi}}, 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 m \geq n, 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,
A= Q[ R\ 0
 ] {{\Pi}}^' P^' {{\Pi}}

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.
  1. Use the QR decomposition of A with column pivoting,
    A= Q[ R\ 0
 ] {{\Pi}}^' 
 = [ Y& Z
 ]
 [ R_1 & R_2 \ 0 & 0
 ] {{\Pi}}^'
    where R= [ R_1 & R_2 
 ] \in{\cal R}^{r x t} is upper trapezoidal, R_1\in{\cal R}^{r x r} is upper triangular and invertible, R_2\in{\cal R}^{r x s}, Q= [Y& Z] is orthogonal, Y\in{\cal R}^{t x r}, Z\in{\cal R}^{t x s}, and {{\Pi}} permutes the columns of A.
  2. Use the transpose L12 of the upper trapezoidal matrix R= [ R_1 & R_2 
 ],
    L_{12} = [ L_1 \ L_2
 ]
 = R^' \in{\cal R}^{t x r}
    with rank(L12) = rank(L1) = r, L_1\in{\cal R}^{r x r} lower triangular, L_2\in{\cal R}^{s x r}. The lower trapezoidal matrix L_{12}\in{\cal R}^{t x r} is premultiplied with r Householder transformations P1, ... ,Pr:
    P_r  ...  P_1 [ L_1 \ L_2
 ]
 = [ L\ 0
 ]
    each zeroing out one of the r columns of L2 and producing the nonsingular lower triangular matrix L\in{\cal R}^{r x r}. Therefore, you obtain
    A= Q[ L^' & 0 \ 0 & 0
 ]{{\Pi}}^' P^' 
 = Y[ L^' & 0
 ] {{\Pi}}^' P^'
    with P= {{\Pi}}P_r ... {P}_1\in{\cal R}^{t x t} and upper triangular L'. This second step is described in Golub and Van Loan (1989, p. 220 and p. 236).

  3. Compute the Moore-Penrose Inverse A- explicitly.
    A^- = P{{\Pi}}[ (L^')^{-1} & 0 \ 0 & 0
 ] Q^' 
 = P{{\Pi}}[ (L^')^{-1} \ 0
 ] Y^'
    (a)
    Obtain Y in Q= [Y& Z
 ] explicitly by applying the r Householder transformations obtained in the first step to [ I_r \ 0 
 ].

    (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 [ (L^')^{-1}Y^' \ 
 0 \ ], which results in the symmetric matrix A^-\in{\cal R}^{t x t}.

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.

Chapter Contents
Chapter Contents
Previous
Previous
Next
Next
Top
Top

Copyright © 1999 by SAS Institute Inc., Cary, NC, USA. All rights reserved.