Джулия: более быстрое вычисление матрицы

#arrays #matrix #optimization #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 amp;amp; dy == 0 amp;amp; dz == 0
                                        mat[ii,jj]  = -(im/2) 
                                    end
                    
                                    if dx == -1 amp;amp; dy == 0 amp;amp; 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 amp;amp; dy == 0 amp;amp; dz ==  1
                                        mat[ii,jj]  = 1/2 
                                    end
                    
                                    if dx == 0 amp;amp; dy == 0 amp;amp; 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) ?

Комментарии:

1. для начала вам следует заменить EYE(Lx,Ly,Lz) на I(Lx*Ly*Lz) . Преобразование UniformScaling матрицы в плотную Matrix в принципе всегда плохая идея.

2. также ваша SINk₁ функция всегда возвращает матрицу, полную нулей. Это намеренно?

3. кроме того, COSk₃ это нулевая матрица, отличная от сдвинутой диагонали 0.5 .

4. Я заменил EYE(Lx, Ly, Lz) на I(Lx Ly Lz), но в данном случае это мало что изменило.

5. SINk₁ это также нулевая матрица, отличная от сдвинутой диагонали 0.5*im , не так ли? Вы можете проверить это для Lx = Ly = Lz = 4.

Ответ №1:

Проблема, с которой вы столкнулись, действительно проста. Вы используете как минимум O (n ^ 12) памяти. Поскольку почти все эти значения равны 0, это огромная трата. Вы почти наверняка должны использовать разреженные массивы. Это должно привести ваше время выполнения к разумным уровням

Комментарии:

1. Поскольку я не знал sparse*dense = dense , sparse*sparse = sparse , я преобразовал только EYE(Lx,Ly,Lz) , SINk₁(Lx,Ly,Lz) , COSk₃(Lx,Ly,Lz) в разреженные массивы раньше. После того, как я получил ваш ответ, я σ₁,σ₂,τ₁,τ₃ тоже попробовал его, и затем время выполнения стало довольно быстрым. Спасибо!

2. Технически вы, вероятно, могли бы получить еще один коэффициент 10, создав определенный тип матрицы, но это было бы намного сложнее.

3. Хм, на данный момент мне кажется, что это сложно. Я подумаю об этом после небольшого изучения.