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?
 A: 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.
A: 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
A: This problem can be solved using singular value decomposition. The original problem can be stated as
$$min\|Ax-b\|_2, \|x\|_2=1$$
By substituting $$A'=[A;-b], x'=w x; w$$ and choosing $w=\frac{\|x\|_2}{\|x;1\|_2}$ we get $\|x'\|_2=w\|x;1\|_2=\|x\|_2$, so the problem becomes
$$min\|A'x'\|_2, \|x'\|_2=1$$
which is same as in this question: Minimizing $\|Ax\|_2$ subject to $\|x\|_2 = 1$
It can be solved by decomposing $A'$ into $U \Sigma V^T$ and then taking the last column of $V$, which corresponds to the smallest singular value of $A'$.
Finally $x$ is calculated by dividing $x'$ by the last element, taking elements of $x'$ except the last one and then normalizing it.
