Frontmatter

If you are publishing this notebook on the web, you can set the parameters below to provide HTML metadata. This is useful for search engines and social media.

Notebook for 2023-03-08

Basic QR iteration

The basic QR iteration is

$$\begin{align*} A^{(k)} &= Q^{(k)} R^{(k)} \\ A^{(k+1)} &= R^{(k)} Q^{(k)} \end{align*}$$
md"""
# Notebook for 2023-03-08

## Basic QR iteration

The basic QR iteration is

$$\begin{align*}
A^{(k)} &= Q^{(k)} R^{(k)} \\
A^{(k+1)} &= R^{(k)} Q^{(k)}
\end{align*}$$
"""
128 ms
using LinearAlgebra
266 μs
using Plots
2.7 s
sympart (generic function with 1 method)
sympart(A) = (A+A')/2
278 μsA
10×10 Matrix{Float64}:
 0.349856  0.52095   0.421871  0.460524   …  0.350273  0.401769  0.337866   0.224796
 0.52095   0.652474  0.377681  0.651064      0.678696  0.678326  0.695169   0.371226
 0.421871  0.377681  0.973693  0.479816      0.545122  0.691404  0.280855   0.245604
 0.460524  0.651064  0.479816  0.177518      0.766478  0.355246  0.0976494  0.591847
 0.526778  0.738552  0.268403  0.392175      0.549213  0.36241   0.238191   0.712749
 0.879979  0.273363  0.897537  0.549474   …  0.934543  0.151518  0.627998   0.512173
 0.350273  0.678696  0.545122  0.766478      0.230484  0.369452  0.471768   0.526793
 0.401769  0.678326  0.691404  0.355246      0.369452  0.700768  0.478094   0.646642
 0.337866  0.695169  0.280855  0.0976494     0.471768  0.478094  0.815433   0.254728
 0.224796  0.371226  0.245604  0.591847      0.526793  0.646642  0.254728   0.615335
A = sympart(rand(10,10))
47.0 μs
naive_qr_iteration (generic function with 1 method)
function naive_qr_iteration(A, nsteps; monitor=(A)->nothing)
for k = 1:nsteps
monitor(A)
Q, R = qr(A)
A = R*Q
end
A
end
1.7 ms

Remarkably enough, the basic iteration does converge! But it takes some time.
md"""
Remarkably enough, the basic iteration does converge! But it takes some time.
"""
116 μs
test_naive_qr (generic function with 1 method)
function test_naive_qr(A, nsteps)
rs = []
naive_qr_iteration(A, nsteps, monitor=(A)->push!(rs, norm(tril(A,-1))))
rs
end
703 μs
plot(test_naive_qr(A, 1000), yscale=:log10)
2.7 s

Hessenberg reduction

The first step for fast QR iteration is reduction to upper Hessenberg form, i.e.

$$H = Q^T A Q$$

where $Q$ is an orthogonal matrix. As described in the notes, we can do this via Householder transformations. It is possible (as with ordinary QR) to save the vectors defining the Householder transformation where we have zeros, but we are not going to bother. Instead, we'll form $Q$ explicitly by accumulating the Householder transformations.
md"""
## Hessenberg reduction

The first step for fast QR iteration is reduction to upper Hessenberg form, i.e.

$$H = Q^T A Q$$

where $Q$ is an orthogonal matrix. As described in the notes, we can do this
via Householder transformations. It is possible (as with ordinary QR) to save
the vectors defining the Householder transformation where we have zeros, but we
are not going to bother. Instead, we'll form $Q$ explicitly by accumulating the
Householder transformations.
"""
204 μs
householder (generic function with 1 method)
function householder(x)
normx = norm(x)
u = copy(x)
u[1] = x[1] + sign(x[1])*normx
u /= norm(u)
u
end
442 μs
my_hessenberg (generic function with 1 method)
function my_hessenberg(A)
H = copy(A)
n = size(A)[1]
Q = Matrix{Float64}(I,n,n)

for j = 1:n-2
u = householder(H[j+1:n,j])
H[j+1:n,j:n] -= 2*u*(u'*H[j+1:n,j:n])
H[j+2:n,j] .= 0.0
H[:,j+1:n] -= 2*(H[:,j+1:n]*u)*u'
Q[:,j+1:n] -= 2*(Q[:,j+1:n]*u)*u'
end

Q, H
end
1.4 ms
begin
B = rand(10,10)
Q, H = my_hessenberg(B)
norm(B*Q-Q*H)/norm(B), norm(tril(H,-2)), norm(Q'*Q-I)
end
225 μs

Bidiagonalization

The bidiagonal reduction involves computing

$$U^T A V = B$$

where $B$ is an upper bidiagonal matrix. We compute it via a sequence of interleaved Householder transformations on the left and right.
md"""
## Bidiagonalization

The bidiagonal reduction involves computing

$$U^T A V = B$$

where $$B$$ is an upper bidiagonal matrix. We compute it via a sequence of interleaved Householder transformations on the left and right.
"""
167 μs
my_bidiagonalization (generic function with 1 method)
function my_bidiagonalization(A)
B = copy(A)
n = size(A)[1]
U = Matrix{Float64}(I,n,n)
V = Matrix{Float64}(I,n,n)

for j = 1:n-1

# Transform from the left to introduce zeros in column j
u = householder(B[j:n,j])
B[j:n,j:n] -= 2*u*(u'*B[j:n,j:n])
B[j+1:n,j] .= 0.0
U[:,j:n] -= 2*(U[:,j:n]*u)*u'

# Transform from the right to introduce zeros in row j
v = householder(B[j,j+1:n])
B[j:n,j+1:n] -= 2*(B[j:n,j+1:n]*v)*v'
B[j,j+2:n] .= 0.0
V[:,j+1:n] -= 2*(V[:,j+1:n]*v)*v'

end

B, U, V
end
1.8 ms
begin
BB, U, V = my_bidiagonalization(B)
norm(U'*B*V-BB)/norm(B), norm(tril(BB,-1)), norm(triu(BB,2)), norm(U'*U-I), norm(V'*V-I)
end
109 μs