# Convex optimization with constraints involving matrix inverse

I have the following convex optimization problem. I would like to ask is there any efficient way to solve it in Python? Can I use CVXOPT package? If so, any detailed instruction? Thanks a lot.

$$\min_{T\in[0, \infty]^{10}}\sum_{i=1}^{10}T_ia_i$$ subject to $$x^{\top}_j(\sum_{i=1}^{10}T_i x_ix_i^{\top})^{-1}x_j\leq b_j, \text{for all} \ j\in.$$ Here, $$\{a_1,\ldots, a_{10}\}\subset \mathbb R$$ and $$\{x_1, \ldots, x_{10}\}\subset\mathbb R^5$$. 10 and 5 are just generic constants. If the inverse does not exist, we could replace it by pseudo inverse.

• You've used $i$ for $x_{i}$ outside of the sum and $x_{i}$ inside of the sum- this is unclear and should be rewritten to use a different index for the sum. – Brian Borchers Aug 17 '19 at 18:12
• Thanks. I have updated it. – Botao Hao Aug 17 '19 at 18:53
• What makes you believe that this is a convex optimization problem? – Brian Borchers Aug 17 '19 at 19:10
• What kind of expression is $(T_ix_ix_i^T)$? If it is matrix-times-vector times vector-transposed, so something like $y_i x_i^T$, then it's a rank-one matrix, but a rank-one matrix does not have an inverse... – Wolfgang Bangerth Aug 17 '19 at 23:08
• Or are you saying that the sum of the ten rank-one matrices has an inverse? – Wolfgang Bangerth Aug 17 '19 at 23:08

Yes this can be done, see the code block below. I'm using the CVXPY library, which is maintained by Boyd's group at Stanford and wraps the CVXOPT solvers (among other things.)

Assuming we've already constructed the input data from the Problem data section, the Construct the problem and Report solution parts take around 550ms in a simple benchmark I just ran. If we increase the problem size by a factor of 2, the solution takes around 3.6s, and 11.5s if we increase by a factor of 3.

import cvxpy as cp
import numpy as np

# Problem data.
# m is matrix dimension, n is number of terms
m, n = 5, 10
X = np.random.randn(m, n)
b = np.abs(np.random.randn(n)) # constraint upper bounds
a = np.abs(np.random.randn(n)) # objective coefficients

# Construct the problem.
T = cp.Variable(n)
A = cp.Variable((m,m))
objective = cp.Minimize(a @ T)
constraints = [A == X @ cp.diag(T) @ X.T] + \
[T >= 0] + \
[cp.matrix_frac(x, A) <= bi for x, bi in zip(X.T, b)]
problem = cp.Problem(objective, constraints)

# Report solution.
result = problem.solve()
print(result)


• Thank you so much! It works well. – Botao Hao Aug 18 '19 at 14:35