# How to speed up this Kronecker Multiplication?

**URL:** <https://discourse.julialang.org/t/how-to-speed-up-this-kronecker-multiplication/112063>\
**Category:** Performance\
**Tags:** question\
**Created:** [March 25, 2024, 5:32am UTC](https://discourse.julialang.org/t/how-to-speed-up-this-kronecker-multiplication/112063 "2024-03-25T05:32:18Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![vikram\_backup](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vikram_backup/32/202851_2.png) [@vikram\_backup](https://discourse.julialang.org/u/vikram_backup)\
**Post date:** [March 25, 2024, 5:32am UTC](https://discourse.julialang.org/t/how-to-speed-up-this-kronecker-multiplication/112063/1 "2024-03-25T05:32:19Z")

</div>

I am writing a code that gives me tensor product of the required Pauli matrices, and I am doing this using `kron` function.

```julia
using StaticArrays

#%% We are definig the Pauli matrices here.

σˣ=SA[0 1.0; 1 0] #\sigma^x

σᶻ=SA[1.0 0; 0 -1] #\sigma^y

σʸ=-im*σᶻ*σˣ #\sigma^z

#println(σʸ,"\n",σᶻ,"\n",σˣ)

σvec=[σˣ,σʸ,σᶻ]

```

the above code stores the Pauli matrices

```julia
function σTensor(A,n::Int)
    include("Pauli_matrices.jl")
    #= Inputs: number of two level systems n (= number of tensor products + 1)
       size of the input array is N and is of the form [1,1,3,2] for N=4 
       1 represents x, 2 represents y and 3 represents z
    =#
    ret=σvec[A[1]]
    for i in 2:n
        ret=kron(ret,σvec[A[i]])
    end
    return ret
end

```

the above code does the required product.

 ![image](https://global.discourse-cdn.com/julialang/original/3X/5/3/5382acbc261a333ecc9c4ee7a6cda9f41308f697.png)  
just for tensor product of 3 matrices it takes 60kib’s  
is there a way to improve this or is it the max we can do?

since the output matrix is sparse is there a way to include that into the code?

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/27854_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [March 25, 2024, 5:41am UTC](https://discourse.julialang.org/t/how-to-speed-up-this-kronecker-multiplication/112063/2 "2024-03-25T05:41:05Z")

</div>

As usual, the question with Kronecker products is: do you really need to build the matrix explicitly? Or can you consider it as a lazy operator in your code?

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [March 25, 2024, 5:53am UTC](https://discourse.julialang.org/t/how-to-speed-up-this-kronecker-multiplication/112063/3 "2024-03-25T05:53:21Z")

</div>

> [@vikram\_backup](#):
>
> `include("Pauli_matrices.jl")`

Don’t structure your code like this, calling `include` from within the body of a function. It’s ugly.

> [@vikram\_backup](#):
>
> `using StaticArrays`

This happens within the `σTensor` function, because of the `include` …

> [@vikram\_backup](#):
>
> just for tensor product of 3 matrices it takes 60kib’s

Provide a complete reproducible example if you want help.

---

<div class="post-metadata">

**Author:** ![vikram\_backup](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vikram_backup/32/202851_2.png) [@vikram\_backup](https://discourse.julialang.org/u/vikram_backup)\
**Post date:** [March 25, 2024, 6:48am UTC](https://discourse.julialang.org/t/how-to-speed-up-this-kronecker-multiplication/112063/4 "2024-03-25T06:48:19Z")

</div>

here `A=[1,2,3]`  
I have changed the code in the following way as you have suggested

```julia
include("PauliMatrices.jl")
function σTensor(A::Array,n::Int)
    #= Inputs: number of two level systems n (= number of tensor products + 1)
       size of the input array is N and is of the form [1,1,3,2] for N=4 
       1 represents x, 2 represents y and 3 represents z
    =#
    ret=σvec[A[1]]
    for i in 2:n
        ret=kron(ret,σvec[A[i]])
    end
    return ret
end

```

and I am getting the following results for different `A`

 ![image](https://global.discourse-cdn.com/julialang/original/3X/e/9/e90cb8be45cac5fc45d684d8b2c48cd814c2e995.png)  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/7/5/7597c6fc8465d8e3cf7645732051021d2626b594.png)  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/f/3/f3f0fabcf64e19cff67280b2c58e74e11baec883.png)

1. is there a way to reduce the number of allocations to 500 or less
2. is there a way to implement sparse matrices here?  
(I am new to `Julia` and also programming in general, so please be noob friendly in your replies)

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [March 25, 2024, 6:54am UTC](https://discourse.julialang.org/t/how-to-speed-up-this-kronecker-multiplication/112063/5 "2024-03-25T06:54:41Z")

</div>

> [@vikram\_backup](#):
>
> since the output matrix is sparse is there a way to include that into the code?

The usual approach is to just use SparseArrays.jl. I don’t think that any further compiler optimization (i.e. StaticArrays.jl, type-level information) is going to help if you consider \> 7 sites.

```julia-repl
julia> using SparseArrays, LinearAlgebra, BenchmarkTools
julia> σx = sparse([0 1; 1.0 0])
2×2 SparseMatrixCSC{Float64, Int64} with 2 stored entries:
  ⋅ 1.0
 1.0 ⋅ 

julia> σy = sparse([0 -im; im 0.0])
2×2 SparseMatrixCSC{ComplexF64, Int64} with 2 stored entries:
     ⋅ 0.0-1.0im
 0.0+1.0im ⋅    

julia> σz = sparse([1 0; 0 -1.0])
2×2 SparseMatrixCSC{Float64, Int64} with 2 stored entries:
 1.0 ⋅ 
  ⋅ -1.0

julia> const σvec = [σx, σy, σz] # note the const! It is needed for performance
3-element Vector{SparseMatrixCSC{Tv, Int64} where Tv}:
 sparse([2, 1], [1, 2], [1.0, 1.0], 2, 2)
 sparse([2, 1], [1, 2], ComplexF64[0.0 + 1.0im, 0.0 - 1.0im], 2, 2)
 sparse([1, 2], [1, 2], [1.0, -1.0], 2, 2)

julia> function σTensor(A)
           n = length(A) # why did you pass it as parameter?
           ret=σvec[A[1]]
           for i in 2:n
               ret=kron(ret,σvec[A[i]])
           end
           return ret
       end
julia> @btime σTensor([1,2,3], 3)
  539.524 ns (13 allocations: 1.08 KiB)
8×8 SparseMatrixCSC{ComplexF64, Int64} with 8 stored entries:
     ⋅ ⋅ ⋅ … ⋅ 0.0-1.0im ⋅    
     ⋅ ⋅ ⋅ ⋅ ⋅ -0.0+1.0im
     ⋅ ⋅ ⋅ ⋅ ⋅ ⋅    
     ⋅ ⋅ ⋅ -0.0-1.0im ⋅ ⋅    
     ⋅ ⋅ 0.0-1.0im ⋅ ⋅ ⋅    
     ⋅ ⋅ ⋅ … ⋅ ⋅ ⋅    
 0.0+1.0im ⋅ ⋅ ⋅ ⋅ ⋅    
     ⋅ -0.0-1.0im ⋅ ⋅ ⋅ ⋅    

```

---

<div class="post-metadata">

**Author:** ![vikram\_backup](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vikram_backup/32/202851_2.png) [@vikram\_backup](https://discourse.julialang.org/u/vikram_backup)\
**Post date:** [March 25, 2024, 6:55am UTC](https://discourse.julialang.org/t/how-to-speed-up-this-kronecker-multiplication/112063/6 "2024-03-25T06:55:35Z")

</div>

The thing is I don’t know and I am new to this lazy matrices. How will I determine this?  
I was using the `Static Arrays` because I read that they are quicker than normal `Array`. My end goal is to implement time evolution of a quantum system, as efficiently as possible.  
the problem with quantum systems if there are `n` qubits in the system, the operators in the system are `2^nx2^n` complex matrix.

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [March 25, 2024, 7:02am UTC](https://discourse.julialang.org/t/how-to-speed-up-this-kronecker-multiplication/112063/7 "2024-03-25T07:02:02Z")

</div>

If you want to do time evolution of quantum systems, then a lazy approach will not be helpful. In the simplest case you want to construct your Hamiltonian and then use a eigen decomposition to compute the matrix exponentials efficiently. This will work for up to ~12-14 spins reasonably fast and easy. Beyond that you will need to employ more and more sophisticated tricks 🙂

Shameless plug: If you want to study simple spin models, you might want to have a look at my library [SpinModels.jl](https://github.com/abraemer/SpinModels.jl) with which you can construct Hamiltonians with single and two-body terms rather efficiently 🙂 (It is not in the registry for now, so you need to add it using `Pkg.add(;url="https://github.com/abraemer/SpinModels.jl")`).

---

<div class="post-metadata">

**Author:** ![vikram\_backup](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vikram_backup/32/202851_2.png) [@vikram\_backup](https://discourse.julialang.org/u/vikram_backup)\
**Post date:** [March 25, 2024, 7:21am UTC](https://discourse.julialang.org/t/how-to-speed-up-this-kronecker-multiplication/112063/8 "2024-03-25T07:21:06Z")

</div>

Thank you for your answer! 😁. I have implemented your code and it works great it brought down the allocations from 1000 to 10. I definitely will check out your SpinModels.jl. Thanks again.