# Implement Hessenberg factorization

**URL:** <https://discourse.julialang.org/t/implement-hessenberg-factorization/40006>\
**Category:** General Usage\
**Tags:** question\
**Created:** [May 23, 2020, 12:37am UTC](https://discourse.julialang.org/t/implement-hessenberg-factorization/40006 "2020-05-23T00:37:25Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![A.U.M](https://avatars.discourse-cdn.com/v4/letter/a/aeb1de/32.png) [@A.U.M](https://discourse.julialang.org/u/A.U.M)\
**Post date:** [May 23, 2020, 12:37am UTC](https://discourse.julialang.org/t/implement-hessenberg-factorization/40006/1 "2020-05-23T00:37:25Z")

</div>

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

```julia
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

```

```julia
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 ?

---

<div class="post-metadata">

**Author:** ![kevbonham](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kevbonham/32/216165_2.png) [@kevbonham](https://discourse.julialang.org/u/kevbonham)\
**Post date:** [May 23, 2020, 12:54am UTC](https://discourse.julialang.org/t/implement-hessenberg-factorization/40006/2 "2020-05-23T00:54:36Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![A.U.M](https://avatars.discourse-cdn.com/v4/letter/a/aeb1de/32.png) [@A.U.M](https://discourse.julialang.org/u/A.U.M)\
**Post date:** [May 23, 2020, 1:08am UTC](https://discourse.julialang.org/t/implement-hessenberg-factorization/40006/3 "2020-05-23T01:08:26Z")

</div>

> [@A.U.M](#):
>
> `Q = Matrix(I,(n,n))`

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

---

<div class="post-metadata">

**Author:** ![ericphanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ericphanson/32/215186_2.png) [@ericphanson](https://discourse.julialang.org/u/ericphanson)\
**Post date:** [May 23, 2020, 1:59am UTC](https://discourse.julialang.org/t/implement-hessenberg-factorization/40006/4 "2020-05-23T01:59:58Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 23, 2020, 1:44pm UTC](https://discourse.julialang.org/t/implement-hessenberg-factorization/40006/5 "2020-05-23T13:44:59Z")

</div>

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

```julia
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.