NSolve cannot find solution to system of trig equations

Question

I was making a program to convert a rotation matrix to a set of euler angles. I wanted to do this for any valid rotation order that the euler angles were made with, and I did not want to do it the hard way, so there is a symbolic component. I was able to complete this with sympy, but it was kinda slow.

I thought I would try Mathematica. Mathematica has always been faster than sympy when I use it, and it would be less painful to use now there is the wolfram engine thingy.

So anyway, I multiply three symbolic rotation matrices together in reverse order (fixed frame) and I get this set of trig functions:

lhs = {{Cos[y]*Cos[z], Sin[x]*Sin[y]*Cos[z] - Sin[z]*Cos[x], Sin[x]*Sin[z] + Sin[y]*Cos[x]*Cos[z]}, {Sin[z]*Cos[y], Sin[x]*Sin[y]*Sin[z] + Cos[x]*Cos[z], -Sin[x]*Cos[z] + Sin[y]*Sin[z]*Cos[x]}, {-Sin[y], Sin[x]*Cos[y], Cos[x]*Cos[y]}} 

Then I grab my handy test numeric rotation matrix:

rhs = {{-0.534895228705377, -0.390488683311612, 0.749270099839694}, {-0.833049961066805, 0.391849908220173, -0.390488683311612}, {-0.141120008059867, -0.833049961066805, -0.534895228705377}}

Then, I try to solve:

NSolve[lhs==rhs, {x, y, z}, Reals]

out[]= {}

The output is an empty list! Oh no, what have I done? Have I forgotten how to use Mathematica? Is there a bug with the online version?

I must go back to Sympy, to check a solution does exist!

euler_ang_soln_constraint = [Eq(L, R) for L, R in zip(list(rot_lhs), list(rotation_matrix))]
    euler_angs = list(nsolve(euler_ang_soln_constraint, (r1, r2, r3), (0, 0, 0), prec=precision)[:, 0])

oh well, continuing...

Using my tricksome python codes, I am able to uncover the following solution:

{x->1-Pi, y->0.141592653589793, z->1-Pi}

(yes, I can simplify it, but that is what sympy gave me, and I am lazy :p)

I proceed to check that mathematica is not broken:

lhsNum = lhs /.{x->1-Pi, y->0.141592653589793, z->1-Pi}

out[]= {{-0.534895,-0.390489,0.74927},{-0.83305,0.39185,-0.390489},{-0.14112,-0.83305,-0.534895}}

rhs

out[]= {{-0.534895,-0.390489,0.74927},{-0.83305,0.39185,-0.390489},{-0.14112,-0.83305,-0.534895}}

lhsNum == rhs

out[]= True

I have lost my way...

So I now turn to you, my fellow citizens of these interwebs, for ridicule and guidance. Does anyone see what I am missing or understand how I can entice the magical wolfram gnomes to give me the answers I need?

Honestly, I have no clue what the problem is. I tried changing the precision, normalizing the numeric columns, changing Method -> {"UseSlicingHyperplanes" -> False}, Reals, WorkingPrecision... Plain Solve[] also reports no solutions, so I hope it is some obvious formatting error or matrix solution peculiarity some more experienced users may know about.

Summary:

Why does mathematica not find any numeric solutions to this system of trig equations with NSolve while sympy can find them with nsolve in 100ms? Does anyone see an obvious problem?

lhs = {{Cos[y]*Cos[z], Sin[x]*Sin[y]*Cos[z] - Sin[z]*Cos[x], Sin[x]*Sin[z] + Sin[y]*Cos[x]*Cos[z]}, {Sin[z]*Cos[y], Sin[x]*Sin[y]*Sin[z] + Cos[x]*Cos[z], -Sin[x]*Cos[z] + Sin[y]*Sin[z]*Cos[x]}, {-Sin[y], Sin[x]*Cos[y], Cos[x]*Cos[y]}} 
rhs = {{-0.534895228705377, -0.390488683311612, 0.749270099839694}, {-0.833049961066805, 0.391849908220173, -0.390488683311612}, {-0.141120008059867, -0.833049961066805, -0.534895228705377}}
NSolve[lhs==rhs, {x, y, z}, Reals]

Out[118]= {}

Answer 1

Try

NMinimize[{1, lhs - rhs == 0}, {x, y, z}]
(*{1., {x -> -2.14159, y -> 0.141593, z -> -2.14159}}*)

or

NMinimize[#.# &[Flatten[lhs - rhs]], {x, y, z}]
(*{8.6012*10^-31, {x -> -2.14159, y -> 0.141593, z -> -2.14159}}*)

Don't know why NSolve fails:

NSolve[Join[Map[# == 0 &,Flatten[lhs - rhs]], {-Pi <= x <= Pi, -Pi <= y <= Pi, -Pi <= z <=Pi}], {x, y, z}, Reals]
(* {} *) 

Answer 2

NSolve appears to work if you only select 3 independent equations for your 3 unknowns. For example, take the anti-diagonal elements.

NSolve[Flatten[rhs][[{3, 5, 7}]] == Flatten[lhs][[{3, 5, 7}]], {x, y, 
   z}] // AbsoluteTiming
(* {0.0279919, {{x -> -2.14159, y -> 0.141593, 
   z -> -2.14159}, {x -> -1., y -> 0.141593, z -> -1.}, {x -> 1., 
   y -> 0.141593, z -> 1.}, {x -> 2.14159, y -> 0.141593, 
   z -> 2.14159}}} *)

Edit

I did some additional testing and it turns out that not all the solutions are valid so I created a test to pick the valid solution. This time with the diagonal elements.

{time, sols} = 
 NSolve[Flatten[rhs][[{1, 5, 9}]] == Flatten[lhs][[{1, 5, 9}]], {x, y,
     z}] // AbsoluteTiming
eps = 0.000001;
totals = ((Total@Flatten@Abs[rhs - (lhs /. #)]) & /@ sols);
sol = Extract[sols, FirstPosition[totals, _?(# < eps &)]];
rhs
(lhs /. sol)
(* {0.035351, {{x -> -2.511, y -> -0.846959, 
   z -> -2.511}, {x -> -2.511, y -> 0.846959, 
   z -> 2.511}, {x -> -2.14159, y -> -0.141593, 
   z -> 2.14159}, {x -> -2.14159, y -> 0.141593, 
   z -> -2.14159}, {x -> -1., y -> -3., z -> 1.}, {x -> -1., y -> 3., 
   z -> -1.}, {x -> -0.630596, y -> -2.29463, 
   z -> -0.630596}, {x -> -0.630596, y -> 2.29463, 
   z -> 0.630596}, {x -> 0.630596, y -> -2.29463, 
   z -> 0.630596}, {x -> 0.630596, y -> 2.29463, 
   z -> -0.630596}, {x -> 1., y -> -3., z -> -1.}, {x -> 1., y -> 3., 
   z -> 1.}, {x -> 2.14159, y -> -0.141593, 
   z -> -2.14159}, {x -> 2.14159, y -> 0.141593, 
   z -> 2.14159}, {x -> 2.511, y -> -0.846959, 
   z -> 2.511}, {x -> 2.511, y -> 0.846959, z -> -2.511}}} *)
(* {{-0.534895, -0.390489, 0.74927}, {-0.83305, 
  0.39185, -0.390489}, {-0.14112, -0.83305, -0.534895}} *)
(* {{-0.534895, -0.390489, 0.74927}, {-0.83305, 
  0.39185, -0.390489}, {-0.14112, -0.83305, -0.534895}} *)