# Minimal example using MOSEK API in python

I want to solve (simplified version) \begin{equation*} \begin{aligned} & \underset{}{\text{find}} & & X\in\mathbb{S}^{n}_{+}, x \in \mathbb{R}^{m}, \nu \in \mathbb{R}, \lambda\geq 0\\ & \text{subject to} & & E^{\top} X E + \lambda A + \nu B \succeq 0, \\ & & & Fx + \lambda a + \nu b = 0, \\ \end{aligned} \end{equation*} where \begin{align*} E &\in \mathbb{R}^{n\times p},\\ A &\in \mathbb{S}^{p}, \\ B &\in \mathbb{S}^{p}, \\ F &\in \mathbb{R}^{q \times m}, \\ a &\in \mathbb{R}^{q}, \\ b &\in \mathbb{R}^{q} \end{align*} are fixed.

I use Python with CVXPY and the MOSEK solver, which I find unreliable. Matlab with CVX and MOSEK is more reliable. I recently found that there is a MOSEK API for Python. However, the tutorial page does not directly cover this use case. I'm looking for a minimal example covering this case.

Edit: The following code passes without errors. However, contrary to Erling's mention below, dualization was not needed. Is this a correct implementation?

from mosek.fusion import *
import mosek.fusion.pythonic
import numpy as np

# Define the dimensions
n = 4
p = 3
m = 2
q = 2

# Generate random fixed data for the problem
np.random.seed()
E = np.random.randn(n, p)
A = np.random.randn(p, p)
A = (A + A.T) / 2  # Make A symmetric
B = np.random.randn(p, p)
B = (B + B.T) / 2  # Make B symmetric
F = np.random.randn(q, m)
a = np.random.randn(q)
b = np.random.randn(q)

# Create a new model
with Model("example") as M:
# Define the variables
X = M.variable("X", Domain.inPSDCone(n))
x = M.variable("x", m, Domain.unbounded())
nu = M.variable("nu", 1, Domain.unbounded())
lam = M.variable("lambda", 1, Domain.greaterThan(0.0))

# PSD constraint: E^T X E + lambda * A + nu * B >= 0
M.constraint("psd_constraint", E.T @ X @ E + lam[0] * A + nu[0] * B, Domain.inPSDCone())

# Linear equality constraint: Fx + lambda * a + nu * b = 0
M.constraint("linear_constraint", F @ x + lam[0] * a + nu[0] * b, Domain.equalsTo(0.0))

# Solve the problem
M.solve()

# Retrieve and print the solution
X_sol = X.level()
x_sol = x.level()
nu_sol = nu.level()
lam_sol = lam.level()

print("Solution X:\n", np.array(X_sol).reshape((n, n)))
print("Solution x:", x_sol)
print("Solution nu:",nu_sol[0])
print("Solution lambda:", lam_sol[0])

• What did you find unreliable in CVXPY/MOSEK? Commented May 23 at 18:18
• @MarkL.Stone The way CVXPY rewrites the feasibility problem to the MOSEK solver leads to worse numerical stability compared to CVX to MOSEK. As I sweep over some data parameters of the feasibility problem, for which I know that feasibility holds by known theorems, (CVX, MOSEK) can recover it to a greater extent than (CVXPY, MOSEK). Of course, the interesting case is the behavior in cases that are not covered by known theory. Commented May 23 at 18:33
• The difference between CVX and CVXPY might come down to order of constraints or minor transformation differences; and if that is enough to sway the results between feasibility and infeasibility, perhaps something is deficient in your model's input data, or the model itself. Have you looked at the beginning of the MOSEK output to see if there are warnings about large or small magnitude numbers? If so, you should rescale the input data (change units) so that all non-zero input data is within a small number of orders of magnitude of 1. Commented May 23 at 18:39
• @MarkL.Stone, I believe the input variables are reasonable. You can have a look at a bad case here: ibb.co/1rbCm2K Commented May 23 at 19:04
• Your code looks good, it is analogous to what you would write in CVXPY and CVX. Regarding dualization a good place to start is docs.mosek.com/latest/faq/… and the links therein. In your case I don't think it will matter since you don't really have an LMI but constraints binding two semidefinite variables ($X$ and the one implicit in $\succeq$.) Commented May 28 at 18:07