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
# This i
function gallery_house(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 = gallery_house(A[i:end,j])
  H = eye(size(A,1),size(A,1))
  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}:
 -1.32486       3.20749   5.55112e-17   0.0       -4.44089e-16
  3.20749      -1.2732    1.06367      -2.27325   -1.61763    
  5.55112e-17   1.06367  -2.63722      -1.41774    2.22075    
  0.0          -2.27325  -1.41774       0.977076   1.18732    
 -4.44089e-16  -1.61763   2.22075       1.18732   -1.89615    
In [3]:
##
B = H*A*H'
H2 = housemat(B,3,2)[1]
display(H2*B*H2')
5×5 Array{Float64,2}:
 -1.32486       3.20749      -2.6036e-16   1.77315e-16  -3.17913e-16
  3.20749      -1.2732       -2.98594      4.44089e-16  -4.44089e-16
 -2.6036e-16   -2.98594       0.566402     2.17101      -0.715965   
  1.77315e-16   3.33067e-16   2.17101     -4.06153       0.683944   
 -3.17913e-16   4.44089e-16  -0.715965     0.683944     -0.0611616  
In [4]:
##
Q = H2*H
display(Q*A*Q')
5×5 Array{Float64,2}:
 -1.32486       3.20749      -1.66533e-16   0.0          -2.22045e-16
  3.20749      -1.2732       -2.98594       1.94289e-16   5.55112e-17
 -1.66533e-16  -2.98594       0.566402      2.17101      -0.715965   
  0.0           2.22045e-16   2.17101      -4.06153       0.683944   
 -2.22045e-16   2.22045e-16  -0.715965      0.683944     -0.0611616  
In [5]:
##
In [6]:
##
using Plots
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] = A[k+1:n,k:n] - 2*vk*(vk'*A[k+1:n,k:n])
    A[1:n,k+1: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(32,32)
A = A + A'
animate_myhess(A)
MethodError: no method matching getindex(::PyPlot.Figure, ::Symbol)
The applicable method may be too new: running in world age 24070, while current world is 24078.

Stacktrace:
 [1] _before_layout_calcs(::Plots.Plot{Plots.PyPlotBackend}) at /Users/dgleich/.julia/v0.6/Plots/src/backends/pyplot.jl:936
 [2] prepare_output(::Plots.Plot{Plots.PyPlotBackend}) at /Users/dgleich/.julia/v0.6/Plots/src/plot.jl:251
 [3] display at /Users/dgleich/.julia/v0.6/Plots/src/output.jl:144 [inlined]
 [4] gui(::Plots.Plot{Plots.PyPlotBackend}) at /Users/dgleich/.julia/v0.6/Plots/src/output.jl:134
 [5] animate_myhess(::Array{Float64,2}) at ./In[6]:18
 [6] include_string(::String, ::String) at ./loading.jl:515
In [7]:
##

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] = A[k+1:n,k:n] - 2*vk*(vk'*A[k+1:n,k:n])
    A[1:n,k+1:n] = A[1:n,k+1:n] - 2*(A[1:n,k+1:n]*vk)*vk'
  end
  return A,V
end
Out[7]:
myhess (generic function with 1 method)
In [8]:
## 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)
5×5 Array{Float64,2}:
  3.38142  -2.28539   -5.55112e-17   4.44089e-16  -1.66533e-16
 -2.28539  -1.04004    0.519925     -1.38778e-16   4.16334e-17
  0.0       0.519925  -0.0288868    -0.984518      1.11022e-16
  0.0       0.0       -0.984518      0.235731      0.0786327  
  0.0       0.0        0.0           0.0786327     0.167332