LinearAlgebra.reflectorApply!

[dlfivefifty]

I'm giving a seminar in a few weeks using QR and just realised that LinearAlgebra.reflectorApply! has issues with if the matrix is not diagonally dominant. Here's a simple demonstration of the problem (only works on Julia v1.8-pre, but can also be seen in Julia v1.6 using ArrayLayouts.reflectorApply!):

julia> LinearAlgebra.reflectorApply!([@interval(1),-1], @interval(1), [[1, @interval(-20,-18.18181818181818)];; ])
2×1 Matrix{Interval{Float64}}:
 [-20, -18.1818]
       [-0.818182, 2.81819]

This actually should just be a permutation, so should have no error:

julia> LinearAlgebra.reflectorApply!([1,-1], 1., [[1, -20.];; ])
2×1 Matrix{Float64}:
 -20.0
   1.0

julia> LinearAlgebra.reflectorApply!([1,-1], 1., [[1, -18.18181818181818];; ])
2×1 Matrix{Float64}:
 -18.18181818181818
   1.0

Do you know of any references on how to implement an interval-arithmetic friendly version of Householder? It's an orthogonal matrix so a change in implementation should fix this.

[dlfivefifty]

Here's a simple demonstration of where the issue is:

julia> y = [[1, @interval(-20,-18.18181818181818)];; ]
2×1 Matrix{Interval{Float64}}:
    [1, 1]
 [-20, -18.1818]

julia> x = [-1,1];  I - x*x'
2×2 Matrix{Int64}:
 0  1
 1  0

julia> (I - x*x')*y # stable to make the matrix first and then apply
2×1 Matrix{Interval{Float64}}:
 [-20, -18.1818]
    [1, 1]

julia> y - x*x'y # multiplying through raises a problem
2×1 Matrix{Interval{Float64}}:
 [-20, -18.1818]
       [-0.818182, 2.81819]

So one option would be to always form the matrix. Here's a quick-and-dirty fix:

function ArrayLayouts.reflectorApply!(x::AbstractVector, τ::Number, y::AbstractVecOrMat{<:Interval})
    v = [1; x[2:end]]
    copyto!(y, (I - τ*v*v')*y)
end

[dlfivefifty]

PS This works now for (almost...) rigorous bounds on roots Bessel functions by simply specifying the 3-term recurrence:

julia> J̃ = (z,_...) -> (SymTridiagonal(-2*(0:∞)/z, ones(∞)) \ [1; zeros(∞)])[1]
#7 (generic function with 1 method)

julia> roots(J̃, 2.4..2.405)
1-element Vector{Root{Interval{Float64}}}:
 Root([2.40482, 2.40483], :unique)

[dpsanders]

This looks like a great example of the classic problem with interval arithmetic, namely the dependency problem. This occurs when you have have the same variable appearing >1 time in a single expression (here, y), and interval arithmetic is "not able to take account" of the fact that those variables are the same.

The simplest example of this is x - x. Instead of computing {x - x: x ∈ X}, it computes {x - y: x ∈ X and y ∈ X}, so that [0, 1] - [0, 1] == [-1, 1] instead of [0, 0].

The solution is indeed the one that you suggest: write it as (I - x*x') * y, so that y occurs only once.

However, usually any iterative algorithm like this will accumulate ever larger errors as you do more steps. I don't know the literature on interval versions of Householder, but I found https://link.springer.com/article/10.1007%2Fs11155-006-2968-5

cc @lucaferranti

[mforets]

Taking products of interval matrices is not even associative, so this problem is (unfortunately) not uncommon. A general rule is to avoid computing the same terms, some articles call this "single-use expressions". For some cases, such as taking matrix powers, we have implemented specialized algorithms in IntervalMatrices.jl.

@lucaferranti's current GSOC is about IntervalLinearAlgebra.jl. Adding a special code for Householder transformations seems like a nice addition to that library.

EDIT: oops, just saw David's comment when I was finishing mine..

[dlfivefifty]

However, usually any iterative algorithm like this will accumulate ever larger errors

My application is asymptotically diagonally dominant which rescues this.

[OlivierHnt]

Closing this issue since linear algebra related features should be reported in IntervalLinearAlgebra.jl.