MATH000 Finding Roots with EigenThu, 16 Nov 2017 04:30:00 +0000

# Finding Roots with Eigen

Eigen is a C++ library for linear algebra. In this episode we'll use it to find the roots of univariate polynomials.

Okay, here's how to do it:

1. Say you have a polynomial $$p$$ for which you need to find all its roots
2. Create the Frobenius companion matrix $$A$$ from the polynomial's coefficients
3. Compute the eigenvalues of $$A$$ using the Eigen library
4. The eigenvalues $$\lambda_1, \lambda_2, \cdots, \lambda_n$$ of $$A$$ are the roots you're looking for

Different authors use different variants of the companion matrix. Here we'll define it as follows:

Let $$p$$ be a univariate polynomial $p(x) = a_n x^n + a_{n-1} x^{n-1} + \cdots + a_1 x + a_0 .$ Its companion matrix is $\mathbf{A} = \begin{pmatrix} -\frac{a_{n-1}}{a_n} & -\frac{a_{n-2}}{a_n} & \cdots & -\frac{a_1}{a_n} & -\frac{a_0}{a_n} \\ 1 & 0 & \cdots & 0 & 0 \\ 0 & 1 & \cdots & 0 & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & \cdots & 1 & 0 \end{pmatrix} .$

## Code Example

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>

Eigen::MatrixXd companion(Eigen::VectorXd const &c)
{
auto const rank = c.rows() - 1;
Eigen::MatrixXd A = Eigen::MatrixXd::Zero(rank, rank);

// compute first row
for (int i = 0; i < rank; ++i)
A(0, i) = - c(i + 1) / c(0);

// set diagonal below first row to ones
for (int i = 0; i < rank - 1; ++i)
A(i + 1, i) = 1.;

return A;
}

int main()
{
// p(x) = x^6 - 2x^5 - 26x^4 + 28x^3 + 145x^2 - 26x - 80
Eigen::VectorXd c(7);
c << 1., -2., -26., 28., 145., -26., -80.; // coefficients of p

auto const A = companion(c);
std::cout << A << std::endl;
/* output:
2   26  -28 -145   26   80
1    0    0    0    0    0
0    1    0    0    0    0
0    0    1    0    0    0
0    0    0    1    0    0
0    0    0    0    1    0
*/

auto const roots = A.eigenvalues();
std::cout << roots << std::endl;
/* values are of type std::complex. but in this example
they're all real (imag = 0). output:
(4.98656,0)
(-3.97823,0)
(3.06868,0)
(-2.16045,0)
(-0.739316,0)
(0.822758,0)
*/
}

## Manual Example

We want to find all roots of $p(x) = - 3x^2 + 4x + 7 \overset{!}{=} 0 .$ Here we can get the roots directly using the quadratic formula. $ax^2 + bx + c = 0$ $x_{1,2} = \frac{ -b \pm \sqrt{b^2 - 4ac} }{ 2a } .$ $\implies x_{1,2} = \frac{ -4 \pm \sqrt{4^2 - \left( 4 \cdot (-3) \cdot 7\right)} }{ 2 \cdot (-3) }$ $\implies x_1 = -1, x_2 = \frac{7}{3}$ Anyway, let's see if the "new trick" leads to the same results: $\mathbf{A} = \begin{pmatrix} -\frac{4}{-3} & -\frac{7}{-3} \\ 1 & 0 \end{pmatrix} = \begin{pmatrix} \frac{4}{3} & \frac{7}{3} \\ 1 & 0 \end{pmatrix}$ Compute the eigenvalues of $$A$$: $\mathbf{I_n} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$ $det\left(\lambda \mathbf{I_n} - \mathbf{A}\right) = det \begin{pmatrix} \lambda - \frac{4}{3} & - \frac{7}{3} \\ - 1 & \lambda \end{pmatrix} \overset{!}{=} 0$ $\implies \left(\lambda - \frac{4}{3}\right) \lambda - \left(- \frac{7}{3}\right) (-1) = 0$ $\iff \lambda^2 - \frac{4}{3}\lambda - \frac{7}{3} = 0$ Quadratic formula again: $\implies \lambda_{1,2} = \frac{ -\left(-\frac{4}{3}\right) \pm \sqrt{\left(-\frac{4}{3}\right)^2 - \left(4 \cdot 1 \cdot \left(-\frac{7}{3}\right)\right)} }{ 2 \cdot 1 }$ $\implies \lambda_1 = \frac{7}{3}, \lambda_2 = -1$ Same results as above, yay.