In [1]:
## 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)]
Out[1]:
10×10 SparseMatrixCSC{Float64,Int64} with 42 stored entries:
  [7 ,  1]  =  -5.0
  [8 ,  1]  =  -5.0
  [9 ,  1]  =  5.0
  [10,  1]  =  -3.0
  [7 ,  2]  =  -3.0
  [8 ,  2]  =  -2.0
  [9 ,  2]  =  -1.0
  [10,  2]  =  1.0
  [7 ,  3]  =  4.0
  [8 ,  3]  =  3.0
  [10,  3]  =  -4.0
  [8 ,  4]  =  2.0
  ⋮
  [4 ,  8]  =  2.0
  [6 ,  8]  =  2.0
  [1 ,  9]  =  5.0
  [2 ,  9]  =  -1.0
  [4 ,  9]  =  -4.0
  [5 ,  9]  =  -3.0
  [6 ,  9]  =  2.0
  [1 , 10]  =  -3.0
  [2 , 10]  =  1.0
  [3 , 10]  =  -4.0
  [4 , 10]  =  3.0
  [5 , 10]  =  -1.0
  [6 , 10]  =  -2.0
In [2]:
##
eigvals(Matrix(B))
Out[2]:
10-element Array{Float64,1}:
 -11.588779551813595     
  -9.213213673152731     
  -4.6189060689107375    
  -0.6946862881029516    
  -1.8209040003065487e-15
   3.597260839706799e-15 
   0.6946862881029459    
   4.618906068910739     
   9.21321367315273      
  11.588779551813591
In [3]:
##
sort(abs.(eigvals(Matrix(B))))
Out[3]:
10-element Array{Float64,1}:
  1.8209040003065487e-15
  3.597260839706799e-15 
  0.6946862881029459    
  0.6946862881029516    
  4.6189060689107375    
  4.618906068910739     
  9.21321367315273      
  9.213213673152731     
 11.588779551813591     
 11.588779551813595
In [4]:
## Final question:
# What are the singular values of A?
svdvals(Matrix(A))
Out[4]:
4-element Array{Float64,1}:
 11.588779551813593 
  9.213213673152733 
  4.6189060689107375
  0.6946862881029496
In [5]:
##
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)
Out[5]:
([0.4082482904638631 0.0 … 0.0 0.4469307559354324; 0.4082482904638631 0.0 … 0.0 0.41026139965122344; … ; 0.0 -0.1373605639486891 … 0.4735401217799603 0.0; 0.0 -0.8241633836921342 … 0.28134549807126485 0.0], [0.0 2.9720924166878353 0.0 0.0; 2.972092416687835 0.0 5.941274454051486 0.0; 0.0 5.941274454051484 0.0 7.378741829546498; 0.0 0.0 7.3787418295465 0.0], 7.346257794251449)
In [6]:
##
U,T,rho = lanczos(B,[ones(6); ones(4)],4)
Out[6]:
([0.31622776601683794 -0.3794823438779049 … -0.4829066931654384 -0.41134130117917816; 0.31622776601683794 -0.18319837290657473 … -0.25972654131998496 0.28696243473928534; … ; 0.31622776601683794 0.07851358838853208 … 0.08639796665289139 -0.4031557832443077; 0.31622776601683794 -0.24862636323035145 … 0.4839588588409059 0.4861879231738003], [-2.2 4.833218389437828 0.0 0.0; 4.833218389437829 -0.4523972602739724 9.254343356795872 0.0; 0.0 9.254343356795875 1.2463575495887615 4.957744664872655; 0.0 0.0 4.957744664872649 -2.7665807783838314], 2.820155990924034)
In [7]:
##
U,T,rho = lanczos(B,[zeros(6); ones(4)],4)
Out[7]:
([0.0 -0.5286548803640718 … -0.2847922219353602 0.0; 0.0 -0.3304093002275449 … 0.012720836536179595 0.0; … ; 0.5 0.0 … 0.0 -0.2667672520449107; 0.5 0.0 … 0.0 -0.15877378758742963], [0.0 7.566372975210778 0.0 0.0; 7.566372975210778 0.0 7.861450358981152 0.0; 0.0 7.861450358981153 0.0 4.798667706348427; 0.0 0.0 4.798667706348427 0.0], 3.8403595754822364)
In [8]:
## 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)
Out[8]:
([-0.5286548803640718 -0.2847922219353601 -0.5543757648117681; -0.3304093002275449 0.01272083653617969 -0.1075932789542194; … ; -0.5947367404095808 0.5701275463325776 0.2619479554450706; 0.46257302031856284 0.07993550318886168 -0.2773741139957646], [0.5 0.7544214703726846 -0.42483575736664697 0.01904794748586141; 0.5 0.14920313202356708 0.850376796998987 -0.06780655219152136; 0.5 -0.4644210165803992 -0.2667672520449106 -0.6805500368046685; 0.5 -0.4392035858158526 -0.15877378758742974 0.729308641510327], 3×3 Bidiagonal{Float64,Array{Float64,1}}:
 diag: 7.566372975210778  4.798667706348426  3.411973753459233
 sub: 7.861450358981152  3.8403595754822364, 8.054989720271308)