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 44 stored entries:
   ⋅     ⋅     ⋅     ⋅    ⋅     ⋅   -1.0  -4.0  -3.0   1.0
   ⋅     ⋅     ⋅     ⋅    ⋅     ⋅   -5.0   2.0  -2.0  -1.0
   ⋅     ⋅     ⋅     ⋅    ⋅     ⋅    4.0  -2.0  -1.0  -5.0
   ⋅     ⋅     ⋅     ⋅    ⋅     ⋅   -5.0  -5.0   4.0    ⋅ 
   ⋅     ⋅     ⋅     ⋅    ⋅     ⋅    2.0    ⋅    3.0   1.0
   ⋅     ⋅     ⋅     ⋅    ⋅     ⋅   -4.0   4.0  -4.0  -5.0
 -1.0  -5.0   4.0  -5.0  2.0  -4.0    ⋅     ⋅     ⋅     ⋅ 
 -4.0   2.0  -2.0  -5.0   ⋅    4.0    ⋅     ⋅     ⋅     ⋅ 
 -3.0  -2.0  -1.0   4.0  3.0  -4.0    ⋅     ⋅     ⋅     ⋅ 
  1.0  -1.0  -5.0    ⋅   1.0  -5.0    ⋅     ⋅     ⋅     ⋅
In [2]:
##
eigvals(Matrix(B))
Out[2]:
10-element Vector{Float64}:
 -10.534200145001812
  -8.971604312666253
  -6.599363520399757
  -4.9989343351967594
  -5.412163639425414e-15
   8.309312122466312e-17
   4.998934335196764
   6.5993635203997565
   8.971604312666262
  10.534200145001808
In [3]:
##
sort(abs.(eigvals(Matrix(B))))
Out[3]:
10-element Vector{Float64}:
  8.309312122466312e-17
  5.412163639425414e-15
  4.9989343351967594
  4.998934335196764
  6.5993635203997565
  6.599363520399757
  8.971604312666253
  8.971604312666262
 10.534200145001808
 10.534200145001812
In [4]:
## Final question:
# What are the singular values of A?
svdvals(Matrix(A))
Out[4]:
4-element Vector{Float64}:
 10.534200145001808
  8.971604312666255
  6.5993635203997565
  4.998934335196757
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.09228859067059805; 0.4082482904638631 0.0 … 0.0 0.4790594329372884; … ; 0.0 -0.2142857142857143 … -0.4028063848238582 0.0; 0.0 -0.6428571428571429 … 0.24893310450716952 0.0], [0.0 5.715476066494084 0.0 0.0; 5.715476066494082 0.0 6.149221509269934 0.0; 0.0 6.149221509269936 0.0 4.774348431569532; 0.0 0.0 4.774348431569531 0.0], 7.495668645088698)
In [6]:
##
U,T,rho = lanczos(B,[ones(6); ones(4)],4)
Out[6]:
([0.31622776601683794 -0.13431339876147755 … -0.00023577165198764634 0.4456257987085685; 0.31622776601683794 -0.059694843893990046 … -0.5324328778917775 -0.031330545924039714; … ; 0.31622776601683794 0.1641608207084726 … 0.1694500436754723 -0.06359966575289865; 0.31622776601683794 -0.28355050849645264 … 0.2679678622640917 0.48968110564002854], [-5.2 4.237924020083418 0.0 0.0; 4.237924020083419 -1.9492204899777275 4.658971809672365 0.0; 0.0 4.658971809672366 0.025168711565927995 7.25588376684862; 0.0 0.0 7.255883766848619 0.1306746216798963], 4.812396112839309)
In [7]:
##
U,T,rho = lanczos(B,[zeros(6); ones(4)],4)
Out[7]:
([0.0 -0.439219063596711 … 0.4786751140636702 0.0; 0.0 -0.37647348308289513 … -0.41016279017397617 0.0; … ; 0.5 0.0 … 0.0 0.8170655276747744; 0.5 0.0 … 0.0 -0.019252780765267714], [0.0 7.968688725254614 0.0 0.0; 7.968688725254614 0.0 3.650585111121684 0.0; 0.0 3.650585111121683 0.0 8.4664317253694; 0.0 0.0 8.466431725369402 0.0], 1.9820581278267)
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.43921906359671103 0.4786751140636701 -0.21289678015603625; -0.3764734830828952 -0.41016279017397617 0.12413082900853112; … ; 0.3764734830828952 0.016320831239238928 0.18263091390769756; -0.5647102246243427 -0.508663242765739 -0.19076767337858463], [0.5 0.6101674226014826 -0.4819143982255567 -0.3813846210569654; 0.5 -0.7820455698131674 -0.3158983486839499 -0.19650180669487746; 0.5 0.06015735152408987 0.8170655276747745 -0.28068668750610576; 0.5 0.11172079568759545 -0.019252780765267707 0.8585731152579458], 3×3 Bidiagonal{Float64, Vector{Float64}}:
 diag: 7.968688725254614  8.466431725369402  7.861572495459137
 sub: 3.6505851111216825  1.9820581278267004, 3.3473933780954677)