Stanford CS205a

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 [ ]: