# Can this random matrix simulation be sped up?

I am a Julia novice and converted the following program from Python (for a simple practice exercise). It computes a N by N random matrix of quaternions, drawn from the Gaussian distribution, and then computes the smallest singular value, and repeats this many times. (I have omitted the graphing of the resulting distribution for simplicity.) We use that quaternions may be represented as 2 by 2 complex matrices.

I have a few basic questions, which elude me even after looking through the manual.

1. The version below hard-codes the constants 100 and 1000. When these are replaced by `n` and `samples` in the `data` function, the code becomes .5 seconds slower (4.8 versus 5.3). Why is this, and is there a way to prevent this slowdown without hard-coding the constants, so that the values are easier to change?

2. This is a literal translation of a Python script. Are there any easy ways to make it faster or more memory efficient? I feel Iâ€™m surely missing a few tricks. In particular, it produces 2.5 GiB of memory allocations, which seems sub-optimal.

3. Can the compilation of the script be cached, so that I do not re-pay the compilation cost each time the script is run without being edited?

4. Are there any Julia packages/modules whose source code is particularly well-written and instructive, which I could read to pick up idiomatic, high-performance Julia? I was thinking of the linear algebra library, but then I realized that itâ€™s probably just calls to OpenBLAS, which is not written in Julia, so probably that would not help. Perhaps the ODE/PDE/SDE solving library?

I apologize for any misunderstandings in these questions â€“ I am still an inexperienced programmer, and an even more inexperienced Julia user. Thank you!

``````using LinearAlgebra

n = 100
samples = 1000

function symplectic(N)
G = Array{Complex{Float64}}(undef, 2N ,2N)

for i in 0:N-1
for j in 0:N-1
a = randn()
b = randn()
c = randn()
d = randn()
@inbounds G[2*i + 1,2*j + 1] = a + b*im
@inbounds G[2*i+1 + 1,2*j+1 + 1] = a - b*im
@inbounds G[2*i+1 + 1,2*j + 1] = -c + d*im
@inbounds G[2*i + 1,2*j+1+ 1] = c + d*im
end
end

M = G * G'
ev = eigvals(M)
return ev[1]*N
end

function data()
for i in 1:1000
println(i," ",symplectic(100))
end
end

@time data()
``````

I wouldnâ€™t recommend my own code here since the ODE/PDE/SDE solving libraries are almost akin to being OpenBLAS sometimes, with its own definitions of broadcast variants, exponentiation, etc. Thereâ€™s some things to learn from it on how to write high-performance generic code, but thereâ€™s also a ton of extra things that I wouldnâ€™t recommend people would do unless their goal was to simply make users as happy as possible.

1 Like

Instead of generating single scalars with rand() running it many times inside the nested loops you could create 100x100 matrices with rand(100,100) before the loops.

Translating code is a non-trivial task - without the original python code it will be hard to find equivalences or mistakes in your julia version.

First, ideally you can pass `n` and `samples` as parameters to `data()`:

`data(n,nsamples)`

avoid using global variables (or, at least, use `const n = 100`, etc).

You could allocate `G` in `data` and pass it as parameter to `sympletic`. And you could also allocate other auxiliary arrays (`M`, for instance), and pass it as well as parameter to `sympletic`.

Something like this (needs checking and still it might be possible to get the eigenvalues into a prealocated array, such that the function should not allocate anything):

``````using LinearAlgebra

function symplectic(N,G,R)

for i in 0:N-1
for j in 0:N-1
@. R = randn()
@inbounds G[2*i + 1,2*j + 1] = R[1] + R[2]*im
@inbounds G[2*i+1 + 1,2*j+1 + 1] = R[1] - R[2]*im
@inbounds G[2*i+1 + 1,2*j + 1] = -R[3] + R[4]*im
@inbounds G[2*i + 1,2*j+1+ 1] = R[3] + R[4]*im
end
end
ev = svdvals!(G)
return ev[2*N]*N
end

function data(n,samples)
G = Array{Complex{Float64}}(undef, 2n ,2n)
R = zeros(4)
for i in 1:samples
symplectic(n,G,R)
end
end

n = 100
samples = 1000
@time data(n,samples)

``````

(of course you can loop from `1 to N` and remove the `+1` from the `G[...] = `.

1 Like

this is equivalent to squaring `svdvals(G)` (albeit sorted in ascending rather than descending order).

2 Likes

Thank you everyone for the suggestions! They are helpful, and I am working through them carefully.

Here are my preliminary observations.

1. Using `const` indeed removes the penalty from using the global variables (compared to hard-coding), which resolves one of my questions.

2. Strangely, using `svdvals(G)` squaring, and taking the element in index `end` instead of computing the eigenvalues of `G * G'` as above almost doubles the running time of the program. Iâ€™m not sure whatâ€™s going on here, but this can be observed both when I run @lmiqâ€™s version, and edit mine to include the SVD.

3. @inbounds does nothing in this script, and should likely be removed from it for clarity.

4. Initializing G in `data()` makes the script about `.1` second slower, but uses about `2 GiB` instead of `2.5 GiB` memory. I donâ€™t see much change initializing the random numbers outside of the function `symplectic()`.

Keep in mind that for a more proper benchmarking you should better use BenchmarkTools. Or, at least, run the code twice, the first to compile, the second to measure the time and allocations.

Thank you for pointing this out. Doing this, I get virtually identical results in time to my original script (4.8 seconds) but with less memory usage (1.9 vs 2.5 GiB), so your suggestion is indeed a win. Sorry for the mistake.

I still wonder if there is a way to cache the compilation between successive runs of the same scriptâ€¦

If two executions occur in the same Julia section functions that did not change (and didnâ€™t receive different types of parameters) wonâ€™t be compiled again.

Check out https://github.com/timholy/Revise.jl which should allow you to revise your code without restarting Julia. That should remove most of the compilation burden, since you will only have to recompile the code youâ€™ve changed.

2 Likes

In my test the allocations are reduced to 250 Mb:

`````` 17.909110 seconds (616.27 k allocations: 256.346 MiB, 1.37% gc time)
``````

This is the code:

``````using LinearAlgebra

function symplectic(N,G,R)

for i in 0:N-1
for j in 0:N-1
@. R = randn()
@inbounds G[2*i + 1,2*j + 1] = R[1] + R[2]*im
@inbounds G[2*i+1 + 1,2*j+1 + 1] = R[1] - R[2]*im
@inbounds G[2*i+1 + 1,2*j + 1] = -R[3] + R[4]*im
@inbounds G[2*i + 1,2*j+1+ 1] = R[3] + R[4]*im
end
end
ev = svdvals!(G)
return ev[2*N]*N
end

function data(n,samples)
G = Array{Complex{Float64}}(undef, 2n ,2n)
R = zeros(4)
for i in 1:samples
println(symplectic(n,G,R))
end
end

n = 100
samples = 1000
data(n,samples)
@time data(n,samples)

``````

And the difference (to 2Gib) is using `svdvals!` (note the exclamation mark)

Edit: I am not sure if this is doing the same as your code, exactly. But appropriate use of `adjoint!` and `mul!` could reduce the allocations even further.

1 Like

Interesting, thank you. I reproduce your drastic decrease in memory usage, to about 220 MiB. However, the time it takes to run the script about doubles when using `svdvals!`.

When I compute the singular values â€śby handâ€ť by first doing `G * G'` and then `eigvals!` (note the exclamation mark), it runs in about half the time as the `svdvals!` version but uses about 5 times as much memory. This is still about half the memory usage of the non-exclamation mark version (the original), very roughly.

I am puzzled why the dedicated svd routine is twice as slow as the direct computation Iâ€™m using.

1 Like

Itâ€™s more accurate. Computing `G * G'` squares the condition number.

4 Likes

Hi,

The correct equivalence between SVD and EIG is when,

M = G * Gâ€™
ev_eig = eigvals(M)
ev_svd = svdvals(M)

And hence, for `i=1,...,N`.

ev_svd[i] == abs(ev_eig[N+1-i])

(Normally you would take the SVD of G, not M, in which case you would take the square roots of `ev_eig`. Also, you might then have to omit some zero eigenvalues for non-square or rank-deficient G.)

2 Likes

Thanks, @stevengj, my apologizes by my misconception