In [1]:
## Reduction to tridiagonal form!
# We can exactly construct a sequence of orthogonal transformations
# to produce a tridiagonal matrix from an arbitrary symmetric matrix.
In [2]:
## How Householder Matrices can be used to zero out
using LinearAlgebra
using SparseArrays
function housevec(x)
  v = copy(x)
  normx = norm(x)
  s = -sign(x[1])*normx
  v[1] = v[1] -s

  y = v*(v'*x)
  beta = -1/(s*v[1])
  return v,beta,s
end

""" Return the householder matrix to zero out below the ith row and jth column """
function housemat(A,i,j)
  v,beta,s = housevec(A[i:end,j])
  H = Matrix(1.0I,size(A)...)
  H[i:end,i:end] = H[i:end,i:end] - beta*v*v'
  return H,v,beta,s
end

A = randn(5,5)
A = A + A'

H = housemat(A,2,1)[1]
display(H*A*H')

5×5 Matrix{Float64}:
  1.22616       2.30041   -1.635e-16  -9.32004e-19  -2.02952e-16
  2.30041      -0.324814  -1.24608    -3.31675       1.66655
 -1.635e-16    -1.24608    1.36951     0.743131     -1.94546
 -9.32004e-19  -3.31675    0.743131   -0.118271      0.545249
 -2.02952e-16   1.66655   -1.94546     0.545249      1.19412
In [3]:
##
B = H*A*H'
H2 = housemat(B,3,2)[1]
display(H2*B*H2')

5×5 Matrix{Float64}:
  1.22616       2.30041      -3.356e-17  8.25656e-17  -2.44906e-16
  2.30041      -0.324814      3.91547    4.96775e-17   9.92785e-18
 -3.356e-17     3.91547       0.804705   1.82952      -0.205638
  8.25656e-17  -1.11637e-16   1.82952    1.51115       1.19528
 -2.44906e-16  -2.0587e-17   -0.205638   1.19528       0.129505
In [4]:
##
Q = H2*H
display(Q*A*Q')

5×5 Matrix{Float64}:
  1.22616       2.30041       2.7223e-17   8.70105e-17  -9.99039e-17
  2.30041      -0.324814      3.91547     -3.29141e-16   1.93421e-16
  2.7223e-17    3.91547       0.804705     1.82952      -0.205638
  8.70105e-17  -1.58227e-16   1.82952      1.51115       1.19528
 -9.99039e-17   2.95691e-16  -0.205638     1.19528       0.129505
In [5]:
##
spones(A::SparseMatrixCSC) = begin fill!(A.nzval, one(eltype(A))); A end
using Plots
pyplot()
using Printf
function animate_myhess(A)
  n = size(A,1)
  V = zeros(n,n) # space to store the householder vectors
  for k = 1:n-2
    x = A[k+1:n,k]
    vk = copy(x)
    vk[1] = x[1] + sign(x[1])*norm(x)
    vk = vk/norm(vk)
    V[k+1:n,k] = vk # store V
    A[k+1:n,k:n] -= 2*vk*(vk'*A[k+1:n,k:n])
    A[1:n,k+1:n] -= 2*(A[1:n,k+1:n]*vk)*vk'
    At = A
    At[abs.(At).<1e-12] .= 0
    spy(-spones(dropzeros!(sparse(At))),markersize = 6.0,
      colorbar=false,axis=false,markercolor=:blue)

    title!(@sprintf("Step %i",k))
    gui()
    sleep(1.5)
  end
  return A,V
end
A = randn(37,37)
A = A + A'
animate_myhess(A)

ArgumentError: Package PyPlot not found in current path.
- Run `import Pkg; Pkg.add("PyPlot")` to install the PyPlot package.

Stacktrace:
  [1] macro expansion
    @ ./loading.jl:1772 [inlined]
  [2] macro expansion
    @ ./lock.jl:267 [inlined]
  [3] __require(into::Module, mod::Symbol)
    @ Base ./loading.jl:1753
  [4] #invoke_in_world#3
    @ ./essentials.jl:926 [inlined]
  [5] invoke_in_world
    @ ./essentials.jl:923 [inlined]
  [6] require(into::Module, mod::Symbol)
    @ Base ./loading.jl:1746
  [7] top-level scope
    @ ~/.julia/packages/Plots/du2dt/src/backends.jl:883
  [8] eval
    @ ./boot.jl:385 [inlined]
  [9] _initialize_backend(pkg::Plots.PyPlotBackend)
    @ Plots ~/.julia/packages/Plots/du2dt/src/backends.jl:882
 [10] backend
    @ ~/.julia/packages/Plots/du2dt/src/backends.jl:245 [inlined]
 [11] pyplot()
    @ Plots ~/.julia/packages/Plots/du2dt/src/backends.jl:86
 [12] top-level scope
    @ In[5]:4
In [6]:
##

function myhess(A)
  # myhess Implement algorithm 26.1 from Trefethen and Bau
  n = size(A,1)
  V = zeros(n,n) # space to store the householder vectors
  for k = 1:n-2
    x = A[k+1:n,k]
    vk = copy(x)
    vk[1] = x[1] + sign(x[1])*norm(x)
    vk = vk/norm(vk)
    V[k+1:n,k] = vk # store V
    A[k+1:n,k:n] -= 2*vk*(vk'*A[k+1:n,k:n])
    A[1:n,k+1:n] -= 2*(A[1:n,k+1:n]*vk)*vk'
  end
  return A,V
end
Out[6]:
myhess (generic function with 1 method)
In [7]:
## show that RQ given QR preserves tridiagonal
A = randn(5,5)
A = A + A'
T = myhess(A)[1]
T = map!(x -> abs(x) > 10*eps(1.0) ? x : 0.0, T, T)
Q,R = qr(T)
display(R*Q)
Q

5×5 Matrix{Float64}:
 -0.679202  5.13236    0.0       -4.44089e-16  -8.29804e-17
  5.13236   1.82358    0.870615   0.0          -1.24471e-16
  0.0       0.870615  -0.671838   1.27967       1.11022e-16
  0.0       0.0        1.27967    1.59845      -0.597202
  0.0       0.0        0.0       -0.597202      1.54224
Out[7]:
5×5 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}