Matrix Inversion

QR Solver - Algorithm 9.1, Gram-Schmidt

Define a numerical zero:

NZERO = 10^(-10);

The function below accomplishes the Gram-Schmidt QR decomposition from Algorithm 9.1 - note that the Q returned has columns (not rows) that are an orthonormal basis for the space spanned by the columns of A.

QRDecompose[A_] := Module[{Q,R,i,j,At,qi},
At = Transpose[A];
R = Table[0,{ii,1,Length[A]},{jj,1,Length[A]}];
Q = {};
For[j = 1, j ≤ Length[A], j = j+1,
qi =N[ At[[j]] ];
For[i = 1, i ≤ j - 1, i = i+1,
R[[i,j]] = Q[[i]].At[[j]];
qi = qi - R[[i,j]] Q[[i]];
];
R[[j,j]] = Sqrt[qi.qi];
If[Abs[R[[j,j]]] ≤ NZERO, Print["Matrix A may be singular"]; Return[{Null,Null}];];
Q = Join[Q,N[{qi/R[[j,j]]}]];
];
Return[{Transpose[Q],R}]
]

Below we solve for x given A and b using the QR decomposition - first we make Q R = A, then perform the back substitution on the upper triangular matrix R.

QRSolve[A_,b_] := Module[{Q,R,x,brhs},
{Q,R} = QRDecompose[A];
If[Q === Null, Print["Matrix A may be singular"]; Return[Null];];
brhs = Transpose[Q].b;
x = UDBacksolve[R,brhs];
Return[x];
]

Back substitution - the pattern comes from (9.18) - (9.20) and is similar in spirit and form to forward subsitution - compare with (9.38).

UDBacksolve[U_,b_] := Module[{xout,i,j,NN},
NN = Length[U];
xout = Table[0,{ii,1,NN}];
xout[[NN]] = {b[[NN]]/U[[NN,NN]]};
For[i = NN-1, i ≥ 1, i = i-1,
xout[[i]] = N[(b[[i]]- Sum[U[[i,jj]]xout[[jj]] ,{jj,i+1,NN}])/U[[i,i]]];
];
xout = Flatten[xout];
Return[xout];
]

Matrix Inversion Tests

Small test

Start with a 4x4 matrix,

A = {{11,2,3,4},{5,60,7,8},{20,21,22,23},{24,25,26,27}};
b = {1,2,3,4};

x = QRSolve[A,b]

the residual has small entries:

A.x - b

Large test

Random matrices with entries drawn on [-1,1] are . . . rarely singular.

ttot = 100;
Ns = 100;
qrres = {};

We’ll perform the inverse on a set of 100 randomly generated matrices, and calculate the maximum absolute entry of the residual vector A x - b:

For[index = 1, index ≤ ttot,index = index+1,
A = Table[RandomReal[{-1,1}],{ii,1,Ns},{jj,1,Ns}];
b = Table[RandomReal[{-1,1}],{ii,1,Ns}];
x = QRSolve[A,b];
qrres = Join[qrres,{Max[Abs[A.x-b]]}];
];

The mean of the maximum (absolute value) element of the residuals is small,

Mean[qrres]

and we can see the scatter over our hundred test matrices below

ListPlot[qrres]

Examples - Figures 9.2 and 9.3

The Vandermonde function takes a set of grid points in x and constructs the Vandermonde matrix appropriate to polynomial approximation of a set of data on the grid (see (9.13) and (9.15)).

VanderMonde[t_] :=Transpose[ Table[t[[kk]]^jj,{jj,0,Length[t]-1},{kk,1,Length[t]}]]

Polynomial Exact

First, we’ll fit a polynomial test function using a polynomial “basis”, as in Section 9.2.2

F[t_] := 1 - 2 t + 3 t^2 - 4 t^3 + 5 t^4 - 6 t^5 + 7 t^6

G0 = Plot[F[t],{t,0,1},AxesOrigin->{0,0},PlotStyle->GrayLevel[0.5]];

The Vanderomonde matrix looks like this

MatrixForm[VanderMonde[{.1,.2,.3,.4,.5,.6,.7}]]

Now set up the problem, in A x = b form, and solve for x - the resulting polynomial is called pF below, and matches the input one

dt = 1./7;
b = Table[F[t],{t,dt,1, dt}];
A = VanderMonde[Table[t,{t,dt,1.,dt}]];
x =QRSolve[A,b];
pF = Sum[x[[jj+1]] t^jj,{jj,0,Length[x]-1}]

G1 = Show[Graphics[Table[Point[{j dt, F[ j dt]}],{j,1,7}]]];
Show[G0,G1]

Polynomial Inexact

In this example, we are missing a term in our Vandermonde matrix - our polynomial is order 6, and our Vandermonde matrix only goes up to order 5.

F[x_] := (x-.1)(x-.4)(x-.5)(x-1)(x+.1)x

G0 = Plot[F[y],{y,.1,.6}]

Again, we set up the problem, and solve, obtaining the polynomial pF below

b = Table[F[x],{x,.1,.6,.1}];
A = VanderMonde[Table[y,{y,.1,.6,.1}]];
x =QRSolve[A,b];
pF = Sum[x[[jj+1]] y^jj,{jj,0,Length[x]-1}]

Compare pF with the original polynomial:

Expand[F[y]]

Beyond the form of the polynomial, we can look at a plot of both the original and the fit:

G1 = Plot[pF,{y,.1,.6}];
Show[G0,G1]

the residual is:

Plot[pF - F[y],{y,.1,.6}]

Oscillatory

Here we try a polynomial approximation to a sinusoidal function as in Section 9.2.2.  The original “data”:

d[t_] := Sin[16 Pi t]

a simple oscillatory signal:

dt = 1.0/7;
G0 = Plot[d[t],{t,0,1+dt},PlotStyle->GrayLevel[0.5]]

Set up the matrix inverse problem using a polynomial “basis”, solve for x and then define the approximating polynomial pF:

b = Table[d[t],{t,dt,1.,dt}];
A = VanderMonde[Table[t,{t,dt,1.,dt}]];
x =QRSolve[A,b];
pF = Sum[x[[jj+1]] t^jj,{jj,0,Length[x]-1}]

G1 = Plot[pF,{t,0,1},PlotStyle->GrayLevel[0.25]];
G2 =Show[Graphics[ Table[Point[{j dt, d[j dt]}],{j,1,7}]]];
Show[{G0,G1,G2}]

Above, we have the original data, the seven sample points used to fit the polynomial, and the polynomial approximation - note that it matches the data at the seven points shown, but is otherwise unlike the data.