Solving a system of polynomial equations with multiple variables

Question

I have a system of equations of the form: $$ l_i^T l_j \cdot m_i^T m_j - m_i^T R l_j \cdot l_i R^T m_j = 0$$ where $R \in \mathbb{R}^{3\times3}$ is an unknown rotation matrix. $l_i, l_j, m_i, m_j \in \mathbb{R}^{3}$ are known vectors.

The equation can also be rewritten: $ r^T Q r = 0 $ where $ r \in \mathbb{R}^9 $ is the rotation matrix written in vector form and $Q \in \mathbb{R}^{9\times9}$ is a known symmetric matrix

My goal is to find $R$ if possible in general form, else numerically.

Since $R$ is a rotation, only 3 variables need to be found and thus 3 equations of the previous form should be enough. Do you think that any software (Mathematica, Maple, Sage, Sympy...) could solve this system in general form?

I have already tried with Sympy and Sage that gave me memory errors. So, in case a general solution could not be found, what do you think are the best optimization techniques to find a numerical solution of this system?

Answer 1

Any rotation matrix can be decomposed into the product of a rotation around the $x$-axis by some angle $\alpha$, a rotation around the $y$-axis by some angle $\beta$, and a rotation around the $z$-axis by some angle $\gamma$. These principal rotation matrices are trigonometric functions of their angles of rotation, so I'm not sure how you would obtain these angles symbolically as a function of the input vectors.

In theory, though, you could use that decomposition in a nonlinear equation solver to find $\alpha, \beta, \gamma \in [0, 2\pi)$; Newton's method (or variants thereof) would probably work fine, assuming your system of equations is well-posed.

You could also opt to use an axis-angle representation of your rotation matrix instead. Then, instead of solving for the principal rotation angles $\alpha, \beta, \gamma$, you have to solve for an angle $\theta \in [0, 2\pi)$ and an axis $u$ such that $u_{x}^{2} + u_{y}^{2} + u_{z}^{2} = 1$. Add the constraint to your system of equations and solve for $u$ and $\theta$; even though you'll be solving for four variables, the presence of the unit norm constraint means you effectively have three degrees of freedom. This representation also involves trigonometric functions, so again, I'm not sure how you would obtain these angles symbolically, but numerically, it would probably work just fine, again, assuming your system of equations is well-posed.

Answer 2

Consider the equation in the form $x^tAx=0$, with $A$ symmetric. Since $A$ is symmetric we can write $A=VDV^t$, where $V$ is orthogonal, and $D$ is the real diagonal matrix of its eigenvalues $\lambda_i$. With the change of variables $y=V^tx$, we have the equation $$ \sum_i \lambda_i y_i^2 = 0. $$

This gives a complete characterization of solutions. In terms of $z_i=y_i^2$, (so that $y_i=\pm\sqrt{z_i}$ and $x=Vy$) all the solutions of the original problem are in one-to-one correspondence with the solutions of $$ \sum_i \lambda_i z_i = 0, \qquad \text{subject to}\ z_i\geq0. $$ The set of all solutions $z$ is convex; it's non-empty whenever at least one $\lambda_i$ is zero or whenever at least two nonzero $\lambda_i$ have different signs. You can also parametrize the set of solutions by separating the indices into those with positive, negative, and zero $\lambda_i$, which makes it easy to deal with the constraint.

In general, you won't be able to find a closed form expression for the eigenvalues $\lambda_i$ in terms of the original matrix $A$, but this can easily be done numerically for any given matrix using standard software.