arrays - Julia：更快的矩阵计算

问题描述

在我的工作中，我必须处理大尺寸矩阵。
例如，我使用以下矩阵。

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

然后，我计算

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

但是，对于较大的 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)

是否有更快的方法来计算克罗内克积、加法或更合适EYE(Lx,Ly,Lz)的SINk₁(Lx,Ly,Lz)定义COSk₃(Lx,Ly,Lz)？

标签： arraysmatrixoptimizationjulia

arrays - Julia：更快的矩阵计算

问题描述

解决方案

推荐阅读