Julia is very slow compared to python/numpy. What am I doing wrong here?

General Usage Performance
,
1

Hi guys,
I wrote a recursive function which construct a matrix in python using Qutip module:
(I just used the tensor function and Pauli matrices definitions from Qutip)

import numpy as np
from qutip import *

def psitest(K):
    if K > 1 and K % 1 == 0:
        return tensor(psitest(K-1), - sigmaz())
    elif K == 1:
        return (1/np.sqrt(2)) * sigmax()

then for example for K=16 I have

Time=0.017930984

I wrote exactly the same function in Julia:

using LinearAlgebra
function ψtest(K)
    if K > 1 && K % 1 == 0
        return kron(ψtest(K-1), [-1. 0.;0. 1.])
    elseif K == 1
        return 1/√(2) * [0. 1.;1. 0.]
    end
end

where for the same K=16 gives

39.406 s (65 allocations: 42.67 GiB)

for K =25 actually python gives it in about 2 second where in Julia I can not even test it.
I know I should not naively transform a code in python to Julia, but I am wondering what is here that makes this huge difference.
Thank you all for your help

1 Like
2

I assume Qutip uses sparse arrays. You can do the same in Julia (like you - I have not tried to optimize code - just added sparse):

julia> using SparseArrays

julia> function ψtest(K)
           if K > 1 && K % 1 == 0
               return kron(ψtest(K-1), sparse([-1. 0.;0. 1.]))
           elseif K == 1
               return 1/√(2) * sparse([0. 1.;1. 0.])
           end
       end
ψtest (generic function with 1 method)

julia> @time ψtest(16);
  0.002063 seconds (182 allocations: 3.012 MiB)

julia> @time ψtest(25);
  1.169190 seconds (290 allocations: 1.500 GiB, 17.61% gc time)

with dense matrix K=25 for sure will not fit in your RAM.

17 Likes
3

Wow, thank you very much :pray:. @bkamins.
I appreciate if you also tell me other kind of optimizations I can do in this kind of computations.
My problem is finally the SYK model in quantum mechanics.

4

I am not a linear algebra expert, but assuming you want to keep the recursive definition then I think this is mostly OK. There might be a way to reduce the amount of allocations by pre-allocating the final matrix and using kron! on views, but it would need to be benchmarked.

4 Likes
5

I should check these. Thanks @bkamins again.

6

You are building a linear operator that you will eventually be applying to one or more vectors to transform them. If so, it will be much more effective to actually compute the result of applying the operator onto the data vector instead of the operator itself. That is, write a function that computes ψtest(k)(v) for any vector v.

For a general methodology to translate recursive operators that involve Kronecker products into code, please allow me to advertise my own work “The Kronecker Product in Approximation and Fast Transform Generation” abstract only or the full PhD pdf.

With Julia’s symbolic manipulation, term rewriting and introspection capabilities, I have been itching to redo the above in a single environment instead of the combination of conditional term rewriting in Mathematica, and the recurrence unrolling, expression simplification and C/FORTRAN/MATLAB code generation in Maple that I used 25 years ago…

8 Likes
7

Thanks @pitsianis for the explanations. I will look into your work for sure.
My problem is coding for many body quantum systems like SYK model. Here for example the matrix finally is the Hamiltonian where I am interested in its eigenvalues. These \psi's are Majorana fermions which satisfy Clifford Algebra. Actually I do not know what I am really to compute :smile:. For now I am just re-deriving what people have obtained.

8

The same advice holds for eigenvalues, since eig(kron(A,B)) = kron(eig(A),eig(B))

The right-hand-side expression is much cheaper to compute than the left.

In your specific recurrence example, it is even easier.

3 Likes
9

Thanks @pitsianis, this formula eig(kron(A,B)) = kron(eig(A),eig(B)) is really cool and I should use it definitely. However in this part of my problem at the end I have H = \psi_1 \psi_2 \psi_3 ... where \psi_i = \sigma_1 \otimes \sigma_2 .... and I want the spectrum of H. Do you have any suggestion in this case? Beside that do you know any cheatsheet like references that I can find these kind of magic formulas?!
thank you very much

10

Start from pages 4-8 in the PhD Thesis that I pointed above. A good reference for me at the time was the book by Willi-Hans Steeb. “Kronecker Product of Matrices and Applications”, Wissenschaftsverlag, 1991.

I also recommend: Charles F. Van Loan, The ubiquitous Kronecker product, Journal of Computational and Applied Mathematics, Volume 123, Issues 1–2, 2000, Pages 85-100

2 Likes
11

Thanks for the references. I start with your thesis.

12

also, your function is not type-stable: it will return nothing when K is less than 1 or non integer

1 Like
13

Thank you. yes, I should take care of that.

14

In case you are not aware QuantumOptics.jl has significant overlap with QuTip.

1 Like
15

Yeah, I just add it and trying to learn it. However I am not quit sure that it is suitable for my kind of problem as I want some rare used algebras like Clifford algebra or others. Maybe it has it. Thanks for the suggestion though.
Sorry a quick question: is numpy.complex128 is the same as Complex{Float64} in Julia. I think they are same but not sure. I am quit surprised that how much Julia is easy and fast by the way. I wrote the code naively in Julia and is 4-5x faster even without using any special package (I thought first that this might be due to precision).

1 Like
16

They are the same.

1 Like
17

That is a relief. THanks

18 
julia> ComplexF64
ComplexF64 (alias for Complex{Float64})

IIRC Julia called this Complex128 , but it was changed before 1.0

3 Likes
19

You might want to check something like this, I don’t know if it’s good for what you want to do

https://github.com/Krastanov/QuantumClifford.jl

1 Like
20

Thanks for sharing this. It is not for my current job. But it is interesting. Maybe in the future.

1 Like