Solving a system of linear equations on the form $a^\top X \, b = c$

Question

I'm implementing a computer vision algorithm for 3D point reconstruction from photos (Tomasi-Kanade). For the last day I've been stuck on a linear algebra problem. The system that I want to solve has $3m$ equations, with $m$ as the amount of photos, for which I want to solve for $L$:

$$a_{i1}^T L a_{i1} = 1$$

$$a_{i2}^T L a_{i2} = 1$$

$$a_{i1}^T L a_{i2} = 0$$

where $a_{ij}$ is a $3$-vector and $L = C C^T$ where $C$ is a $3 \times 3$ matrix. To solve it, it should be in the form of $Ax = b$. Had $a_{ij}$ been an invertible matrix, the equations could be rewritten as follows

$$a_{ij}^T L a_{ij} a_{ij}^{-1} = c a_{ij}^{-1} \rightarrow a_{ij}^T L = c a_{ij}^{-1}$$

but this can't be done since they are vectors.

I'm sure there is a simple step to solve this, but my linear algebra is a lot worse than it used to be and I just can't wrap my head around it.

Answer 1

$L$ can be written in detail as $$ L = \begin{bmatrix} L_{11} & L_{12} & L_{13} \\ L_{12} & L_{22} & L_{23} \\ L_{13} & L_{23} & L_{33} \end{bmatrix}. $$ If $u = \begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix}$ then \begin{align} u^T L u &= u_1^2 L_{11} + u_2^2 L_{22} + u_3^2 L_{33} + 2 u_1 u_2 L_{12} + 2 u_1 u_3 L_{13} + 2 u_2 u_3 L_{23} \\ &= \begin{bmatrix} u_1^2 & u_2^2 & u_3^2 & 2 u_1 u_2 & 2 u_1 u_3 & 2 u_2 u_3 \end{bmatrix} \begin{bmatrix} L_{11} \\ L_{22} \\ L_{33} \\ L_{12} \\ L_{13} \\ L_{23} \end{bmatrix}. \end{align} This shows how to convert equations 1-3 to matrix notation.

Once you have computed $L$ (by finding a least squares solution to your system of equations), you can then get $C$ by computing the Cholesky factorization of $L$.

Answer 2

Given vectors $\mathrm a_1, \mathrm a_2, \dots, \mathrm a_n \in \mathbb R^3$ and scalars $b_1, b_2, \dots, b_n \in \{0,1\}$, we would like to solve the following system of $n$ linear equations in symmetric and positive semidefinite matrix $\mathrm X \in \mathbb R^{3 \times 3}$

$$\begin{aligned} \mathrm a_1^\top \mathrm X \,\mathrm a_1 &= b_1\\ \mathrm a_2^\top \mathrm X \,\mathrm a_2 &= b_2\\ &\vdots\\ \mathrm a_n^\top \mathrm X \,\mathrm a_n &= b_n\end{aligned}$$

One option would be to solve the following least-squares problem

$$\begin{array}{ll} \text{minimize} & \displaystyle\sum_{k=1}^n \left( \mathrm a_k^\top \mathrm X \,\mathrm a_k - b_k \right)^2\\ \text{subject to} & \mathrm X \succeq \mathrm O_3\end{array}$$

Let us define

$$\mathrm A := \begin{bmatrix} — \left( \mathrm a_1 \otimes \mathrm a_1 \right)^\top —\\ — \left( \mathrm a_2 \otimes \mathrm a_2 \right)^\top —\\ \vdots\\ — \left( \mathrm a_n \otimes \mathrm a_n \right)^\top —\end{bmatrix} \qquad \qquad \mathrm b := \begin{bmatrix} b_1\\ b_2\\ \vdots\\ b_n\end{bmatrix}$$

where $\otimes$ is the Kronecker product (in NumPy: kron). Hence, the optimization problem above can be rewritten in the more succinct form

$$\begin{array}{ll} \text{minimize} & \| \mathrm A \,\mbox{vec} (\mathrm X) - \mathrm b \|_2^2\\ \text{subject to} & \mathrm X \succeq \mathrm O_3\end{array}$$

where $\mbox{vec}$ is the vectorization operator. In CVXPY, given NumPy arrays A and b:

X = Semidef(3)

# create optimization problem
objective = Minimize( norm(A * vec(X) - b, 2) )
prob = Problem(objective, [])

# solve optimization problem
prob.solve()