Optimization of a Trace Operator with Symmetric, PSD and Unit Trace Constraints

Question

I'm trying to solve a simple optimisation problem, we want to have a complex valued Hermitian matrix as its variable (topic is quantum mechanics).

using Convex  #load the optimization solvers
using SCS

# define pauli-y+ projector
# by construction a positive operator valued hermitian matrix
y_plus = [1,im]/sqrt(2)
My0 = y_plus*y_plus'

# define the variable; a 2x2 density matrix
rho = Variable(2, 2)
problem.constraints += [rho == rho']     # hermitian
problem.constraints += [trace(rho) == 1] # unit trace
problem.constraints += [rho in :SDP]     # positive definite
    

# define the objective
problem = maximize(trace(rho*My0))

# solve
solve!(problem,SCSSolver(verbose=false))

problem.optval

The trouble is, Julia/JuMP/Convex.jl all give errors when it comes to

maximize(trace(rho*My0))

Since the trace of rho*My0 can in principle be complex, however we should be assured that rho*My0 is real given the constraints on rho and My0.

How to deal with these issues? There might be a simple solution. Worst case we probably have to split the real and imaginary parts.

Claudeclaudel

Given that adding complex support is an open project for summer of code applicants, you're at the edge of supported behavior. You could try something like maximize(trace(convert(Array{Float64}, rho*My0))) or maximize(trace(real(rho*My0))), but I'd be surprised if that works. – Kobe

Why not just maximize the abs(trace(rho*My0)))? If the trace is real, it's the same. – Bacchant

@DanGetz, What if the maximum value is negative. Won't abs() mean the minimum will be obtained? In this case the My0 and rho are PSD matrices hence the trace of their multiplication always positive (math.stackexchange.com/questions/888677) so this will work in this case. But not for any 2 matrices under those conditions. – Aton

Moreover, I think it won't hold the DCP conditions. – Aton

Answer 1

You may write the problem as:

$$ \begin{align*} \arg \min_{\boldsymbol{A}} \quad & \operatorname{Tr} \left( \boldsymbol{A} \boldsymbol{B} \right), ; \boldsymbol{B} \in \mathcal{S} \ \text{subject to} \quad & \begin{aligned} \boldsymbol{A} & \in \mathcal{S}_{+} \ \operatorname{Tr} \left( \boldsymbol{A} \right ) & = 1 \end{aligned} \end{align*} $$

Where the set is the set of symmetric positive semidefinite matrices.

The problems are equivalent as you may always use $\boldsymbol{B} = -\boldsymbol{M}$ from you problem above.

The projection steps are 3:

1. Projection to Symmetric Matrices Set.
2. Projection to PSD Matrices Set.
3. Projection to Matrices with Unit Trace.

Since the projection to each set is known, we can solve the problem pretty easily.
Yet, since the set to project onto is not a sub space, we can't just project on each iteratively, we must use a solver for the projection (See Orthogonal Projection onto the Intersection of Convex Sets).

The way I see it, there are 3 options:

1. Projected Gradient Descent
   Gradient step for the objective function and then solving the projection by the methods in Orthogonal Projection onto the Intersection of Convex Sets.
2. 3 Terms ADMM
   Unite the projection into SPSD into a single function. Then you have 3 terms ADMM (Where on of the projections will require inner iterations). You may use PD3O or just the classic 3 Operators ADMM (See A Three Operator Splitting Scheme and its Optimization Applications). Specifically, for this case, A Three Operator Splitting Perspective of a Three Block ADMM for Convex Quadratic Semi Definite Programming and Extensions.
3. 4 Terms ADMM
   I haven't found a paper which shows the guarantees for convergence in this case. But it is worth a try.

I implemented Projected Gradient Descent (1) and PD3O for this problem in Julia:

The reference was the solution with Convex.jl.
In the case of @WillemHekman, the matrix $ \boldsymbol{B} $ is a PSD matrix. Hence the trace $ \opertaorname{Tr}\left( \boldsymbol{A} \boldsymbol{B} \right) $ is guaranteed to be non negative (See The Non Negative of Matrix Multiplication of 2 PSD Matrices).
Which means I could minimize the abs() of the objective function and not only the actual value. This means it will work for complex matrices.

The implementations does not rely on that and can solve the problem for arbitrary matrix $ \boldsymbol{B} $ as long as it is square.

The Julia code is available at my StackOverflow Q35813091 GitHub Repository (Look at the StackOverflow\Q35813091 folder).

Remark: The motivation of this answer is a solver which is independent of general convex solvers as Convex.jl and JuMP.jl.

Aton