Quadratic optimisation with quadratic equality constraints

Question

I would like to solve the following optimisation problem:

$$\text{minimize} \quad x'Ax \qquad \qquad \text{subject to} \quad x'Bx = x'Cx = 1$$

Where $A$ is symmetric and $B$ and $C$ are diagonal.

Does anyone have a suggestion for an efficient way of solving this?

Thank you.

Answer 1

This is a partial answer - you can solve for $$ \min x'Ax \\ s.t. x'Qx = 1 $$ for $Q=B$ or $Q=C$ precisely through strong duality. Similarly, you could also solve for: $$ \min x'Ax \\ s.t. x'x = 2 $$ where I have used the fact that $B$ and $C$ are diagonal to have that $x'Bx + x'Cx = x'(B+C)x$.

I'm not sure though how you could combine these various solutions to ensure that you achieve a minimum for your original problem.

Answer 2

You can have your own opinion as to whether this is efficient in some sense. But if you want to solve the problem, then finding a solution can be considered more efficient than not finding a solution.

Here is the code to solve this using the global optimizer BARON under YALMIP.

x=sdpvar(length(A),1)
optimize([x'*B*x==1,x'*C*x==1],x'*A*x,sdpsettings('solver','baron'))

Depending on the values of B and C, this can solve quite quickly to find the global minimum.

Given that B and C have non-overlapping non-zeros, per the OP's comment, the constraints are the product of two Riemannian manifolds, so manifold optimization could be used, as for instance with MANOPT https://manopt.org/ . I have no idea how well it will do vs. BARON, but it will not be guaranteed to find a global minimum.

Answer 3

Echoing the answers of MotiNK, since there are only two quadratic constraints you can expect to find a global optimum via a semidefinite relaxation (SDP).

In Matlab with the CVX toolbox installed, you could proceed as follows. As long as the optimal matrix $X$ has rank 1 (which is typically the case, at least when $B$ and $C$ are positive semidefinite and such that $B+C$ is positive definite), then the vector $x$ constructed below is a global optimum.

n = 6;
A = randsym(n);
B = diag([rand(n/2, 1) ; zeros(n/2, 1)]);
B = B/trace(B);
C = diag([zeros(n/2, 1) ; rand(n/2, 1)]);
C = C/trace(C);

cvx_begin
    variable X(n, n)
    minimize trace(A*X)
    subject to
        trace(B*X) == 1;
        trace(C*X) == 1;
        X == semidefinite(n);
cvx_end

[V, D] = eig(X);
% Confirm that D(end, end) is the only positive entry.
diag(D)'
% If so, then x is optimal.
x = sqrt(D(end, end))*V(:, end);

Alternatively, as Mark L. Stone pointed out, given the special structure of $B$ and $C$, you could also consider using optimization on two spheres.

Indeed, it seems that $x$ is really a concatenation (perhaps with some re-ordering of entries) of two vectors, as $x = [x_1 ; x_2]$ with $x_1 = \tilde B^{-1/2}y_1$ and $x_2 = \tilde C^{-1/2}y_2$ -- here, $\tilde B$ denotes the top-left diagonal block of $B$, and $\tilde C$ denotes the bottom-right diagonal block of $C$.

Essentially, we get $x^\top B x = x_1^\top \tilde B x_1 = y_1^\top y_1 = 1$ if $y_1$ is unit norm.

So they idea is that we can minimize $x^\top A x$, which is really a function of $y_1, y_2$, both unit-norm vectors. That's a homogeneous quadratic cost function to minimize over a product of two spheres.

The non-convexity is not too bad and reasonably well understood. See for example Corollary 5.8 in this paper for guarantees if we relax the problem just slightly.