Language Reference


APPCORT Call

CALL APPCORT (prqb, lindep, a, b, <, sing> );

If $\bA $ is rank-deficient, then the least squares problem $\min _{x} \|  \bA \mb{x} - \mb{b} \| ^2_2$ 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 $\mb{b}$.

The input arguments to the APPCORT subroutine are as follows:

a

is an $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

is a $m \times p$ matrix, $\bB $.

sing

is an optional scalar that specifies a singularity criterion.

The APPCORT subroutine returns the following values:

prqb

is an $n \times p$ matrix product

\[ \bP \bPi \left[ \begin{array}{cc} (\bL ^{\prime })^{-1} & \mb{0} \\ \mb{0} & \mb{0} \end{array} \right] \bQ ^{\prime } \bB \]

which is the minimum Euclidean-norm solution of the rank-deficient least squares problem $\|  \bA \mb{x} - \mb{b} |^2_2$.

lindep

is the number of linearly dependent columns in the matrix $\bA $ that are detected by applying the r Householder transformations. That is, lindep is $n-r$, where r is the numerical rank of $\bA $.

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;

Figure 25.39: Solution to a Rank-Deficient Least Squares Problem

x
0.3
0.6



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 25.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