Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2004
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

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).

Let's look a bit more closely at your example. Start with

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

In[2]:= InputForm[Timing[es = Eigensystem[mat]]]
Out[2]//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.

In[3]:= InputForm[esrad = ToRadicals[es]]
Out[3]//InputForm=
{{-5/3 - ((7/2)^(2/3)*(1 + I*Sqrt[3]))/(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
    ((1 - I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6,
   -5/3 - ((7/2)^(2/3)*(1 - I*Sqrt[3]))/(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
    ((1 + I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6,
   -5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt[3])/2)^(1/3)) +
    ((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3)/3, 0},
  {{-1 - 3*(-5/3 - ((7/2)^(2/3)*(1 + I*Sqrt[3]))/
        (3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
       ((1 - I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6) -
     (-5/3 - ((7/2)^(2/3)*(1 + I*Sqrt[3]))/(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
       ((1 - I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6)^2,
    2*(-5/3 - ((7/2)^(2/3)*(1 + I*Sqrt[3]))/(3*(-1 + 
(3*I)*Sqrt[3])^(1/3)) -
       ((1 - I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6) +
     (-5/3 - ((7/2)^(2/3)*(1 + I*Sqrt[3]))/(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
       ((1 - I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6)^2,
    -5/3 - ((7/2)^(2/3)*(1 + I*Sqrt[3]))/(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
     ((1 - I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6, 1},
   {-1 - 3*(-5/3 - ((7/2)^(2/3)*(1 - I*Sqrt[3]))/
        (3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
       ((1 + I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6) -
     (-5/3 - ((7/2)^(2/3)*(1 - I*Sqrt[3]))/(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
       ((1 + I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6)^2,
    2*(-5/3 - ((7/2)^(2/3)*(1 - I*Sqrt[3]))/(3*(-1 + 
(3*I)*Sqrt[3])^(1/3)) -
       ((1 + I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6) +
     (-5/3 - ((7/2)^(2/3)*(1 - I*Sqrt[3]))/(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
       ((1 + I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6)^2,
    -5/3 - ((7/2)^(2/3)*(1 - I*Sqrt[3]))/(3*(-1 + (3*I)*Sqrt[3])^(1/3)) -
     ((1 + I*Sqrt[3])*((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3))/6, 1},
   {-1 - 3*(-5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt[3])/2)^(1/3)) +
       ((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3)/3) -
     (-5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt[3])/2)^(1/3)) +
       ((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3)/3)^2,
    2*(-5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt[3])/2)^(1/3)) +
       ((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3)/3) +
     (-5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt[3])/2)^(1/3)) +
       ((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3)/3)^2,
    -5/3 + 7^(2/3)/(3*((-1 + (3*I)*Sqrt[3])/2)^(1/3)) +
     ((7*(-1 + (3*I)*Sqrt[3]))/2)^(1/3)/3, 1}, {0, 0, 0, 1}}}

In[4]:= LeafCount[es]
Out[4]= 414

In[5]:= LeafCount[esrad]
Out[5]= 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[6]:= InputForm[N[es]]
Out[6]//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.}}}

In[7]:= InputForm[N[esrad]]
Out[7]//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:

In[8]:= Timing[N[esrad,100];]
Out[8]= {0.02 Second, Null}

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

To answer your question, one can avoid the Root form by using ToRadicals 
but for virtually all purposes I can think of you'd be best off leaving 
it as is.

Daniel Lichtblau
Wolfram Research


  • Prev by Date: Re: Re: Linear Programming
  • Next by Date: suggestion for frontend
  • Previous by thread: Root function
  • Next by thread: Re: Root function