using Convex
n = 20
P = randn(n,n) + im*randn(n,n)
P = P*P'
Q = randn(n,n) + im*randn(n,n)
Q = Q*Q'
Z = ComplexVariable(n,n)
objective = 0.5*real(trace(Z+Z'))
constraint = [P Z;Z' Q] ⪰ 0
problem = maximize(objective,constraint)
solve!(problem)
computed_fidelity = evaluate(objective)

# Verify that computer fidelity is equal to actual fidelity
P1,P2 = eig(P)
sqP = P2 * diagm([p1^0.5 for p1 in P1]) * P2'
Q1,Q2 = eig(Q)
sqQ = Q2 * diagm([q1^0.5 for q1 in Q1]) * Q2'

actual_fidelity = sum(svd(sqP * sqQ)[2])