# Least-squares solution to an over-defined probem.

This problem arose in my stereo vision project.

I have two matrix equations:

$$\left( \begin{array}{ccc} x_1.w_1 \\ y_1.w_1\\ w_1 \end{array} \right) = \left( \begin{array}{ccc} A_{11} & A_{12} & A_{13} & A_{14} \\ A_{21} & A_{22} & A_{23} & A_{24} \\ A_{31} & A_{32} & A_{33} & A_{34} \end{array} \right) * \left( \begin{array}{ccc} X \\ Y \\ Z \\ 1 \end{array} \right)$$

$$\left( \begin{array}{ccc} x_2.w_2 \\ y_2.w_2\\ w_2 \end{array} \right) = \left( \begin{array}{ccc} B_{11} & B_{12} & B_{13} & B_{14} \\ B_{21} & B_{22} & B_{23} & B_{24} \\ B_{31} & B_{32} & B_{33} & B_{34} \end{array} \right) * \left( \begin{array}{ccc} X \\ Y \\ Z \\ 1 \end{array} \right)$$

$x_{1}$,$y_{1}$,$x_{2}$,$y_{2}$ and all $A_{ij}$ and $B_{ij}$ are know scalars. Find the best-fit to $X$,$Y$,$Z$ (in a least squares sense).

I should point out that I don't know the values of $w_1$ and $w_2$, otherwise this would be the same question as in my other question (I accidentally poorly-defined my problem in that one).

The answer to the above problem is in here (page 6), but in the way described in the link it would be necessary to compute the matrix inversions every time, when $X$,$Y$ and $Z$ changes. I would like to know how to write the solution to the above problem in the following way:

$$\left( \begin{array}{ccc} X \\ Y \\ Z \\ 1 \end{array} \right) = U*\left( \begin{array}{ccc} x_{1} \\ y{1} \\ 1 \\ x_{2}\\ y_{2} \\ 1 \end{array} \right)$$ (or in some similar form)

Where the matrix U does not depend on the values of $x_{1}$,$y_{1}$,$x_{2}$,$y_{2}$.

-

$$\left( \begin{array}{ccc} x_1.w_1 \\ y_1.w_1\\ w_1 \end{array} \right) = \left( \begin{array}{ccc} A_{11} & A_{12} & A_{13} & A_{14} \\ A_{21} & A_{22} & A_{23} & A_{24} \\ A_{31} & A_{32} & A_{33} & A_{34} \end{array} \right) * \left( \begin{array}{ccc} X \\ Y \\ Z \\ 1 \end{array} \right)$$ This is the same as $$\left( \begin{array}{ccc} w_1 \\ w_1\\ w_1 \end{array} \right) = \left( \begin{array}{ccc} A_{11}/x_1 & A_{12}/x_1 & A_{13}/x_1 & A_{14}/x_1 \\ A_{21}/y_1 & A_{22}/y_1 & A_{23}/y_1 & A_{24}/y_1 \\ A_{31} & A_{32} & A_{33} & A_{34} \end{array} \right) * \left( \begin{array}{ccc} X \\ Y \\ Z \\ 1 \end{array} \right)$$ I am supposing that this scaling is what leads you (and maybe the authors of the paper) to believe that an inverse need be calculated every time. This is because I am supposing that your $A$ and $B$ are constant but that $x$ and $y$ are not. I believe you stated about as much.

So instead of $X=MP$ (to use terms from your previous question), you now have $$X=DMP$$ where $D$ is a diagonal matrix. So instead of $[M^\top M]^{-1}M^\top$ we now have (using $D^\top = D$) $$[M^\top D^2M]^{-1}DM^\top$$ And indeed the scaling alters things. Note that the $D$ matrix (the $x$ and $y$ values) can not be commuted with $M$, so your desire to write it as a matrix left multiply will not be met.

If you are willing to "go under the hood" with the least squares, once you have the in-span left inverse of $M$, it is possible to update things to account for the matrix $D$ in a more efficient way. But it is a bit involved, and will only save computation if you need to do this often. Also, it would save more computation as dimension is increased, so may not be worth while here. To describe the way to do it would require much more than I am willing to write right now (sorry), and I am not certain it is something that would easily be found in literature.

(The short story is this: scale your least squares inverse appropriately to keep it a left inverse, then re-orthogonalize it to have it as the in-span left inverse. I have not done this specifically so there may be some "gotcha" details, but I believe this would be more efficient than re-calculating each time. But as I said, involved and fairly complicated.)

-
I didn't fully understood your answer. Is it possible to write: $\left( \begin{array}{ccc} X \\ Y \\ Z \\ 1\end{array} \right)=U*\left( \begin{array}{ccc} x_1 \\ y_1 \\ 1 \\ x_2 \\ y_2 \\ 1\end{array} \right) ?$ –  JLagana Feb 4 '13 at 12:03
No, the main reason is that $D=\pmatrix{x_1 & & \\ & y_1 & \\ & & 1}$ does not commute to have $x_1$ and $y_1$ (and $x_2$, $y_2$) as a right multiplication with $U$. I would be willing to write up how to use least squares for $P$ in $X=MP$ to efficiently get $P$ in $X=DMP$ with $D$ diagonal. If that is indeed what you want. Maybe pose it as another question "How to efficiently go from $X=MP$ least squares to $X=DMP$ least squares solution of $P$, with $D$ diagonal and $M$ a $6\times 4$ matrix?" I could pose the question myselr, but then would take me a day to be able to answer it... –  adam W Feb 4 '13 at 16:19
That would be really helpful, if you have the time I would be grateful. –  JLagana Feb 6 '13 at 1:05
It is easier done than explained, though I think from reading from your pdf link that it is not $DMX = P$ with psuedo-inverse of $DM$ desired, but something different. I am not certain there is any possibility of an efficient update in the case of your project. –  adam W Feb 6 '13 at 18:51
If there is such a possibility, I would look at the inverse of $\left[\matrix{M & Q}\right]$, where $Q$ is the normalized null column space of $M$, ie $Q^\top M =0$ and $Q^\top Q = I$. Q may be found in the QR factorization of $M$, it would be the null portion of the Q factor there. The inverse of $\left[\matrix{M & Q}\right]$ contains the least squares pseudo-inverse, and changes to $M$ could efficiently be applied to this inverse to keep track of the changes to the pseudo-inverse. As I stated, easier done than explained. If you formulate your needs exactly in the next question... –  adam W Feb 6 '13 at 18:54