LU instability for non-float matrices due to lack of pivoting

[stevengj]

As was pointed out recently on discourse, the lu function suffers from numerical instability for non-floating-point matrices, where it tries to do the factorization without pivoting first.

(The example on the mailing list was det([1 1 0 0 1 0 0 0; 1 0 1 0 0 1 0 0; 1 0 0 1 0 0 1 0; 0 1 1 1 0 0 0 0; 0 1 0 0 0 0 1 1; 0 0 1 0 1 0 0 1; 0 0 0 1 1 1 0 0; 0 0 0 0 1 1 0 1]).)

The problem seems to stem from PR JuliaLang/julia#10215 by @andreasnoack, and indeed the abovementioned determinant started being wrong in Julia 0.4.

As @tpapp suggested, we should default to pivoting at least whenever lutype(T) <: AbstractFloat, but we also want to handle cases like lutype(Complex{Int}) == Complex{Float64}. Perhaps the right criterion is to pivot whenever abs(zero(lutype(T))) isa AbstractFloat since pivoting is based on the abs function.

(Not sure whether we should try to support types that don't have abs via try-catch?)

[andreasnoack]

Yes. Something like abs(zero(lutype(T))) isa AbstractFloat would be a good check. I would prefer nok to introduce try-catch and requiring abd doesn't seem very restrictive.

[barpit20]

Hi,
I am a first time contributor. I tried fixing this bug by adding the following lines of code before computing the LU factorisation without pivoting-

if(abs(zero(S)) isa AbstractFloat)
    return lu!(AA, Val(true); check = check)
end

This seemed to fix the bug but make testall failed. The testset "lu with type whose sum is another type" is giving an error because the abs function is not defined for the structs in the trickyarithmetic.jl file. I was able to define the abs function for the structs A,B and C but the struct D always gives me a no method matching error. I tried different method signatures but none seem to work. Can I please get some help with this?

Struct D for reference:

https://github.com/JuliaLang/julia/blob/589b96dd6c47aa497e5be9f7ce683b222b76bd8c/stdlib/LinearAlgebra/test/trickyarithmetic.jl#L28-L35

[KlausC]

You could chang the type check to

S <: Union{AbstractFloat,Complex{<:AbstractFloat}}

[stevengj]

@KlausC, it's not enough to check for real or complex, since we also want to handle lots of other types (e.g. user-defined types for quaternions, dimensionful quantities, and so on).

Probably it should be abs(one(...)) and not abs(zero(...)) to correctly handle dimensionful types (for which one will be dimensionless).

We could do something like:

onex = one(lutype(T))
if applicable(abs, onex) && abs(onex) isa AbstractFloat
    return lu!(....)
end

to handle the case where abs does not work.

[stevengj]

As long as we are calling applicable, it seems like a more reliable method is

if applicable(float, onex) && float(onex) === onex

since that also works for types like Measurement from Measurement.jl and Dual from ForwardDiff.jl, whereas the abs method doesn't work for those types.

[barpit20]

On changing the condition to if applicable(float, onex) && float(onex) === onex, one test fails due to the error- MethodError: no method matching AbstractFloat(::TrickyArithmetic.D{TrickyArithmetic.C,TrickyArithmetic.C}).

However, for the condition if applicable(abs, onex) && abs(onex) isa AbstractFloat, all tests pass. Which one should we stick with?

[stevengj]

I think we should get rid of the float(x) = AbstractFloat(x) fallback. Producing an AbstractFloat type is simply not what float does in general, so this fallback is wrong.

Once/if we do this, the applicable(float, onex) method will work.

[barpit20]

Removing or modifying the callback causes a LoadError for mathconstants.jl at line 14. Initially, I thought that this was due to the float() function not being specifically defined for types like BigInt and BigFloat but they are defined in mpfr.jl. Can I get some help on how to proceed with the issue?

[stevengj]

AbstractIrrational should have a float method, maybe