# SchurDecomposition doesn't work (in a weird way)

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 ?
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?
• If the matrix is symmetric (after all, you are dealing with quadratic forms), then you shouldn't need SchurDecomposition[]. Apr 21, 2020 at 13:16
• @J.M. that was another ironic issue. Sometimes it seems neither SVD nor eigenvectors give the desired result for a symmetric matrix. Apr 21, 2020 at 13:18
• I just "plugged" the polynomial 1 - 3 x + 3 x^2 - 9 y + 9 x y + 27 y^2 - 6 z + 12 x z + 18 y z + 12 z^2 into canonicalForm[], and lo and behold!, I got 1. What does this mean? (I'm somewhat of a novice in this area.) See my preceding questions math.stackexchange.com/questions/3660652/… and mathoverflow.net/questions/359459/…. I used Q_1=x, Q_2= y, Q_3 =z. I'll now try the code of J. M. in his answer. May 6, 2020 at 0:00

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.

• Your method is indeed simpler and more efficient than mine. But still, the issue with SchurDecomposition is unaddressed, and furthermore, it still doesn't work on the third equation. Apr 21, 2020 at 14:33
• I just tried it on the third equation; you need to do MapAt[FullSimplify @* Orthogonalize, Eigensystem[coefs[[3]]], {2}] instead, and quad /. Thread[{x, y, z} -> Transpose[vecs].{xt, yt, zt}] ought to give a result that has Root[] objects in the coefficients, but is otherwise free of cross-terms. My point anyway was that you don't need to use SchurDecomposition[] here, unless your coefficients are already inexact numbers. Apr 21, 2020 at 14:39
• Per my just posted comment to the question itself, I took quad =1 - 3 x + 3 x^2 - 9 y + 9 x y + 27 y^2 - 6 z + 12 x z + 18 y z + 12 z^2, and fed it through the three indicated steps and got 1 - ((61 + 13 Sqrt[61]) Sqrt[3538 + 272 Sqrt[61]] xt)/1098 + 3/2 (14 + Sqrt[61]) xt^2 + ( Sqrt[3538 - 272 Sqrt[61]] (61 - 13 Sqrt[61]) yt)/1098 - 3/2 (-14 + Sqrt[61]) yt^2 and Sqrt[2/3] xt + 6 xt^2 + (2 yt)/Sqrt[3] + 3 yt^2 + zt (Sqrt[2] + 2 zt) Clear? Signifigance? May 6, 2020 at 0:08
• This is why I said writing a robust routine is hard work, @Paul. A good routine ought to be able to detect degenerate cases like yours (cf. the criteria listed here and many other references). May 6, 2020 at 0:41
• Thanks--and thanks for the quadric surfaces link also, J. M. But maybe you could further "spoon-feed" me as to the nature of the degeneracy? What conclusions can one make as to the nature of quad =1 - 3 x + 3 x^2 - 9 y + 9 x y + 27 y^2 - 6 z + 12 x z + 18 y z + 12 z^2? What transformation can be employed to reduce it to a presumably simpler form? Actually, as the previous links of mine indicate, I'm primarily interested in the inequality (!) constraint $x^2+3 x y +\left(3 y+z\right){}^2<3 y+2 x z$. . So, what can I conclude as to the volume that the quadric surface encloses. May 6, 2020 at 1:44