COMPORT Call :: SAS/IML(R) 13.1 User's Guide

COMPORT Call

CALL COMPORT (q, r, p, piv, lindep, a <, b> <, sing> ) ;

The COMPORT subroutine provides the complete orthogonal decomposition by Householder transformations of a matrix $\bA$ .

The subroutine returns the following values:

q: is a matrix. If $b$ is not specified, $q$ is the $m \times m$ orthogonal matrix $\bQ$ that is the product of the $\min (m,n)$ separate Householder transformations. If $b$ is specified, $q$ is the $m \times p$ matrix ${\bQ }^{\prime } \bB$ that has the transposed Householder transformations ${\bQ }^{\prime }$ applied to the $p$ columns of the argument matrix $\bB$ .
r: is the $n \times n$ upper triangular matrix $\bR$ that contains the $r \times r$ nonsingular upper triangular matrix ${\bL }^{\prime }$ of the complete orthogonal decomposition, where $r\leq n$ is the rank of $\bA$ . The full $m \times n$ upper triangular matrix $\bR$ of the orthogonal decomposition of matrix $\bA$ can be obtained by vertical concatenation of the $(m-n) \times n$ zero matrix to the result r.
p: is an $n \times n$ matrix that is the product $\bP \bPi$ of a permutation matrix $\bPi$ with an orthogonal matrix $\bP$ . The permutation matrix is determined by the vector piv.
piv: is an $n \times 1$ vector of permutations of the columns of $\bA$ . That is, the QR decomposition is computed, not of $\bA$ , but of the matrix with columns $[\bA _{piv[1]}\ldots \bA _{piv[n]}]$ . The vector piv corresponds to an $n \times n$ permutation matrix, $\bPi$ , of the pivoted QR decomposition in the first step of orthogonal decomposition.
lindep: specifies the number of linearly dependent columns in the matrix $\bA$ detected by applying the $r$ Householder transformation in the order specified by the argument piv. That is, lindep is $n-r$ .

The input arguments to the COMPORT subroutine are as follows:

a: specifies the $m \times n$ matrix $\bA$ , with $m\geq n$ , which is to be decomposed into the product of the $m \times m$ orthogonal matrix $\bQ$ , the $n \times n$ upper triangular matrix $\bR$ , and the $n \times n$ orthogonal matrix $\bP$ ,

$\bA = \bQ \left[ \begin{array}{c} \bR \\ \mb {0} \end{array} \right] \bPi ^{\prime } \bP ^{\prime } \bPi$
b: specifies an optional $m \times p$ matrix $\bB$ that is to be left-multiplied by the transposed $m \times m$ matrix $\bQ ^{\prime }$ .
sing: is an optional scalar that specifies a singularity criterion.

The complete orthogonal decomposition of the singular matrix $\bA$ can be used to compute the Moore-Penrose inverse ${\bA }^-$ , $r = \mr {rank}(\bA ) < n$ , or to compute the minimum Euclidean-norm solution of the rank-deficient least squares problem $\| \bA \bm {x}- \bm {b}\| ^2_2$ .

Use the QR decomposition of ${\bA }$ with column pivoting,

$\bA = \bQ \left[ \begin{array}{c} \bR \\ \mb {0} \end{array} \right] \bPi ^{\prime } = \left[ \begin{array}{cc} \bY & \bZ \end{array} \right] \left[ \begin{array}{cc} \bR _1 & \bR _2 \\ \mb {0} & \mb {0} \end{array} \right] \bPi ^{\prime }$

where $\bR = [ \begin{array}{cc} {\bR }_1 & {\bR }_2 \end{array}] \in {\mathcal R}^{r \times t}$ is upper trapezoidal, ${\bR }_1\in {\mathcal R}^{r \times r}$ is upper triangular and invertible, ${\bR }_2\in {\mathcal R}^{r \times s}$ , ${\bQ } = [\begin{array}{cc} {\bY } & {\bZ } \end{array}]$ is orthogonal, ${\bY }\in {\mathcal R}^{t \times r}$ , ${\bZ }\in {\mathcal R}^{t \times s}$ , and ${\bPi }$ permutes the columns of ${\bA }$ .
Use the transpose ${\bL }_{12}$ of the upper trapezoidal matrix ${\bR } = \left[ \begin{array}{cc} {\bR }_1 & {\bR }_2 \end{array} \right]$ ,

${\bL }_{12} = \left[ \begin{array}{c} {\bL }_1 \\ {\bL }_2 \end{array} \right] = {\bR }^{\prime } \in {\mathcal R}^{t \times r}$

with $\mr {rank}({\bL }_{12}) = \mr {rank}({\bL }_1) = r$ , ${\bL }_1\in {\mathcal R}^{r \times r}$ lower triangular, ${\bL }_2\in {\mathcal R}^{s \times r}$ . The lower trapezoidal matrix ${\bL }_{12}\in {\mathcal R}^{t \times r}$ is premultiplied with $r$ Householder transformations ${\bP }_1,\ldots ,{\bP }_ r$ ,

$\bP _ r \ldots \bP _1 \left[ \begin{array}{c} \bL _1 \\ \bL _2 \end{array} \right] = \left[ \begin{array}{c} \bL \\ \mb {0} \end{array} \right]$

each zeroing out one of the $r$ columns of ${\bL }_2$ and producing the nonsingular lower triangular matrix ${\bL }\in {\mathcal R}^{r \times r}$ . Therefore, you obtain

$\bA = \bQ \left[ \begin{array}{cc} \bL ^{\prime } & \mb {0} \\ \mb {0} & \mb {0} \end{array} \right]\bPi ^{\prime } \bP ^{\prime } = \bY \left[ \begin{array}{cc} \bL ^{\prime } & \mb {0} \end{array} \right] \bPi ^{\prime } \bP ^{\prime }$

with $\bP = {\bPi }{\bP }_ r\ldots {\bP }_1\in {\mathcal R}^{t \times t}$ and upper triangular ${\bL }^{\prime }$ . This second step is described in Golub and Van Loan (1989).
Compute the Moore-Penrose inverse ${\bA }^-$ explicitly:

$\bA ^- = \bP \bPi \left[ \begin{array}{cc} (\bL ^{\prime })^{-1} & \mb {0} \\ \mb {0} & \mb {0} \end{array} \right] \bQ ^{\prime } = \bP \bPi \left[ \begin{array}{c} (\bL ^{\prime })^{-1} \\ \mb {0} \end{array} \right] \bY ^{\prime }$
- Obtain ${\bY }$ in ${\bQ } = \left[\begin{array}{cc} {\bY } & {\bZ } \end{array}\right]$ explicitly by applying the $r$ Householder transformations obtained in the first step to $\left[ \begin{array}{c} {\bI }_ r \\ \mb {0} \end{array} \right]$ .
- Solve the $r \times r$ lower triangular system $({\bL }^{\prime })^{-1}{\bY }^{\prime }$ with $t$ right-hand sides by using backward substitution, which yields an $r \times t$ intermediate matrix.
- Left-apply the $r$ Householder transformations in ${\bP }$ on the $r \times t$ intermediate matrix $\left[ \begin{array}{c} ({\bL }^{\prime })^{-1}{\bY }^{\prime } \\ \mb {0} \\ \end{array} \right]$ , which results in the symmetric matrix ${\bA }^-\in {\mathcal R}^{t \times t}$ .

The GINV function computes the Moore-Penrose inverse ${\bA }^-$ by using the singular value decomposition of ${\bA }$ . Using complete orthogonal decomposition to compute ${\bA }^-$ 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 ${\bA }$ , than the singular value decomposition.

The following example demonstrates some uses of the COMPORT subroutine:

/* 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
   msg = "Pseudoinverse conditions not satisfied";
else
   msg = "Pseudoinverse conditions satisfied";
print msg;

Figure 24.75: Results from a Complete Orthogonal Factorization

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

msg
Pseudoinverse conditions satisfied