Householder routines reflector! and reflectorApply! in LinearAlgebra

Specific Domains Numerics
1

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]
2

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.

3

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

4

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
5

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