detx: don't catch exceptions and retry with wider type

[nsajko]

This puts detx in line with the common Julia conventions (the limited width of the input is the caller's responsibility), while improving performance and type stability.

[scheinerman]

I'm not seeing how this works. Here's what I tried:

julia> A = rand(Int,20,20) .% 1000;

julia> det(A//true)
ERROR: OverflowError: 481953496777 * -94142821213 overflowed for type Int64
Stacktrace:
  [1] throw_overflowerr_binaryop(op::Symbol, x::Int64, y::Int64)
    @ Base.Checked ./checked.jl:163
....

julia> detx(A)
11376312482736357124775885872759468327049946982629425088942433216

I'm not comfortable accepting this PR.

[nsajko]

Tried to explain everything above already. Basically, there's a very important principle in Julia called "type stability". The idea is that we want the compiler to be able to infer the exact concrete type of the return value of most functions (especially public API, like this one) given the function and the types of the arguments.

detx currently fails the above test:

julia> using LinearAlgebraX, Test

julia> A = rand(Int,20,20);

julia> @inferred detx(A)
ERROR: return type BigInt does not match inferred return type Union{Int64, BigInt}

There's basically three things you can do:

1. Nothing. Keep the current type unstable behavior. This makes your package largely unsuitable for use by others, especially as a dependency of other packages.
2. Accept this PR.
3. Promote to big unconditionally, instead of the try-catch. This is a better choice than (1), but generally goes against the good practice mentioned above, "the limited width of the input is the caller's responsibility".

[scheinerman]

Thanks @nsajko. Your explanation makes a lot of sense. While I might not fully grasp the importance of type stability, I get that it's important and understand that I break it with my code. Your suggestion is this:

function detx(A::AbstractMatrix{T}) where {T<:IntegerX}
    det(A // true)
end

For an integer matrix, though, the resulting determinant would be Rational and not some flavor of Integer. Would it be better to do this instead:

    T.(det(A//true))

What do you think?

[nsajko]

I might not fully grasp the importance of type stability

The problem with type unstable code is that instabilities often propagate through to other code, causing inefficiency and run time dispatch. For example, imagine some code like this: f(g(x)). If the g(x) call is not type stable, Julia can't infer the type of g(x), so it doesn't know which method of f should be called after g is done. Then a run time dispatch is necessary, which is much slower than static dispatch and inhibits compiler optimization.

For an integer matrix, though, the resulting determinant would be Rational and not some flavor of Integer. Would it be better to do this instead

You're absolutely right, except without broadcasting, I think, so T(det(A//true)).

Now that I think about it, it might make sense to instead promote narrow integers to BigInt, like here (from JuliaAlgebra/MultivariatePolynomials.jl#282):

https://github.com/JuliaAlgebra/MultivariatePolynomials.jl/blob/1d31b075fdc815baf6e507eac1d33d2beae7f5be/src/det.jl#L46-L72

Or maybe fully embrace option (3) and promote rationals to big rationals, too. This package is specifically meant for exact computation so it could make sense to do that.

[scheinerman]

Thanks! I think option (3) makes the most sense and is consistent with how I mean for this package to work. I'll try to make those changes soon.

[nsajko]

OK, then feel free to close this if you don't need it.

[scheinerman]

@nsajko I just submitted version 0.2.0 to be registered. I think I fixed the type stability issues at the cost of always returning big values (where it made sense). Thanks for explaining the importance of this. If you find any problems, please let me know. In 0.2.1 I hope to include exact complex numbers, but for now those are not implemented.

[nsajko]

Nice!