# Householder routines reflector! and reflectorApply! in LinearAlgebra

Hello,

I would like to use the routines `reflector!` and `reflectorApply!` of LinearAlgebra (julia/generic.jl at master · JuliaLang/julia · GitHub) to efficiently apply an Householder reflection to a vector a.

The setup is the following:

Step 1: Given a vector b, compute the associated Householder vector v, such that Hb = (I-2 vv^\top)b = -sign(b[1]) ||b||_2 \boldsymbol{e}_1. This can be done using the routine `reflector!`. The code is based on http://www.netlib.org/lapack/lawnspdf/lawn72.pdf

Step 2: Apply the Householder reflection to a vector a, i.e. compute (I - 2 v v^\top)a. This can be done using the routine `reflectorApply!`.

In LinearAlgebra, `reflectorApply!(x::AbstractArray{T,1} where T, tau::Number, A::StridedArray{T, 2} where T)` is only defined for a `StridedArray{T,2}` `A`, what is the best way to apply this reflector when `A` is an `AbstractVector`? I am using `repeat(a,1,1)` to create an `Array{Float64,2}` from `a`. Can we avoid this conversion?

``````using LinearAlgebra

b = randn(5)
@show norm(b)
v = deepcopy(b)
bcopy = deepcopy(b)
a = randn(5)

tau = deepcopy(LinearAlgebra.reflector!(v))

Hb = LinearAlgebra.reflectorApply!(v, tau, repeat(b, 1, 1))
@show Hb
@show norm(bcopy)
@show -norm(bcopy)*sign(bcopy[1])

#  What is the best way to apply the reflection to a vector a?
@show Ha = LinearAlgebra.reflectorApply!(v, tau, repeat(a,1,1))
``````

Output:

``````norm(b) = 1.9799448125238643
Hb = [-1.9799448125238646; -5.551115123125783e-17; 2.220446049250313e-16; 2.220446049250313e-16; -2.220446049250313e-16]
norm(bcopy) = 1.9799448125238643
-(norm(bcopy)) * sign(bcopy[1]) = -1.9799448125238643
Ha = LinearAlgebra.reflectorApply!(v, tau, repeat(a, 1, 1)) = [-0.5062333940143671; -1.1850776857048684; -1.295711281614973; -0.16447887645751813; 1.065353231871041]
``````

These routines were intended for DenseArrays. The LAPACK approach, which was the inspiration behind the linear algebra routines in Julia, is to reuse as much of the input memory as possible so mutation is ubiquitous in numerical linear algebra. Why is your vector not a `StridedVector`? Sometimes the problem can be fixed by simply making a custom array a subtype of `DenseArray` instead of `AbstractArray` but, indeed, that might not always be feasible.

Thank you for the explanation. My post was a bit confusing about the `x` vector. I have modified it.
My question is how can we use the function `reflectorApply!(x::AbstractArray{T,1} where T, tau::Number, A::StridedArray{T, 2} where T)` when `A` is an `Array{Float64,1}`? Can we avoid the conversion to an `Array{Float64,2}`?

I see. I made the signature more restrictive than necessary because I wrote these with the QR factorization in mind. It should be fine to make it `StridedVecOrMat` so please open a PR and ping me if you think the functionality would be useful. It will, of course, take a while before you can benefit from the change so meanwhile, you can use `reshape` since it will just alias the original array but with a different shape so in that way, you can make it a `Matrix` without a copy.

1 Like

Thank you for the advice. I forgot that reshape is acting like a view.