When computing the Hilbert-Schmidt inner product of two matrices a,b it is much faster to use dot(a,b) than the definition tr(a*b), as this avoids a matrix multiplication. I've benchmarked it numerically to confirm, and the difference is orders of magnitude. However, it doesn't work with JuMP variables, for some reason dot(a,b) is the slower one!

I've tried replacing dot(a,b) with dot(vec(a),vec(b)), and it does makes things better. For small dimensions it is still slower than tr(a*b), but as we get to dimensions where the speed matters (like 200) it gets orders of magnitude faster. So there is a workaround already, but such a low-level hack shouldn't be necessary, dot should just do the right thing.

Here is the code I used for benchmarking:

using JuMP
using LinearAlgebra

function random_matrix(d)
    x = randn(Float64, (d,d))
    return Symmetric(x+x')
end

function timejump(d)
    model = JuMP.Model()
    a = random_matrix(d)
    b = JuMP.@variable(model, [1:d,1:d] in JuMP.PSDCone())
    display(@elapsed dot(a,b))
    display(@elapsed tr(a*b))
end

function timejumpvec(d)
    model = JuMP.Model()
    a = random_matrix(d)
    b = JuMP.@variable(model, [1:d,1:d] in JuMP.PSDCone())
    display(@elapsed dot(vec(a),vec(b)))
    display(@elapsed tr(a*b))
end