com.opengamma.strata.math.impl.rootfinding.EigenvaluePolynomialRootFinder Maven / Gradle / Ivy

Go to download

Show more of this group Show more artifacts with this name
Show all versions of strata-math Show documentation

Mathematic support for Strata

There is a newer version: 2.12.46

/*
 * Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies
 *
 * Please see distribution for license.
 */
package com.opengamma.strata.math.impl.rootfinding;

import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.EigenDecomposition;
import org.apache.commons.math3.linear.RealMatrix;

import com.opengamma.strata.collect.ArgChecker;
import com.opengamma.strata.math.impl.function.RealPolynomialFunction1D;

/**
 * The eigenvalues of a matrix $\mathbf{A}$ are the roots of the characteristic
 * polynomial $P(x) = \mathrm{det}[\mathbf{A} - x\mathbb{1}]$. For a 
 * polynomial 
 * $$
 * \begin{align*}
 * P(x) = \sum_{i=0}^n a_i x^i
 * \end{align*} 
 * $$
 * an equivalent polynomial can be constructed from the characteristic polynomial of the matrix
 * $$
 * \begin{align*}
 * A = 
 * \begin{pmatrix}
 * -\frac{a_{m-1}}{a_m}  & -\frac{a_{m-2}}{a_m} & \cdots & -\frac{a_{1}}{a_m} & -\frac{a_{0}}{a_m} \\
 * 1                      & 0                     & \cdots & 0                   & 0                   \\
 * 0                      & 1                     & \cdots & 0                   & 0                   \\
 * \vdots                &                       & \cdots &                     & \vdots             \\
 * 0                      & 0                     & \cdots & 1                   & 0                   
 * \end{pmatrix}
 * \end{align*}
 * $$
 * and so the roots are found by calculating the eigenvalues of this matrix.
 */
public class EigenvaluePolynomialRootFinder implements Polynomial1DRootFinder {

  @Override
  public Double[] getRoots(RealPolynomialFunction1D function) {
    ArgChecker.notNull(function, "function");
    double[] coeffs = function.getCoefficients();
    int l = coeffs.length - 1;
    double[][] hessianDeref = new double[l][l];
    for (int i = 0; i < l; i++) {
      hessianDeref[0][i] = -coeffs[l - i - 1] / coeffs[l];
      for (int j = 1; j < l; j++) {
        hessianDeref[j][i] = 0;
        if (i != l - 1) {
          hessianDeref[i + 1][i] = 1;
        }
      }
    }
    RealMatrix hessian = new Array2DRowRealMatrix(hessianDeref);
    double[] d = new EigenDecomposition(hessian).getRealEigenvalues();
    Double[] result = new Double[d.length];
    for (int i = 0; i < d.length; i++) {
      result[i] = d[i];
    }
    return result;
  }

}