# Exact linear algebra

Hey all!

I’ve been using Julia for algebra computations (think polynomial rings, ideals, Groebner basis, etc). Linear algebra is often fundamental to that, and Julia seems very suited to that!

There’s a few choices in the current `linalg` package that make it a bit challenging at the moment. Here’s a two simple examples:

``````julia> M = BigInt[1 1; 2 2]
julia> typeof(det(M)) # expecting a BigInt obtained from "a*c - b*d" (or a more advanced, but still exact, algorithm)
BigFloat
julia> nullspace(M) # expecting the space spanned by [-1, 1]
ERROR: MethodError: no method matching svdfact!(::Array{BigFloat,2}; thin=false)
``````

What are community opinions on this? Would there be any way to change this behaviour in future versions of Julia? If so, is this desirable? If yes, how can I help?

1 Like

For the `nullspace` thing: you could use https://github.com/JuliaMath/GenericSVD.jl, which according to https://github.com/JuliaLang/julia/issues/5429 should make it into `Base.LinAlg` at some point. You could help pulling that into Base.LinAlg by submitting a pull request to https://github.com/JuliaLang/julia (perhaps after checking with the author of GenericSVD.jl and the people who replied to the issue, and making sure that everything is sufficiently tested).

Your `det` example goes through `det(A::AbstractArray{T,2}) where T in Base.LinAlg at linalg/generic.jl:1205` on 0.6, which calls `lufact`, which entails doing divisions, and `typeof(BigInt(1) / BigInt(2))` is `BigFloat`. I think this can be called a bug, and you should report it to Julia. Again, if you code up an algorithm that does it, please do consider contributing it back to `Base.LinAlg`.

In general, I think the Julia developers are fully on board with supporting generic linear algebra algorithms for arbitrary scalar types, but there might be some gaps here and there.

1 Like

`det` for generic matrices works via `lufact`. There are methods for linear algebra with integers, but they would need to be implemented.

The vast majority of linear algebra algorithms involve division, so try `Rational{BigInt}`.

As you can see, the generic `nullspace()` uses an SVD algorithm which is iterative and so would get out of hand with rationals. Hence there is a forced conversion to the associated `Float` type. (Since the algorithm is not exact this seems appropriate.)

``````julia> svd(Matrix{Rational{BigInt}}([1 1; 2 2]))
(BigFloat[-4.472135954999579392818347337462552470881236719223051448541794490821041851275639e-01 8.944271909999158785636694674925104941762473438446102897083588981642083702551192e-01; -8.944271909999158785636694674925104941762473438446102897083588981642083702551192e-01 -4.472135954999579392818347337462552470881236719223051448541794490821041851275639e-01], BigFloat[3.162277660168379331998893544432718533719555139325216826857504852792594438639231, 0.000000000000000000000000000000000000000000000000000000000000000000000000000000], BigFloat[-7.071067811865475244008443621048490392848359376884740365883398689953662392310596e-01 7.07106781186547524400844362104849039284835937688474036588339868995366239231051e-01; -7.07106781186547524400844362104849039284835937688474036588339868995366239231051e-01 -7.071067811865475244008443621048490392848359376884740365883398689953662392310596e-01])

julia> svd(Matrix{Rational{Int}}([1 1; 2 2]))
(Float32[-0.447214 -0.894427; -0.894427 0.447214], Float32[3.16228, 8.42937f-8], Float32[-0.707107 -0.707107; -0.707107 0.707107])

``````

@andreasnoack Why is the associated float for `Rational{Int}` only `Float32`?

1 Like

This is not a bug. It was decided early on that `/` for two integers in Julia should produce a floating-point value. If you want a `Rational` result, use `Rational` inputs.

It might be interesting to have a library for exact linear-algebra calculations on integer matrices, perhaps along the lines of LinBox. But since the algorithms appear to be quite different, it seems like it should be a separate package from “ordinary” generic linear algebra routines that rely on division.

2 Likes

Sorry, I was imprecise. I meant that the fact that `det(::AbstractMatrix{BigInt})` returns a `BigFloat` could be called a bug (as, theoretically, no divisions should be needed to compute the determinant). I just wanted to show that `typeof(BigInt(1) / BigInt(2))` being `BigFloat` is the reason for the current behavior.

It is the consequence of `promote_rule(Rational{Int},Float32)` which like `promote_rule(Int,Float32)` returns `Float32`. The idea is to promote to the smallest BLAS type possible.

Thanks for the pointers everyone! In particular, thanks to @tkoolen for pointing to the `GenericSVD` package and for the sanity check about what types we should ideally expect.

I’ll be looking into whether this solves my use case now.

1 Like

You may be away of this, but in Nemo we have algorithms for this. While previously you had to use the integers provided by Nemo (`fmpz` from flint), due to the efforts @wbhart you can now use all of our generic matrix functionality with `BigInt`:

``````julia> using Nemo

julia> Matrix(JuliaZZ, 2, 2, BigInt[3, 3, 1, 2])
[3 3]
[1 2]

julia> det(ans)
3

julia> typeof(ans)
BigInt
``````

(Current master branch.)

We are also working on splitting Nemo into two packages, so that you can use these algorithms without installing any C libraries. But this is still work in progress.

3 Likes

What is the difference in performance between using the Nemo integers and the `BigInt`?

The difference will be quite big.

The algorithms which will be called for `BigInt` are very generic. They basically only use that you can add and multiply two elements from the underlying ring. Over the integers they are prone to coefficient explosion. We could actually implement a fast modular algorithm for the determinant of matrices of `BigInts`, but this is currently not on our agenda.

If you use the Nemo integers, the computations will be delagated to flint, which has specialized code for matrices over integers and it will be quite fast.

Here is a quick benchmark:

``````julia> function perf_bigint(B, n)
A = rand(-B:B, n, n)
M = Matrix(JuliaZZ, n, n, A)
return det(M)
end

julia> function perf_fmpz(B, n)
A = rand(-B:B, n, n)
M = Matrix(FlintZZ, n, n, A)
return det(M)
end

julia> @time perf_fmpz(1000, 100);
0.023958 seconds (15.53 k allocations: 1.403 MiB)

julia> @time perf_bigint(1000, 100);
0.319814 seconds (1.49 M allocations: 58.265 MiB, 6.19% gc time)

julia> @time perf_fmpz(2^60, 100);
0.056312 seconds (30.57 k allocations: 1.921 MiB, 19.65% gc time)

julia> @time perf_bigint(2^60, 100);
1.362713 seconds (1.55 M allocations: 212.311 MiB, 3.04% gc time)
``````
3 Likes