Searching for Roots

The default procedure for finding roots is a three step process:

  1. Searching for sensible upper and lower bounds bracketing all roots
  2. Bisecting the bracket above to obtain bracket with a single root each
  3. Using a standard root-finding method to search for the root within that bracket

The first step is carried out by iteratively computing the eigenvalues of \(H(s)\) and setting the candidate lower(upper) boundary below(above) the smallest(largest) eigenvalue. Additionally, the upper boundary is set to a value such that \(\mathop{det}H(s)\) is strictly positive. There are two convenience functions to encapsulate this functionality, find_upper_bound_for_roots() and find_lower_bound_for_roots().

python:
from numpy import all
from dcprogs.likelihood import eig
from dcprogs.likelihood import find_upper_bound_for_roots, find_lower_bound_for_roots,        \
                               find_root_intervals, brentq, find_roots, QMatrix, DeterminantEq

qmatrix = QMatrix([ [-3050,        50,  3000,      0,    0], 
                    [2./3., -1502./3.,     0,    500,    0], 
                    [   15,         0, -2065,     50, 2000], 
                    [    0,     15000,  4000, -19000,    0], 
                    [    0,         0,    10,      0,  -10] ], 2)
det = DeterminantEq(qmatrix, 1e-4);


upper_bound = find_upper_bound_for_roots(det);
lower_bound = find_lower_bound_for_roots(det);

get_eigenvalues = lambda s: eig(det.H(s))[0].T
assert all(get_eigenvalues(lower_bound) > lower_bound) 
assert all(get_eigenvalues(upper_bound) < upper_bound)

Note

This package can be compiled to use 128bit reals. However, numpy does not provide all of its linear algebra utilities for such type. As consequence, this package exposes some of the functionality that it needs for its reals.
c++11:
#include <iostream>

#include <likelihood/determinant_equation.h>
#include <likelihood/root_finder.h>
#include <likelihood/brentq.h>
 
int main() {

  // Define parameters.
  DCProgs::t_rmatrix matrix(5 ,5);
  matrix << -3050,        50,  3000,      0,    0, 
            2./3., -1502./3.,     0,    500,    0,  
               15,         0, -2065,     50, 2000,  
                0,     15000,  4000, -19000,    0,  
                0,         0,    10,      0,  -10;
  DCProgs::QMatrix qmatrix(matrix, /*nopen=*/2);
  DCProgs::DeterminantEq det(qmatrix, 1e-4);

  // Find upper and lower bound
  DCProgs::t_real upper_bound = DCProgs::find_upper_bound_for_roots(det);
  DCProgs::t_real lower_bound = DCProgs::find_lower_bound_for_roots(det);


  // computes eigenvalues of H(s) for given s
  auto get_eigenvalues = [&det](DCProgs::t_real _s) -> DCProgs::t_rvector {
    return Eigen::EigenSolver<DCProgs::t_rmatrix>(det.H(_s)).eigenvalues().real();
  };

  // Checks bounds are correct.
  if((get_eigenvalues(lower_bound).array() < lower_bound).any()) 
    throw DCProgs::errors::Runtime("Incorrect lower bound.");
  if((get_eigenvalues(upper_bound).array() > upper_bound).any()) 
    throw DCProgs::errors::Runtime("Incorrect upper bound.");

  return 0;
}

The snippets above check that upper and lower bounds are indeed upper and lower bound, as advertised. It is possible that overflow errors make it difficult to find an upper or lower bound. It is possible in most function and classes to pass actual values for the upper and lower bounds (rather than the default, None).

The second step of the process is to bisect the input bracket until intervals are found which contain a single root (e.g. a single eigenvalue of H(s)).

python:
intervals = find_root_intervals(det, lower_bound, upper_bound)
c++11:
  std::vector<DCProgs::RootInterval> intervals
    = DCProgs::find_root_intervals(det, lower_bound, upper_bound);

The third step is performed by calling (by default) the brentq() subroutine. This method is copied straight from Scipy, with some modifications to allow it to cope with long doubles, if need be.

python:
for (start, end), multiplicity in intervals:
  root, iterations, function_calls = brentq(det, start, end)
  print("  * Root interval: [{0}, {1}]\n"
        "    Corresponding root: {2}\n".format(start, end, root))
c++11:
  for(DCProgs::RootInterval const& interval: intervals) {
    auto brentq_result = DCProgs::brentq(det, interval.start, interval.end);
    std::cout << "  * Root interval: [" << interval.start << ", " << interval.end << "]\n" 
              << "    Corresponding root: " << std::get<0>(brentq_result) << "\n\n";
  }

The whole procedure is encapsulated within the function find_roots(). On top of the parameters in the snippets below, find_roots() can take variety of parameters to control the root-finding procedure. Most notably, it accepts lower_bound and upper_bound keywords, allowing users to by-pass the first step if need be.

python:
roots = find_roots(det);
print("  * All roots: {0}\n".format([root for root, multiplicity in roots]))
c++11:
  std::vector<DCProgs::Root> roots = DCProgs::find_roots(det);
  std::cout <<  "  * All roots: ";
  for(DCProgs::Root const &root: roots) std::cout << root.root << " ";
  std::cout << "\n";