## Reduction to tridiagonal form!
# We can exactly construct a sequence of orthogonal transformations
# to produce a tridiagonal matrix from an arbitrary symmetric matrix.

## 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 Array{Float64,2}:
 -0.786842      2.58457  -2.115e-16   5.14036e-17  -7.79959e-17
  2.58457      -1.56744   1.01906    -2.91169       1.95535    
 -2.115e-16     1.01906   1.84522    -2.75286      -0.32045    
  5.14036e-17  -2.91169  -2.75286     0.963228     -0.359506   
 -7.79959e-17   1.95535  -0.32045    -0.359506      0.736494   

##
B = H*A*H'
H2 = housemat(B,3,2)[1]
display(H2*B*H2')

5×5 Array{Float64,2}:
 -0.786842      2.58457       1.41747e-16  -1.68774e-16   6.98651e-17
  2.58457      -1.56744      -3.65238      -3.0166e-16   -3.36155e-16
  1.41747e-16  -3.65238       2.40268      -1.351         1.47652    
 -1.68774e-16  -7.96157e-17  -1.351        -1.0009       -0.747327   
  6.98651e-17  -1.1411e-16    1.47652      -0.747327      2.14316    

##
Q = H2*H
display(Q*A*Q')

5×5 Array{Float64,2}:
 -0.786842      2.58457       1.43493e-16  -1.77876e-16   1.23098e-16
  2.58457      -1.56744      -3.65238      -2.04182e-16  -1.55982e-16
  1.43493e-16  -3.65238       2.40268      -1.351         1.47652    
 -1.77876e-16  -2.27584e-16  -1.351        -1.0009       -0.747327   
  1.23098e-16   9.32264e-17   1.47652      -0.747327      2.14316    

##
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)

([1.6368617153752258 -8.808317733951808 … 0.0 0.0; -8.808317733951808 -4.055803925667554 … 0.0 0.0; … ; 0.0 0.0 … -1.3932924109386884 1.3801333741205086; 0.0 0.0 … 1.3801333741205104 0.5479952964627801], [0.0 0.0 … 0.0 0.0; 0.7211559710182545 0.0 … 0.0 0.0; … ; -0.017553019339705166 -0.14612954495066524 … 0.0 0.0; -0.17172241809637437 0.06134032652868235 … 0.0 0.0])

##

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

myhess (generic function with 1 method)

## 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 Array{Float64,2}:
 2.54959    0.323006   2.22045e-16  -2.82291e-16   1.53821e-16
 0.323006  -1.22303   -2.9905       -4.44089e-16   0.0        
 0.0       -2.9905     0.690147      2.00179      -2.22045e-16
 0.0        0.0        2.00179      -1.10299      -1.90437    
 0.0        0.0        0.0          -1.90437      -3.1724     

5×5 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
 -0.866551  -0.247104   0.00746166  -0.371725    0.223144  
  0.499088  -0.42904    0.0129555   -0.645415    0.387437  
  0.0       -0.868829  -0.00851976   0.424437   -0.254786  
  0.0        0.0       -0.999852    -0.0147536   0.00885648
  0.0        0.0        0.0          0.51468     0.857382