# Solving a system of polynomial equations with multiple variables

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?

• I don't see how you compute the value of $Q^{-1}$ but I think that you use the fact that Q is positive semidefinite which is false, in general (at least for the sqrt of diag($Q$)) – Yohann Salaun Sep 11 '14 at 17:05
• I dont see the purpose in inverting it and I dont see why you can affirm that $Q^{-1} = rr^T$ – Yohann Salaun Sep 11 '14 at 17:39
• Woah, sorry about that, I read your original statement incorrectly. – Aurelius Sep 11 '14 at 17:54
• Is the equation supposed to read $l_i^T l_j \cdot m_i^{T} R R^T m_j - m_i^{T} R l_j \cdot l_i^{T} R^T m_j = 0$? – Geoff Oxberry Sep 11 '14 at 22:06
• It may help to note that $R R^T = I$ whenever $R$ is an orthogonal matrix. – tmyklebu Sep 12 '14 at 16:22

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.

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.

• Suppose you had to solve $x_{1}^{2} - x_{2}^{2} = 0$. Your method transforms the equation to $z_{1} - z_{2} = 0$, $z_{1}, z_{2} \geq 0$, and argues that the solution then would be $z_{1} = z_{2}, z_{1} \geq 0$. In that case, $\sqrt{y_{1}} = \sqrt{y_{2}}$, so $y_{1} = y_{2}, y_{1} \geq 0$. However, the set of solutions is the union of the lines $y_{1} - y_{2} = 0$ and $y_{1} + y_{2} = 0$ for all $y_{1} \in \mathbb{R}$. The problem is that your mapping from $y_{i}$ to $z_{i}$ is not one-to-one. – Geoff Oxberry Sep 14 '14 at 4:16
• @GeoffOxberry You're right! I inserted a $\pm$ sign into the answer, which should be all right now. My point in transforming to $z_i$ was make it clear that the equations are linear, I forgot about signs. – Kirill Sep 14 '14 at 4:21