22) Low Rank#
Last time#
Condition number and SVD
SVD for solving systems
Costs of decompositions
Today#
Reflection on algorithm choices
Low-rank structure
Primer on interpolation
using LinearAlgebra
using Plots
default(linewidth=4, legendfontsize=12)
function vander(x, k=nothing)
if isnothing(k)
k = length(x)
end
m = length(x)
V = ones(m, k)
for j in 2:k
V[:, j] = V[:, j-1] .* x
end
V
end
function gram_schmidt_classical(A)
m, n = size(A)
Q = zeros(m, n)
R = zeros(n, n)
for j in 1:n
v = A[:,j]
R[1:j-1,j] = Q[:,1:j-1]' * v
v -= Q[:,1:j-1] * R[1:j-1,j]
R[j,j] = norm(v)
Q[:,j] = v / R[j,j]
end
Q, R
end
function qr_householder(A)
m, n = size(A)
R = copy(A)
V = [] # list of reflectors
for j in 1:n
v = copy(R[j:end, j])
v[1] += sign(v[1]) * norm(v) # <--- here we pick the sign of v so that moves it the largest distance
v = normalize(v)
R[j:end,j:end] -= 2 * v * v' * R[j:end,j:end]
push!(V, v)
end
V, R
end
function qr_chol(A)
R = cholesky(A' * A).U
Q = A / R
Q, R
end
function qr_chol2(A)
Q, R = qr_chol(A)
Q, R1 = qr_chol(Q)
Q, R1 * R
end
function peanut()
theta = LinRange(0, 2*pi, 50)
r = 1 .+ .4*sin.(3*theta) + .6*sin.(2*theta)
r' .* [cos.(theta) sin.(theta)]'
end
function circle()
theta = LinRange(0, 2*pi, 50)
[cos.(theta) sin.(theta)]'
end
function Aplot(A)
"Plot a transformation from X to Y"
X = peanut()
Y = A * X
p = scatter(X[1,:], X[2,:], label="in")
scatter!(p, Y[1,:], Y[2,:], label="out")
X = circle()
Y = A * X
q = scatter(X[1,:], X[2,:], label="in")
scatter!(q, Y[1,:], Y[2,:], label="out")
plot(p, q, layout=2, aspect_ratio=:equal)
end
Aplot (generic function with 1 method)
1. Reflection on algorithm choices#
Recap on Condition number via SVD#
\[\begin{split} U \overbrace{\begin{bmatrix} \sigma_{\max} && \\ & \ddots & \\ && \sigma_{\min} \end{bmatrix}}^{\Sigma} V^T = A \end{split}\]
(8)#\[\begin{align}
\lVert A \rVert &= \sigma_{\max} & \textrm{, and } \;
\kappa(A) &= \frac{\sigma_{\max}}{\sigma_{\min}} = \texttt{cond}(A)
\end{align}\]
A = randn(2,2) # nonsymmetric
A = A + A' # make it symmetric
2×2 Matrix{Float64}:
-1.31938 -0.0909455
-0.0909455 -2.10179
@show svdvals(A) # Julia built-in
U, S, V = svd(A) # Julia built-in
@show U - U' # U is symmetric
Aplot(A)
svdvals(A) = [2.112222153284016, 1.308949301207618]
U - U' =
[0.0 -2.220446049250313e-16; 2.220446049250313e-16 0.0]
Real-world example: autonomous vehicles#
Need to solve least squares problems in real time
Weight/cost/size increase with compute
What algorithm to choose?
What precision to use?
Factors to consider:
How many objects?
Speed (of robot and objects)
Aerial, wheeled, walking
Fog, light – longer memory?
Tolerences (how accurate does the solution need to be?)
Consequences of being wrong, who bears those consequences?
A = rand(5000, 500)
A_32 = Float32.(A)
@show cond(A)
@time qr(A); # Householder; backward stable
@time qr_chol(A); # Unstable
@time qr(A_32); # Julia built-in; best in terms of memory allocations; Run twice!
cond(A) = 56.36412841584856
0.084675 seconds (20.78 k allocations: 20.649 MiB, 24.21% compilation time)
0.812533 seconds (790.06 k allocations: 76.448 MiB, 21.22% gc time, 76.82% compilation time)
0.120157 seconds (56.96 k allocations: 13.393 MiB, 48.07% compilation time)
V = vander(LinRange(-1, 1, 20))
@show cond(V)
Q, R = qr(Float32.(V)) # Julia built-in, but with single-precision Floats
@show norm(Q' * Q - I)
Q, R = qr_chol(V) # Unstable; really not orthogonal
@show norm(Q' * Q - I)
cond(V) = 2.7224082312417406e8
norm(Q' * Q - I) =
1.6641898f-6
norm(Q' * Q - I) = 0.1749736012761826
0.1749736012761826
2. Low-rank structure#
Best low rank approximation#
The SVD can be truncated to yield the best rank-\(k\) approximation of a matrix.
n, k = 2, 1
A = randn(n, n)
Aplot(A)
@show U, S, V = svd(A)
(U, S, V) = svd(A) = SVD{Float64, Float64, Matrix{Float64}, Vector{Float64}}([-0.966377620054708 0.2571270025831574; 0.2571270025831574 0.9663776200547077], [1.6540697594694416, 0.36740419488930603], [-0.6307829286820065 0.7759593397101104; -0.7759593397101104 -0.6307829286820065])
SVD{Float64, Float64, Matrix{Float64}, Vector{Float64}}
U factor:
2×2 Matrix{Float64}:
-0.966378 0.257127
0.257127 0.966378
singular values:
2-element Vector{Float64}:
1.6540697594694416
0.36740419488930603
Vt factor:
2×2 Matrix{Float64}:
-0.630783 0.775959
-0.775959 -0.630783
@show Uhat = U[:, 1:k] # Uhat is now 2x1
@show Shat = S[1:k] # truncate to first k singular values, in this case 1
@show Vhat = V[:, 1:k] # Vhat is now 2x1
@show Ahat = Uhat * diagm(Shat) * Vhat'
@show norm(Ahat)
Aplot(Ahat - A) # we have squished every point onto a line
Uhat = U[:, 1:k] = [-0.966377620054708; 0.2571270025831574;;]
Shat = S[1:k] = [1.6540697594694416]
Vhat = V[:, 1:k] = [-0.6307829286820065; 0.7759593397101104;;]
Ahat = Uhat * diagm(Shat) * Vhat' =
[1.008278755510557 -1.2403368604227443; -0.2682757638344614 0.3300201624038536]
norm(Ahat) = 1.654069759469442
Example: Galaxies#
Suppose we have two galaxies of size \(n_1 = 100\) and \(n_2 = 200\), each randomly distributed around their respective centers.
galaxy(center, sigma, n) = reshape(center, 1, 3) .+ sigma*randn(n, 3)
g1 = galaxy([0 0 0], 1, 100)
g2 = galaxy([10 0 0], 1, 100)
scatter(g1[:,1], g1[:,2], aspect_ratio=:equal)
scatter!(g2[:,1], g2[:,2])
Forces between stars#
Consider Newton’s law of universal gravitation between two bodies with spatial extent (i.e., not point masses), then we can write the gravitational force from a star at position \(x_2\) acting on a star at position \(x_1\),
\[ F_{1,2} = G \frac{m_1 m_2}{\lVert \mathbf{x}_2 - \mathbf{x}_1 \rVert^3} (\mathbf{x}_2 - \mathbf{x}_1) \]
where \(m_1\) and \(m_2\) are the masses of each star, respectively.
function gravitational_force(g1, g2)
m = size(g1, 1)
n = size(g2, 1)
F = zeros(3*m, n)
for i in 0:m-1
for j in 1:n
r = g2[j,:] - g1[1+i,:]
F[1+3*i:3*(i+1),j] = r / norm(r)^3
end
end
F
end
# Let's apply it to our two galaxies
gravitational_force(g1, g2)
300×100 Matrix{Float64}:
0.00809386 0.00897236 0.00844164 … 0.00943165 0.00775605
-0.00181779 -0.00180624 -0.0013674 -0.000702235 -0.00200539
0.000982464 -0.000237697 -0.00093514 -0.00135195 0.000931452
0.0122178 0.0130191 0.0115109 0.0123042 0.0117607
-0.00118271 -0.000830598 -0.000259041 0.00110427 -0.00162736
0.000493909 -0.00180677 -0.00269537 … -0.00344481 0.00046878
0.00992049 0.01032 0.00919458 0.00974749 0.00959912
-0.000844702 -0.000569455 -0.000168943 0.000800277 -0.00117753
-2.80498e-5 -0.00170067 -0.00229494 -0.00283945 -2.7614e-5
0.0145824 0.0162305 0.0143837 0.0156672 0.0139343
-0.00183898 -0.00147391 -0.000633043 … 0.00124611 -0.00239165
0.00156283 -0.00144075 -0.00282418 -0.00383751 0.00147254
0.00927423 0.0103313 0.00958244 0.0104896 0.00894124
⋮ ⋱
0.0105814 0.0112595 0.0100109 0.0102884 0.0103802
0.000449664 0.000882565 0.00111823 0.00228796 2.47183e-5
0.000717063 -0.00109081 -0.00186992 … -0.00233153 0.000694871
0.007004 0.00700916 0.0063048 0.00663281 0.00680573
-0.00070345 -0.000530222 -0.000276511 0.000262548 -0.000899509
-0.000686559 -0.00166534 -0.00193805 -0.0022879 -0.000661561
0.00883871 0.00980147 0.0090399 0.00948996 0.0086603
0.000104009 0.000429515 0.000685197 … 0.0016719 -0.000217438
0.00138061 7.12902e-5 -0.000717486 -0.00105272 0.00133833
0.00771815 0.00864538 0.008086 0.00834497 0.00763244
0.000668245 0.00103405 0.00119456 0.00205664 0.000386986
0.00161927 0.000604731 -0.000117539 -0.000346841 0.00158522
Spectrum#
g1 = galaxy([0 0 0], 1, 500)
g2 = galaxy([10 0 0], 1, 500)
F = gravitational_force(g1, g2)
@show size(F)
U, S, V = svd(F) # U is 1500 x 500, S is 500 long, and V is 500 x 500
scatter(S, yscale=:log10, ylims=(1e-10, 10), xlims=(0, 200))
size(F) = (1500, 500)
k = 10 # let's truncate at the first 10 singular values
Uhat = U[:,1:k] # Uhat is now 1500 x 10
Shat = S[1:k] # Shat is now 10 long
Vhat = V[:,1:k] # Vhat is now 500 x 10
Fhat = Uhat * diagm(Shat) * Vhat' # Fhat is still 1500 x 500
@show norm(F)
@show norm(F - Fhat) # Fhat is the best rank-10 approximation of F and it is not too far off from F indeed
norm(F) = 5.451123280381913
norm(F - Fhat) = 0.005749113583614124
0.005749113583614124