## Lecture 29 - Bipartite Matrices
using LinearAlgebra
using SparseArrays
using Random
Random.seed!(0)
A = rand(-5.0:5, 6,4)
B = [spzeros(6,6) A; A' spzeros(4,4)]
##
eigvals(Matrix(B))
##
sort(abs.(eigvals(Matrix(B))))
## Final question:
# What are the singular values of A?
svdvals(Matrix(A))
##
function lanczos(A,b,k)
n = size(A,1)
V = zeros(n,k+1)
T = Tridiagonal(zeros(k-1), zeros(k), zeros(k-1))
rho = 0.0
V[:,1] = b/norm(b)
for j=1:k
y = A*V[:,j]
for i=max(j-1,1):j
T[i,j] = V[:,i]'*y
y -= T[i,j]*V[:,i]
end
rho = norm(y)
V[:,j+1] = y/rho
if j < k
T[j+1,j] = rho
end
end
return V,T,rho
end
U,T,rho = lanczos(B,[ones(6); zeros(4)],4)
##
U,T,rho = lanczos(B,[ones(6); ones(4)],4)
##
U,T,rho = lanczos(B,[zeros(6); ones(4)],4)
## The Golub-Kahan process
# This "simplifies" Lanczos on a bipartite matrix
function golubkahan(A,v,k)
B = Bidiagonal(zeros(k),zeros(k-1),:L)
m,n = size(A)
V = zeros(n,k+1)
v = v/norm(v)
V[:,1] = v
U = zeros(m,k)
u = zeros(m)
beta = 0.0
for i=1:k
if i>1
B[i,i-1] = beta
end
Av = A*v - beta*u
B[i,i] = beta = norm(Av)
u = U[:,i] = normalize!(Av)
v = A'*Av - beta*v
beta = norm(v)
v = V[:,i+1] = normalize!(v)
end
return U,V,B,beta
end
U,V,BD,beta = golubkahan(A,ones(4),3)