Possible regression in A \ b operations on 4x4 matrices since commit 9b6427fd

[gwater]

One of the test cases for IntervalRootFinding.jl recently broke and we traced the issue to this recent commit in StaticArrays.jl. Currently, it is unclear how exactly the error emerges but previous to that commit our test case works fine. Unfortunately, that means we cannot suggest a possible fix or provide a minimal example right now.

As far as I can tell, the code in question is called from here. Our test case uses 4-element SVectors and static 4x4 matrices which explains why the commit only affects that specific test case.

[gwater]

I can now provide an example where StaticArray behaviour deviates from Base.Array

using IntervalArithmetic
using StaticArrays

m = IntervalBox(2.7912..2.79121, 2.7912..2.79121, 1.16882..1.16883, 1.16882..1.16883)
jac = [-18.7299..12.0483 -13.6028..13.2686 -1.61778..3.89494 -1.25794..2.11343; -11.7012..10.1435 -11.88..8.31026 -1.25794..2.11343 -0.889139..2.35177; -10.7905..10.0659 -9.6589..8.27525 -2.69591..1.26513 -0.903859..1.51855; -5.78849..4.3997 -4.51621..4.90624 -0.628891..1.05658 -1.72913.. -0.108871]

@show jac \ m
#4-element Array{IntervalArithmetic.Interval{Float64},1}:
# [-∞, ∞]
# [-∞, ∞]
# [-∞, ∞]
# [-∞, ∞]

@show SMatrix{4,4}(jac) \ SVector{4}(m)
#4-element StaticArrays.SArray{Tuple{4},IntervalArithmetic.Interval{Float64},1,4}:
# [-∞, ∞]                
# [-∞, ∞]                
# [-∞, ∞]                
# [-10.736, -0.675958]

As you can see the new left divide method produces some error in this case. I assume this somehow related to the values contained in the matrix and how they interact through IntervalArithmetic.jl.

[tkoolen]

The LU decomposition seems to be different:

julia> sL, sU, sp = lu(SMatrix{4,4}(jac));

julia> L
4×4 Array{IntervalArithmetic.Interval{Float64},2}:
  [1, 1]   [0, 0]   [0, 0]  [0, 0]
 [-∞, ∞]   [1, 1]   [0, 0]  [0, 0]
 [-∞, ∞]  [-∞, ∞]   [1, 1]  [0, 0]
 [-∞, ∞]  [-∞, ∞]  [-∞, ∞]  [1, 1]

julia> sL
4×4 LowerTriangular{IntervalArithmetic.Interval{Float64},StaticArrays.SArray{Tuple{4,4},IntervalArithmetic.Interval{Float64},2,16}}:
 [1, 1]    ⋅       ⋅       ⋅   
 [0, 0]  [1, 1]    ⋅       ⋅   
 [0, 0]  [0, 0]  [1, 1]    ⋅   
 [0, 0]  [0, 0]  [0, 0]  [1, 1]

julia> U
4×4 Array{IntervalArithmetic.Interval{Float64},2}:
   [-18.73, 12.0484]     [-13.6029, 13.2687]      [-1.61779, 3.89495]      [-1.25795, 2.11344]
 [0, 0]               [-∞, ∞]                 [-∞, ∞]                  [-∞, ∞]                
 [0, 0]                [0, 0]                 [-∞, ∞]                  [-∞, ∞]                
 [0, 0]                [0, 0]                  [0, 0]                  [-∞, ∞]                

julia> sU
4×4 UpperTriangular{IntervalArithmetic.Interval{Float64},StaticArrays.SArray{Tuple{4,4},IntervalArithmetic.Interval{Float64},2,16}}:
   [-18.73, 12.0484]    [-13.6029, 13.2687]     [-1.61779, 3.89495]  [-1.25795, 2.11344] 
   ⋅                    [-11.8801, 8.31027]     [-1.25795, 2.11344]  [-0.88914, 2.35178] 
   ⋅                    ⋅                       [-2.69592, 1.26514]  [-0.90386, 1.51856] 
   ⋅                    ⋅                      ⋅                     [-1.72914, -0.10887]

@martinholters probably knows the most about the LU code.

[martinholters]

Obviously, neither L*U nor sL*sU give jac. But I suppose that jac should be within (w.r.t. to the intervals) sL*sU, right? (It is for L*U.) My guess would be that I incorrectly assume x-x == zero(x) or x/x == one(x) somewhere, but I'll take a closer look.

[martinholters]

No, it's the check whether pivotization was successful. This should fix it:

diff --git a/src/lu.jl b/src/lu.jl
index cbd89d0..6666211 100644
--- a/src/lu.jl
+++ b/src/lu.jl
@@ -101,7 +101,7 @@ function __lu(A::StaticMatrix{M,N,T}, ::Type{Val{Pivot}}) where {M,N,T,Pivot}
         # Scale first column
         Akkinv = inv(A[kp,1])
         Ls = A[ps,1] * Akkinv
-        if !isfinite(Akkinv)
+        if iszero(A[kp,1])
             Ls = zeros(Ls)
         end

The ideal solution probably would be to use isinf(Akkinv) and overload isinf in IntervalArithmetic to return true only for [-∞, -∞] and [∞, ∞], i.e. only if the value has to be Inf.

[gwater]

I can confirm that with this change the original test case passes without error.

Maybe @dpsanders can comment on the proposal to change isinf() in IntervalArithmetic. I suspect its behavior might be prescribed by IEEE 1788-2015.

EDIT: I looks like isinf(::Interval) was redefined very recently.

[dpsanders]

It's subtle, since [-∞, -∞] and [∞, ∞] are not valid intervals.

I'm a bit worried about this change, since an interval may contain 0, without being equal to 0.

The solution for intervals is probably to move to using interval-specific versions of \.

[c42f]

So @dpsanders I'm assuming we shouldn't go ahead with Martin's one line solution then?

Is there any action we should be taking in StaticArrays itself to make this work better for you?

[dpsanders]

@gwater Do you have a suggestion as to how to fix LU to make it work generically here?
I see that with NumberIntervals it fails because it encounters a missing. This is a nice solution, but it still doesn't give the correct answer!

[gwater]

Hmm, hard to say. The trail ends at a copysign operation which fails because at least one interval contains positive and negative numbers. Furthermore, if we look at the determinant of this matrix we find:

julia> det(SMatrix{4,4}(NumberInterval.(jac)))
x ∈ [-10294, 10667.8]

which to me indicates that the set of matrices represented by this interval matrix contains some matrices which are singular. I don't know enough about LU decomposition to say if it could nevertheless work. I think, it makes more sense to divert to some other procedure and avoid LU if the interval matrix contains singular matrices.

[gwater]

Also, and I just remembered this, returning a vector that contained the entire real line in each component was fine. The problems started when the interval in the fourth component was wrongly constrained too much. So, IntervalRootFinding could check if the determinant interval contains zero and instead of calling ldiv it could instantly return the aforemention vector of real lines. Not elegant - but a generic solution to this specific bug.