# How do iterative methods applied to the companion matrix of a polynomial $p(\lambda)$ relate to $p$ itself?

A few days ago, I had a vague question in my mind about "matrix methods" for finding roots of a polynomial. Now I can ask at least a semi-precise question, thanks to the post How to calculate complex roots of a polynomial. I do not have a reference for this, but the idea behind the stated method is as follows. First create the companion matrix $C$ of the given polynomial $p$, whose characteristic equation is $p$ itself. Now, the eigenvalues of $C$ will give the roots of $p$.

Take a favorite iterative procedure for computing eigenvalues; for concreteness, I'll restrict to the QR algorithm, or the power iteration for the dominant eigenvalue. Does applying this procedure to $C$ correspond to a natural root-finding algorithm for $p$?

Surely such methods can be translated into methods to be applied directly to the polynomial when the matrix is just the companion matrix. Of course, they can be rather involved in description in that case.

I was quite surprised by this, since I was actually expecting a simple and intuitive translation. Is there one?

-
For the QR algorithm (at least, as first designed by Francis); it's very unlikely since the Frobenius companion matrix starts out sparse and is usually dense (only above the subdiagonal, of course) after the first similarity transformation. That is, the entries of the matrix are rather complicated functions of the polynomial's coefficients. On the other hand, there is now a new version of the QR algorithm that is specially designed to exploit the sparsity of the Frobenius companion matrix, but I've no experience with that... – J. M. Jul 25 '11 at 5:13

Take a favorite iterative procedure for computing eigenvalues; for concreteness, I'll restrict to ... the power iteration for the dominant eigenvalue. Does applying this procedure to $C$ correspond to a natural root-finding algorithm for $p$?

I'll restrict myself to power iteration proper as well, lest this answer get too long. Most of the variants of power iteration translate readily to a corresponding operation on polynomials, which I might discuss later in a separate answer if there's interest. (I've already said something on the QR algorithm in the comments.)

I'll consider the following form of the (upper Hessenberg) Frobenius companion matrix for $p(x)=c_n x^n +\cdots +c_1 x+c_0$:

$$\mathbf C=\begin{pmatrix}-\frac{c_{n-1}}{c_n}&\cdots&-\frac{c_1}{c_n}&-\frac{c_0}{c_n}\\1&&&\\&\ddots&&\\&&1&\end{pmatrix}$$

which is similar to another form you might be accustomed to:

$$\mathbf T\mathbf C\mathbf T^{-1}=\begin{pmatrix}&&&-\frac{c_0}{c_n}\\1&&&-\frac{c_1}{c_n}\\&\ddots&&\vdots\\&&1&-\frac{c_{n-1}}{c_n}\end{pmatrix}$$

with the similarity transformation $\mathbf T$ given by the unit upper triangular Toeplitz matrix

$$\mathbf T=\begin{pmatrix}1&\frac{c_{n-1}}{c_n}&\cdots&\frac{c_1}{c_n}\\&\ddots&\ddots&\vdots\\&&\ddots&\frac{c_{n-1}}{c_n}\\&&&1\end{pmatrix}$$

Now, the key is that there is an intimate relationship between polynomials, difference equations, and the Frobenius companion matrix; to wit, the difference equation

$$c_n u_{k+n}+\dots+c_1 u_{k+1}+c_0 u_k=0$$

has $p(x)$ as its characteristic polynomial, and its solutions $u_k$ take the form

$$u_k=w_1 \mu_1^k+\dots+w_n \mu_n^k$$

where the $\mu_j$ satisfy $p(\mu_j)=0$ with $|\mu_1|\geq \cdots \geq |\mu_n|$, and the $w_j$ are arbitrary constants. For the rest of this answer, I assume that there is only one dominant root $\mu_1$ (thus, $\mu_1$ is necessarily real and simple).

As $k\to \infty$, the ratio $\dfrac{u_{k+1}}{u_k}$ eventually tends to $\mu_1$, with the convergence rate depending on how far away the dominant root is from the other roots. Thus, one method to find the dominant root $\mu_1$ of a polynomial $p(x)$ would consist of iterating the recurrence relation for the $u_k$ (with starting values chosen such that $w_1$ in the expression for $u_k$ is nonzero) and forming successive ratios. This method is called the Bernoulli iteration.

The connection with companion matrices can be seen if we note that

$$\begin{pmatrix}u_{k+1}\\u_k\\\vdots\\u_{k-n+2}\end{pmatrix}=\begin{pmatrix}-\frac{c_{n-1}}{c_n}&\cdots&-\frac{c_1}{c_n}&-\frac{c_0}{c_n}\\1&&&\\&\ddots&&\\&&1&\end{pmatrix}\cdot\begin{pmatrix}u_k\\u_{k-1}\\\vdots\\u_{k-n+1}\end{pmatrix}$$

(see for instance this related question); we can thus consider Bernoulli's method as equivalent to repeatedly multiplying some arbitrary starting vector $\mathbf u_0$ with $\mathbf C$ ($\mathbf u_j=\mathbf C\mathbf u_{j-1}$), after which the dominant eigenvalue is estimated as the modified Rayleigh quotient

$$\frac{\mathbf e_1^\top\mathbf C\mathbf u_j}{\mathbf e_1^\top\mathbf u_j}$$

where $\mathbf e_1$ is the first column of the identity.

This leaves the problem of starting up the iteration. The customary way is to derive the starting values from the solutions of the equation

$$\mathbf T\mathbf u_0=\mathbf d$$

where $\mathbf d=(c_1\quad 2c_2\quad \cdots\quad n c_n)^\top$ represents the coefficients of $p^\prime(x)$.

As I mentioned earlier, most of the variants of power iteration, when specialized to the Frobenius companion matrix case, give rise to recursive methods for estimating the roots of a polynomial that directly act on the polynomial's coefficients. In particular, the Jenkins-Traub method, the modern descendant of the Bernoulli iteration, can be thought of either as a recursion on polynomials or on companion matrices. Some of the details are discussed in this paper.

Here's some sundry Mathematica code for playing around with Bernoulli's method:

n = 5; p = Expand[Times @@ (x - Range[n])];

cmat[poly_, x_] := Module[
{cl = CoefficientList[poly, x], n = Exponent[poly, x]},
SparseArray[{{1, i_} :> -cl[[n - i + 1]]/Last[cl],
Band[{2, 1}] -> 1}, {n, n}]]

start = LinearSolve[ToeplitzMatrix[UnitVector[n, 1],
(Reverse[Rest[#]/Last[#]]) &[CoefficientList[p, x]]],
CoefficientList[D[p, x], x]];

With[{prec = 20, its = 20},
Ratios[Drop[LinearRecurrence[
-(Reverse[Most[#]/Last[#]]) &[CoefficientList[p, x]],
Reverse[N[start, prec]], its + n + 1], n - 1]]]

With[{prec = 20, its = 20},
Table[((UnitVector[n, 1].cmat[p, x].#)/(UnitVector[n, 1].#)) &[
MatrixPower[cmat[p, x], k, N[start, prec]]], {k, 0, its}]]


You can change p here to any polynomial whose largest (in magnitude) root is real (of course, you must then set n to be the degree of p).

-
 Thanks J.M. It's quite a nice coincidence that I was talking about this question to a friend just 2 days back. I think I need some more background to understand your answer. So I'll take some more days to read it carefully. – Srivatsan Sep 29 '11 at 20:09 No problem. It took me days to figure out how to write a short answer for this question, as my first version was too long and rambling. ;) Just ping me if you've queries... – J. M. Sep 29 '11 at 23:30 Thanks @J.M for the nice answer. The relationship seems a lot clearer now. (I finally sat down and read the answer. :)) – Srivatsan Oct 10 '11 at 20:50