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

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 2

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 3

Here is an alternative way of solving - let $K$ be sufficiently large, and consider the problem

$$ \min x'(A + K(B-C))x \\ s.t. x'(B-C)x \geq 0 \\ x'x \geq 2 $$ where $K$ is set as a penalty to ensure that the first constraint is tight. The second constraint will necessarily be tight at the minimum.

This problem has two quadratic constraints and may be solvable through strong duality - in particular see Theorem 3.5. Note that you need to solve the dual and then you can check whether there is no duality gap by that Theorem.

Answer 4

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.