7
$\begingroup$

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.

$\endgroup$
  • 4
    $\begingroup$ Bel, thanks for bringing the question here and welcome to the site! In the future, please do reference other places you have asked this question as well as answers you've received, particularly when they are from the principal authors of the package in question :) $\endgroup$ – Aron Ahmadia Mar 10 '13 at 23:06
  • $\begingroup$ The link provided in the above comment by @Aron Ahmadia has moved to a new location, which is ask.cvxr.com/t/trace-of-s-2-in-cvx/132 $\endgroup$ – Mark L. Stone Feb 5 '18 at 14:52
8
$\begingroup$

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))
$\endgroup$
  • $\begingroup$ what do you mean by for a given S? My objective function to minimize is $f(S)=\mathrm{trace}(S)+\mathrm{trace}(S^{-2})$ where S is the unknown that we are looking to find. What you have suggested is equivalent to have $X\geq S^{-1}$ and $X^TZ^{-1}X\geq I$. So this constraints would impose an upper bound for $S^{-1}$ and an another lower bound for the quadratic form of $Z^{-1}$. Do you mean that by adding these two constraints I can avoid writing $Z=S^2$ which cannot be accepted in a convex problem because it's not an affine equality? $\endgroup$ – user2987 Mar 12 '13 at 2:22
  • $\begingroup$ The case with a fixed S was just to create a simple test to check correctness. If it wouldn't work for a fixed S, it would obviously not work when you optimize over S (i.e., when you have S linearly parameterized in some decision variabables). You never add an explicit constraint $Z=S^2$ or something like that. It would kill convexity. The epigraph-based code above is complete, except that you have some more code to define S, what ever it might be. Note that it only works if you minimize the $trace(S^{-2})$ term, as it is based on an upper, in optimality tight, bound. $\endgroup$ – Johan Löfberg Mar 12 '13 at 7:14
  • $\begingroup$ Actually I tried your code using cvx for not fixed S and I found that $S\ne Z^2$. $\endgroup$ – user2987 Mar 12 '13 at 7:25
  • $\begingroup$ The BoundSquare is not equivalent to $X^TZ^{-1}X\geq I$, you have reversed the inequality. $\endgroup$ – Johan Löfberg Mar 12 '13 at 7:35
  • 2
    $\begingroup$ Well, then your problem formulation does not capture what you really want to compute, i.e., it is really not an issue about the implementation of your problem, but the relevance of the model to begin with. $\endgroup$ – Johan Löfberg Mar 13 '13 at 12:10
9
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Welcome to scicomp, Dr. Grant. Feel free to get in touch if you have any questions about the site, and please consider weighing in here about how we can improve the experience of CVX users on the site. $\endgroup$ – Aron Ahmadia Mar 18 '13 at 12:20
  • $\begingroup$ Michael Grant wrote "I am unaware of an AD/SD engine that handles matrix expressions such as f(S)." The MATLAB automatic differentiator ADiMat can compute the gradient (and I suspect Hessian) of $trace(S^{-2})$ $\endgroup$ – Mark L. Stone Sep 13 '15 at 16:30
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Yes you're right and I have already mentioned this in my question. I think that this is a limitation in cvx because it's considering my problem in its own way. But this should work in cvx without any constraint since $f$ is convex and you can check the proof here: math.stackexchange.com/questions/325875/… $\endgroup$ – user2987 Mar 10 '13 at 19:41
  • $\begingroup$ Is there any other solver that I can use in matlab to minimize this function. I tried fmincon but I have a nonlinear constraint which contains a constant matrix from the main program that @mycon function doesn't recognize.. $\endgroup$ – user2987 Mar 10 '13 at 19:56
  • $\begingroup$ You write in your question "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", but my point is that such a constraint is not a convex equality. I have edited my answer above. $\endgroup$ – Geoff Oxberry Mar 10 '13 at 20:16
  • $\begingroup$ which YALMIP version do you think that would solve my problem? Is it the free version SDPT3? Thanks! $\endgroup$ – user2987 Mar 10 '13 at 20:49
  • $\begingroup$ YALMIP is a free interface. Some of the solvers it interfaces with are not free. You can see which solvers are free and which are not here. I believe you must install the solvers separately, and you have many choices. Also, some of the commercial solvers have free academic licenses (I know that Gurobi and CPLEX have free academic licenses, for instance). $\endgroup$ – Geoff Oxberry Mar 10 '13 at 21:06
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ It's good to have you around again Prof. Neumaier :) $\endgroup$ – Aron Ahmadia Mar 11 '13 at 10:18
  • $\begingroup$ Prof. Neumaier many thanks for your answer. Do you mean that the cholesky decomposition could replace the psd constraint? $\endgroup$ – user2987 Mar 11 '13 at 10:57
  • $\begingroup$ @Bel: $C$ psd is equivalent to $C=R^TR$. But the problem in $R$ is more nonlinear than the problem in $C$. But since your problem is nonlinear anyway, it won't hurt much. However, in terms of $R$ the problem is also nonconvex, which means that there might be additional stationary points in the modified problem. So a convex solver able to handle the problem would be likely to be more robust. $\endgroup$ – Arnold Neumaier Mar 11 '13 at 17:22
  • $\begingroup$ @Bel: What are the constraints in your problem? (Without constraints, the solution is not unique.) $\endgroup$ – Arnold Neumaier Mar 11 '13 at 17:24
  • $\begingroup$ @ArnoldNeumaier: My constraint is $\|S-C\|_1\leq \alpha$ where $C\in \mathcal{M}_{m,m}$ and $\alpha$ are given and of course my second constraint is to have $S$ spd otherwise $f$ is not anymore convex. $\endgroup$ – user2987 Mar 11 '13 at 20:43
1
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.