How to solve the following SDP with Python?

Question

Supposing that $\{B_{ij}\}_{i,j}$ are all Hermitian matrices and $\{c_{ijk}\}_{i,j,k}$ are all real numbers, the corresponding SDP(Semidefinite Programming) problem is as follows: $$ \begin{aligned} &\!\!\min_{\{X_{ij} \;\!\}}\quad \mathrm{tr}\sum_{i,j}X_{ij}B_{ij} \\ &\mathrm{s.t.}\quad\sum_{i,j} c_{ijk}X_{ij}\succeq0,\;\forall k \\ &\qquad\;\, \mathrm{tr}\sum_{i,j,k}c_{ijk}X_{ij}=1 \end{aligned} $$ where $\{X_{ij}\}_{i,j}$ are all Hermitian matrices.

How should I solve it by programming?

I know that semi-definite programming is a kind of convex optimization problem, and Python has two commonly used libraries for solving convex optimization: cvxpy and cvxopt.

I also checked the corresponding documents, but there is no such form of SDP in the document, so how do I program it to solve it?

Can you give me some ideas?

---

Update:

Now, I will restate my question with the original notation in the paper:

Assuming that $x$ can take values $0,1,\cdots,m-1$, and $a$ can take values $0,1,\cdots,n-1$, now we define a function:

$$ \lambda(\cdot): \{0,1,\cdots,m-1\}\to\{0,1,\cdots,n-1\},\quad x\mapsto a $$ Thus the function $\lambda$ can be uniquely determined by the tuple $(a_{x=0},a_{x=1},\cdots,a_{x=m-1})$. It can be seen from the definition that there are $n^m$ such functions.

For example, let $m=3,n=2$, $\lambda=(1,0,1)$, it means that $$ \lambda(0)=1,\quad\lambda(1)=0,\quad \lambda(2)=1 $$ Then we define a function again: $$ D(a|x,\lambda)= \delta_{a,\lambda(x)}= \begin{cases} 1,&\lambda(x)=a \\ 0,&\lambda(x)\neq a\\ \end{cases} $$ Using the example mentioned above, we know $$ D(1|0,\lambda)=1,\quad D(0|1,\lambda)=1,\quad D(1|2,\lambda)=1 \\ D(a|x,\lambda)=0,\quad \text{others} $$ Now, I assume you have understood the definition of $D(a|x,\lambda)$. Next, I will introduce the SDP that appears in the paper.

Let $\{F_{a|x}\}_{a,x}$ be a series of Hermitian matrices, given a series of matrices $\{\sigma_{a|x}\}_{a,x}$ and a series of distributions $\{D(a|x,\lambda)\}_{\lambda}$, the corresponding SDP is described as follows: $$ \begin{aligned} &\!\!\min_{\{F_{a|x} \;\!\}}\quad \mathrm{tr}\sum_{a,x}F_{a|x}\sigma_{a|x} \\ &\mathrm{s.t.}\quad\sum_{a,x} D(a|x,\lambda)F_{a|x}\succeq0,\;\forall \lambda\\ &\qquad\;\, \mathrm{tr}\sum_{a,x,\lambda}D(a|x,\lambda)F_{a|x}=1 \end{aligned} $$

What is done between $F_{a|x}$ and $\sigma_{a|x}$ is matrix multiplication.

Note that the matrices mentioned above are complex matrices.

Answer 1

Edit: I looked at the paper you linked in the comment. I'm no expert in quantum computing, but it seems like a hot mess. The notation is certainly not clear to me, but if you can clear it up then I'll edit this formulation. The following is my best guess.

---

I believe you have a number of typos. First, I believe you are conflating equivalent expressions for the Frobenius inner product.

• The first expression is $\text{tr}(B^{T} X)$

• The second expression is $\sum_{i,j = 1}^{n} B_{i,j} X_{i,j}$, where $X_{i,j}$ is the $(i,j)$ entry of $X$.

Second, I think you mean to optimize over $X$, not $B$. Otherwise your constraints have no variables.

If you can edit the question to correct these issues then I will edit my answer. Otherwise I'll assume that by $\text{tr}(\sum_{i,j=1}^{n} B_{i,j} X_{i,j})$ you mean Frobenius inner product, which seems the most reasonable way to interpret your notation. I'll also assume that you mean the hadamard product in your PSD constraints, which also seems like the most reasonable interpretation. That is, I'm solving

$$\min_{X \succeq 0} \; \text{tr}(B^{T} X)$$ $$\text{s.t.} \quad C_{k} \odot X \succeq 0 \quad \forall k = 1,...,m$$ $$ \quad \quad \text{tr}(C_{k}^{T} X) = 1 \quad \forall k = 1,...,m$$

Here's my code to solve.

import cvxpy as cvx

def solve_Problem(B, Cs):
    '''
    B is an n x n numpy array
    Cs is a list such that each Cs[i] is an n x n numpy array 
    each of which appears in one of the constraints'''
   
    n = B.shape[0]
    X = cvx.Variable((n,n), hermitian=True)
    obj = cvx.sum(cvx.multiply(B, X)) #this is Frob. product
    cons = [cvx.sum(cvx.multiply(Cs[i], X)) == 1 for i in range(len(Cs))]
    cons += [cvx.multiply(Cs[i], X)) >> 0 for i in range(len(Cs))] #psd constraints
    prob = cvx.Problem(cvx.Minimize(obj), cons)
    prob.solve()
    return X.value

2nd Edit: In response to the clarification of the problem, here's how one could approach it with cvxpy. I'll assume here that you have access to python lists Lambda, A, and X which enumerate the $a$ values, $x$ values, and $\lambda$ values respectively. I'll assume Sigma is a 2D list of $n\times n$ numpy arrays. This list is of size len(A) x len(X). Similarly I'll assume D is a 3D list of scalars. This list is of size len(A) x len(X) x len(Lambda).

import cvxpy as cvx

n = Sigma[0][0].shape[0] #get size of Sigmas
Fs = [[cvx.Variable((n,n), hermitian=True) for _ in range(len(X))] for _ in range(len(A))]
obj = cvx.trace(sum(Fs[i][j] @ Sigmas[i][j] for i in range(len(A)) for j in range(len(X))))
#trace is linear
cons = [sum(cvx.trace(D[i][j][k]*sum(Fs[i][j] for k in range(len(Lambda))) for i in range(len(A)) for j in range(len(X))) == 1]
for k in range(len(Lambda)):
    cons += [sum(D[i][j]*Fs[i][j] for i in range(len(A)) for j in range(len(X))) >> 0]
prob = cvx.Problem(cvx.Minimize(obj), cons)
prob.solve()

#now access your solution with Fs[i][j].value for i in range(len(A))
#and j in range(len(X))

Note that, despite being a convex program, this is not a tractable problem since the number of constraints depends exponentially on len(A) and len(X). That is, you shouldn't expect to solve this problem in a polynomial number of floating point operations, where the polynomial is in variables len(A), len(X).