Language Reference: COMPORT Call :: SAS/IML(R) 9.2 User's Guide

Language Reference

COMPORT Call

provides complete orthogonal decomposition by Householder transformations

CALL COMPORT( , , , piv, lindep, $a\lt$ , b<, sing>);

The COMPORT subroutine returns the following values:

: If is not specified, is the orthogonal matrix that is the product of the $\min(m,n)$ separate Householder transformations. If is specified, is the matrix ${q}^' {b}$ that has the transposed Householder transformations ${q}^'$ applied on the columns of the argument matrix .
: is the upper triangular matrix that contains the nonsingular upper triangular matrix ${l}^'$ of the complete orthogonal decomposition, where $r\leq n$ is the rank of . The full upper triangular matrix of the orthogonal decomposition of matrix can be obtained by vertical concatenation (IML operator //) of the zero matrix to the result .
p: is an matrix that is the product ${p}{{{{\pi}}}}$ of a permutation matrix ${{{{\pi}}}}$ with an orthogonal matrix . The permutation matrix is determined by the vector piv.
piv: is an vector of permutations of the columns of . That is, the QR decomposition is computed, not of , but of the matrix with columns $[{a}_{piv[1]} ... {a}_{piv[n]}]$ . The vector piv corresponds to an 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 detected by applying the Householder transformation in the order specified by the argument piv. That is, lindep.

The inputs to the COMPORT subroutine are as follows:

specifies the

matrix

, with $m \geq n$ , which is to be decomposed into the product of the

orthogonal matrix

, the

upper triangular matrix

, and the

orthogonal matrix

${a}= {q}[ {{r}}\ 0 ] {{{{\pi}}}}^' {p}^' {{{{\pi}}}}$

specifies an optional

matrix

that is to be left multiplied by the transposed

matrix ${q}^'$ .

sing

is an optional scalar specifying a singularity criterion.

The complete orthogonal decomposition of the singular matrix

can be used to compute the Moore-Penrose inverse ${a}^-$ , $r = {\rm rank}({a}) \lt n$ , or to compute the minimum 2-norm solution of the (rank deficient) least squares problem $\vert {a}{x}- {b}\vert^2_2$ .

Use the QR decomposition of 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}$ , is orthogonal, ${y}\in{\cal r}^{t x r}$ , ${z}\in{\cal r}^{t x s}$ , and ${{{{\pi}}}}$ permutes the columns of .
Use the transpose ${l}_{12}$ of the upper trapezoidal matrix ${{r}}= [ {{r}}_1 & {{r}}_2 ]$ ,
${l}_{12} = [ {l}_1 \ {l}_2 ] = {{r}}^' \in{\cal r}^{t x r}$
with ${\rm rank}({l}_{12}) = {\rm rank}({l}_1) = 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 Householder transformations ${p}_1, ... ,{p}_r$ :
${p}_r ... {p}_1 [ {l}_1 \ {l}_2 ] = [ {l}\ 0 ]$
each zeroing out one of the columns of ${l}_2$ 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).
Compute the Moore-Penrose inverse ${a}^-$ explicitly:
${a}^- = {p}{{{{\pi}}}}[ ({l}^')^{-1} & 0 \ 0 & 0 ] {q}^' = {p}{{{{\pi}}}}[ ({l}^')^{-1} \ 0 ] {y}^'$

(a)
Obtain in explicitly by applying the Householder transformations obtained in the first step to $[ {i}_r \ 0 ]$ .

(b)
Solve the lower triangular system $({l}^')^{-1}{y}^'$ with right-hand sides by using backward substitution, which yields an intermediate matrix.

(c)
Left-apply the Householder transformations in on the 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}^-$ by using the singular value decomposition of . Using complete orthogonal decomposition to compute ${a}^-$ usually requires far fewer floating-point operations. However, it can be slightly more sensitive to rounding errors, which can disturb the detection of the true rank of , than singular value decomposition.

The following example demonstrates the use of the APPCORT subroutine, as well as the resulting output:

  
    /* Only four linearly independent columns */ 
    A = {1 0 1 0 0, 
         1 0 0 1 0, 
         1 0 0 0 1, 
         0 1 1 0 0, 
         0 1 0 1 0, 
         0 1 0 0 1 }; 
    m = nrow(A); 
    n = ncol(A); 
  
    call comport(q,r,p,piv,lindep,A); 
    fullR = r // j(m-n, n, 0); 
    perm = i(n); 
    perm[piv,] = i(n); 
  
    /* recover A from factorization */ 
    A2 = q*fullR*p`*perm`; 
    reset fuzz; 
    print A2; 
  
    /* compute Moore-Penrose generalized inverse */ 
    rankA = n - lindep; 
    subR = R[1:rankA, 1:rankA]; 
    fullRinv = j(n, n, 0); 
    fullRinv[1:rankA, 1:rankA] = inv(subR); 
    Ainv = perm*p*fullRinv*q[,1:n]`; 
    print Ainv; 
  
    /* verify generalized inverse */ 
    eps = 1e-12; 
    if any(A*Ainv*A-A > eps) | 
       any(Ainv*A*Ainv-Ainv > eps) | 
       any((A*Ainv)`-A*Ainv > eps) | 
       any((Ainv*A)`-Ainv*A > eps) then 
           print "Pseudoinverse conditions not satisfied"; 
       else 
          print "Pseudoinverse conditions satisfied"; 
  
                              A2 
  
              1         0         1         0         0 
              1         0         0         1         0 
              1         0         0         0         1 
              0         1         1         0         0 
              0         1         0         1         0 
              0         1         0         0         1 
  
                                  AINV 
  
      0.2666667 0.2666667 0.2666667 -0.066667 -0.066667 -0.066667 
      -0.066667 -0.066667 -0.066667 0.2666667 0.2666667 0.2666667 
            0.4      -0.1      -0.1       0.4      -0.1      -0.1 
           -0.1       0.4      -0.1      -0.1       0.4      -0.1 
           -0.1      -0.1       0.4      -0.1      -0.1       0.4 
  
                       Pseudoinverse conditions satisfied

Top of Page