#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. Хм, на данный момент мне кажется, что это сложно. Я подумаю об этом после небольшого изучения.