1

I got some code from a friend which reads like this:

import numpy as np
import cvxpy as cp

def tomoSolve(m, P,prin=True):
    # get d from projector length
    d = int(np.sqrt(P.shape[1]))

    # n.b. due to cvxpy's implementation of variable special properties,
    # we must define 2 variables with the positive semi-definite and
    # complex properties respectively and constrain them to be equal

    # initialise target variable with Hermitian constraint
    X = cp.Variable((d, d), hermitian=True)

    # create target var with complex value constraint
    x = cp.Variable((d, d), complex=True)

    # define objective
    obj = cp.Minimize(cp.norm((P @ cp.reshape(X, (d**2, 1)))-m))

    # herm&complex, PSD,  unit trace constraint
    const = [X == x, X >> 0, cp.trace(X)-1 == 0]

    # construct problem in cvxpy
    prob = cp.Problem(obj, const)

    # solve problem and update variable fields
    # prob.solve(verbose=True) # for statistics output
    prob.solve(verbose=False) # for no print
    # print exit status
    if prin:
        print(f'Status: {prob.status}')

    # return solution
    return X.value.T

Now I would like to implement another constraint which limits the resulting matrix to Rank 1. I tried cp.real(cp.trace(cp.power(X, 2))) == 1 but it doesn't work cause ValueError: Arguments to power cannot be complex.. Anyone got any idea how to implement this?

Thanks a lot, Markus

Markus
  • 37
  • 6
  • Such questions should be at Math SE or SciComp SE. Stack Overflow does not allow MathJax and your question is about mathematics, not coding. – Rodrigo de Azevedo Jun 24 '20 at 06:28
  • The desired equality constraint is non-convex. CVXPY is for convex problems. The rank constraint is easy to handle by replacing the matrix variable by the outer product of the same vector variable. – Rodrigo de Azevedo Jun 24 '20 at 06:30
  • I know that the way I wrote it is mathematically correct, but not in terms of what I can use as a constraint in cvxpy, that's why I put it here. Sorry if that is not correct. The entity I calculate is called purity in physics, and if the purity is 1 the rank is 1. – Markus Jun 24 '20 at 06:39
  • will try to implement it as a vector but I'm totally new to this so idk if your answer makes it easier for me. Thanks anyway, was worth a shot :) – Markus Jun 24 '20 at 06:41
  • If you posted your question in an SE site that allows MathJax, perhaps you would be able to communicate your ideas succinctly and not have to post Python code. – Rodrigo de Azevedo Jun 24 '20 at 06:44

0 Answers0