 # Implement Hessenberg factorization

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 ?

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

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

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

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