Solving a system of linear equations on the form $a^\top X \, b = c$ 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.
 A: $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$.
A: 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()

