Minimal example using MOSEK API in python

Question

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])

Answer 1

First observe problems that are not robust feasible or robust infeasible by by definition are nasty. Since you are testing feasibility of your system, then it is likely to be the case sometimes, because that might be baked into what you are doing.

If you have a nasty (read illposed) problem, then tiny changes in the problem can change the outcome. As far as I recall CVX has an automatic dualizer which is likely to be invoked in this case and maybe it helps for some reason. For this and other reasons CVX and CVXPY are likely to feed different formulations to Mosek and hence produce different outcomes on nasty problems.

Btw if you want to use MOSEK directly, then using Fusion is likely to be easier. However, it still requires some work. SDP are just more complicated to do.

Since your problem includes an LMI, it is normally better to dualize the problem when you use Mosek directly since Mosek does not have an automatic dualizer yet.