Message Boards

WOLFRAM COMMUNITY

5276 Views

7 Replies

6 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Mathematics Algebra Wolfram Language

Calculate the Rational Canonical Form of a square matrix

Hongyi Zhao

Posted 2 years ago

Based on the code snippet here, I try to calculate the Rational Canonical Form of a square matrix as follows: In[84]:= CompanionMatrix[p_, x_] := Module[ {n, w = CoefficientList[p, x]}, w = -w/Last[w]; n = Length[w] - 1; SparseArray[{{i_, n} :> w[[i]], {i_, j_} /; i == j + 1 -> 1}, {n, n}]] MatrixMinimalPolynomial[a_List?MatrixQ, x_] := Module[{i, n = 1, qu = {}, mnm = {Flatten[IdentityMatrix[Length[a]]]}}, While[Length[qu] == 0, AppendTo[mnm, Flatten[MatrixPower[a, n]]]; qu = NullSpace[Transpose[mnm]]; n++]; First[qu] . Table[x^i, {i, 0, n - 1}]] rationalMatrixForm[a_?(MatrixQ[#, NumericQ] &)] := Module[(version 8/24/13 2PM) {p, q, min, c = {}, moreFactors = True, z, x}, p = CharacteristicPolynomial[a, x]; min = MatrixMinimalPolynomial[a, x]; While[moreFactors, q = PolynomialQuotient[p, min, x]; If[q === 0, moreFactors = False; If[Not[FreeQ[p, x]], z = CompanionMatrix[p, x]; AppendTo[c, z] ], z = CompanionMatrix[min, x]; AppendTo[c, z]; p = q ] (* if ) ]; (end WHILE more factorization needed) SparseArray[Band[{1, 1}] -> c] ] In[87]:= x={{0,1,1,1},{1,0,1,1},{1,1,0,1},{1,1,1,0}}; A=x//rationalMatrixForm//Normal Inverse[A] . X . A Out[88]= {{0, 3, 0, 0}, {1, 2, 0, 0}, {0, 0, 0, 3}, {0, 0, 1, 2}} Out[89]= {{-(2/3), 5/3, 1/3, 5/3}, {1/3, 2/3, 1/3, 5/3}, {1/3, 5/ 3, -(2/3), 5/3}, {1/3, 5/3, 1/3, 2/3}} But the above result is inconsistent with the one given by GAP, as shown below: gap> x:=[[0, 1, 1, 1], > [1, 0, 1, 1], > [1, 1, 0, 1], > [1, 1, 1, 0]]; [ [ 0, 1, 1, 1 ], [ 1, 0, 1, 1 ], [ 1, 1, 0, 1 ], [ 1, 1, 1, 0 ] ] gap> A:=RationalCanonicalFormTransform(x); [ [ 0, 0, 0, 1 ], [ 1, 0, 0, 1 ], [ -1, 1, 0, 1 ], [ 0, -1, 1, 0 ] ] gap> A^-1 x * A; [ [ -1, 0, 0, 0 ], [ 0, -1, 0, 0 ], [ 0, 0, 0, 3 ], [ 0, 0, 1, 2 ] ] Any tips for fixing this problem will be appreciated. Regards, Zhao

Based on the code snippet here, I try to calculate the Rational Canonical Form of a square matrix as follows:

In[84]:= CompanionMatrix[p_, x_] := Module[
  {n, w = CoefficientList[p, x]},
  w = -w/Last[w];
  n = Length[w] - 1;
  SparseArray[{{i_, n} :> w[[i]], {i_, j_} /; i == j + 1 -> 1}, {n, 
    n}]]

MatrixMinimalPolynomial[a_List?MatrixQ, x_] := 
 Module[{i, n = 1, qu = {}, 
   mnm = {Flatten[IdentityMatrix[Length[a]]]}}, 
  While[Length[qu] == 0, AppendTo[mnm, Flatten[MatrixPower[a, n]]];
   qu = NullSpace[Transpose[mnm]];
   n++];
  First[qu] . Table[x^i, {i, 0, n - 1}]]

rationalMatrixForm[a_?(MatrixQ[#, NumericQ] &)] := Module[(*version 8/24/13 2PM*)
  {p, q, min, c = {}, moreFactors = True, z, x},

  p = CharacteristicPolynomial[a, x];
  min = MatrixMinimalPolynomial[a, x];

  While[moreFactors,
   q = PolynomialQuotient[p, min, x];

   If[q === 0,

    moreFactors = False;
    If[Not[FreeQ[p, x]],
     z = CompanionMatrix[p, x];
     AppendTo[c, z]
     ],

    z = CompanionMatrix[min, x];
    AppendTo[c, z];
    p = q

    ] (* if *)

   ]; (*end WHILE more factorization needed*)

  SparseArray[Band[{1, 1}] -> c]
  ]

In[87]:= x={{0,1,1,1},{1,0,1,1},{1,1,0,1},{1,1,1,0}}; 

A=x//rationalMatrixForm//Normal
Inverse[A] . X . A

Out[88]= {{0, 3, 0, 0}, {1, 2, 0, 0}, {0, 0, 0, 3}, {0, 0, 1, 2}}

Out[89]= {{-(2/3), 5/3, 1/3, 5/3}, {1/3, 2/3, 1/3, 5/3}, {1/3, 5/
  3, -(2/3), 5/3}, {1/3, 5/3, 1/3, 2/3}}

But the above result is inconsistent with the one given by GAP, as shown below:

gap> x:=[[0, 1, 1, 1],
> [1, 0, 1, 1],
> [1, 1, 0, 1],
> [1, 1, 1, 0]];
[ [ 0, 1, 1, 1 ], [ 1, 0, 1, 1 ], [ 1, 1, 0, 1 ], [ 1, 1, 1, 0 ] ]
gap> A:=RationalCanonicalFormTransform(x);
[ [ 0, 0, 0, 1 ], [ 1, 0, 0, 1 ], [ -1, 1, 0, 1 ], [ 0, -1, 1, 0 ] ]
gap> A^-1 * x * A;
[ [ -1, 0, 0, 0 ], [ 0, -1, 0, 0 ], [ 0, 0, 0, 3 ], [ 0, 0, 1, 2 ] ]

Any tips for fixing this problem will be appreciated.

Regards, Zhao

POSTED BY: Hongyi Zhao

7 Replies

Sort By:

Daniel Lichtblau

Daniel Lichtblau, Wolfram Research

Posted 2 years ago

Can be done like so, using the resource function `PolynomialSmithDecomposition`. I changed slightly the definition of `companionMatrix` to evade `SparseArray` complaints. mat = {{0, 1, 1, 1}, {1, 0, 1, 1}, {1, 1, 0, 1}, {1, 1, 1, 0}}; companionMatrix[p_, x_] := Module[{n, w = CoefficientList[p, x]}, w = -w/Last[w]; n = Length[w] - 1; Normal@ SparseArray[{{i_, n} :> w[[i]], {i_, j_} /; i == j + 1 -> 1}, {n, n}]] {uu, ss, vv} = ResourceFunction["PolynomialSmithDecomposition"][ mat - tIdentityMatrix[Length[mat]], t]; diag = Diagonal[ss]; compmats = Map[companionMatrix[#, t] &, diag]; compmat = Normal[SparseArray[Band[{1, 1}] -> (compmats /. {} -> Nothing)]] ( Out[17]= {{-1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, 0, 3}, {0, 0, 1, 2}} *)

Can be done like so, using the resource function PolynomialSmithDecomposition. I changed slightly the definition of companionMatrix to evade SparseArray complaints.

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

companionMatrix[p_, x_] := 
 Module[{n, w = CoefficientList[p, x]}, w = -w/Last[w];
  n = Length[w] - 1;
  Normal@
   SparseArray[{{i_, n} :> w[[i]], {i_, j_} /; i == j + 1 -> 1}, {n, 
     n}]]

{uu, ss, vv} = 
  ResourceFunction["PolynomialSmithDecomposition"][
   mat - t*IdentityMatrix[Length[mat]], t];
diag = Diagonal[ss];
compmats = Map[companionMatrix[#, t] &, diag];

compmat = 
 Normal[SparseArray[Band[{1, 1}] -> (compmats /. {} -> Nothing)]]

(* Out[17]= {{-1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, 0, 3}, {0, 0, 1, 2}} *)

POSTED BY: Daniel Lichtblau

Bob Sandheinrich

Bob Sandheinrich, Wolfram Research

Posted 2 years ago

This looks like a cloud object paclet loading issue that I've seen before. Try it again after <<CloudObject` Or Quit

POSTED BY: Bob Sandheinrich

Hongyi Zhao

Posted 2 years ago

Thank you for pointing out this trick. But I met the following failure: In[47]:= mat = {{0, 1, 1, 1}, {1, 0, 1, 1}, {1, 1, 0, 1}, {1, 1, 1, 0}}; companionMatrix[p_, x_] := Module[{n, w = CoefficientList[p, x]}, w = -w/Last[w]; n = Length[w] - 1; Normal@ SparseArray[{{i_, n} :> w[[i]], {i_, j_} /; i == j + 1 -> 1}, {n, n}]] {uu, ss, vv} = ResourceFunction["PolynomialSmithDecomposition"][ mat - t*IdentityMatrix[Length[mat]], t]; diag = Diagonal[ss]; compmats = Map[companionMatrix[#, t] &, diag]; compmat = Normal[SparseArray[Band[{1, 1}] -> (compmats /. {} -> Nothing)]] During evaluation of In[47]:= URLFetch::invurl: CloudObject[https://www.wolframcloud.com/obj/resourcesystem/api/1.0/AcquireResource] is not a valid URL During evaluation of In[47]:= ResourceAcquire::rsunavail: Failed to receive data from the resource system. During evaluation of In[47]:= Throw::nocatch: Uncaught Throw[$Failed] returned to top level. Out[49]= Hold[Throw[$Failed]] During evaluation of In[47]:= Diagonal::list: List expected at position 1 in Diagonal[ss]. Out[52]= {{}}

Thank you for pointing out this trick. But I met the following failure:

In[47]:= mat = {{0, 1, 1, 1}, {1, 0, 1, 1}, {1, 1, 0, 1}, {1, 1, 1, 0}};

companionMatrix[p_, x_] := 
 Module[{n, w = CoefficientList[p, x]}, w = -w/Last[w];
  n = Length[w] - 1;
  Normal@
   SparseArray[{{i_, n} :> w[[i]], {i_, j_} /; i == j + 1 -> 1}, {n, 
     n}]]

{uu, ss, vv} = ResourceFunction["PolynomialSmithDecomposition"][ mat - t*IdentityMatrix[Length[mat]], t];
diag = Diagonal[ss];
compmats = Map[companionMatrix[#, t] &, diag];

compmat = Normal[SparseArray[Band[{1, 1}] -> (compmats /. {} -> Nothing)]]

During evaluation of In[47]:= URLFetch::invurl: CloudObject[https://www.wolframcloud.com/obj/resourcesystem/api/1.0/AcquireResource] is not a valid URL

During evaluation of In[47]:= ResourceAcquire::rsunavail: Failed to receive data from the resource system.

During evaluation of In[47]:= Throw::nocatch: Uncaught Throw[$Failed] returned to top level.

Out[49]= Hold[Throw[$Failed]]

During evaluation of In[47]:= Diagonal::list: List expected at position 1 in Diagonal[ss].

Out[52]= {{}}

POSTED BY: Hongyi Zhao

Hongyi Zhao

Posted 2 years ago

I tried again and it works. But this time the check of divisibility between polynomials failed to simplify: In[80]:= mat = {{0, 1, 1, 1}, {1, 0, 1, 1}, {1, 1, 0, 1}, {1, 1, 1, 0}}; companionMatrix[p_, x_] := Module[{n, w = CoefficientList[p, x]}, w = -w/Last[w]; n = Length[w] - 1; Normal@ SparseArray[{{i_, n} :> w[[i]], {i_, j_} /; i == j + 1 -> 1}, {n, n}]] {u, s, v} = ResourceFunction["PolynomialSmithDecomposition"][ mat - tIdentityMatrix[Length[mat]], t]; diag = Diagonal[s]; compmats = Map[companionMatrix[#, t] &, diag]; compmat = Normal[SparseArray[Band[{1, 1}] -> (compmats /. {} -> Nothing)]] Out[85]= {{-1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, 0, 3}, {0, 0, 1, 2}} In[86]:= (Check the matrix equation:) Expand[u . mat . v - s] (Check that u and v are unimodular; that is, they are constants:) {Det[u], Det[v]} (Check that the diagonal entries of s divide into subsequent entries:*) Apply[Divisible, Reverse[#]] & /@ Partition[Diagonal[s], 2, 1]//Simplify Out[86]= {{0, 0, 0, t}, {0, -t, 0, 0}, {0, 0, -t, 0}, {-t, 0, 0, 2 t - 2 t^2}} Out[87]= {1, 1} Out[88]= {Divisible[1 + t, 1], True, Divisible[-3 - 2 t + t^2, 1 + t]}

I tried again and it works. But this time the check of divisibility between polynomials failed to simplify:

In[80]:= mat = {{0, 1, 1, 1}, {1, 0, 1, 1}, {1, 1, 0, 1}, {1, 1, 1, 0}};

companionMatrix[p_, x_] := 
 Module[{n, w = CoefficientList[p, x]}, w = -w/Last[w];
  n = Length[w] - 1;
  Normal@
   SparseArray[{{i_, n} :> w[[i]], {i_, j_} /; i == j + 1 -> 1}, {n, 
     n}]]

{u, s, v} = ResourceFunction["PolynomialSmithDecomposition"][ mat - t*IdentityMatrix[Length[mat]], t];
diag = Diagonal[s];
compmats = Map[companionMatrix[#, t] &, diag];

compmat = Normal[SparseArray[Band[{1, 1}] -> (compmats /. {} -> Nothing)]]

Out[85]= {{-1, 0, 0, 0}, {0, -1, 0, 0}, {0, 0, 0, 3}, {0, 0, 1, 2}}

In[86]:= (*Check the matrix equation:*)
Expand[u . mat . v - s]
(*Check that u and v are unimodular; that is, they are constants:*)
{Det[u], Det[v]}
(*Check that the diagonal entries of s divide into subsequent entries:*)
Apply[Divisible, Reverse[#]] & /@ Partition[Diagonal[s], 2, 1]//Simplify

Out[86]= {{0, 0, 0, t}, {0, -t, 0, 0}, {0, 0, -t, 0}, {-t, 0, 0, 
  2 t - 2 t^2}}

Out[87]= {1, 1}

Out[88]= {Divisible[1 + t, 1], True, Divisible[-3 - 2 t + t^2, 1 + t]}

POSTED BY: Hongyi Zhao

Daniel Lichtblau

Daniel Lichtblau, Wolfram Research

Posted 2 years ago

Your first check appears to be incorrect. The second one is using a method from the WFR reference page that will not work in general. I blame the author of the function for that.

POSTED BY: Daniel Lichtblau

Hongyi Zhao

Posted 2 years ago

I changed it to the following version, and it seems to work: In[691]:= (Check the matrix equation:) Expand[u . mat . v - s] (Check that u and v are unimodular; that is, they are constants:) {Det[u], Det[v]} (Check that the diagonal entries of s divide into subsequent entries:) PolynomialRemainder[Sequence@@Reverse[#],t]& /@ Partition[Diagonal[s], 2, 1]//Simplify Out[691]= {{0, 0, 0, t}, {0, -t, 0, 0}, {0, 0, -t, 0}, {-t, 0, 0, 2 t - 2 t^2}} Out[692]= {1, 1} Out[693]= {0, 0, 0}

I changed it to the following version, and it seems to work:

In[691]:= (*Check the matrix equation:*)
Expand[u . mat . v - s]
(*Check that u and v are unimodular; that is, they are constants:*)
{Det[u], Det[v]}
(*Check that the diagonal entries of s divide into subsequent entries:*)
PolynomialRemainder[Sequence@@Reverse[#],t]& /@ Partition[Diagonal[s], 2, 1]//Simplify

Out[691]= {{0, 0, 0, t}, {0, -t, 0, 0}, {0, 0, -t, 0}, {-t, 0, 0, 
  2 t - 2 t^2}}

Out[692]= {1, 1}

Out[693]= {0, 0, 0}

POSTED BY: Hongyi Zhao

Hongyi Zhao

Posted 2 years ago

See the related discussion here.

POSTED BY: Hongyi Zhao

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback