spaugment

Purpose

Create an augmented matrix associated with least squares problems.

Synopsis

S = spaugment(A,c)
S = spaugment(A)

Description

S = spaugment(A,c) creates the sparse, square, symmetric, and indefinite matrix. 

S = [c*I  A; A'  0]
The augmented matrix can compute the solution to the least squares problem.

If r = b - Ax is the residual and c is a residual scaling factor, the least squares problem is equivalent to the system of linear equations

If A is m-by-n, the augmented matrix S is (m+n)-by-(m+n). The augmented linear system can be written 

S * [r/c; x] = [b; z]
where z = zeros(n,1).

S = spaugment(A) uses a value of max(max(abs(A)))/1000 for c. The optimum value of the residual scaling factor, c, involves min(svd(A)) and norm(r), which are usually too expensive to compute. You can use the command spparms to change the default.

The sparse backslash operation x = A\b automatically forms the same augmented matrix, with the default value of c. spparms can set other values of c for use with the backslash operation.

Examples

Here is a simple system of three equations in two unknowns.

Since it is such a small system, it is ordinarily solved in MATLAB with

A = [1  2; 3  0; 0  4];
b = [5  6  7]';
x = A\b
This produces the least squares solution

x =
    1.9592
    1.7041
The sparse augmented systems approach uses

S = spaugment(A,1)
which gives

full(S) =
    1    0    0    1    2
    0    1    0    3    0
    0    0    1    0    4
    1    3    0    0    0
    2    0    4    0    0
This example used c = 1 and printed full(S) to improve readability. The default value of c for this example is 0.0040. Even in this small example, more than half of the elements of S are zero.

The augmented systems also use an extended right-hand side

y = [b; 0; 0];
Now

z = S\y
involves a square matrix, so MATLAB factors that matrix by Gaussian elimination instead of orthogonalization. For sparse matrices, this usually involves less fill-in and requires less storage.

The computed solution is

z =
    -0.3673
     0.1224
     0.1837
     1.9592
     1.7041
The first three components are the residual r and the last two components are the same solution x, as was computed by the first approach.

See Also

\, diag, spparms

(c) Copyright 1994 by The MathWorks, Inc.