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.94196 0.937797
0.937797 -1.03437
@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.2128001268646975, 1.3052093673307865]
U - U' =
[0.0 0.0; 0.0 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.24692410056267
0.078043 seconds (20.78 k allocations: 20.649 MiB, 25.59% compilation time)
0.555163 seconds (790.06 k allocations: 76.448 MiB, 3.66% gc time, 91.29% compilation time)
0.098542 seconds (56.96 k allocations: 13.393 MiB, 46.42% 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.7707367195355678 0.6371537562939981; 0.6371537562939981 0.7707367195355677], [1.5796579122907455, 0.23586033139977886], [0.4752532682751691 0.8798490387525408; 0.8798490387525408 -0.4752532682751691])
SVD{Float64, Float64, Matrix{Float64}, Vector{Float64}}
U factor:
2×2 Matrix{Float64}:
-0.770737 0.637154
0.637154 0.770737
singular values:
2-element Vector{Float64}:
1.5796579122907455
0.23586033139977886
Vt factor:
2×2 Matrix{Float64}:
0.475253 0.879849
0.879849 -0.475253
@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.7707367195355678; 0.6371537562939981;;]
Shat = S[1:k] = [1.5796579122907455]
Vhat = V[:, 1:k] = [0.4752532682751691; 0.8798490387525408;;]
Ahat = Uhat * diagm(Shat) * Vhat' =
[-0.5786210239365152 -1.071216519057767; 0.4783352726388646 0.8855548355515195]
norm(Ahat) = 1.5796579122907453
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.00779732 0.00598369 0.0055781 … 0.00606115 0.00702459
-2.33319e-5 0.000585375 -5.20352e-5 0.000735937 0.000317306
-0.000329506 0.000157679 -0.000602207 -0.00048417 -0.00070004
0.0123498 0.00845431 0.00851145 0.00904661 0.0111406
0.000542236 0.0013661 0.000235249 0.00173673 0.00113673
0.00165815 0.0016333 0.000185053 … 0.000583738 0.000581704
0.0120206 0.00882434 0.00780622 0.00850269 0.0101807
0.000546239 0.00143767 0.000225543 0.00164393 0.00105369
-0.00230512 -0.000751094 -0.00191406 -0.00186066 -0.00260529
0.0169918 0.0109691 0.0108046 0.0114855 0.0146648
0.00122855 0.00220795 0.00051977 … 0.00271882 0.00201674
0.00101539 0.00152574 -0.000563136 -9.8384e-5 -0.000441728
0.00876278 0.00632551 0.00621327 0.00633521 0.00765724
⋮ ⋱
0.00925116 0.00729166 0.00641937 0.00740612 0.00840886
-0.00140195 -0.000158721 -0.000863168 1.63836e-5 -0.000766602
-0.000633065 7.40801e-5 -0.000872881 … -0.000790136 -0.00109811
0.0113108 0.00834224 0.0074734 0.00818074 0.00974111
0.000247049 0.0011543 7.47429e-5 0.00136811 0.000767714
-0.00162444 -0.000397471 -0.00152163 -0.00143779 -0.00202345
0.00771857 0.00576913 0.00574163 0.00614775 0.00718699
-2.44457e-5 0.000568106 -5.46935e-5 … 0.000751643 0.000326419
0.00111072 0.00110094 0.000286817 0.00052969 0.000559298
0.0107101 0.00763113 0.00746111 0.00803095 0.0096708
0.000191761 0.00100158 5.29373e-5 0.00127741 0.000706712
0.00088057 0.00110119 -0.000106163 0.000189357 8.56869e-5
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.309830732097176
norm(F - Fhat) = 0.00482056930385322
0.00482056930385322