# Build a Companion Matrix of a Polynomial?

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[
~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?

• Have you seen the docs and MathWorld? – J. M. will be back soon Apr 8 '17 at 20:16
• Ah, thank you @J.M. Always nice when useful stuff like this gets hidden away. Apparently my Google-Fu is not up to snuff. – erfink Apr 8 '17 at 20:19

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 =
NRootsCompanionMatrix[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

• Awesome, thanks. I hadn't seen the Fiedler Companion matrix before! – erfink Apr 8 '17 at 23:13
• This is the Fiedler eigenvector Fiedler, I assume? – Igor Rivin Apr 9 '17 at 2:59
• The very same Fiedler, @Igor. – J. M. will be back soon Apr 9 '17 at 3:25

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 Characteristic Polynomial 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`.