# Solve a linear equation system, while enforcing a unit vector solution

I have an equation system of the form Aix + Biy + Ciz = Di, where (x,y,z) is a unit vector, and (Ai, Bi, Ci, Di) are sets of measurements from a noisy system (with typically 3-5 independant readings).

My first intuition to solve this problem was to pose this as an overdetermined linear equation system AX = B where X = (x,y,z), and to solve for X. However, with that approach, I have no way to enforce that the solution for vector X is a unit vector.

Is there an elegant (or standard) solution to that problem, or should I simply dive into non-linear equation solving solutions?

-
Suggestion: If you have $n$ measurement sets, then you have $\binom{n}{3}=% \frac{n(n-1)(n-2)}{3!}$ ways of choosing $3$ of them and for each one solve the corresponding linear equation system $AX=B$. You could determine $% X^{\ast }=(x^{\ast },y^{\ast },z^{\ast })$ so that $\left\vert 1-\left( x^{\ast }\right) ^{2}-\left( y^{\ast }\right) ^{2}-\left( z^{\ast }\right) ^{2}\right\vert$ is minimum. Finally you would compute $X^{\ast \ast }=(x^{\ast \ast },y^{\ast \ast },z^{\ast \ast })=\frac{X^{\ast }}{\left\vert X^{\ast }\right\vert }$. –  Américo Tavares Nov 5 '10 at 22:31

Unfortunately, the problem is almost certainly innately non-linear, since your solution space is inherently so. Pragmatically, the best approach is probably to treat your overdetermined linear system as a minimization-problem, least-squares style; this gives you a (positive definite) function to be minimized. Find the (exact) minimum value in $\mathbb{R}^3$, project that down to $S^2$ (i.e., normalize your vector X) to find the closest point on the sphere, and use that as a starting point for some descent-based method (i.e., find the component of the gradient of your goodness-of-fit function tangent to the sphere, add that to your current point on $S^2$ and renormalize down to the sphere again, and keep doing that until you meet whatever convergence criterion you want to use. There's no guarantee that this will work if your data is particularly ill-conditioned, but in practice it's likely to be the most straightforward method.

-

Your objective is to solve $Ax=b$ subject to the constraint $||x||=1$. One way to do this is to recognize that the solution to $Ax=b$ is an affine space (translated subspace). This means that out of all the solutions to $Ax=b$, there is a unique "smallest" one, i.e. where $||x||$ is as small as possible. To find this solution, you can solve the minimium-norm optimization problem, whose solution is: $$x_0 = \lim_{\lambda \to 0^+} (A^TA+\lambda I)^{-1}A^Tb$$ In the case where $A$ is skinny and full-rank, this simplifies to $(A^TA)^{-1}A^Tb$. In general, the solution can be found by taking the SVD of $A$.

Once you have verified that $x_0$ is indeed a solution, the next thing to do is check that $||x_0|| \leq 1$, for if that were not the case, then you could confidently assert that there is no solution to $Ax=b$ with norm 1.

The next step is to find a basis for the nullspace of $A$, i.e. the set of vectors $x$ such that $Ax=0$. This can be done using the SVD again. If $N$ is a matrix whose columns form a basis for the nullspace of $A$, then the set of all solutions to $Ax=b$ is: $$x = x_0 + Ny$$ where $y$ is an arbitrary vector of appropriate dimension. The key point here is that $x_0$ is orthogonal to all vectors in the nullspace of $A$ (because $x_0$ is the solution point closest to the origin). Thus, you may compute the norm of $x$ using the Pythagorean theorem: $$||x||^2 = ||x_0||^2 + ||Ny||^2$$ Now we see that the if we want $||x||=1$, we simply need to choose $y$ so that $$||Ny|| = \sqrt{1 - ||x_0||^2}$$ We may always choose the columns of $N$ to be orthonormal, so let's assume we have done that. Then the set of desired $y$ is precisely the sphere with radius $\sqrt{1 - ||x_0||^2}$. By choosing any $y$ on this sphere, we parametrize all solutions to your problem. If you only care about a particular solution, then you can choose $y$ of the form $[\alpha,0,...,0]$, where $\alpha = \sqrt{1 - ||x_0||^2}$.f

-
The 2-norm may or may not be appropriate here, I think; the OP never said anything about how noisy his data is. –  Guess who it is. Nov 5 '10 at 22:33
The 2-norm is the appropriate norm if the noise is IID Gaussian, which is as good an assumption as any. "how noisy" it is makes no difference -- the only thing that matters is the pdf. –  Laurent Lessard Nov 6 '10 at 0:21
Yes, I'm just giving a note that one should be careful, since not all noisy data have Gaussian noise. –  Guess who it is. Nov 6 '10 at 0:26
The solution to $Ax = b$ would be a subspace if the system were underdetermined, wouldn't it? The OP stated that their system was overdetermined. –  Rahul Nov 6 '10 at 1:04
indeed. I assumed in my write-up that the problem actually had a solution -- my method is then exact. In the case where it does not ($Ax_0 \neq b$), then $x_0$ is the least-squares estimate, which minimizes the mean squared error across all measurements. Similarly, if $A$ has no nullspace, we can let $N$ be the right-singular-vector corresponding to the smallest singular value of $A$. –  Laurent Lessard Nov 8 '10 at 3:51