Least Squares and the Eigenvalue Problem

Least Squares

Polynomial data - Figure 10.1

Here we fit polynomial data with sine functions (the opposite of the fitting example from Chapter 9) - this example is in Section 10.2.  First define the data function

d[t_] :=t (t-1/5)(t-2/5)(t-3/5)(t-4/5)(t-1)

now set up the matrices for the residual minimization defined by (10.3) and solved from (10.9):

dt = 1/100.0;
Fmat = Table[Sin[k Pi t],{t,dt,1,dt},{k,1,10}];
b = Table[d[t],{t,dt,1,dt}];

We’ll use the built-in QRDecomposition function (as opposed to the one in the previous, Chapter 9, notebook) - note that Mathematica returns the transpose of what we call Q (due to its obsession with row-based matrices), hence the “Transpose” below:

{Q,R} = QRDecomposition[Fmat];
Q = Transpose[Q];

Again, using a built-in to perform the matrix inversion, we have the solution, x, to (10.9), and can write the finite sine-series sum in terms of those coefficients

x = LinearSolve[R, Transpose[Q].b];
func = Sum[x[[k]] Sin[k Pi t],{k,1,Length[x]}]

A plot of the two shows that the approximation is not too bad:

G0 = Plot[d[t],{t,0,1},PlotStyle->GrayLevel[0.5]];
G1 = Plot[func,{t,0,1},PlotStyle->GrayLevel[.25]];
Show[G0,G1]

The Eigenvalue Problem

Quantum Mechanics Examples

The function below makes the matrix from (10.27) with potential V inserted (dimensionless, and defined as a function of a single variable q = 0 -> Nn dq).  Note that in the construction of this matrix,
we assume that the wave function vanishes at zero and the right-hand end-point, Nn dq - we have immersed the problem in an infinite square well (or, alternatively, brought spatial infinity in to a pair of finite values).

FDSchrodingerMat[V_, Nn_, dq_] := Module[{retmat,index},
retmat = Table[0.0,{i,1,Nn},{j,1,Nn}];
retmat[[1,1]] = -2.0 - dq^2 V[dq];
retmat[[1,2]] = 1.0;
For[index = 2, index ≤ Nn-1, index=index+1,
retmat[[index,index-1]] = 1.0;
retmat[[index,index]] =-2.0 - dq^2 V[index dq];
retmat[[index,index+1]] = 1.0;
];
retmat[[Nn,Nn-1]] = 1.0;
retmat[[Nn,Nn]] = -2.0 - dq^2 V[Nn dq];
retmat = -1.0/dq^2 retmat;
Return[retmat];
]

Square Well - Figure 10.2

Our first example is the infinite square well, with V = 0

V[q_] := 0

We’ll set 102 grid points extending from 0 to 1 - the action occurs on the interior 100 gridpoints (excluding 0 and 1)

Nn = 100;
dq = 1.0/(Nn+1);
Mat = FDSchrodingerMat[V,Nn,dq];

For now, we can use the built-in function Eigenvalues to compute the eigenvalues

vals = Eigenvalues[Mat];

Take just the first 20 for comparison (we know the spectrum here, just (((k π))^2) for integer k):

firstfew = Sort[vals][[1;;20]];
G0 = ListPlot[firstfew,PlotStyle->GrayLevel[0.0]];
G1 = ListPlot[Table[k^2 Pi^2,{k,1,20}],PlotStyle->GrayLevel[0.5]];
Gd = Show[G1,G0]

Using the built-in “Eigensystem” command, we have access to the eigenvectors as well as the eigenvalues, and can establish that the correct wave functions
emerge numerically:

{vals,vecs} = Eigensystem[Mat];
vec1 = vecs[[Nn]];
vec2 = vecs[[Nn-1]];
vec3 = vecs[[Nn-2]];
vec1 = Table[{j dq, vec1[[j]]},{j,1,Length[vec1]}];
vec2 = Table[{j dq, vec2[[j]]},{j,1,Length[vec2]}];vec3 = Table[{j dq, vec3[[j]]},{j,1,Length[vec3]}];

Ga = Show[{ListPlot[vec1,PlotStyle->GrayLevel[0.25]],Plot[Sqrt[2 dq] Sin[Pi x],{x,0,1},PlotStyle->GrayLevel[0.5]]}];
Gb = Show[{ListPlot[vec2,PlotStyle->GrayLevel[0.25]],Plot[Sqrt[2 dq] Sin[2Pi x],{x,0,1},PlotStyle->GrayLevel[0.5]]}];
Gc = Show[{ListPlot[vec3,PlotStyle->GrayLevel[0.25]],Plot[Sqrt[2 dq] Sin[3 Pi x],{x,0,1},PlotStyle->GrayLevel[0.5]]}];

GraphicsGrid[{{Ga,Gb},{Gc,Gd}}]

Step in Well - Figure 10.5

We can put a step inside the well, in this case, extending from the middle to the right-hand side - this example appears in Section 10.7.1, where we can compare the numerical
result with the perturbative one provided the step is small (it’s not, for this first case)

V[q_] := If[q ≥ .5, 475.,0.0]

Nn = 100;
dq = 1.0/(Nn+1);
Mat = FDSchrodingerMat[V,Nn,dq];

Take just the first 10 eigenvalues to get a sense of the growth:

vals = Sort[Eigenvalues[Mat]][[1;;10]]

We can get the eigenvectors, too:

{vals,vecs} = Eigensystem[Mat];
vec1 = vecs[[Nn]];
vec2 = vecs[[Nn-1]];
vec3 = vecs[[Nn-2]];
vec4 = vecs[[Nn-3]];
vec1 = Table[{j dq, vec1[[j]]^2},{j,1,Length[vec1]}];
vec2 = Table[{j dq, vec2[[j]]^2},{j,1,Length[vec2]}];
vec3 = Table[{j dq, vec3[[j]]^2},{j,1,Length[vec3]}];
vec4 = Table[{j dq, vec4[[j]]^2},{j,1,Length[vec4]}];
Ga = Show[{ListPlot[vec1,PlotStyle->GrayLevel[0.25]],Plot[(Sqrt[2 dq] Sin[Pi x])^2,{x,0,1},PlotStyle->GrayLevel[0.5]]}];
Gb = Show[{ListPlot[vec2,PlotStyle->GrayLevel[0.25]],Plot[(Sqrt[2 dq] Sin[2Pi x])^2,{x,0,1},PlotStyle->GrayLevel[0.5]]}];
Gc = Show[{ListPlot[vec3,PlotStyle->GrayLevel[0.25]],Plot[(Sqrt[2 dq] Sin[3 Pi x])^2,{x,0,1},PlotStyle->GrayLevel[0.5]]}];
Gd = Show[{ListPlot[vec4,PlotStyle->GrayLevel[0.25]],Plot[(Sqrt[2 dq] Sin[4 Pi x])^2,{x,0,1},PlotStyle->GrayLevel[0.5]]}];
GraphicsGrid[{{Ga,Gb},{Gc,Gd}}]

If we instead use a small value for the step, then we can compare with time-independent perturbation theory:

V[q_] := If[q >= .5-dq, 1.,0.0]

Nn = 2000;
dq = 1.0/(Nn+1);
Mat = FDSchrodingerMat[V,Nn,dq];

vals = Eigenvalues[Mat];

Sort[vals][[1;;10]]

and we can compare with the perturbative result:

N[Table[n^2 Pi^2 + 1/2,{n,1,10}]]

Bump in Well

We can also put a bump in the middle of the well, modeling a repulsive force at the center - one can check this setup against a delta spike.

V[q_] := If[q ≥ .49 && q ≤ .51, 100.,0.0]

Nn = 100;
dq = 1.0/(Nn+1);
Mat = FDSchrodingerMat[V,Nn,dq];

{vals,vecs} = Eigensystem[Mat];

Here are the first few eigenvectors, obtained numerically:

G0 = ListLinePlot[vecs[[100]], PlotRange-> {-.15,.15}];
G1 = ListLinePlot[vecs[[99]], PlotStyle-> RGBColor[0,0,1]];
G2 = ListLinePlot[vecs[[98]], PlotStyle -> RGBColor[1,0,0]];
G3 = ListLinePlot[vecs[[97]],  PlotStyle -> RGBColor[0,1,0]];
Show[G0,G1,G2,G3]

The Power Method

The method, described in Section 10.5, amounts to multiplying by the matrix A over and over.  Start with a 3x3 symmetric matrix, and an arbitrary vector x:

A = N[{{1,2,3},{2,4,5},{3,5,6}}];
x = {1.,0.,0.};
MatrixForm[A]

Multiply by A repeatedly:

w =A.A. A.A.A.A.A.A.A.x

Normalize the resulting vector so that it has length 1:

what =N[ w/Sqrt[w.w]]

The eigenvalue, if what (read double-u-hat) is an eigenvector, is

L1 = what.A.what

Let’s check

A.what - L1 what

that looks pretty good.

Image Compression - Figure 10.3, 10.4

Read in the jpeg data, trim it a little, and then we’ll take the RGB colors and project them onto the the gray-scale line, zeroing small entries (the magic numbers below, like Sqrt[3] - .1, are just to clean up the image):

data = Import["http://people.reed.edu/~jfrankli/Book/downloads/files/Smiley.jpg","Data"];
data = data[[100;;420,100;;420]];
sdata = Table[({255,255,255} - data[[i,j]]).{1,1,1}/Sqrt[3 255^2],{i,1,Length[data]},{j,1,Length[data]}];
sdata = Table[If[sdata[[i,j]] < Sqrt[3] - .1,0,sdata[[i,j]]],{i,1,Length[data]},{j,1,Length[data]}];

The image is already pretty symmetric, and we can ensure symmetry, after our surgery above, by taking an equal admixture of sdata and its transpose:

sdata =(1/2)(sdata + Transpose[sdata]);

MatrixPlot[sdata]

Find the eigenvalues and eigenvectors of this matrix:

{S,V} = Eigensystem[N[sdata]];
SM = DiagonalMatrix[S];

The eigenvalues decay in a way that will allow compression:

ListPlot[Table[SM[[j,j]],{j,1,Length[SM]}],PlotRange->{-70,70},PlotStyle->GrayLevel[0.5]]

The error we make in truncating can be estimated using (10.51) - if we just take the first 128 eigenvectors, then the norm of A - (OverTilde[A]) is less than:

Sqrt[Sum[S[[j]]^2,{j,129,Length[S]}]]

compare that with the norm of the matrix itself:

Norm[N[sdata],2]

We'll take the first Nlen (10 below) eigenvectors to approximate the matrix

Nlen = 10;
Vs = Take[V,{1,Nlen}];
SM = DiagonalMatrix[Take[S,{1,Nlen}]];
PseudoG0 =Chop[Transpose[Vs].SM.Vs];
MatrixPlot[PseudoG0]

Pretty good, let's see how the compression works at different values of Nlen.

lens = {2,4,8,16,32,64,128,256,Length[sdata]};
gout = {};
For[index = 1, index ≤ Length[lens], index = index+1,
Nlen = lens[[index]];
Vs = Take[V,{1,Nlen}];
SM = DiagonalMatrix[Take[S,{1,Nlen}]];
PseudoG0 =Chop[Transpose[Vs].SM.Vs];
 gout = Join[gout,{MatrixPlot[PseudoG0,ColorFunction->"GrayTones",Frame->False]}];
];
GG = GraphicsGrid[Partition[gout,3]]

At 64 eigenvectors, our compression is:

N[64 Length[sdata]/(Length[sdata]^2)]