I am trying to solve the following optimization problem: \begin{align} &\min\limits_{s} \rm{tr}\left(S^T S\right) + \mathrm{tr}\left(\left(S^T S\right)^{-2}\right)\\ &\text{subject to }\rm\left\|S^TS-Z^TZ\right\|_{norm}\leq \epsilon\end{align} Where $\rm \|\cdot\|_{norm}$ could be either the Frobenius norm or the $\ell_1$-norm.

That is, $\rm S^TS$ belongs to the ball with radius $\epsilon$ and center $\rm Z^TZ$. Actually, I want to solve my problem with $\rm Q = S^TS$ which is a symmetric positive semi-definite matrix, thus, a convex problem to solve but I couldn't find another way to impose that property on my solution.

So I have tried first with minFunc toolbox by using the lagrange expression, and the Frobenius norm for the constraint $\rm \|\cdot\|_{norm}$ i.e.,

$$\min\limits_{s} \rm \mathrm{tr}\left(S^T S\right) + \mathrm{tr}\left(\left(S^T S\right)^{-2}\right) + \lambda\cdot\mathrm{tr}\left(\left(S^T S - Z^TZ\right)\left(S^T S - Z^TZ\right)^T\right)$$

But it stops after few iterations with a very high error about $10^{2}$ for a small $S\in\mathrm{R}^{10\times10}$.

Here is my code:

options.Method = 'cg';%'lbfgs'; 
options.maxIter = 100000000;
options.MaxFunEvals = 20000000000000;
lambda = 0.001;
x0 = rand(100,1);
Z = rand(15,10);
funObj = @(x)myfunc(x, lambda, Z);
[sol, f_val] = minFunc(funObj,x0,options);

function [f,g] = myfunc(x, lambda, Z)
s = reshape(x, [sqrt(numel(x)),sqrt(numel(x))]);
f = trace(s'*s) + trace(inv(s'*s)^2) + lambda* trace((s'*s - Z'*Z)*(s'*s - Z'*Z)');
g = 2*s -4*s*(inv(s'*s))^3 + lambda*(2*s*s'*s + 2*s'*s*s -2*s*Z'*Z -2*Z'*Z*s);
g = g(:);

I have tried also with minConf_PQN using the $\ell_1$-norm constraint but returns this error:

close to singular or badly scaled. Results may be inaccurate. RCOND = NaN

And here is the code that I have used:

funObj = @(x)myfunc(x);

Z = rand(15,10);
S0 = load('S_true.mat')
tau = sum(abs(L(:))); %true tau
r = reshape(Z'*Z, [size(Z,2)^2,1]);
funProj = @(w)sign(w).*projectRandom2C(abs(w'w-r),tau);
sol = minConf_PQN(funObj,x0,funProj,options);

function [f,g] = myfunc(x)

s = reshape(x, [sqrt(numel(x)),sqrt(numel(x))]);

f = trace(s'*s) + trace(inv(s'*s)^2);

g = 2*s -4*s*(inv(s'*s))^3;
g = g(:);


Any help solving my problem would be really appreciated.

  • 2
    $\begingroup$ ..."which is a convex problem". Well, the problem you have stated is not convex. Just because the norm is convex doesn't mean that the composition of a quadratic and a norm is convex. It is a convex set in $Q=S^TS$, but not in $S$ $\endgroup$ – Johan Löfberg Feb 27 '16 at 18:52
  • $\begingroup$ $S^TS$ itself is a symmetric psd matrix and the objective function has been proven to be convex in this post plus the constraint is a convex space. But I am not sure if it's convex in terms of only $S$ and not $S^TS$. $\endgroup$ – user2987 Feb 27 '16 at 18:58
  • $\begingroup$ I would like to solve the problem for $Q$ symmetric using minConf or minFunc so that is the only way that I found to impose that constraint. $\endgroup$ – user2987 Feb 27 '16 at 19:02
  • $\begingroup$ I just edited the question and removed the convexity claim. $\endgroup$ – user2987 Feb 27 '16 at 19:03

Not a direct answer to your title question, but I think you are better off attacking this problem from the semidefinite domain instead. Trivial approach is to linearize the objective at some initial guess, solve the linearized problem, perform a line-search along the computed direction, and repeat until the objective doesn't improve. The code below does this in an implementation using YALMIP (disclaimer: developed by me). It terminates in a second or so, using an SDP solver such as Mosek or SDPT3.

Z = rand(15,10);

% Take a step D from current solution Qi, Q = Qi+D
D = sdpvar(10);
% Initial guess Q = Qi = Z'*Z
Qi = Z'*Z;
% Values to test in brute-force line-search
alpha = logspace(-4,3,100);
% Run some iterations
for i = 1:10
    % Linearized objective
    Objective = trace(Qi + D) + trace(Qi^-2 - Qi^-2*D*Qi^-1 - Qi^-1*D*Qi^-2)
    % Solve linear SDP
    optimize([Qi+D>=0, norm(Qi + D - Z'*Z,'fro') <= .01],Objective)
    % Perform naive line-search
    for j = 1:length(alpha)
        Qtest = Qi + alpha(j)*value(D);
        if min(eig(Qtest))>0  && norm(Qtest - Z'*Z,'fro')<=0.01
            Objtest(j) = trace(Qtest) + trace(Qtest^-2);
            Objtest(j) = inf;
    if all(Objtest >= trace(Qi) + trace(Qi^-2))
    % Pick best step
    [~,best] = min(Objtest);
    % Update solution
    Qi = Qi + alpha(best)*value(D);

Adding another answer, as I just realized that the problem is easily solved as a linear SDP. Let $Q=S^TS$ and you have the objective $\mathrm{trace}~Q + \mathrm{trace}~Q^{-2}$. Introduce an upper bound $X\succeq Q^{-1}$ and minimize $\mathrm{trace}~Q + \mathrm{trace}~X^{2}$. At optimality you will have $X= Q^{-1}$. The constraints $X\succeq Q^{-1}$ and $Q\succeq 0$ is converted to an LMI using a Schur complement, and the objective is convex quadratic.

Also here implemented in YALMIP

Q = sdpvar(n);
X = sdpvar(n);
Model = [[X eye(n);eye(n) Q]>=0, norm(Q-Z'*Z)<= .01]
Objective = trace(Q) + trace(X*X);
  • $\begingroup$ Thanks a lot for your answer. But I am trying to solve with minFunc because my matrix $Q$ is about $1000\times1000$ or even larger whereas YALMIP can't solve it for a larger scale than 100 I guess. $\endgroup$ – user2987 Mar 1 '16 at 16:46
  • 2
    $\begingroup$ For 1000x1000 case, you will have to develop your own solver from scratch, seriously exploiting structure in the problem. That is roughly 500000 variables, and most solvers will choke on that. $\endgroup$ – Johan Löfberg Mar 1 '16 at 17:48

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.