Build a Companion Matrix of a Polynomial?

Question

The eigenvalues of a square matrix $A$ are the roots of its characteristic polynomial $\chi(\lambda)$. Conversely, if we have a monic polynomial $p(\lambda)=a_0 + a_1 \lambda + \cdots + a_{n-1}\lambda^{n-1} + \lambda^n,$ we can define a companion matrix of the polynomial $p$:

$C=\begin{bmatrix} 0 & 0 & \dots & 0 & -a_0 \\ 1 & 0 & \dots & 0 & -a_1 \\ 0 & 1 & \dots & 0 & -a_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \dots & 1 & -a_{n-1} \end{bmatrix}.$

$C$ will have the property that its characteristic polynomial will recover $p$; i.e., $\chi(\lambda) = p(\lambda)$.

Question #1: Is there a built-in function to produce the companion matrix for a polynomial? I don't seem to be finding one.

In the meantime, I have cobbled together the following:

p[x_] := With[{maxDegree = 10000}, 
         Sum[a[n] x^n, {n, 0, maxDegree}] + x^maxDegree]

companionMatrix[poly_]:= With[
    {coeff = -Drop[CoefficientList[poly[x], x],-1]}, 
    Transpose[
    PadLeft[IdentityMatrix[Length@coeff -1],{Length@coeff-1,Length@coeff} ]
     ~Join~ {coeff}
    ]]

companionMatrix[p]//MatrixForm 

Question #2: Is this code reasonable? Is there a better way?

My code seems to run reasonably fast ($\approx 0.03$ seconds when I bump $p$ up to a degree 1000 monic polynomial, $\approx 4.8$ seconds for degree 10,000). However, there are a lot of zeroes in this matrix. I presume that this would be much much more efficient as a SparseArray?

Answer 1

Here is a slight variation of your method for generating a (Frobenius) companion matrix. This version also yields an upper Hessenberg matrix, but has the (monicized) coefficients appear at the top instead of at the rightmost part of the matrix:

frobeniusCompanion[poly_, x_] /; PolynomialQ[poly, x] := 
         Module[{n = Exponent[poly, x], coef},
                coef = CoefficientList[poly, x]; coef = -Most[coef]/Last[coef];
                SparseArray[{{1, j_} :> coef[[-j]], Band[{2, 1}] -> 1}, {n, n}]]

An example:

MatrixForm[mat = frobeniusCompanion[5 + 4 x + 3 x^2 + 2 x^3 + x^4 + x^5, x]]

$$\begin{pmatrix} -1 & -2 & -3 & -4 & -5 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ \end{pmatrix}$$

CharacteristicPolynomial[mat, x]
   -5 - 4 x - 3 x^2 - 2 x^3 - x^4 - x^5

---

As it turns out, there is a built-in, but undocumented function for generating the Frobenius companion matrix from a polynomial (in a slightly different format):

MatrixForm[mat =
           NRoots`CompanionMatrix[CoefficientList[5 + 4 x + 3 x^2 + 2 x^3 + x^4 + x^5, x]]

$$\begin{pmatrix} 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ -5 & -4 & -3 & -2 & -1 \\ \end{pmatrix}$$

Note that this yields a lower Hessenberg matrix, with the monicized coefficients appearing at the bottom. Transposing this matrix will yield the format given in the OP.

---

It should be noted that the Frobenius companion matrix is not the only possible companion matrix for a polynomial. There is a whole class of "congenial matrices" (see this as well). One interesting recent development is a pentadiagonal companion matrix due to Miroslav Fiedler:

fiedlerCompanion[poly_, x_] /; PolynomialQ[poly, x] := 
       Module[{n = Exponent[poly, x], coef},
              coef = CoefficientList[poly, x]; coef = -Most[coef]/Last[coef];
              SparseArray[{{1, 1} -> coef[[-1]], 
                           Band[{1, 2}] -> Riffle[Take[coef, {-2, 1, -2}], 0], 
                           Band[{2, 1}] -> Prepend[Riffle[Take[coef, {-3, 1, -2}], 0], 1],
                           {j_, k_} /; (j - k == 2 && EvenQ[j]) ||
                                       (k - j == 2 && OddQ[k]) :> 1}, {n, n}]]

Using the same polynomial as above,

MatrixForm[mat = fiedlerCompanion[5 + 4 x + 3 x^2 + 2 x^3 + x^4 + x^5, x]]

$$\begin{pmatrix} -1 & -2 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & -3 & 0 & -4 & 1 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & -5 & 0 \\ \end{pmatrix}$$

CharacteristicPolynomial[mat, x]
   -5 - 4 x - 3 x^2 - 2 x^3 - x^4 - x^5

Answer 2

Following @J.M.'s comment, we have the following functions.

Turning the code from MathWorld into a function gives:

mathWorldCompanionMatrix[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}]]

Note too that the mathWorldCompanionMatrix function will normalize an arbitrary polynomial to be monic.

Turning the code from the documentation for CharacteristicPolynomial into a function gives:

docsCompanionMatrix[p_, x_] := Module[{n, cl = CoefficientList[p, x]},
     n = Length[cl] - 1;
     SparseArray[{{i_, n} :> -cl[[i]], Band[{2, 1}] -> 1}, {n, n}]]

While all of these give the same result (and both of these functions are much cleaner), I get the following timing results:

p = RandomInteger[9, 10000].x^Range[0, 9999] + x^10000; (*test function*)

companionMatrix[p]; // Timing (*function from the original question*)
   {0.625, Null}

docsCompanionMatrix[p,x];//Timing
   {47.375, Null}

mathWorldCompanionMatrix[p,x];//Timing
   {104.125, Null}

A perhaps unreasonably large test polynomial, but revealing nonetheless. I'm not an expert, but presumably the functions in MathWorld and the Docs waste a lot of time on the pattern matching / Mathematica is fast at running IdentityMatrix. Both of these new functions also return sparse arrays, which can also be achieved by

SparseArray[companionMatrix[p]];//Timing
    {0.641, Null}

Presumably an even better version could be built using ArrayFlatten.

Answer 3

I've decided to write another answer, as the following generalizes the companion matrices featured in my other answer.

De Terán, Dopico, and Mackey, in their paper, show a method to construct a generalized Fiedler companion matrix, which has the Frobenius and Fiedler matrices in my other answer as special cases. Here is a Mathematica implementation of their algorithm:

generalizedFiedler[perm_, poly_, x_] /; PolynomialQ[poly, x] &&
                                        Sort[perm] == Range[Exponent[poly, x]] := 
           Module[{n = Exponent[poly, x], coef, diff, idx},
                  coef = CoefficientList[poly, x]; coef = -Most[coef]/Last[coef];
                  diff = Sign[Differences[InversePermutation[perm]]];
                  idx = Thread[If[diff[[1]] > 0, Identity, Reverse]
                               [{{2, 1}, {1, 1}, {1, 2}}] -> Append[Take[coef, 2], 1]];
                  Do[idx = If[diff[[k + 1]] > 0,
                     Join[Thread[{{1, 1}, {1, 2}} -> {coef[[k + 2]], 1}], 
                          MapAt[# + {1, Boole[Last[#] > 1]} &, #, 1] & /@ idx],
                     Join[Thread[{{1, 1}, {2, 1}} -> {coef[[k + 2]], 1}], 
                          MapAt[# + {Boole[First[#] > 1], 1} &, #, 1] & /@ idx]],
                     {k, n - 2}];
                  SparseArray[idx, {n, n}]]

perm can be any permutation of Range[n], where n is the polynomial degree.

Using a cubic polynomial as an example:

MatrixForm[generalizedFiedler[#, z^3 + C[2] z^2 + C[1] z + C[0], z]] & /@ 
Permutations[Range[3]]

$$\small \left\{\begin{pmatrix} -c_2 & 1 & 0 \\ -c_1 & 0 & 1 \\ -c_0 & 0 & 0 \\ \end{pmatrix}, \begin{pmatrix} -c_2 & -c_1 & 1 \\ 1 & 0 & 0 \\ 0 & -c_0 & 0 \\ \end{pmatrix}, \begin{pmatrix} -c_2 & 1 & 0 \\ -c_1 & 0 & -c_0 \\ 1 & 0 & 0 \\ \end{pmatrix}, \begin{pmatrix} -c_2 & 1 & 0 \\ -c_1 & 0 & -c_0 \\ 1 & 0 & 0 \\ \end{pmatrix}, \begin{pmatrix} -c_2 & -c_1 & 1 \\ 1 & 0 & 0 \\ 0 & -c_0 & 0 \\ \end{pmatrix}, \begin{pmatrix} -c_2 & -c_1 & -c_0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ \end{pmatrix} \right\}$$

Note that the first and last entries are Frobenius matrices, while the second and the fourth are the usual Fiedler matrices.

Different permutations can yield the same generalized Fiedler companion matrix. To demonstrate, the following snippet groups permutations of $(1\;2\;3\;4)$ by whether they yield the same companion matrix for a quartic polynomial:

With[{n = 4}, 
     Values[Sort /@ GroupBy[{generalizedFiedler[#, z^n + Sum[C[j] z^j, {j, 0, n - 1}], z],
                             #} & /@ Permutations[Range[n]], First -> Last]]]
{{{1, 2, 3, 4}},
 {{1, 2, 4, 3}, {1, 4, 2, 3}, {4, 1, 2, 3}},
 {{1, 3, 2, 4}, {1, 3, 4, 2}, {3, 1, 2, 4}, {3, 1, 4, 2}, {3, 4, 1, 2}},
 {{1, 4, 3, 2}, {4, 1, 3, 2}, {4, 3, 1, 2}},
 {{2, 1, 3, 4}, {2, 3, 1, 4}, {2, 3, 4, 1}},
 {{2, 1, 4, 3}, {2, 4, 1, 3}, {2, 4, 3, 1}, {4, 2, 1, 3}, {4, 2, 3, 1}},
 {{3, 2, 1, 4}, {3, 2, 4, 1}, {3, 4, 2, 1}},
 {{4, 3, 2, 1}}}

---

In general, for a degree-$n$ polynomial p[x], we have the identities

frobeniusCompanion[p[x], x] == generalizedFiedler[Range[n, 1, -1], p[x], x]

and

perm = Flatten[If[OddQ[n], Identity, Reverse][GatherBy[Range[n], OddQ]]];
fiedlerCompanion[p[x], x] == generalizedFiedler[perm, p[x], x]