# Julia slow with Hilbert matrices

#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

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).

#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:

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

Couldn't be good to postpone 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)
``````

#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.

#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.