       Re: Root function

• To: mathgroup at smc.vnet.net
• Subject: [mg51033] Re: [mg51009] Root function
• From: Daniel Lichtblau <danl at wolfram.com>
• Date: Sat, 2 Oct 2004 03:18:07 -0400 (EDT)
• References: <200410010849.EAA11569@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

```Ben wrote:
> When I try to diagonalize this matrix in mathematica 5:
> Eigensystem[{{-1, 1, 0, 0}, {1, -2, 1, 0}, {0, 1, -2, 0}, {0, 0, 1, 0}}]
> I get this ugly mess
>
> \!\({{Root[1 + 6\ #1 + 5\ #1\^2 + #1\^3 &, 1], Root[1 + 6\ #1 +
>              5\ #1\^2 + #1\^3 &, 2], Root[
>        1 + 6\ #1 + 5\ #1\^2 + #1\^3 &, 3], 0}, {{\(-1\) - 3\
>              Root[1 + 6\ #1 + 5\ #1\^2 + #1\^3 &, 1] - Root[1 + 6\ #1 +
>   etc etc etc
>
> Mathematica 4 gives me the right answer.
> How do I get v5 to not spit this out and what does it mean?

The Root form is a concise way of expressing algebraic numbers via the
minimal polynomial they satisfy, along with a canonical ordering in the
complex plane (hence those numbers 1-3 in the second field).

mat = {{-1, 1, 0, 0}, {1, -2, 1, 0}, {0, 1, -2, 0}, {0, 0, 1, 0}};

In:= InputForm[Timing[es = Eigensystem[mat]]]
Out//InputForm=
{0.05*Second, {{Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 1, 0],
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 2, 0], Root[1 + 6*#1 + 5*#1^2 +
#1^3 & ,
3, 0], 0}, {{-1 - 3*Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 1, 0] -
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 1, 0]^2,
2*Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 1, 0] +
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 1, 0]^2,
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 1, 0], 1},
{-1 - 3*Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 2, 0] -
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 2, 0]^2,
2*Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 2, 0] +
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 2, 0]^2,
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 2, 0], 1},
{-1 - 3*Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 3, 0] -
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 3, 0]^2,
2*Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 3, 0] +
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 3, 0]^2,
Root[1 + 6*#1 + 5*#1^2 + #1^3 & , 3, 0], 1}, {0, 0, 0, 1}}}}

Fast, and concise. To get a radical form, you can use ToRadicals.

Out//InputForm=
{{-5/3 - ((7/2)^(2/3)*(1 + I*Sqrt))/(3*(-1 + (3*I)*Sqrt)^(1/3)) -
((1 - I*Sqrt)*((7*(-1 + (3*I)*Sqrt))/2)^(1/3))/6,
-5/3 - ((7/2)^(2/3)*(1 - I*Sqrt))/(3*(-1 + (3*I)*Sqrt)^(1/3)) -
((1 + I*Sqrt)*((7*(-1 + (3*I)*Sqrt))/2)^(1/3))/6,
-5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt)/2)^(1/3)) +
((7*(-1 + (3*I)*Sqrt))/2)^(1/3)/3, 0},
{{-1 - 3*(-5/3 - ((7/2)^(2/3)*(1 + I*Sqrt))/
(3*(-1 + (3*I)*Sqrt)^(1/3)) -
((1 - I*Sqrt)*((7*(-1 + (3*I)*Sqrt))/2)^(1/3))/6) -
(-5/3 - ((7/2)^(2/3)*(1 + I*Sqrt))/(3*(-1 + (3*I)*Sqrt)^(1/3)) -
((1 - I*Sqrt)*((7*(-1 + (3*I)*Sqrt))/2)^(1/3))/6)^2,
2*(-5/3 - ((7/2)^(2/3)*(1 + I*Sqrt))/(3*(-1 +
(3*I)*Sqrt)^(1/3)) -
((1 - I*Sqrt)*((7*(-1 + (3*I)*Sqrt))/2)^(1/3))/6) +
(-5/3 - ((7/2)^(2/3)*(1 + I*Sqrt))/(3*(-1 + (3*I)*Sqrt)^(1/3)) -
((1 - I*Sqrt)*((7*(-1 + (3*I)*Sqrt))/2)^(1/3))/6)^2,
-5/3 - ((7/2)^(2/3)*(1 + I*Sqrt))/(3*(-1 + (3*I)*Sqrt)^(1/3)) -
((1 - I*Sqrt)*((7*(-1 + (3*I)*Sqrt))/2)^(1/3))/6, 1},
{-1 - 3*(-5/3 - ((7/2)^(2/3)*(1 - I*Sqrt))/
(3*(-1 + (3*I)*Sqrt)^(1/3)) -
((1 + I*Sqrt)*((7*(-1 + (3*I)*Sqrt))/2)^(1/3))/6) -
(-5/3 - ((7/2)^(2/3)*(1 - I*Sqrt))/(3*(-1 + (3*I)*Sqrt)^(1/3)) -
((1 + I*Sqrt)*((7*(-1 + (3*I)*Sqrt))/2)^(1/3))/6)^2,
2*(-5/3 - ((7/2)^(2/3)*(1 - I*Sqrt))/(3*(-1 +
(3*I)*Sqrt)^(1/3)) -
((1 + I*Sqrt)*((7*(-1 + (3*I)*Sqrt))/2)^(1/3))/6) +
(-5/3 - ((7/2)^(2/3)*(1 - I*Sqrt))/(3*(-1 + (3*I)*Sqrt)^(1/3)) -
((1 + I*Sqrt)*((7*(-1 + (3*I)*Sqrt))/2)^(1/3))/6)^2,
-5/3 - ((7/2)^(2/3)*(1 - I*Sqrt))/(3*(-1 + (3*I)*Sqrt)^(1/3)) -
((1 + I*Sqrt)*((7*(-1 + (3*I)*Sqrt))/2)^(1/3))/6, 1},
{-1 - 3*(-5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt)/2)^(1/3)) +
((7*(-1 + (3*I)*Sqrt))/2)^(1/3)/3) -
(-5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt)/2)^(1/3)) +
((7*(-1 + (3*I)*Sqrt))/2)^(1/3)/3)^2,
2*(-5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt)/2)^(1/3)) +
((7*(-1 + (3*I)*Sqrt))/2)^(1/3)/3) +
(-5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt)/2)^(1/3)) +
((7*(-1 + (3*I)*Sqrt))/2)^(1/3)/3)^2,
-5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt)/2)^(1/3)) +
((7*(-1 + (3*I)*Sqrt))/2)^(1/3)/3, 1}, {0, 0, 0, 1}}}

In:= LeafCount[es]
Out= 414

Out= 1284

This, by the way, is actually slightly smaller than the radical version
4 result.

In addition to size, there are other reasons to prefer the Root[] form.

(i) It is typically faster to obtain. In this example, which is fairly
easy, I see a factor of three or so gain in version 5 vs 4.2.

(ii) It is numerically more stable to evaluate. In general radical
formulations are prone to numeric problems. Root[] objects do not have
this liability.

(iii) When the roots of an irreducible cubic are all real, but not
rational, the so-called "casus irreducibilus" shows that they still must
be expressed in terms of Sqrt[-1]. This means that numeric evaluation
will give small imaginary parts unless, by happenstance, they exactly
cancel. Small numeric error from round off makes this unlikely. For example:

In:= InputForm[N[es]]
Out//InputForm=
{{-3.246979603717467, -1.554958132087371, -0.19806226419516174, 0.},
{{-1.801937735804838, 4.048917339522306, -3.246979603717467, 1.},
{1.2469796037174667, -0.6920214716300959, -1.554958132087371, 1.},
{-0.4450418679126288, -0.3568958678922094, -0.19806226419516174, 1.},
{0., 0., 0., 1.}}}

Out//InputForm=
{{-3.246979603717467 + 0.*I, -1.5549581320873713 -
2.220446049250313*^-16*I,
-0.19806226419516182 + 1.1102230246251565*^-16*I, 0.},
{{-1.801937735804838 + 0.*I, 4.048917339522306 + 0.*I,
-3.246979603717467 + 0.*I, 1.},
{1.2469796037174672 - 2.4406313453516104*^-17*I,
-0.6920214716300959 + 2.464509183785474*^-16*I,
-1.5549581320873713 - 2.220446049250313*^-16*I, 1.},
{-0.4450418679126286 - 2.8908825018377506*^-16*I,
-0.35689586789220956 + 1.7806594772125943*^-16*I,
-0.19806226419516182 + 1.1102230246251565*^-16*I, 1.}, {0., 0., 0.,
1.}}}

(iv) For sufficiently complicated algabraics, it is often faster to
evaluate the Root form numerically, at least at high precision. For example:

Out= {0.02 Second, Null}

In:= Timing[N[es,100];]
Out= {0. Second, Null}