Julia slow with Hilbert matrices

General Usage Performance
1

Hello,
I am new to Julia, which I learn in the hope to port some mathematical
software (group theory and arithmetic) to. I did a simple test with Hilbert
matrices, and to my disappointment I found that on this test Julia is
slower than all the other CAS systems I have access to.

Here is the Julia code I used:

 hilb=n->[1//(BigInt(i)+j-1) for i in 1:n, j in 1:n]
 function test(n) hilb(n)^-1 end
 julia> @time test(80);
 5.158600 seconds (42.41 M allocations: 1.072 GiB, 31.96% gc time)

Here is the code in GAP

hilb:=n->List([1..n],i->List([1..n],j->1/(i+j-1)));
hilb(80)^-1;;time;
706   (0.706 second)

Now the Magma version:

hilb:=function(n) return Matrix(n,n,[1/(i+j-1) : i in [1..n], j in[1..n]]);
end function;
time hilb(80)^-1;
140 (0.140 second)

Now the Maple version

with(linalg);
m:=Matrix(80,(i,j)->1/(i+j-1));
inverse(m); # takes 1.5 sec

It is also slower than in the languages Ruby and Python:

In Ruby

require "benchmark"
require "matrix"
m=Matrix.build(80,80){|i,j| Rational(1,i+j+1)}
Benchmark.measure{m.inv}
=> 2.4 sec

In Python via SageMath:

sage: a=Matrix([[1/(i+j-1) for j in (1..80)]for i in (1..80)])
sage: timeit('a^-1')
5 loops, best of 3: 1.03 s per loop

So what I am doing wrong in Julia? Is there another type than BigInt which
guarantees that hilb(80)*hilb(80)^-1 is the identity matrix? Is there a way
in Julia which is not 35 times slower than Magma?

Best regards,

2

Did you run your Julia code twice? The first time you run a bit of code includes the compilation time which can slow execution (the first time) by an order of magnitude or more. Subsequent runs will use the cached compiled version and run much faster. Consequently it is standard practice to first run the function to compile the time it. (Or run it twice and discard the first timing result).

I’m new here so apologies if this is off base.

3

The problem is that Rational numbers do a reduction step after each operation. @JeffreySarnoff I think was working on changing that? Also, is using BigInt necessary here? They are slower than machine integers (of course). Are the other languages defaulting to machine integers or “infinite sized” integers?

But most importantly, here you’re using a generic matrix inversion algorithm instead of one specifically for matrices of rational numbers (i.e. it’s being sent through an algorithm which was developed with floating point numbers in mind). I assume that the CAS’s are calling some smarter algorithm in this case. Specialized linear algebra for rational numbers isn’t part of Julia’s Base (it just has generic fallbacks to standard numerical methods which work because of how the Rational type works, but isn’t necessarily the best thing to do), but I think you may find it in a package like Nemo.jl (Nemo.jl is essentially a CAS for Julia, like SageMath is for Python).

2 Likes
4

I did run my code several times…

5

All the other languages default to “infinite sized” integers, which is necessary
to have hilb(80)*hilb(80)^-1 equal to the identity matrix.

6

Big{Int,Float} are on the slower side currently:

https://github.com/JuliaLang/julia/issues/11202

(improvements there should be a non-breaking change as far as I know, and can be done in a point release after 1.0)

7

To second @ChrisRackauckas suggestion, you really should visit http://nemocas.org given your domain of interest, and check their benchmarks. Here is a quick try at your problem:

julia> using Nemo
Welcome to Nemo version 0.7.0
Nemo comes with absolutely no warranty whatsoever
julia> S=MatrixSpace(QQ,80,80)
Matrix Space of 80 rows and 80 columns over Rational Field
julia> hilb=n->[one(QQ)//(QQ(i)+j-1) for i in 1:n, j in 1:n]
(::#1) (generic function with 1 method)
julia> function test(n) inv(S(hilb(n))); end
test (generic function with 1 method)
#[skip first run for compilation]
julia> @time test(80);
  0.786288 seconds (293.70 k allocations: 48.261 MiB, 0.84% gc time)
3 Likes
8

Nemo does support this (I get ~5.5x speedup);

using Nemo
M = MatrixSpace(QQ, 80, 80)
h = hilbert(M)
@elapsed inv(h)

There may well be a way to do it faster in Nemo, idk.

I did mess with Rational – got no response, so I trashed it.

3 Likes
9

Tank you for your answer. Do you know whether the timing difference of Nemo is mainly due to a better handling of Rationals?

By the way, there seems to be some compatiblity problems with Nemo. I succeeded in installing it only on one of 3 computers.

10

Nemo uses packages that are best-in-class for Number Theorists. The timing difference is not about a “better handling of Rationals” — Nemo handles Integers, Rationals … as Number Fields, so it is more accurate to say Nemo uses a “different handling of Rationals” that is well suited to the specific problem you posed.

Nemo should install on any of the official v0.6.2 downloads, with Pkg.add("Nemo"). If you have installed ArbFloats, first
Pkg.rm("ArbFloats"); Pkg.rm("ArbFloats");. And [then] do
Pkg.rm("Nemo");Pkg.rm("Nemo");Pkg.add(Nemo); quit();
(yes, twice each)
now using Nemo should work (maybe with some warnings, there should not be errors). Then re-add ArbFloats, if you had it before.

2 Likes