Build a Companion Matrix of a Polynomial?

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

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.


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]