## Designing and Analyzing Linear Systems¶

Supplemental Julia Notebook

Prof. Doug James

## Example: Tikhonov regularization¶

In [29]:
A = [1 1;
1 1.00001]   ## An ill-conditioned matrix
cond(A)

Out[29]:
400002.0000028382
In [31]:
b = [1; 0.99];

In [32]:
x = A \ b

Out[32]:
2-element Array{Float64,1}:
1001.0
-1000.0
In [33]:
norm(A*x-b)/norm(b)

Out[33]:
6.46964603432798e-15
In [35]:
norm(x)

Out[35]:
1414.920845833358

## Regularized solve¶

In [24]:
α = 0.00000000001
B = A'*A + α*eye(2)  ## Regularized

Out[24]:
2×2 Array{Float64,2}:
2.0      2.00001
2.00001  2.00002
In [25]:
x = B \ (A'*b)  ## Solve Normal equations

Out[25]:
2-element Array{Float64,1}:
1000.99
-999.995
In [26]:
relErr = norm(A*x-b)/norm(b)

Out[26]:
2.552917567500873e-8
In [27]:
norm(x)

Out[27]:
1414.9136610830515

## Cholesky Factorization¶

In [1]:
?(chol)

search: chol cholfact cholfact! searchsortedlast CachingPool chop chown chomp


Out[1]:
chol(A) -> U

Compute the Cholesky factorization of a positive definite matrix A and return the UpperTriangular matrix U such that A = U'U.

chol(x::Number) -> y

Compute the square root of a non-negative number x.

In [245]:
chol(rand(3,3))

ArgumentError: matrix is not symmetric/Hermitian. This error can be avoided by calling chol(Hermitian(A)) which will ignore either the upper or lower triangle of the matrix.

in chol(::Array{Float64,2}) at ./linalg/cholesky.jl:166
In [7]:
U = chol(eye(3))

Out[7]:
3×3 UpperTriangular{Float64,Array{Float64,2}}:
1.0  0.0  0.0
⋅   1.0  0.0
⋅    ⋅   1.0
In [5]:
U'*U

Out[5]:
3×3 Array{Float64,2}:
1.0  0.0  0.0
0.0  1.0  0.0
0.0  0.0  1.0
In [247]:
R = rand(3,3);
A = R' * R

Out[247]:
3×3 Array{Float64,2}:
1.03716   0.485069  0.700163
0.485069  0.273739  0.295844
0.700163  0.295844  0.765217
In [248]:
L = chol(A)'

Out[248]:
3×3 LowerTriangular{Float64,Array{Float64,2}}:
1.01841     ⋅         ⋅
0.476301   0.216509   ⋅
0.687507  -0.14603   0.520793
In [249]:
b=ones(3);

In [250]:
A \ b

Out[250]:
3-element Array{Float64,1}:
-6.84456
13.0564
2.5217 
In [251]:
y = L \ b

Out[251]:
3-element Array{Float64,1}:
0.981924
2.45859
1.31328 
In [252]:
L' \ y

Out[252]:
3-element Array{Float64,1}:
-6.84456
13.0564
2.5217 

## Cholesky on Sparse Matrices¶

In [253]:
using PyPlot
N = 15;
R = sprand(N, N, 0.1);
A = R*R' + 1.e-7 * speye(N)

Out[253]:
15×15 sparse matrix with 29 Float64 nonzero entries:
[1 ,  1]  =  1.0e-7
[2 ,  2]  =  0.959033
[4 ,  2]  =  0.0536958
[5 ,  2]  =  0.568696
[9 ,  2]  =  0.188129
[3 ,  3]  =  1.11993
[5 ,  3]  =  0.1869
[2 ,  4]  =  0.0536958
[4 ,  4]  =  0.131619
[9 ,  4]  =  0.237006
⋮
[4 ,  9]  =  0.237006
[9 ,  9]  =  2.02699
[10, 10]  =  0.807924
[12, 10]  =  0.349218
[11, 11]  =  1.0e-7
[10, 12]  =  0.349218
[12, 12]  =  0.150947
[6 , 13]  =  0.211823
[13, 13]  =  0.212114
[14, 14]  =  1.0e-7
[15, 15]  =  1.0e-7
In [254]:
spy(A)

Out[254]:
PyObject <matplotlib.image.AxesImage object at 0x7fb20c898d50>
In [256]:
svdvals(full(A))'

Out[256]:
1×15 Array{Float64,2}:
2.09612  1.34565  1.09208  0.958871  …  1.0e-7  1.0e-7  1.0e-7  1.0e-7
In [257]:
chol(A)

Use cholfact() instead of chol() for sparse matrices.

in chol(::SparseMatrixCSC{Float64,Int64}) at ./sparse/linalg.jl:918
In [258]:
F = cholfact(A)

Out[258]:
Base.SparseArrays.CHOLMOD.Factor{Float64}
type:          LLt
method: simplicial
maxnnz:         22
nnz:            22

In [259]:
b = ones(N);
xF = F \ b

Out[259]:
15-element Array{Float64,1}:
1.0e7
-2.92591
-0.0943949
9.39129
5.91608
6867.08
2.21476
1.0e7
-0.333173
-2.06776e6
1.0e7
4.78381e6
-6852.95
1.0e7
1.0e7    
In [260]:
xA = A \ b   ## "\" is smart, and actually uses sparse Cholesky to solve system(!)

Out[260]:
15-element Array{Float64,1}:
1.0e7
-2.92591
-0.0943949
9.39129
5.91608
6867.08
2.21476
1.0e7
-0.333173
-2.06776e6
1.0e7
4.78381e6
-6852.95
1.0e7
1.0e7    
In [261]:
norm(xA-xF)/norm(xA)

Out[261]:
2.5888015819409084e-11
In [ ]:



## 1D Laplacian matrix¶

In [174]:
N=100;
A = spzeros(N,N)  ## sparse zero matrix

Out[174]:
100×100 sparse matrix with 0 Float64 nonzero entries
In [175]:
## Create 1D Laplacian matrix (h=1)
for i=1:N
A[i,i] = -2.;
if (i>1) A[i,i-1] = 1; end
if (i<N) A[i,i+1] = 1; end
end
A

Out[175]:
100×100 sparse matrix with 298 Float64 nonzero entries:
[1  ,   1]  =  -2.0
[2  ,   1]  =  1.0
[1  ,   2]  =  1.0
[2  ,   2]  =  -2.0
[3  ,   2]  =  1.0
[2  ,   3]  =  1.0
[3  ,   3]  =  -2.0
[4  ,   3]  =  1.0
[3  ,   4]  =  1.0
[4  ,   4]  =  -2.0
⋮
[96 ,  97]  =  1.0
[97 ,  97]  =  -2.0
[98 ,  97]  =  1.0
[97 ,  98]  =  1.0
[98 ,  98]  =  -2.0
[99 ,  98]  =  1.0
[98 ,  99]  =  1.0
[99 ,  99]  =  -2.0
[100,  99]  =  1.0
[99 , 100]  =  1.0
[100, 100]  =  -2.0
In [176]:
nnz(A)

Out[176]:
298
In [177]:
nnz(A)/(N*N)

Out[177]:
0.0298
In [ ]:


In [178]:
spy(A)

Out[178]:
PyObject <matplotlib.image.AxesImage object at 0x7fb20e6e0990>
In [179]:
# SOLVE Laplaces equation for point charge:
#     A * potential = chargeDensity
chargeDensity = zeros(N); chargeDensity[div(N,2)] = 1.;  # point charge at N/2 position
plot(chargeDensity)

Out[179]:
1-element Array{Any,1}:
PyObject <matplotlib.lines.Line2D object at 0x7fb20e610310>
In [185]:
plot(-A \ chargeDensity) # electric potential = solution to 1D Laplace's equation (-gradient --> electric field)
# Matrix A is actually only positive semidefinite (!)

Out[185]:
1-element Array{Any,1}:
PyObject <matplotlib.lines.Line2D object at 0x7fb20e4851d0>

## Visualizing "Fill In"¶

In [212]:
N = 15;
L = speye(N, N)
L[1:N,1] = ones(N);

In [213]:
spy(L)

Out[213]:
PyObject <matplotlib.image.AxesImage object at 0x7fb20db9ae90>
In [218]:
nnz(L)

Out[218]:
29
In [219]:
nnz(L*L')

Out[219]:
225
In [220]:
N*N

Out[220]:
225
In [241]:
## Random sparse lower triangular matrix
N = 225
L = tril(sprand(N,N,0.1) + speye(N))
print("nnz(L)/(N*N) = ",(nnz(L)/(N*N)))
spy(L)

nnz(L)/(N*N) = 0.052997530864197534
Out[241]:
PyObject <matplotlib.image.AxesImage object at 0x7fb20ca34b10>
In [242]:
print("nnz(L*L')/(N*N) = ",(nnz(L*L')/(N*N)))
spy(L*L')

nnz(L*L')/(N*N) = 0.5019456790123457
Out[242]:
PyObject <matplotlib.image.AxesImage object at 0x7fb20c966810>
In [ ]: