Take the 2-minute tour ×
MathOverflow is a question and answer site for professional mathematicians. It's 100% free, no registration required.

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.

share|improve this question
1  
$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
add comment

4 Answers

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.

share|improve this answer
add comment

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.

share|improve this answer
    
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
add comment

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))
}
share|improve this answer
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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.