kbarros
commented
1 year ago Here is some code to construct the spin operators as tensor products, and to use SVD to compress operators.
using Sunny, LinearAlgebra
# Produces a matrix representation of a tensor product of operators, C = A⊗B.
# Like built-in `kron` but with permutation. Returns C_{acbd} = A_{ab} B_{cd}.
function kron_operator(A::AbstractMatrix, B::AbstractMatrix)
TS = promote_type(eltype(A), eltype(B))
C = zeros(TS, size(A,1), size(B,1), size(A,2), size(B,2))
for ci in CartesianIndices(C)
a, c, b, d = Tuple(ci)
C[ci] = A[a,b] * B[c,d]
end
return reshape(C, size(A,1)*size(B,1), size(A,2)*size(B,2))
end
# Use SVD to find the decomposition D = ∑ₖ Aₖ ⊗ Bₖ, where Aₖ is N₁×N₁ and Bₖ is
# N₂×N₂. Returns the list of matrices [(A₁, B₁), (A₂, B₂), ...].
function svd_tensor_expansion(D::Matrix{T}, N1, N2) where T
@assert size(D, 1) == size(D, 2) == N1*N2
D = permutedims(reshape(D, N1, N2, N1, N2), (1,3,2,4))
(; S, U, V) = svd(reshape(D, N1*N1, N2*N2))
ret = []
for (k, σ) in enumerate(S)
if abs(σ) > 1e-12
u = reshape(U[:, k], N1, N1)
v = reshape(V[:, k], N2, N2)
push!(ret, (σ*u, conj(v)))
end
end
return ret
end
# Returns the spin operators for two sites, with Hilbert space dimensions N₁ and
# N₂, respectively.
function spin_pair(N1, N2)
S1 = Sunny.spin_matrices(N1)
S2 = Sunny.spin_matrices(N2)
I1 = Ref(Matrix(I, N1, N1))
I2 = Ref(Matrix(I, N2, N2))
return (kron_operator.(S1, I2), kron_operator.(I1, S2))
end
N1 = 3
N2 = 3
# Check Kronecker product of operators
A1 = randn(N1,N1)
A2 = randn(N1,N1)
B1 = randn(N2,N2)
B2 = randn(N2,N2)
@assert kron_operator(A1, B1) * kron_operator(A2, B2) ≈ kron_operator(A1*A2, B1*B2)
# Check SVD decomposition
S1, S2 = spin_pair(N1, N2)
B = (S1' * S2)^2 # biquadratic interaction
D = svd_tensor_expansion(B, N1, N2) # a sum of 9 tensor products
@assert sum(kron_operator(d...) for d in D) ≈ B # consistency check
Currently we support scalar biquadratic interactions in
:dipolemode. Eventually we should allow this also for:SUNmode. The latter would traditionally require thinking about a basis for SU(N) generators that explicitly defines spin operators $\hat S$ as well as operators $\hat Q$ quadratic in spin.Rather than hand-coding these generators, it might be better to work entirely with polynomials of the spin operators, and their tensor products, e.g. $S_j = S \otimes I$ and $S_k = I \otimes S$ to represent spins for sites $j$ and $k$ in the context of the bond $(j,k)$. This approach would generalize to any order exchange interaction, as well as to entangled units.