perform embedding of an operator in a joint hilbert space

Author: alexander papageorge

src/operators.jl

@@ -6,8 +6,10 @@ export AbstractOperator, DataOperator, length, basis, dagger, ishermitian, tenso
 
 import Base: ==, +, -, *, /, ^, length, one, exp, conj, conj!, transpose
 import LinearAlgebra: tr, ishermitian
+import SparseArrays: sparse
 import ..bases: basis, tensor, ptrace, permutesystems,
             samebases, check_samebases, multiplicable
+# import ..operators_sparse: SparseOperator
 import ..states: dagger, normalize, normalize!
 
 using ..sortedindices, ..bases, ..states
@@ -131,11 +133,72 @@ function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
 
     return tensor(out...)
 end
+
+
+"""
+    embed(basis1[, basis2], indices::Vector, operators::Vector)
+
+Embed operator acting on a joint Hilbert space where missing indices are filled up with identity operators.
+"""
+function embed(basis_l::CompositeBasis, basis_r::CompositeBasis,
+               indices::Vector{Int}, op::T) where T<:AbstractOperator
+    N = length(basis_l.bases)
+    @assert length(basis_r.bases) == N
+
+    @assert reduce(tensor, basis_l.bases[indices]) == op.basis_l
+    @assert reduce(tensor, basis_r.bases[indices]) == op.basis_r
+
+    index_order = [idx for idx in 1:length(basis_l.bases) if idx ∉ indices]
+    all_operators = AbstractOperator[identityoperator(T, basis_l.bases[i], basis_r.bases[i]) for i in index_order]
+
+    for idx in indices[end:-1:1]
+        pushfirst!(index_order, idx)
+    end
+    pushfirst!(all_operators, op)
+
+    sortedindices.check_indices(N, indices)
+
+    # Create the operator.
+    permuted_op = tensor(all_operators...)
+
+    # Create a copy to fill with correctly-ordered objects (basis and data).
+    unpermuted_op = copy(permuted_op)
+
+    # Create the correctly ordered basis.
+    unpermuted_op.basis_l = permutesystems(permuted_op.basis_l, index_order)
+    unpermuted_op.basis_r = permutesystems(permuted_op.basis_r, index_order)
+
+    # Reorient the matrix to act in the correctly ordered basis.
+    # Get the dimensions necessary for index permuting.
+    dims_l = [b.shape[1] for b in basis_l.bases]
+    dims_r = [b.shape[1] for b in basis_r.bases]
+
+    # Get the order of indices to use in the first reshape. Julia indices go in
+    # reverse order.
+    expand_order = index_order[end:-1:1]
+    # Get the dimensions associated with those indices.
+    expand_dims_l = dims_l[expand_order]
+    expand_dims_r = dims_r[expand_order]
+
+    # Prepare the permutation to the correctly ordered basis.
+    perm_order_l = [indexin(idx, length(dims_l):-1:1)[1] for idx in expand_order]
+    perm_order_r = [indexin(idx, length(dims_r):-1:1)[1] for idx in expand_order]
+
+    # Perform the index expansion, the permutation, and the index collapse.
+    M = (reshape(permuted_op.data, tuple([expand_dims_l; expand_dims_r]...)) |>
+         x -> permutedims(x, [perm_order_l; perm_order_r .+ length(dims_l)]) |>
+         x -> sparse(reshape(x, (prod(dims_l), prod(dims_r)))))
+
+    unpermuted_op.data = M
+
+    return unpermuted_op
+end
 embed(basis_l::CompositeBasis, basis_r::CompositeBasis, index::Int, op::AbstractOperator) = embed(basis_l, basis_r, Int[index], [op])
 embed(basis_l::CompositeBasis, basis_r::CompositeBasis, index::Vector{Int}, op::AbstractOperator) = embed(basis_l, basis_r, index, [op])
 embed(basis::CompositeBasis, index::Int, op::AbstractOperator) = embed(basis, basis, Int[index], [op])
 embed(basis::CompositeBasis, indices::Vector{Int}, operators::Vector{T}) where {T<:AbstractOperator} = embed(basis, basis, indices, operators)
 embed(basis::CompositeBasis, index::Vector{Int}, op::AbstractOperator) = embed(basis, basis, index, [op])
+embed(basis::CompositeBasis, indices::Vector{Int}, op::AbstractOperator) = embed(basis, basis, indices, op)
 
 """
     embed(basis1[, basis2], operators::Dict)
@@ -202,7 +265,31 @@ normalize!(op::AbstractOperator) = throw(ArgumentError("normalize! is not define
 
 """
     expect(op, state)
+index_order = [idx for idx in 1:length(dims) if idx ∉ indices]
+
+eyedim = prod(dims[index_order])
+
+for idx in indices[end:-1:1]
+    pushfirst!(index_order, idx)
+end
+
+alldim = prod(dims[index_order])
+
+index_order
+
+expand_order = [dims[idx] for idx in index_order][end:-1:1]
+
+perm_order = [indexin(idx, expand_order)[1] for idx in dims[end:-1:1]]
+
+println(index_order)
+println(expand_order)
+println(perm_order)
+
+my_matrix = kron(m, eye(eltype(m), eyedim))
 
+M = (reshape(my_matrix, tuple([expand_order; expand_order]...)) |>
+     x -> permutedims(x, [perm_order; perm_order .+ length(dims)]) |>
+     x -> sparse(reshape(x, (prod(dims), prod(dims)))))
 Expectation value of the given operator `op` for the specified `state`.
 
 `state` can either be a (density) operator or a ket.

test/test_operators.jl

@@ -62,6 +62,7 @@ op_test3 = test_operators(b1 ⊗ b2, b2 ⊗ b1, randoperator(b, b).data)
 
 @test embed(b, b, [1,2], op) == embed(b, [1,2], op)
 @test embed(b, Dict{Vector{Int}, SparseOperator}()) == identityoperator(b)
+@test_throws AssertionError embed(b1⊗b2, [2], [op1])
 
 b_comp = b⊗b
 @test embed(b_comp, [1,3,4], [op1,op]) == dense(op1 ⊗ one(b2) ⊗ op)