# Householder routines reflector! and reflectorApply! in LinearAlgebra

**URL:** <https://discourse.julialang.org/t/householder-routines-reflector-and-reflectorapply-in-linearalgebra/56766>\
**Category:** Numerics\
**Tags:** linearalgebra\
**Created:** [March 8, 2021, 10:13pm UTC](https://discourse.julialang.org/t/householder-routines-reflector-and-reflectorapply-in-linearalgebra/56766 "2021-03-08T22:13:34Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [March 8, 2021, 10:13pm UTC](https://discourse.julialang.org/t/householder-routines-reflector-and-reflectorapply-in-linearalgebra/56766/1 "2021-03-08T22:13:34Z")

</div>

Hello,

I would like to use the routines `reflector!` and `reflectorApply!` of LinearAlgebra ([julia/generic.jl at master · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/generic.jl)) 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](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?

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

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

```

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [March 9, 2021, 8:55am UTC](https://discourse.julialang.org/t/householder-routines-reflector-and-reflectorapply-in-linearalgebra/56766/2 "2021-03-09T08:55:31Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [March 9, 2021, 4:05pm UTC](https://discourse.julialang.org/t/householder-routines-reflector-and-reflectorapply-in-linearalgebra/56766/3 "2021-03-09T16:05:53Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![andreasnoack](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/andreasnoack/32/27_2.png) [@andreasnoack](https://discourse.julialang.org/u/andreasnoack)\
**Post date:** [March 9, 2021, 8:36pm UTC](https://discourse.julialang.org/t/householder-routines-reflector-and-reflectorapply-in-linearalgebra/56766/4 "2021-03-09T20:36:40Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![mleprovost](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mleprovost/32/7166_2.png) [@mleprovost](https://discourse.julialang.org/u/mleprovost)\
**Post date:** [March 10, 2021, 3:15am UTC](https://discourse.julialang.org/t/householder-routines-reflector-and-reflectorapply-in-linearalgebra/56766/5 "2021-03-10T03:15:33Z")

</div>

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