21) 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.77458 -0.293769
-0.293769 0.0297292
@show svdvals(A) # Julia built-in
U, S, V = svd(A) # Julia built-in
@show U - U' # U is symmetric
Aplot(A)
svdvals(A) = [1.8227160219614589, 0.018403077457637293]
U - U' =
[0.0 -2.7755575615628914e-17; 2.7755575615628914e-17 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.57956085979877
0.091929 seconds (20.78 k allocations: 20.649 MiB, 7.00% gc time, 16.07% compilation time)
0.281116 seconds (748.90 k allocations: 73.596 MiB, 8.64% gc time, 87.07% compilation time)
0.276054 seconds (56.96 k allocations: 13.407 MiB, 50.88% gc time, 12.04% 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.47870606451107356 0.8779752296052096; 0.8779752296052096 0.47870606451107345], [1.8757788549148868, 0.12067438349951641], [0.8576906491355568 -0.5141660727677657; 0.5141660727677657 0.8576906491355568])
SVD{Float64, Float64, Matrix{Float64}, Vector{Float64}}
U factor:
2×2 Matrix{Float64}:
-0.478706 0.877975
0.877975 0.478706
singular values:
2-element Vector{Float64}:
1.8757788549148868
0.12067438349951641
Vt factor:
2×2 Matrix{Float64}:
0.857691 -0.514166
0.514166 0.857691
@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.47870606451107356; 0.8779752296052096;;]
Shat = S[1:k] = [1.8757788549148868]
Vhat = V[:, 1:k] = [0.8576906491355568; -0.5141660727677657;;]
Ahat = Uhat * diagm(Shat) * Vhat' = [-0.7701604996161654 0.4616937352501302; 1.4125198981424731 -0.8467736117517749]
norm(Ahat) = 1.8757788549148877
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.0104259 0.0150896 0.00809692 … 0.0103273 0.0113001
0.000374816 -0.000139851 0.000601325 0.000282173 -0.000193855
-0.00236521 -0.00691677 -0.000592256 -0.000829134 -0.000141851
0.00995598 0.0173012 0.00705352 0.0088178 0.00908384
0.000680328 0.000629145 0.000708303 0.000520267 0.000164425
0.000717091 -0.000866648 0.00125079 … 0.0017686 0.00250321
0.00939907 0.0172575 0.00641756 0.00779966 0.00770018
0.000760255 0.000894178 0.000718442 0.000552615 0.000225544
0.00268441 0.00393771 0.00233023 0.00319644 0.00381307
0.0127977 0.0248749 0.00854716 0.010859 0.0109872
0.000695526 0.000271012 0.000790334 … 0.000481064 -3.39784e-5
0.00168555 0.000198254 0.00205494 0.00299105 0.00400438
0.0104549 0.0179927 0.00744119 0.00923621 0.00940944
⋮ ⋱
0.00912179 0.0155053 0.00653848 0.00800563 0.00811096
-0.000105633 -0.00101553 0.000196483 -0.000142754 -0.00048588
0.00101799 0.000118701 0.00135925 … 0.00187113 0.0024891
0.00769977 0.0126248 0.00562202 0.00688829 0.00710804
0.00101249 0.00152705 0.000857893 0.000834576 0.000617544
0.00075603 8.59404e-7 0.00106892 0.00146226 0.00198641
0.0086289 0.013894 0.00637882 0.00778118 0.00796141
-0.000444365 -0.00158955 -3.76035e-5 … -0.000437582 -0.000782425
0.000350424 -0.00109739 0.000914164 0.00125376 0.00182505
0.0111115 0.0204468 0.00761576 0.00938789 0.00932516
-0.000450175 -0.00231049 7.32083e-5 -0.000435491 -0.000882865
0.00192783 0.00151591 0.00205979 0.00288493 0.00365784
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.379612303943579
norm(F - Fhat) = 0.004660806587794126
0.004660806587794126