Minimizing $\mathrm{trace}(S)+\mathrm{trace}(S^{-2})$ using CVX

Question

In Matlab, I would like to minimize the function

$$f(S)=\mathrm{trace}(S)+\mathrm{trace}(S^{-2})$$

where $S \in \mathcal{M}_{m,m}$ is symmetric and positive definite, which is definitely a convex function.

I tried the following code:

cvx_solver sdpt3
cvx_begin quiet
variable S(m,m) symmetric;
S == semidefinite(m);
minimize (trace(S)+trace_inv(square(S)));
cvx_end

After running this, I got the following error:

??? Error using ==> cvx.trace_inv at 9
Input must be affine.

Actually this means that what is written in the code is equivalent to have a quadratic equality in the constraint. That is, if we replace in the code $S^2$ by another $Q$ matrix, we have to add as constraint $Q=S^2$ which is convex equality since it's quadratic whereas in order to have a convex optimization problem we must have as constraints affine equalities and/or a convex inequalities. However, in reality, this is not true since my function is convex and I'm not imposing any constraint but cvx is making its own reformulation of the problem.

I even tried to solve the problem in different manner as follow:

cvx_solver sdpt3
cvx_begin quiet
variable Q(m,m) symmetric;
Q == semidefinite(m);
minimize (trace_sqrtm(Q) + trace_inv(Q));
cvx_end
S=sqrtm(Q);

but I got the following error:

??? Error using ==> cvx.plus at 83
Disciplined convex programming error:
   Illegal operation: {concave} + {convex}

which means that the function $\mathrm{trace}(S^{1/2})$ is a concave function.

So I'm wondering how can I write my problem in a different manner in order to be accepted by cvx given that my function $f(S)=\mathrm{trace}(S)+\mathrm{trace}(S^{-2})$ is convex. If there is another option rather than cvx that would work too for me.

Answer 1

I am pretty sure you can model the trace of the squared inverse by $trace(Z)$ where $Z$ and $X$ satisifes $\begin{bmatrix} X & I\\I & S\end{bmatrix} \succeq 0$ and $\begin{bmatrix}Z& X\\X&I\end{bmatrix} \succeq 0$ (Derived using Schur complements, $X$ is used to model/upper bound the inverse of $S$ and $Z$ is used for the square of $X$). Hence, solvable using standard LMI methods.

The following test computes the value for a fixed matrix $S$ (in YALMIP but the CVX version would be similar)

S = randn(3);S = S*S';
X = sdpvar(3);
Z = sdpvar(3);
BoundInverse = [X eye(3);eye(3) S]>=0;
BoundSquare = [Z X;X eye(3)] >= 0;
solvesdp([BoundInverse, BoundSquare],trace(Z))

% Is the objective equal to trace(S^-2)
trace(inv(S)^2)
double(trace(Z))
% Is Z equal to S^-2 at optimality
norm(double(Z) - S^-2)

...and to follow up on the comment, here is a model which solves a problem where $S$ is a decision variable ($S$ should be psd, but also close in 1-norm to a given matrix $C$, and the objective is $trace(S)+trace(S^{-2})$)

S = sdpvar(3);
C = randn(3);
X = sdpvar(3);
Z = sdpvar(3);
Model = [S >= 0, norm(S-C,1) <= 4];
BoundInverse = [X eye(3);eye(3) S]>=0;
BoundSquare = [Z X;X eye(3)] >= 0;
optimize([Model, BoundInverse, BoundSquare],trace(S) + trace(Z))

Answer 2

Johan has certainly done a more complete job of answering Bel's question here than I did on the CVX Forum. But I do think there is a conceptual difficulty at play here that is worth discussing.

In my response on the CVX Forum, I said that "you simply cannot expect CVX to accept any expression that you want even if you prove it is convex." Johan offers a solution that involves rewriting the objective $f(S)=\mathop{\textrm{trace}}(S)+\mathop{\textrm{trace}}(S^{-2})$ considerably. In other words, YALMIP doesn't directly handle the function either, although he does successfully show that it is possible. To my knowledge, no modeling framework does, yet.

In order for any modeling framework to "handle" a particular expression, be it $f(S)$ or any other, it needs to be able to construct a computational description of that expression that is compatible with the underlying solvers.

For instance, many solvers require code to compute the values and derivatives (typically the first two) of the expressions involved. You could supply those derivatives by hand---but for $f(S)$, this is by no means straightforward. In most cases, automatic or symbolic differentiation does the difficult work---but this only works if the expression is composed of the functions and operations that the AD/SD system already knows. I am unaware of an AD/SD engine that handles matrix expressions such as $f(S)$.

CVX and YALMIP have a different approach at their disposal: they can construct graph representations of expressions. Roughly speaking, a graph representation uses equations and inequalities involving "known" expressions to describe the geometry of a "new" inequality. A simple example is a constraint involving the absolute value: $$|x| \leq y$$ Smooth solvers will have difficulty with the non-differentiability at $x=0$. But the solution is simple and well-known: replace the constraint with two linear inequalities: $$-y \leq x \leq y$$ Graph representations let CVX and YALMIP generalize and automate these kinds of transformations. These modeling frameworks provide a variety of functions that are not available in a strict function/derivative approach. But the libraries are still finite, which means that users must construct their models using only the functions provided.

What Johan is doing in his post is showing how to construct a graph representation of $f(S)$ by hand, so that it can be solved using YALMIP or CVX. I'm reasonably certain that his derivation is correct.

Modeling frameworks hide and automate the complexity of translating models from their natural form to a machine-solvable form. This is almost always a good thing. But occasionally, it is helpful to understand how the proverbial sausage is made. With a convex optimization framework, it is not enough simply to prove that your model is convex. One must still consider the problem of translating to solvable form. CVX and YALMIP automate this process if you are willing to follow the conventions they require.

EDIT: I wrote a paper on graph representations if anyone is interested.

Answer 3

Nonlinear equality constraints are nonconvex, which is why CVX won't accept your problem.

Here is a scalar example:

Consider the equality constraint $x^{2} = y$, and consider its associated feasible set $T = \{(x,y): x^{2} - y = 0\}$. Both $(0,0)$ and $(1,1)$ are in $T$. If $T$ is a convex set, then the line $(z,z)$ for $z \in [0,1]$ should be in $T$. However, $(1/2, 1/2)$ is not in $T$, so $T$ is not convex.

In any case, the source of your problem could be that $S^{-n}$, for $n$ positive integer, and $S$ positive semidefinite, is not an atomic expression recognized by CVX as convex. If you had to prove it yourself, I doubt it's in there, and I couldn't find it in this reference. If you can add $S^{-n}$ as an atomic expression that is convex, you should be okay.

If you want to use fmincon, you should theoretically be able to, but I would think that the main problem would be enforcing the constraint that $S$ is positive semidefinite. You would likely be better off using something like YALMIP, which is a frontend for many solvers, including SDP solvers.

Answer 4

You might wish to ask Michael Grant, the author of CVX, to discipline your function in the next release of CVX.

In the mean time, try IPOPT form COIN-OR, which is a fairly robust general purpose solver for nonlinear programming problems. To make it work, you need to remove the psd constraint by replacing each occurrence of $C$ by $R^TR$, where $R$ is an upper triangular matrix with nonnegative diagonal entries.

Since IPOPT expects smooth constraints, you also need to replace the nonsmooth L1 constraint $\|S−C\|_1\le \alpha$ (assuming the 1-norm to be the column sum norm) by the smooth constraints $$-X\le R^TR-C \le X,~~~ e^TX\le \alpha e,$$ where $e$ is the all-one vector, and $X$ is a matrix of new variables.

Answer 5

As Michael Grant (principal author of CVX) responds in his answer here, you're trying to shove a general convex programming problem into CVX, which is designed to handle disciplined convex programming problems.

This is why you are not getting anything reasonable out of CVX, you must first try to reformulate your problem into disciplined convex form.

As Geoff suggests, you might have luck with more general solvers, particularly if you know your problem is convex and you are not solving a large problem.