Myrkgard logo
Series logo MATH000 Finding Roots with Eigen
Thu, 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)
    */
}
plot

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.