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

General Usage Performance
,
21

Hi just saw this post, and just want to show off Yao.jl here, you can use Yao’s operator system to construct matrix (it is not just a circuit simulator!), the syntax is also more readable than plain matrices

julia> using Yao

julia> function yao(K)
           if K > 1 && K % 1 == 0
               kron(yao(K-1), -Z)
           elseif K == 1
               1/sqrt(2) * X
           end
       end
yao (generic function with 1 method)

julia> @time mat(yao(25)); # with compile time
  7.625354 seconds (23.85 M allocations: 2.792 GiB, 6.06% gc time, 90.30% compilation time)

julia> @time mat(yao(25)); # without compile time
  0.563007 seconds (1.63 k allocations: 1.500 GiB, 10.98% gc time)

and on my machine the qutip version

In [1]: 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()
   ...:

In [2]:

In [2]: %timeit psitest(25)
1.6 s ± 3.69 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

with Julia SparseArrays

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.021314 seconds (36.51 k allocations: 5.085 MiB, 88.99% compilation time)

julia> @time ψtest(25);
  1.408559 seconds (290 allocations: 1.500 GiB, 17.22% gc time)
7 Likes
22

Hi @Roger-luo. Thank you very much :pray::pray::pray:. You guys are amazing :star_struck:.
I should work with Yao then. I think I should start with the documentation.

23

Sorry, Just a quick question. What is the identity matrix (for general dim) in Yao?
There is IdentityGate but I can not take kron(Identity{2},X).
And can I work with LuxurySparse.jl and get the same result?
I need this because I am trying to construct a Majorana fermion:

function ψ(K, i)
    if K > 1 && K % 1 == 0
        if i%1 == 0 && 1 <= i <= (2*K)-2
            return kron(ψ(K-1,i), sparse([-1. 0.;0. 1.]))
        elseif i == (2*K) -1
            return kron(sparse(I,2^(K-1),2^(K-1)), 1/√2 * sparse([0 1;1 0]))
        elseif i == (2*K)
            return kron(sparse(I,2^(K-1),2^(K-1)), 1/√2 * sparse([0 -im;im 0]))
        else
            print("i = $i out of range!")
        end
    elseif K == 1
        if i == 1
            return   sparse(1/√2 * [0 1;1 0])
        elseif i == 2
            return   sparse(1/√2 * [0 -im;im 0])
        else
            print("i = $i out of range!")
        end
    else
        print("K = $K is out of range!")
    end
end
24

Just to report back. I tried LuxurySparse and there is an improvement. But Yao might be better (if only I can find the identity :smile:).

@btime [ψLuxspa(10,i) for i=1:2*10];        149.987 μs (1859 allocations: 839.81 KiB)
@btime [ψ(10,i) for i=1:2*10];          450.979 μs (4016 allocations: 1.30 MiB)
25

You want to construct an operator not using the type,

julia> IdentityGate{2}(10)|>nqubits
10

here 2 means it’s 2-level system operator, 10 stands for 10 sites/qubits

2 Likes
26

Great. Thank you very much. so dim here is 2^10. :+1:

1 Like