INVUPDT Function :: SAS/IML(R) 13.1 User's Guide

INVUPDT Function

INVUPDT (matrix, vector<, scalar>) ;

The INVUPDT function updates a matrix inverse.

The arguments to the INVUPDT function are as follows:

matrix: is an $n \times n$ nonsingular matrix. In most applications matrix is symmetric positive definite.
vector: is an $n \times 1$ or $1 \times n$ vector.
scalar: is a numeric scalar.

The Sherman-Morrison-Woodbury formula is

$(\mb {A} + \mb {U}\mb {V}^{\prime })^{-1} = \mb {A}^{-1} - \mb {A}^{-1} \mb {U} (\mb {I} + \mb {V}^{\prime }\mb {A}^{-1}\mb {U})^{-1} \mb {V}^{\prime }\mb {A}^{-1}$

where $\mb {A}$ is an $n \times n$ nonsingular matrix and $\mb {U}$ and $\mb {V}$ are $n\times k$ . The formula shows that a rank $k$ update to $\mb {A}$ corresponds to a rank $k$ update of $\mb {A}^{-1}$ .

The INVUPDT function implements the Sherman-Morrison-Woodbury formula for rank-one updates with $\mb {U} = w \mb {X}$ and $\mb {V} = \mb {X}$ , where $\mb {X}$ is an $n \times 1$ vector and $w$ is a scalar.

If $\mb {M} = \mb {A}^{-1}$ , then you can call the INVUPDT function as follows:

R = invupdt(M, X, w);

This statement computes the following matrix:

$\mb {R} = \mb {M} - w\mb {MX} (\mb {I} + w\mb {X}^{\prime }\mb {MX})^{-1} \mb {X}^{\prime }\mb {M}$

The matrix $\mb {R}$ is equivalent to $(\mb {A} + w\mb {XX}^{\prime })^{-1}$ . If $\mb {A}$ is symmetric positive definite, then so is $\mb {R}$ .

If $w$ is not specified, then it is given a default value of $1$ .

A common use of the INVUPDT function is in linear regression. If $\mb {Z}$ is a design matrix, $\mb {M} = (\mb {Z}^{\prime }\mb {Z})^{-1}$ is the associated inverse crossproduct matrix, and $\mb {v}$ is a new observation to be used in estimating the parameters of a linear model, then the inverse crossproducts matrix that includes the new observation can be updated from $\mb {M}$ by using the following statement:

M2 = invupdt(M, v);

If $w$ is $1$ , the function adds an observation to the inverse; if $w$ is $-1$ , the function removes an observation from the inverse. If weighting is used, $w$ is the weight.

To perform the computation, the INVUPDT function uses about $2n^2$ multiplications and additions, where $n$ is the row dimension of the positive definite argument matrix.

The following program demonstrates adding or removing observations from a linear fit and updating the inverse crossproduct matrix:

X = {0, 1, 1, 1, 2, 2, 3, 4, 4};
Y = {1, 1, 2, 6, 2, 3, 3, 3, 4};

/* find linear fit */
Z = j(nrow(X), 1, 1) || X;           /* design matrix */
M = inv(Z`*Z);

b = M*Z`*Y;                          /* LS estimate */
resid = Y - Z*b;                     /* residuals */
print "Original Fit", b resid;

/* residual for observation (1,6) seems too large.
   Take obs number 4 out of data set and refit. */
v = z[4,];
M = invupdt(M, v, -1);               /* update inverse crossprod */

keepObs = (1:3) || (5:nrow(X));
Z = Z[keepObs, ];
Y = Y[keepObs, ];
b = M*Z`*Y;                          /* new LS estimate */
print "After deleting observation 4", b;

/* Add a new obs (x,y) = (0,2) and refit. */
obs = {0 2};
v = 1 || obs[1];                     /* new row in design matrix */
M = invupdt(M, v);

Z = Z // v;
Y = Y // obs[2];
b = M*Z`*Y;                          /* new LS estimate */
print "After adding observation (0,2)", b;

Figure 24.169: Refitting Linear Regression Models

Original Fit

b	resid
2.0277778	-1.027778
0.375	-1.402778
	-0.402778
	3.5972222
	-0.777778
	0.2222222
	-0.152778
	-0.527778
	0.4722222

After deleting observation 4

b
1
0.6470588

After adding observation (0,2)

b
1.3
0.5470588