Implement Hessenberg factorization

General Usage
1

I am trying to Implement Hessenberg factorization in the following code.

function GivensRotation(a::Float64, b::Float64)
    # Calculate the Given's rotation that rotates [a;b] to [r;0]:
c = 0.; s = 0.; r = 0.
if ( b == 0. )
    c = sign(a)
    s = 0.
    r = abs(a)
elseif ( a == 0. )
    c = 0.
    s = -sign(b)
    r = abs(b)
elseif ( abs(a) .> abs(b) )
    t = b/a
    u = sign(a)*abs(sqrt(1+t^2))
    c = 1/u
    s = -c*t
    r = a*u
else
    t = a/b
    u = sign(b)*abs(sqrt(1+t^2))
    s = -1/u
    c = -s*t
    r = b*u
end
    return (c, s, r)
end

function HessenbergQR( R::Matrix )
    # Compute the QR factorization of an upper-Hessenberg matrix: 
    
    n = size(R, 1)
    Q = Matrix(I,(n,n))
    # Convert R from upper-hessenberg form to upper-triangular form using n-1 Givens rotations:
    for k = 1:n-1
        (c, s, r) = GivensRotation(R[k,k], R[k+1,k])
        # Compute G*R[k:k+1,:] and Q[:,k:k+1]*G', where G = [c -s ; s c]
        for j = 1:n
            newrow = c*R[k,j] - s*R[k+1,j]
            R[k+1,j] = s*R[k,j] + c*R[k+1,j]
            R[k,j] = newrow
            newcol = c*Q[j,k] - s*Q[j,k+1]
            Q[j,k+1] = s*Q[j,k] + c*Q[j,k+1]
            Q[j,k] = newcol
        end
    end
    return (Q, R)
end

H = triu(rand(3,3),-1)
HessenbergQR( H )

I got this error
InexactError: Bool(-0.7006778795038133)

any tips ?

2

Not sure I can help, but you might have better luck if you change the title to be something a bit more descriptive.

3

I feel like the problem in this line ( Q = Matrix(I,(n,n)) ) but i don’t know who to fix it

4

I think changing it to Matrix{Float64}(I, n, n) might work then. (Didn’t test because I’m on my phone). The issue being that it creates a boolean matrix there unless you specify the type.

3 Likes
5

Note that the LinearAlgebra standard library does support Hessenberg factorization via hessenberg(A):

julia> H = hessenberg(rand(5,5)).H
5×5 UpperHessenberg{Float64,Array{Float64,2}}:
  0.455787  -0.531629   0.318809    0.409985   -0.213419
 -1.59706    2.04202   -0.353893   -0.275307    0.308906
   ⋅         0.128048  -0.0584893   0.0974416   0.130365
   ⋅          ⋅         0.146707   -0.482426    0.160377
   ⋅          ⋅          ⋅          0.404674    0.598712

However, this built-in UpperHessenberg type does not currently have an optimized QR factorization routine (it just falls back to the generic dense QR).

I’m curious, in what application do you have an upper-Hessenberg matrix and need its QR factors?

Adding an efficient algorithm here would be a nice addition to the standard library. However, an efficient implementation would not compute the Q matrix explicitly. Instead, you’d store the list of Givens rotations in a new subtype of AbstractQ, so that you can multiply a vector by Q or Q' in O(N) time and storage.

3 Likes