SchurDecomposition doesn't work (in a weird way)

Question

I have written a function to get the canonical form of a 3D quadric surface. This involves finding a unitary transformation of the coordinates which eliminates the cross-terms such as $xy$ or $xz$. For example, let $q(x,y,z)=2xy-x+y+z-1$. Then by a rotation of axes, this equation is converted to $y^2- x^2+ \sqrt{2}x+z-1$.

The code is as follows

canonicalForm[polynomial_] := Module[{a, b, c, svd, p, v, transform},
  v = Normal@ CoefficientArrays[polynomial];
  a = Last[v];
  b = v[[2]];
  c = First[v]; (** form of x\[Transpose].a.x + b.x + c = 0 **)
  a = (a + Transpose[a])/2;           (** make it symmetric **)
  svd = SingularValueDecomposition[a];
  If[Equal @@ Drop[svd, {2}],
    p = First[svd],
    p = Normalize/@ Last@ Eigensystem[a]; (* THIS IS LINE 10 *)
    If[Expand[p.Transpose@ p] != IdentityMatrix[3],
      p = Chop@ First@ SchurDecomposition@ N@ a]
  ];
  transform = p.{{x}, {y}, {z}};
  a = Simplify@Flatten[Transpose[transform].a.transform];
  b = Together[b.transform];
  Chop@ First[a + b + c]]

Three equations are used to test the code: $$\begin{align} q_1&=3 x^2-2 x y+2 x z+5 y^2-2 y z+3 z^2+2 z\\ q_2&=2 x y-\sqrt{2} x-\sqrt{2} y+z-1\\ q_3&=y z-x y+z^2+z-2 \end{align}$$ The code works perfectly for $q_1$ and $q_2$, but for $q_3$ it takes a long time and then gives a very lengthy output in terms of the roots of some polynomials.

canonicalForm[3x^2 -2x y +2x z +5y^2 -2y z +3z^2 +2z]

$6 x^2+3 y^2+2 z^2+\frac{1}{3} \left(\sqrt{6} x+2 \sqrt{3} y+3 \sqrt{2} z\right)$

canonicalForm[2x y -Sqrt[2]x -Sqrt[2]y +z -1]

$-x^2+y^2-2 y+z-1$

canonicalForm[y z -x y +z^2 +z -2]

a vary lengthy output which is also mathematically wrong

For $q_1$ singular value decomposition results in a unitary transform. For $q_2$ the normalized eigenvectors give the desired unitary transform. But for $q_3$ none of these works and the code switches to SchurDecomposition, which only accepts approximate numerical values.

Now, if I skip calculating the normalized eigenvectors and replace the 10th, 11th and 12th lines with this:

p = Chop@ First@ SchurDecomposition@ N@ a

The code will produce the desired correct results, but in approximate numerical form.

canonicalForm[y z -x y +z^2 +z -2]

$-0.585043 x^2+0.233192 x+0.344446 y^2+0.397113 y+1.2406 z^2+0.88765 z-2$

So it seems calculating the eigenvectors before using SchurDecomposition somehow affects its performance. I wasn't able to suppress this error by defining temporary variables and any other method that I could think of. The only way is deleting those three lines and directly using SchurDecomposition.

Questions

1. Is there any way other than circumventing eigenvectors, that I could make this work? Is this a bug?
2. Is there any other method that results in a unitary transform and gives the exact values, rather than approximate numerical ones given by SchurDecomposition?

Answer 1

Writing a robust routine for dealing with quadric surfaces (and conic sections for that matter) is a lot of work, so I'll give a skeleton from which I hope the OP (or other more motivated users) can build on further.

I'll take the first quadric as an example:

quad = 3 x^2 - 2 x y + 2 x z + 5 y^2 - 2 y z + 3 z^2 + 2 z;

The first important step is to ensure that CoefficientArrays[] returns a symmetric matrix:

coefs = Normal[CoefficientArrays[quad, Variables[quad], "Symmetric" -> True]]
   {0, {0, 0, 2}, {{3, -1, 1}, {-1, 5, -1}, {1, -1, 3}}}

Then, one should also recall that Eigensystem does not yield orthonormal eigenvectors by default, so one has to do some extra work:

{vals, vecs} = MapAt[Orthogonalize, Eigensystem[coefs[[3]]], {2}]
   {{6, 3, 2}, {{1/Sqrt[6], -Sqrt[(2/3)], 1/Sqrt[6]},
                {1/Sqrt[3], 1/Sqrt[3], 1/Sqrt[3]},
                {-(1/Sqrt[2]), 0, 1/Sqrt[2]}}}

Now that you have an orthogonal matrix, you can do this:

quad /. Thread[{x, y, z} -> Transpose[vecs].{xt, yt, zt}] // Simplify
   Sqrt[2/3] xt + 6 xt^2 + (2 yt)/Sqrt[3] + 3 yt^2 + zt (Sqrt[2] + 2 zt)

which should now yield a result that is amenable to completing the square.