2

I would like to use linear least squares to solve for $x \in \mathbb{R}^5$ where

$$ Ax = b \rightarrow x = (A^TA)^{-1}A^Tb $$

but would like to include the constraint

$$ x_1^2 + x_2^2 = 1 $$

I looked around and found a few examples of how to solve using linear constraints of the form

$$ Cx = d$$

but unsure how to continue with quadratic constraints. Thanks!

Rodrigo de Azevedo
  • 23,223
  • 2
    So effectively you are saying you want to minimize $(Ax-b)^T(Ax-b)$ under the constraint $x^Tx=1$. – Johan Löfberg Nov 06 '14 at 19:13
  • @Johan Yes, but I should have mentioned $ x \in \mathbb{R}^5 $. Updated question. –  Nov 06 '14 at 19:29
  • $x^\top \begin{pmatrix}1&0&0&0&0\0&1&0&0&0\0&0&0&0&0\0&0&0&0&0\0&0&0&0&0\end{pmatrix}x=1$ – davcha Nov 06 '14 at 19:50
  • Using Royi's answer https://math.stackexchange.com/a/2402360/125757 should yield a good approach: by relaxing $x^Tx = 1 \to x^Tx \le 1$ and then finding an optimal multiplier for the constraint such that it is tight. – Guillermo Angeris Aug 22 '17 at 15:47

1 Answers1

1

Quadratic optimization over a single quadratic constraint is a classical problem. One way to address it in a global fashion is to employ a semidefinite relaxation, which can be shown to be tight in this special case (also called a lossless application of the S-procedure)

One reference out of many http://arxiv.org/abs/1312.1400

Here is some MATLAB reference code to setup and solve the semidefinite relaxation (called moment relaxation sometimes). It assumues you have the modelling toolbox YALMIP installed (disclaimer, developed by me), and a semidefinite solver, such as sedumi, sdpt3, mosek.

A = randn(10,5);
b = rand(10,1);
x = sdpvar(5,1);

p = (A*x-b)'*(A*x-b);
h = 1-x(1)^2-x(2)^2;

[diagnostics,xhat] = solvemoment(h == 0,p)
xhat{1}

The semidefinite relaxation boils down to a generalized eigenvalue problem in this case, so a general semidefinite programming approach might be a bit of an overkill, but it is convenient and fast enough (~0.1s) unless you have to solve many of these.

Having said that, this problem is also small enough to be easily solved using a standard global nonlinear solver, such as baron, scip or bmibnb in YALMIP.

 optimize(h == 0,p,sdpsettings('solver','bmibnb'))
 value(x)
Johan Löfberg
  • 9,737
  • 1
  • 16
  • 15