Exact linear algebra

Specific Domains Numerics
1

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
2

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
3

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

4

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
5

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
6

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.

7

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.

8

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
9

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
10

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

11

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