things.jl

things.jl

using Base.LinAlg: givensAlgorithm
import Base: convert

struct Householder{T}
    u::T
end

"""
Returns a reflector s.t. Q * x = norm(x) * e1.
"""
function householder(x)
    u = copy(x)
    u[1] -= norm(x)
    u ./= norm(u)
    Householder(u)
end

convert(::Type{Matrix}, Q::Householder) = I - 2 * Q.u * Q.u'

function run_qr(H_original, debug = true)
    @assert size(H_original, 2) == size(H_original, 1)

    n = size(H_original, 1)
    indices = collect("₁₂₃₄₅₆₇₈₉")
    H = copy(H_original)
    Q = eye(n)

    if debug
        println("Initial H")
        display("text/plain", H)
        println()
    end

    # α, β = H[n-1, n-1], H[n-1, n]
    # γ, δ = H[n, n-1],   H[n, n]
    # tr = α + δ
    # det = α * δ - γ * β
    # D = tr^2 - 4det

    # if D < 0
    #     return
    # end

    # μ₁ = tr + √D

    μ₁ = 1.10
    μ₂ = 1.4

    # Compute the entries of (H - μ₂)(H - μ₁)e₁.
    # Please don't do it like this :p, we should not store the household reflection as a
    # vector -- if we would wanna store it as a vector, then we should use SVector from 
    # StaticArrays.jl; but yeah, 2 givens rotations work OK as well
    fst = H[1:3, 1] # H * e₁
    fst[1] -= μ₁ # - μ₁ * e₁
    fst .= H[1:3, 1:2] * fst[1:2] - μ₂ * fst # well...

    # Sanity check
    @assert fst ≈ ((H - μ₂ * I) * ((H - μ₁ * I)[:, 1]))[1:3]

    G = eye(n)
    G[1:3,1:3] .= convert(Matrix, householder(fst))
    H = G * H
    
    if debug
        println("Compute Householder reflection that maps G₁ * (H - μ₂)(H - μ₁)e₁ = r₁₁e₁")
        display("text/plain", H)
        println()
    end

    H *= G' # and apply the rotation from the right.
    Q *= G' # callback, somehow

    if debug
        println("Apply the Housholder reflection from the rhs: G₁ * H * G₁'")
        display("text/plain", H)
        str = "G₁ * H * G₁'"
    end

    # Restore to Hessenberg form
    # Bulge is initially in H[1:3,1:3]
    for i = 2 : n - 1

        # Restore...
        range = i : min(i + 2, n) # when i = n - 1, this is just a given's rotation
        G = eye(n)
        G[range,range] .= convert(Matrix, householder(H[range, i - 1]))
        H = G * H

        if debug
            println("\n------\n\n")
            str = "G$(indices[i]) * " * str
            println("Restore Hessenberg form: $str")
            display("text/plain", H)
            println()
        end

        # Destroy
        H *= G'
        Q *= G'  # callback, somehow

        if debug
            str = str * " * G$(indices[i])'"
            println("Apply from the rhs: $str")
            display("text/plain", H)
            println()
        end
    end

    return H_original, H, Q
end

function example_2()
    H = triu(rand(10, 10), -1)
    run_qr(H, true)
end