# Least square given constraint on subcomponents

Hi all,

I have to find a set of parameters that fit to a set of data with constraint to a subset of the parameters. In summary, I want to solve $\min ||A[x_1 \, x_2]^T -b ||$ given $||x_2|| = g$.

I thought the problem is trivial but it turns out that it is not trivial at all. Lagrange multiplier approach gives a very complicated matrix equation (still not solvable for me). Any idea, or numerical/analytical solution, is greatly appreciated.

Thanks.

-
$x_1$ and $x_2$ are vectors? –  J. M. Nov 23 '10 at 1:16
Hi all, I haven't understand all your answers. But I will look into them. Thank you so much for your replies. And just to clarify, $x_1$ and $x_2$ are vectors. –  user11175 Nov 29 '10 at 19:50

Writing $A = [A_1 A_2]$ and taking without loss of generality $g=1$ (scale $A_2$ appropriately), your problem is equivalent to $$\min_{x_1,x_2}\ ||A_1 x_1 + A_2 x_2 - b||^2 \text{ s.t. } ||x_2||=1.$$ Notice that the solution of $\min_{x_1} ||A_1 x_1 - c||^2$ can be obtained in closed form (assuming column independence in $A_1$): it is equal to $c^T P_1 c$ with $P_1 = I - A_1 (A_1^T A_1)^{-1}A_1^T$ (note $P_1$is positive semidefinite). Now you simply have to solve $$\min_{x_2}\ (b-A_2 x_2)^T P_1 (b-A_2 x_2) \text{ s.t. } ||x_2||=1$$ which can be done with standard techniques. I am not sure a closed-form solution can be obtained, but you can for example obtain a scalar equation in the Lagrange multiplier, which you solve numerically, and then obtain $x_2$. See also this link, where the problem is reduced to a quadratic eigenvalue problem.

-

You can express this problem as a quadratically constrained quadratic program (QCQP). Unfortunately, because of the equality constraint, the QCQP will be nonconvex. However, there is some discussion on the Wikipedia page for handling nonconvex QCQP's, and a Google search should turn up more. This paper "Relaxations and Randomized Methods for Nonconvex QCQP's" might help, too; Example 1.2.1 in the paper is very similar to your problem.

-
Why the reason for the downvote? If someone will explain what's wrong with my answer I'll gladly delete it. –  Mike Spivey Nov 24 '10 at 0:09

Have you tryed a simple newton algorithm (with the constraint added to the algo)?

Let $(\alpha_{k})$ be defined as: $\alpha_k=1/k^2$

Initialisation :

$x^0=[0,...,0]$

Compute $H=(A^* A)^{-1}$

Loop for $k$ in $1:m$

$x^k=x^{k-1}-\alpha_k H A^*(Ax^{k-1}-b)$

$x^{k}=g*x^{k}/\|x_2^k\|$

end for loop

Obviously, there are more adaptive way of choosing $\alpha_k$... but maybe you don't need such sofistication to solve a norm minimization problem. If $A^* A$ has very small eigen values you can use $H_k=(A^* A+\epsilon_k)^{-1}$ instead of $H$ ($\epsilon_k$ decreasing to zero)...

Note that this type of code is relatively general when you want to find the saddle point of a lagrangian and you know how to find the maxima with respect to Lagrange multipliers (in the dual space) (second step of the loop) but you need a gradient descent (or here Newton algo) for finding the minima in the principal space.

Here is the corresponding R code:

A=t(array(1:1000,c(10,100)))
m=100; b=1:10; g=3; l=5; p=10;
alpha=1:m
alpha=1/alpha^2
x=array(0,c(m,p))
H=t(A)%*%A
svdH=svd(H)
H=svdH$v%*%diag(1/svdH$d)%*%t(svdH$u) for (k in 2:m) { x[k,]=x[k-1,]-alpha[k]*H%*%(t(A)%*%(A%*%x[k-1,]-b)) x[k,]=g*x[k,]/sqrt(sum(x[k,(l+1):p]^2)) print(sum((A%*%x[k,]-b)^2)) }  - add comment First eliminate$x_1$by solving an ordinary least squares, and then you need to essentially solve a problem of the form:$\min x_2^TMx_2$s.t.$\|x_2\|=g$, for appropriate$M$. This problem is the famous trust-region subproblem, aka, TRS. Please have a look at the following references (and references therein), which provide algorithms and discussion on how to solve such problems; perhaps you can simplify or adapt one of their methods: 1. LSTRS: http://ta.twi.tudelft.nl/wagm/users/rojas/lstrs-paper.pdf 2. Moré-Sorensen TRS algorithm (in the book on Trust-region subproblems) 3. http://www.optimization-online.org/DB_HTML/2002/09/530.html Depending on how large$A\$ is, or what kind of structure it has, different TRS methods may be preferred. Example, for small matrices, where you can afford to do Cholesky, the More-Sorensen method is usually very hard to beat. If your matrix is however large and sparse, then you might prefer the LSTRS method instead.

-