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.
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.
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))
}
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:
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.
