# Julia slow with Hilbert matrices

**URL:** https://discourse.julialang.org/t/julia-slow-with-hilbert-matrices/7905
**Category:** Performance
**Created:** [December 21, 2017, 10:46pm UTC](https://discourse.julialang.org/t/julia-slow-with-hilbert-matrices/7905 "2017-12-21T22:46:46Z")
**Posts on this page:** 10
**Page:** 1

<div class="post-metadata">

### Author: ![jmichel](https://avatars.discourse-cdn.com/v4/letter/j/b3f665/32.png) [@jmichel](https://discourse.julialang.org/u/jmichel)
#### Post date: [December 21, 2017, 10:46pm UTC](https://discourse.julialang.org/t/julia-slow-with-hilbert-matrices/7905/1 "2017-12-21T22:46:46Z")

</div>

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,

---

<div class="post-metadata">

### Author: ![Jordan\_Cluts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jordan_cluts/32/13753_2.png) [@Jordan\_Cluts](https://discourse.julialang.org/u/Jordan_Cluts)
#### Post date: [December 21, 2017, 11:39pm UTC](https://discourse.julialang.org/t/julia-slow-with-hilbert-matrices/7905/2 "2017-12-21T23:39:43Z")

</div>

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.

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 21, 2017, 11:56pm UTC](https://discourse.julialang.org/t/julia-slow-with-hilbert-matrices/7905/3 "2017-12-21T23:56:21Z")

</div>

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

---

<div class="post-metadata">

### Author: ![jmichel](https://avatars.discourse-cdn.com/v4/letter/j/b3f665/32.png) [@jmichel](https://discourse.julialang.org/u/jmichel)
#### Post date: [December 22, 2017, 1:14am UTC](https://discourse.julialang.org/t/julia-slow-with-hilbert-matrices/7905/4 "2017-12-22T01:14:11Z")

</div>

I did run my code several times…

---

<div class="post-metadata">

### Author: ![jmichel](https://avatars.discourse-cdn.com/v4/letter/j/b3f665/32.png) [@jmichel](https://discourse.julialang.org/u/jmichel)
#### Post date: [December 22, 2017, 1:14am UTC](https://discourse.julialang.org/t/julia-slow-with-hilbert-matrices/7905/5 "2017-12-22T01:14:35Z")

</div>

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

---

<div class="post-metadata">

### Author: ![ihnorton](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ihnorton/32/26_2.png) [@ihnorton](https://discourse.julialang.org/u/ihnorton)
#### Post date: [December 22, 2017, 1:19am UTC](https://discourse.julialang.org/t/julia-slow-with-hilbert-matrices/7905/6 "2017-12-22T01:19:36Z")

</div>

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

[https://github.com/JuliaLang/julia/issues/11202](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)

---

<div class="post-metadata">

### Author: ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)
#### Post date: [December 22, 2017, 3:18am UTC](https://discourse.julialang.org/t/julia-slow-with-hilbert-matrices/7905/7 "2017-12-22T03:18:20Z")

</div>

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

```julia
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)

```

---

<div class="post-metadata">

### Author: ![JeffreySarnoff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeffreysarnoff/32/1980_2.png) [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)
#### Post date: [December 22, 2017, 3:20am UTC](https://discourse.julialang.org/t/julia-slow-with-hilbert-matrices/7905/8 "2017-12-22T03:20:52Z")

</div>

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

```julia
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.

---

<div class="post-metadata">

### Author: ![jmichel](https://avatars.discourse-cdn.com/v4/letter/j/b3f665/32.png) [@jmichel](https://discourse.julialang.org/u/jmichel)
#### Post date: [December 22, 2017, 7:04am UTC](https://discourse.julialang.org/t/julia-slow-with-hilbert-matrices/7905/9 "2017-12-22T07:04:15Z")

</div>

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.

---

<div class="post-metadata">

### Author: ![JeffreySarnoff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeffreysarnoff/32/1980_2.png) [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)
#### Post date: [December 22, 2017, 7:17am UTC](https://discourse.julialang.org/t/julia-slow-with-hilbert-matrices/7905/10 "2017-12-22T07:17:05Z")

</div>

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](https://julialang.org/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.
