## 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 Matrix{Float64}:
-0.779083 3.77777 -1.58747e-17 -1.21998e-16 -9.8519e-17
3.77777 0.175024 -0.106515 -0.40171 -0.634336
-1.58747e-17 -0.106515 0.477158 1.43702 2.49054
-1.21998e-16 -0.40171 1.43702 -0.93729 0.626665
-9.8519e-17 -0.634336 2.49054 0.626665 0.838653
##
B = H*A*H'
H2 = housemat(B,3,2)[1]
display(H2*B*H2')
5×5 Matrix{Float64}:
-0.779083 3.77777 1.49262e-16 -4.52965e-17 2.26e-17
3.77777 0.175024 0.758352 9.67017e-17 2.18976e-16
1.49262e-16 0.758352 1.68758 1.90662 1.9899
-4.52965e-17 -1.43206e-17 1.90662 -0.76218 0.32621
2.26e-17 2.46872e-17 1.9899 0.32621 -0.546874
##
Q = H2*H
display(Q*A*Q')
5×5 Matrix{Float64}:
-0.779083 3.77777 -9.04361e-17 -1.58106e-17 8.29753e-17
3.77777 0.175024 0.758352 1.17925e-16 2.44333e-16
-9.04361e-17 0.758352 1.68758 1.90662 1.9899
-1.58106e-17 1.40027e-17 1.90662 -0.76218 0.32621
8.29753e-17 2.60237e-16 1.9899 0.32621 -0.546874
##
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:2296 [inlined]
[2] macro expansion
@ ./lock.jl:273 [inlined]
[3] __require(into::Module, mod::Symbol)
@ Base ./loading.jl:2271
[4] #invoke_in_world#3
@ ./essentials.jl:1089 [inlined]
[5] invoke_in_world
@ ./essentials.jl:1086 [inlined]
[6] require(into::Module, mod::Symbol)
@ Base ./loading.jl:2260
[7] top-level scope
@ ~/.julia/packages/Plots/FFuQi/src/backends.jl:883
[8] eval
@ ./boot.jl:430 [inlined]
[9] _initialize_backend(pkg::Plots.PyPlotBackend)
@ Plots ~/.julia/packages/Plots/FFuQi/src/backends.jl:882
[10] backend(pkg::Plots.PyPlotBackend)
@ Plots ~/.julia/packages/Plots/FFuQi/src/backends.jl:245
[11] pyplot()
@ Plots ~/.julia/packages/Plots/FFuQi/src/backends.jl:86
[12] top-level scope
@ In[5]:4
##
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 Matrix{Float64}:
-0.528236 -2.89719 -4.44089e-16 1.18401e-15 -1.01855e-15
-2.89719 0.993984 2.48371 4.44089e-16 -7.63916e-16
0.0 2.48371 0.718501 2.81394 4.44089e-16
0.0 0.0 2.81394 0.189853 1.46831
0.0 0.0 0.0 1.46831 0.859402
5×5 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}