symmetric_projector for arbitrary types

Project.toml

@@ -8,7 +8,6 @@ GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 Hypatia = "b99e6be6-89ff-11e8-14f8-45c827f4f8f2"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

src/Ket.jl

@@ -6,7 +6,7 @@ import LinearAlgebra as LA
 import SparseArrays as SA
 import Nemo
 import JuMP
-import MathOptInterface as MOI
+const MOI = JuMP.MOI
 import Hypatia, Hypatia.Cones
 
 include("basic.jl")

src/basic.jl

@@ -313,22 +313,28 @@ export orthonormal_range
     symmetric_projection(dim::Integer, n::Integer; partial::Bool=true)
 
 Orthogonal projection onto the symmetric subspace of `n` copies of a `dim`-dimensional space. By default (`partial=true`)
-it returns amatrix (say, `P`) with columns forming an orthogonal basis for the symmetric subspace. If `partial=false`, then it
-returns the actual projection `P * P'`.
+it returns an isometry (say, `V`) encoding the symmetric subspace. If `partial=false`, then it
+returns the actual projection `V * V'`.
 
 Reference: [Watrous' book](https://cs.uwaterloo.ca/~watrous/TQI/), Sec. 7.1.1
 """
-function symmetric_projection(dim::Integer, n::Integer; partial::Bool = true)
-    P = SA.spzeros(dim^n, dim^n)
+function symmetric_projection(::Type{T}, dim::Integer, n::Integer; partial::Bool = true) where {T}
+    if T <: SA.CHOLMOD.VTypes #sparse qr decomposition fails for anything other than Float64 or ComplexF64
+        P = SA.spzeros(T, dim^n, dim^n)
+    else
+        P = zeros(T, dim^n, dim^n)
+    end
     perms = Combinatorics.permutations(1:n)
     for perm in perms
-        P += permutation_matrix(dim, perm; sp=true)
+        P .+= permutation_matrix(dim, perm; sp=true)
     end
-    P /= length(perms)
+    P ./= length(perms)
     if partial
-        P = orthonormal_range(P)
-        size(P, 2) != binomial(n + dim - 1, dim - 1) && throw(AssertionError("Rank computation failed"))
+        V = orthonormal_range(P)
+        size(V, 2) != binomial(n + dim - 1, dim - 1) && throw(AssertionError("Rank computation failed"))
+        return V
     end
-    P
+    return P
 end
-export symmetric_projection
+export symmetric_projection
+symmetric_projection(dim::Integer, n::Integer; partial::Bool = true) = symmetric_projection(Float64, dim, n; partial)

src/entanglement.jl

@@ -210,11 +210,11 @@ function entanglement_dps(
 
     # Dimension of the extension space w/ bosonic symmetries: AA' dim. + `n` copies of BB'
     sym_dim = d1 * binomial(n + d2 - 1, d2 - 1)
-    P = kron(LA.I(d1), symmetric_projection(d2, n; partial=true)) # Bosonic subspace projector
+    V = kron(LA.I(d1), symmetric_projection(T, d2, n; partial=true)) # Bosonic subspace isometry
 
     model = JuMP.GenericModel{_solver_type(T)}()
     JuMP.@variable(model, Q[1:sym_dim, 1:sym_dim] in JuMP.HermitianPSDCone())
-    JuMP.@expression(model, lifted, P * Q * P')
+    JuMP.@expression(model, lifted, V * Q * V')
     JuMP.@expression(model, reduced, partial_trace(lifted, collect(3:n+1), dims))
 
     if !witness
@@ -256,4 +256,4 @@ function entanglement_dps(
     end
     return obj
 end
-export entanglement_dps
+export entanglement_dps