How to solve a constrained optimization problem using minFunc or minConf

Question

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(:);
end

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')
L=S0'*S0-Z'*Z; 
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(:);

end

Any help solving my problem would be really appreciated.

Answer 1

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);
        else
            Objtest(j) = inf;
        end
    end
    if all(Objtest >= trace(Qi) + trace(Qi^-2))
        break
    end
    plot(alpha,Objtest)
    % Pick best step
    [~,best] = min(Objtest);
    % Update solution
    Qi = Qi + alpha(best)*value(D);
end

Answer 2

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);
optimize(Model,Objective)