Julia: faster matrix calculation

Question

In my work, I have to deal with large size matrices.
For example, I use the following matrices.

using LinearAlgebra

#Pauli matrices
σ₁ = [0 1; 1 0]
σ₂ = [0 -im; im 0]
τ₁ = [0 1; 1 0]
τ₃ = [1 0; 0 -1]

#Trigonometric functions in real space
function EYE(Lx,Ly,Lz)
 　　　
    N = Lx*Ly*Lz
   
    mat = Matrix{Complex{Float64}}(I, N, N)
    
    return mat
end


function SINk₁(Lx,Ly,Lz)
 　　　
    N = Lx*Ly*Lz
   
    mat = zeros(Complex{Float64},N,N) 
    
    for ix = 1:Lx        
        for iy = 1:Ly 
            for iz = 1:Lz 
                for dx in -1:1             
                 jx = ix + dx 
                 jx += ifelse(jx > Lx,-Lx,0) 
                 jx += ifelse(jx < 1,Lx,0)
                    
                    for dy in -1:1
                    jy = iy + dy
                    jy += ifelse(jy > Ly,-Ly,0) 
                    jy += ifelse(jy < 1,Ly,0)   
                        
                        for dz in -1:1
                        jz = iz + dz
                        
                        
                            ii = (iz-1)*Lx*Ly + (ix-1)*Ly + iy 
                            jj = (jz-1)*Lx*Ly + (jx-1)*Ly + jy

                                if 1 <= jz <= Lz
                                                         
                                    if dx == +1 && dy == 0 && dz == 0
                                        mat[ii,jj] += -(im/2)　
                                    end
                    
                                    if dx == -1 && dy == 0 && dz == 0
                                        mat[ii,jj] += im/2
                                    end
                                       
                                end    
                    
                        end
                    end
                end
            end
        end
    end
    
    return mat
end


function COSk₃(Lx,Ly,Lz)
 　　　
    N = Lx*Ly*Lz
   
    mat = zeros(Complex{Float64},N,N) 
    
     for ix = 1:Lx
        for iy = 1:Ly 
            for iz = 1:Lz 
                for dx in -1:1             
                 jx = ix + dx 
                 jx += ifelse(jx > Lx,-Lx,0) 
                 jx += ifelse(jx < 1,Lx,0)
                    
                    for dy in -1:1
                    jy = iy + dy
                    jy += ifelse(jy > Ly,-Ly,0) 
                    jy += ifelse(jy < 1,Ly,0)   
                        
                        for dz in -1:1
                        jz = iz + dz
                        
                        
                            ii = (iz-1)*Lx*Ly + (ix-1)*Ly + iy 
                            jj = (jz-1)*Lx*Ly + (jx-1)*Ly + jy

                                if 1 <= jz <= Lz
                                                         
                                    if dx == 0 && dy == 0 && dz == +1
                                        mat[ii,jj] += 1/2　
                                    end
                    
                                    if dx == 0 && dy == 0 && dz == -1
                                        mat[ii,jj] += 1/2
                                    end
                                       
                                end    
                    
                        end
                    end
                end
            end
        end
    end
    
    return mat
end

Then, I calculate

kron(SINk₁(Lx,Ly,Lz),kron(σ₁,τ₁)) + kron(EYE(Lx,Ly,Lz) + COSk₃(Lx,Ly,Lz),kron(σ₂,τ₃)) 

This calculation, however, takes much times for large Lx,Ly,Lz;

Lx = Ly = Lz = 15
@time kron(SINk₁(Lx,Ly,Lz),kron(σ₁,τ₁)) + kron(EYE(Lx,Ly,Lz) + COSk₃(Lx,Ly,Lz),kron(σ₂,τ₃)) 

4.692591 seconds (20 allocations: 8.826 GiB, 6.53% gc time)

Lx = Ly = Lz = 20
@time kron(SINk₁(Lx,Ly,Lz),kron(σ₁,τ₁)) + kron(EYE(Lx,Ly,Lz) + COSk₃(Lx,Ly,Lz),kron(σ₂,τ₃)) 

52.687861 seconds (20 allocations: 49.591 GiB, 2.69% gc time)

Are there faster ways to calculate Kronecker product, Addition or more proper definitions of EYE(Lx,Ly,Lz), SINk₁(Lx,Ly,Lz), COSk₃(Lx,Ly,Lz)?

Answer 1

The problem you're having is really simple. You are using at least O(n^12) memory. Since almost all of these values are 0, this is a huge waste. You should almost certainly be using Sparse arrays. This should bring your runtime to reasonable levels