If
is rank-deficient, then the least squares problem
has infinitely many solutions (Golub and Van Loan, 1989, p. 241). However, there is a unique solution which has the smallest Euclidean norm. The APPCORT subroutine computes the
minimum Euclidean-norm solution of the (rank-deficient) least squares problem by applyng a complete orthogonal decomposition
by Householder transformations to the vector
.
The input arguments to the APPCORT subroutine are as follows:
is an
matrix
, with
, which is to be decomposed into the product of the
orthogonal matrix
, the
upper triangular matrix
, and the
orthogonal matrix
,
is a
matrix,
.
is an optional scalar that specifies a singularity criterion.
The APPCORT subroutine returns the following values:
is an
matrix product
which is the minimum Euclidean-norm solution of the rank-deficient least squares problem
.
is the number of linearly dependent columns in the matrix
that are detected by applying the
Householder transformations. That is, lindep is
, where
is the numerical rank of
.
See the section COMPORT Call for information about complete orthogonal decomposition.
The following example uses the APPCORT call to solve a rank-deficient least squares problem:
/* compute solution for rank-deficient least squares problem:
min |Ax-b|^2
The range of A is a line; b is a point not on the line. */
A = {1 2,
2 4,
-1 -2};
b = {1, 3, -2};
call appcort(x,lindep,A,b);
print x;
The argument b can also be a matrix. If b is an identity matrix, then you can use the APPCORT subroutine to form a generalized inverse, as shown in the following example:
/* A has 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 };
/* compute Moore-Penrose generalized inverse */
b = i(nrow(A)); /* identity matrix */
call appcort(Ainv, lindep, A, b);
print Ainv;
/* verify generalized inverse conditions (Golub & Van Loan, p. 243) */
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.40: Generalized Inverse
| 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 |