4

This problem arose in my stereo vision project. $$ P_{1c} = A*P $$

$$ P_{2c} = B*P $$ where: $P_{1c}$ and $P_{2c}$ are $3\times1$ vectors, $A$ and $B$ are $3 \times 4$ matrices and $P$ is a $4\times1$ vector.

Or, more verbosely:

$$ \left( \begin{array}{ccc} a_1 \\ b_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} a_2 \\ b_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) $$ $P_{1c}, P_{2c}, A$ and $B$ are known, find $P$.

An explicit solution is a nice start, and in fact I already have it:

(Explicit Solution - page 6)

but from the nature of the problem, I can tell that it will not always be possible to find $P$ using the explicit solution (since I only have estimates for $P_{1c},P_{2c},A$ and $B$), so I was trying to tackle the problem like this:

$f(P) = (A*P-P_{1c})^2+(B*P-P_{2c})^2$

Minimize $f(P)$ with respect to $P$.

I believe this notation is not as rigorous as one would expect, since $P$ is a vector. So, to put it in another way:

Knowing $P_{1c},P_{2c},A$ and $B$, find the vector $P$ wich gives the "best-fit" to $$ P_{1c} = A*P $$

$$ P_{2c} = B*P $$ (in a least squares sense).

Any help would be great, and please feel free to approach this problem with a different strategy than the one I'm trying, but I would appreciate if someone could answer this without "opening up" the matrices (instead of using the $A_{ij}$'s and working out with the explicit equations, only use the matrices, and matrix operations like $A^T, A^{-1}$), if that's even possible.

Hope I made myself clear, thank's in advance!

Sigur
  • 6,416
  • 3
  • 25
  • 45
JLagana
  • 462

1 Answers1

2

You have one matrix equation:

$$\pmatrix{P_{1c} \\ P_{2c}} = \pmatrix{A \\ B}\cdot P \tag{1}$$

It is an overdetermined system (unless there are too many dependencies in the data).

Here is what you want for the least squares solution to it:

$$P=\left[\pmatrix{A^\top & B^\top}\pmatrix{A \\ B}\right]^{-1}\pmatrix{A^\top & B^\top}\pmatrix{P_{1c} \\ P_{2c}}$$

There are methods that are more efficient than this explicit formula here, that avoids the extra calculation of squaring then inverting. You may notice that I got this using what is called a left inverse of $\pmatrix{A \\ B}$ on the equation (1). There are possibly more than one left inverse in general, but this one is the specific left inverse that gives the least squares solution.

Let $M = \pmatrix{A \\ B}$. Using user7530's helpful comment, it may be factored as $M=QR$, which is a standard form available in any computer package, simply called the QR factorization. I will show how to use that form later, but first the derivation of the formula I gave you.

To solve for $P$ in the equation $$X= MP \tag{2}$$ where $M$ has more rows than columns, consider a left inverse $\hat{M}$ such that $\hat{M}M=I$. At first glance, assuming the existence of $\hat{M}$, this seems to be the solution since it gives $$\hat{M}X = \hat{M}MP = IP = P$$ and $$P=\hat{M}X$$ gives the solution. But the left inverse may not be unique. So the question is which left inverse do we want. Consider equation (2). We are looking for some $P$ to solve it, but it may not be possible. So instead we look at the closest possible solution. Any $P$ when right multiplied with $M$ will necessarily give a result that is in the column span of $M$, since a right multiply can only mix the columns of $M$. This along with one-sided inverses turns out to be all the considerations necessary to find the unique solution.

Since the result of $\hat{M}X$ should have the same column-span as $M$, this gives a unique $\hat{M}$ and it is

$$\hat{M} = \left[M^\top M\right]^{-1}M^\top$$

To use the QR factorization, we have $M=QR$. Then \begin{align} P&=\hat{M} X \\ &= \left[M^\top M\right]^{-1}M^\top X\\ & = \left[R^\top Q^\top Q R\right]^{-1} M^\top X\\ & = \left[R^\top R\right]^{-1} M^\top X\\ \end{align}

Where the last line there is easier to compute since $R$ is triangular with zero rows.

adam W
  • 5,565
  • +1 One good "more efficient" approach is to compute the QR decomposition of $\left(\begin{array}{c}A\B\end{array}\right)$. – user7530 Jan 23 '13 at 17:03
  • Could you post a derivation for your answer? Or maybe cite a reference that would guide me to a better understanding of how you got there? Thanks! – JLagana Jan 23 '13 at 17:04
  • I am more of a self study type of person and am not too good at references, and I personally find the usual Wikipedia links to be difficult; the articles never explain terribly well in my opinion. But I would be happy to expand on my answer, and will add to it now. – adam W Jan 23 '13 at 17:16
  • Just to be clear, I think I understood how you got to P=[(A⊤B⊤)(AB)]−1(A⊤B⊤)(P1cP2c) (multiply both sides by the transpose and then by the inverse of the product between the transpose and AB). But why is this the least squares answer? – JLagana Jan 23 '13 at 17:22
  • It has to do with the row-span: the least squares solution lies in the row-span of $\pmatrix{A \ B}$ – adam W Jan 23 '13 at 17:38
  • I edited the answer. And I think that yet still there is a more efficient way of doing least squares. I programmed a "homemade" version that costs the same number of calculations as a regular inverse, ie that of multiplication. As I stated, I am not too good on references, and explaining how I made mine is a chapter in a book itself probably. – adam W Jan 23 '13 at 18:00
  • I googled for "least squares inverse" and found something to get me started! Thank you for the consideration =] – JLagana Jan 23 '13 at 18:12
  • Sorry about the earlier comment about row-span, I meant column-span. See the new edit. – adam W Jan 24 '13 at 18:27
  • @user7530 Thanks for that. You might be interested to know I am writing up my algorithm for computing (in exact arithmetic) the psuedo-inverse. I had thought the process to be known, but it seems apparently it is not (in publication I can find anyway). It requires twice the operations of the normal inverse, but benefits from the non-square nature, and may also be easily made parallel. – adam W Feb 26 '13 at 21:45
  • @user7530 (and h3now) I have published my procedure – adam W Feb 28 '13 at 18:22