# Optimizing a quadratic restricted to the sphere

Let $A$ be $p\times p$ symmetric positive definite with distinct eigenvalues and $x_p\in \mathbb{R}^p$ and consider the problem

Minimize $x'Ax + b'x$

Subject to $x'x=1$

Most of the information I've found is is either very general/theoretical or specific to linear constraints, although I'm largely flitting around optimization texts and crossing my fingers since I don't know exactly what I'm looking for.

Anyway, my first pass was to use a Lagrange multiplier;

$f(x, \lambda) = x'Ax + b'x + \lambda(x'x-1)$

Taking derivatives and setting to zero gives

$x = -\frac{1}{2} (A+\lambda I)^{-1}b$

$\frac{1}{4} b'(A+\lambda I)^{-2}b = 1$

I've got my system in $p+1$ equations and I can go about solving them. Analytically I haven't gotten anywhere, except simplifying things a little with the eigendecomposition of $A$. When $b=0$ the solution is trivially $x=e_1$, the first eigenvector of $A$, so let's ignore that case. So my first question: is there an analytical solution that I'm too mathematically challenged to see? If not, what is the best way to solve this problem?

(To help quantify "best", I have potentially many such problems to solve for smallish $p$, say 5-10, and $A$ is the same but $b$ changes. An approximate solution is OK, in fact an approximate solution near the correct global solution is better than an exact local one.)

-
Not sure this keyword is of any help to you, this seems to be a quite particular QCQP (quad. const. quad. prog.)en.wikipedia.org/wiki/… Unfortunately, this is more or less were my knowledge on this stops. Minor remark: it seems to me there is typo in your def of f(x,lambda), a superflous lambda. – user9072 Apr 10 '11 at 23:53
On second thought, I am not so sure anymore the QCQP suggestion is so good, equality vs inequality. Sorry, for the noise. – user9072 Apr 11 '11 at 0:14
Fixed the error, thanks. QCQP was part of my offhand remark on very general/theoretical info :) Although the remark at the bottom that interior point methods could be applied might be a lead. – JMS Apr 11 '11 at 0:14
Here is the ref. to the classic paper on solving such trust-region subproblems: Computing a Trust Region Step. Jorge J. Moré and D. C. Sorensen. SIAM J. Sci. and Stat. Comput. 4, pp. 553-572 (20 pages) – Suvrit Apr 11 '11 at 7:20

Your problem has been studied extensively in the context of trust region methods for optimization, and there are a number of algorithms that have been developed.

See for example:

W. W. Hager, Minimizing a quadratic over a sphere. SIAM Journal on Optimization, 12:188-208, 2001.

Hager's paper gives a lemma that characterizes the solutions to your problem and a solution in terms of the eigenvalues and eigenvectors of A as well as algorithms for solving the problem. Since then there have been several other papers written on this topic, with particular attention to algorithms for solving instances where A is large and sparse- this isn't a particular issue for you.

-
The problem is actually much more classic---and I guess the most famous algorithm for solving it is the Moré-Sorensen Newton algorithm. For the large-sparse case, one can download the LSTRS software from: ta.twi.tudelft.nl/wagm/users/rojas/lstrs.html – Suvrit Apr 11 '11 at 7:25
@Suvrit- I agree that the TRS problem has been studied for a long time. I didn't mean to imply that this was the first paper on the topic, but rather cited this paper because I think it gives the simplest answer to this question. The original poster has many instances of a very scale version of the problem with a fixed $A$ matrix and many different $b$ vectors, so simply using the eigenvalue-eigenvector decomposition of $A$ is the way to go here rather than using something like LSTRS. – Brian Borchers Apr 11 '11 at 13:09
A copy of the Hager paper can be found at: math.ufl.edu/~hager/papers/Regular/sphere.pdf – Brian Borchers Apr 11 '11 at 13:11
Exactly what I was looking for, thanks! Kind of embarassing I couldn't turn it up in an afternoon of googling... – JMS Apr 11 '11 at 16:31

This is not an answer, but it was getting too long to be a comment. Also, I'm not familiar with such problems, so hopefully someone will be able to give you a complete and more satisfying answer and some references. Here are some general remarks. It's a bit messy and mostly guesswork :S

If $|A^{-1}b|<1$ there might be several solutions, i.e. several $x$ minimizing $\varphi(x):=\frac{1}{2}\langle x|Ax \rangle+\langle b|x \rangle$ (I added $\frac{1}{2}$ for convenience only), precisely, if you call $\lambda_1>\dots>\lambda_p>0$ the distinct eigenvalues of $A$ and you have $e_1,\dots, e_p$ an associated orthonormal frame and $b$ such that $$|A^{-1}b|<1~\mathrm{and}~b\in\mathrm{Span}(e_1,\dots,e_r)\setminus\mathrm{Span}(e_1,\dots,e_{r-1})$$

I would expect the set of minimizing $x$ to be a sphere of dimension $p-r-1$ and a point when $r=p$

If $|A^{-1}b|>1$, a picture convinces me there is only one solution to your problem (and the "projection onto a convex closed set"-theorem in inner product spaces proves it). If $x_0$ is a solution, then we have to have that the differential of $\varphi$ vanishes at $x_0$ which translates into $\langle Ax_0|- \rangle+\langle b|- \rangle$ vanishing on the orthogonal of $x_0$ i.e. there exists some $\mu\in\mathbb{R}$ such that $Ax_0+b=\mu x_0$, and again, a picture suggests there are preciseyl two solutions and that we have to look for the solutions with $\mu < 0$. Then there is indeed only one solution, because if we look for $x_0\in\mathbb{S}^{n-1}$ satisying above equation, we can write down $x_0=-(A-\mu)^{-1}b$ because the matrix is invertible. The norm of above vector is stricly increasing with $\mu$ and equals $|A^{-1}b|>1$ at $\mu=0$ and tends to $0$ as $\mu$ tends to $-\infty$. So there'll be only one $\mu<0$ such that $-(A-\mu)^{-1}b$ has unit length. But finding that given $\mu$ requires solving the equation $\frac{b_1^2}{(\lambda_1-\mu)^2}+\dots+\frac{b_p^2}{(\lambda_p-\mu)^2}=1$

So if $|A^{-1}b|>1$, you can find the solution as long as you can find the negative $\mu$ that satisfies the above equation...

-