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

The subroutine returns the following values:

q

is a matrix. If is not specified, is the orthogonal matrix that is the product of the separate Householder transformations. If is specified, is the matrix that has the transposed Householder transformations applied to the columns of the argument matrix .

r

is the upper triangular matrix that contains the nonsingular upper triangular matrix of the complete orthogonal decomposition, where is the rank of . The full upper triangular matrix of the orthogonal decomposition of matrix can be obtained by vertical concatenation of the zero matrix to the result r.

p

is an matrix that is the product of a permutation matrix 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 . The vector piv corresponds to an 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 detected by applying the Householder transformation in the order specified by the argument piv. That is, lindep is .

The input arguments to the COMPORT subroutine are as follows:

a

specifies the matrix , with , which is to be decomposed into the product of the orthogonal matrix , the upper triangular matrix , and the orthogonal matrix ,

b

specifies an optional matrix that is to be left-multiplied by the transposed matrix .

sing

is an optional scalar that specifies a singularity criterion.

The complete orthogonal decomposition of the singular matrix can be used to compute the Moore-Penrose inverse , , or to compute the minimum Euclidean-norm solution of the rank-deficient least squares problem .

1. Use the QR decomposition of with column pivoting,

   where is upper trapezoidal, is upper triangular and invertible, , is orthogonal, , , and permutes the columns of .

2. Use the transpose of the upper trapezoidal matrix ,

   with , lower triangular, . The lower trapezoidal matrix is premultiplied with Householder transformations ,

   each zeroing out one of the columns of and producing the nonsingular lower triangular matrix . Therefore, you obtain

   with and upper triangular . This second step is described in Golub and Van Loan (1989).

3. Compute the Moore-Penrose inverse explicitly:

   • Obtain in explicitly by applying the Householder transformations obtained in the first step to .

   • Solve the lower triangular system with right-hand sides by using backward substitution, which yields an intermediate matrix.

   • Left-apply the Householder transformations in on the intermediate matrix , which results in the symmetric matrix .

The GINV function computes the Moore-Penrose inverse by using the singular value decomposition of . Using complete orthogonal decomposition to compute 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 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;