Finding the roots of a univariate function is a basic numerical procedure that nearly all engineers and scientists encounters at some point, because the root finding question addresses the basic equation,

by finding the roots to,

Algorithms

I will demonstrate three well known univariate root finding algorithms, Newton-Raphson, Ridders, and Secant. All three methods require a delegate to the function and a bracketed range where a root is expected to reside, additionally, Newton-Raphson requires a continuous, non-zero derivative of f(x)in the region of the root. None of these methods are recommended for finding roots of polynomials in general, although Newton-Raphson will work well if the polynomial isn’t ill-conditioned by containing an even number of duplicated roots (Bracketing such roots isn’t possible because the function doesn’t change sign in the region of these roots.), or by containing root pairs extremely close together.

C# Example with NMath

All of these roots finders have a similar interface and properties in NMath, so exchanging one with the other is easy. As seen in the example below, both the Ridders and Secant root finders implement the IOneVariableRootFinder interface, and the Newton-Raphson root finder implements the slightly different IOneVariableDRootFinder, where the ‘D’ stands for derivative.

Using Wolfram|Alpha’s online services, I’ve computed and graphed the roots of the following non-linear univariate function.

The code block below demonstrates the use of each of these root finders, and shows the how to define a function delegate using a lamda expression. I find lamda expressions are convenient for defining delegates to mathematical expressions and that they create very readable code. (As an aside, lamda expressions are a C# language feature, and are not part of the .NET framework. Therefore, applications targeting the .NET 2.0 framework can leverage lamda expressions as long as the compiler supports them, as VS9 and VS10 both do.).

using CenterSpace.NMath.Analysis;
 
void RootFindingExample()
    {
 
     // Create the function delegate with a lamda expression, 
     // and the interval which contains a root.
      OneVariableFunction f = 
        new OneVariableFunction(x => 9 * Math.Log(x) * Math.Sin(x) - x * x * x);
      Interval bracket = new Interval(1, 1.5, Interval.Type.Closed);
 
      IOneVariableDRootFinder nrRoots = new NewtonRaphsonRootFinder(3);
      Double rootnewton = nrRoots.Find(f, f.Derivative(), bracket);
      Console.WriteLine(String.Format("Newton-Raphson root is = {0}", rootnewton));
 
      IOneVariableRootFinder rootFinder = new RiddersRootFinder(3);
      Double rootridder = rootFinder.Find(f, bracket);
      Console.WriteLine(String.Format("Ridder's root is = {0}", rootridder));
 
      rootFinder = new SecantRootFinder(3);
      Double rootsecant = rootFinder.Find(f, bracket);
      Console.WriteLine(String.Format("Secant root is = {0}", rootsecant));
 
    }

With each root finder limited to just 3 iterations, this program creates the output,

Newton-Raphson root is = 1.26760302048967
Ridder's root is = 1.26760302063032
Secant root is = 1.3334994219432

Where Newton-Raphson found the first root most accurately with it’s knowledge of derivative information. Besides specifying the maximum number of iterations, a root tolerance can also be used to control convergence effort.

A second somewhat more difficult example which can cause problems for both the Secant and Newton-Raphson methods.

I gave the Secant method 1000 iterations, and still it’s root estimate was very poor – in fact with an initial bracket of [-5, -0.01], the Secant method will never converge on the root. There is a stable orbit around which the Secant method spins forever, preventing convergence. Because the Newton-Raphson requires a continuous derivative any starting bracket that contains the origin will cause an exception to be thrown because of the discontinuous derivative at zero.

Newton-Raphson root is = -1.05539938820662 (with 8 iterations)
Ridder's root is = -1.05539938548287 (with 8 iterations)
Secant root is = -0.0912954387917323 (with 1000 iterations)

As a general recommendation, I would start my root search using Ridder’s method and then finish up by polishing with Newton-Raphson.

-Happy Computing,
Paul

We are sometimes asked how to reproduce the various Excel Trendline types in NMath, including printing out the form of the equation and the R2 value (coefficient of determination). Excel offers these trend types:

• Linear
• Logarithmic
• Exponential
• Power
• Polynomial
• Moving Average

These can all be easily computed using NMath. Let’s see how. First, let’s start with a simple data series:

DoubleVector x = new DoubleVector(11, 15, 18, 23, 26, 31, 39, 44, 54, 64, 74);
DoubleVector y = new DoubleVector(0.00476, 0.0105, 0.0207, 0.0619, 0.337, 0.74, 1.7, 2.45, 3.5, 4.5, 5.09);

These data represent the evolution of an algal bloom in the Adriatic Sea. The x-values are time expressed in days, and the y-values are the size of the surface of the bloom in mm².

Linear

The Linear trendline fits a line, y = a*x + b, to the data. This is easily computed using the LinearRegression class in NMath Stats, as shown in the following C# code:

bool addIntercept = true;
LinearRegression lr = new LinearRegression(new DoubleMatrix(x), y, addIntercept);
LinearRegressionAnova lrAnova = new LinearRegressionAnova(lr);
double a = lr.Parameters[0];
double b = lr.Parameters[1];
double r2 = lrAnova.RSquared;
 
Console.WriteLine("y = {0}x + {1}", a, b);
Console.WriteLine("r2 = {0}", r2);

Output

y = -1.63908651994092x + 0.0913403802489977
r2 = 0.967828167696715

Note that to compute the coefficient of determination (R2), we construct a LinearRegressionAnova object from the LinearRegression instance. This class tests overall model significance for regressions computed by LinearRegression.

Logarithmic

The Logarithmic trendline fits a line to ln(x), y–that is, y = a*ln(x) + b. Again, we can use the LinearRegression class for this:

DoubleVector logX = NMathFunctions.Log(x);
lr = new LinearRegression(new DoubleMatrix(logX), y, addIntercept);
a = lr.Parameters[0];
b = lr.Parameters[1];
r2 = new LinearRegressionAnova(lr).RSquared;
 
Console.WriteLine("y = {0}ln(x) + {1}", a, b);
Console.WriteLine("r2 = {0}", r2);

Output:

y = -8.17805531083818ln(x) + 2.87283105614145
r2= 0.840393353070466

Exponential

The Exponential trendline fits the function y = a * e ^(b * x). This can also be computed using LinearRegression by fitting a line to x, ln(y). Taking the log of both sides of the function gives:

ln(y) = ln(a * e ^(b * x)) = ln(a) + bx

The following C# code does this:

DoubleVector logY = NMathFunctions.Log(y);
lr = new LinearRegression(new DoubleMatrix(x), logY, addIntercept);
a = Math.Exp(lr.Parameters[0]);
b = lr.Parameters[1];
r2 = new LinearRegressionAnova(lr).RSquared;
 
Console.WriteLine("y = {0}e^{1}x", a, b);
Console.WriteLine("r2 = {0}", r2);

Output:

y = 0.00552689525130073e^0.112876522212724x
r2 = 0.812366433832176

Note that because the intercept of the fitted line is the natural log of the “a” parameter in the exponential function, we need to take the exponential of the found intercept (using Math.Exp) to recover the value of “a”.

Power

The Power trendline fits the function y = a * x^b. This can be computed using LinearRegression by fitting a line to ln(x), ln(y). Taking the log of both sides of the equation gives:

ln(y) = ln(a * x^b) = ln(a) + b * ln(x)

lr = new LinearRegression(new DoubleMatrix(logX), logY, addIntercept);
lrAnova = new LinearRegressionAnova(lr);
a = Math.Exp(lr.Parameters[0]);
b = lr.Parameters[1];
r2 = new LinearRegressionAnova(lr).RSquared;
 
Console.WriteLine("y = {0}x^{1}", a, b);
Console.WriteLine("r2 = {0}", r2);

Output:

y = 2.46993343563889E-07x^4.11443630332377
r2 = 0.947447653331871

Again, we recover the value of parameter “a” by taking the exponential of the found intercept.

Polynomial

NMath provides class PolynomialLeastSquares for fitting a polynomial of the specified degree to paired vectors of x- and y-values. For example, this code fits a cubic to our data:

int degree = 3;
PolynomialLeastSquares pls = new PolynomialLeastSquares(degree, x, y);
 
// compute R2
DoubleVector predictions = pls.FittedPolynomial.Evaluate(x);
double regressionSumOfSquares = StatsFunctions.SumOfSquares(predictions - StatsFunctions.Mean(y));
double residualSumOfSquares = pls.LeastSquaresSolution.Residuals.TwoNormSquared();
r2 = regressionSumOfSquares / (regressionSumOfSquares + residualSumOfSquares);
 
Console.WriteLine("y = " + pls.FittedPolynomial);
Console.WriteLine("r2 = {0}", r2);

Output:

y = -4.68278158094222E-05x^3 + 0.00640408381023593x^2 -
      0.1643720340709x + 1.12703300920657
r2 = 0.998634459376868

Note that PolynomialLeastSquares does not currently provide the R2 value, so in the code above we compute it directly.

Moving Average

Unlike the trendlines we’ve examined so far, a moving average does not fit a functional form to data. Rather, it filters data in order to smooth out noise. NMath provides class MovingWindowFilter for this purpose.

Class MovingWindowFilter replaces data points with a linear combination of the data points immediately to the left and right of each point, based on a given set of coefficients to use in the linear combination. Static class methods are provided for generating coefficients to implement a moving average filter and a Savitzky-Golay smoothing filter.

For example, the following C# code filters the data using a window of width 3 (the “period” parameter in Excel):

int numberLeft = 1;
int numberRight = 1;
DoubleVector movingAvgCoefficents = MovingWindowFilter.MovingAverageCoefficients(numberLeft, numberRight);
MovingWindowFilter filter = new MovingWindowFilter(numberLeft, numberRight, movingAvgCoefficents);
DoubleVector yfiltered = filter.Filter(y, MovingWindowFilter.BoundaryOption.DoNotFilterBoundaryPoints);
Console.WriteLine("yfiltered = " + yfiltered);

Output:

yfiltered = [ 0.00476 0.0119866666666667 0.0310333333333333
              0.139866666666667 0.379633333333333 0.925666666666667 1.63
              2.55 3.48333333333333 4.36333333333333 5.09 ]

Advanced Curve Fitting

That covers the simple trendlines produced by Excel. For advanced curve fitting, NMath provides class OneVariableFunctionFitter which fits arbitrary one variable functions to a set of points. (You must supply at least as many data points to fit as your function has parameters.) In a future post, we’ll take a look at this class. In the meantime, see this page for an example demonstrating the use of OneVariableFunctionFitter, including how to define your own function.

Ken

Application Note

Power Spectral Density Calculation

I’ve had several customers ask about computing the PSD in C# with NMath, so I though it was time for a post on the subject. The power spectral density provides an estimate of the power present within each slice of spectrum, and is presented as graph of the signal power versus frequency. It’s a common signal processing calculation across many fields from acoustics to chemistry, and can provide insight into periodicities contained within a time domain signal.

For stationary, square summable signals, the PSD is expressed as,

where F(w) is the Fourier transform of the time domain signal f(t), and T is the width of the time domain sampled signal. Naturally we can never sample the entire signal, so calculations of the PSD are all estimates. Techniques for estimating the PSD can be divided into two classes, parametric (model based), and nonparametic (non-model based). We will discuss only nonparametic techniques here. For discrete signals that are not square summable (i.e. non-stationary signals – and so the Fourier transform does not exist), estimates of the power spectral density can be derived from,

Which is the Fourier transform of the autocorrelation as the correlation width of the sampled signal tends to infinity. For a concise derivation of both of these formulas read these short lecture notes.

Computing the PSD using NMath

Concentrating on stationary (periodic) signals, the PSD is most efficiently computed by applying smoothing to discrete periodograms .

Where n is the number of signal samples. Each point in the periodogram represents the relative contribution to the variance in the time domain signal at that frequency. (Visualization provided courtesy of Infragistics.)

An example periodogram of sunspot data. This is smoothed in some fashion to estimate the PSD.

Another estimation technique involves computing multiple windowed periodograms and then combining these together to get a progressively more accurate estimate (Welch’s Method, similarly MTM with Slepian windows). The time domain signal should be detrended before any of these operations.
(more…)

A set of classes for working with complex numbers is not including in the .NET framework. These classes are frequently hand rolled by programmers to fill an immediate need, but this forces the developer into an on-going task developing compatible numeric algorithms with these custom classes. CenterSpace’s NMath libraries solve this issue by providing a framework with an extensive set of numeric classes that work natively with complex numbers.

Complex Number Classes

The CenterSpace NMath framework contains two classes to represent complex numbers at two different levels of precision. Additionally, classes are included for supporting vectors and matrices of complex numbers which are integrated into NMath’s extensive linear algebra classes.

FloatComplex
DoubleComplex
FloatComplexVector
DoubleComplexVector
FloatComplexMatrix
DoubleComplexMatrix

Each class supports a natural set of overloaded operations (more…)

A customer recently asked how to reproduce the results of the MATLAB polyfit function in NMath. This is easily done using NMath in C#.

For example, here’s some MATLAB code for fitting a polynomial of degree 6 to some data, and then evaluating the fitted polynomial at a point (x = 0.25).

 x = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 
        1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 2.2 2.3 2.4 2.5]
 
 y = [0 0.112462916018285 0.222702589210478 0.328626759459127 0.428392355046668 0.520499877813047 0.603856090847926 0.677801193837419 0.742100964707661 0.796908212422832 0.842700792949715 0.880205069574082 0.910313978229635 0.934007944940652 0.952285119762649 0.966105146475311 0.976348383344644 0.983790458590775 0.989090501635731 0.992790429235257 0.995322265018953 0.997020533343667 0.998137153702018 0.998856823402643 0.999311486103355 0.999593047982555]
 
p = polyfit(x,y,6)
p = 0.0084   -0.0983    0.4217   -0.7435    0.1471    1.1064    0.0004
f = polyval(p, 0.25); 
f = 0.2762

Here’s the equivalent C# NMath code:

#using CenterSpace.NMath.Analysis
public void FitData()
{
  DoubleVector x = new DoubleVector("[0 0.1 0.2 0.3 0.4 0.5 
    0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2 2.1 
    2.2 2.3 2.4 2.5]");
 
  DoubleVector y = new DoubleVector("[0 0.112462916018285 0.222702589210478 0.328626759459127   0.428392355046668 0.520499877813047 0.603856090847926 0.677801193837419 0.742100964707661 0.796908212422832 0.842700792949715 0.880205069574082 0.910313978229635 0.934007944940652 0.952285119762649 0.966105146475311 0.976348383344644 0.983790458590775 0.989090501635731 0.992790429235257 0.995322265018953 0.997020533343667 0.998137153702018 0.998856823402643 0.999311486103355 0.999593047982555]");
 
  PolynomialLeastSquares pls = new PolynomialLeastSquares( 6, x, y );
  Console.WriteLine( pls.Coefficients );
  Console.WriteLine();
  Console.WriteLine( pls.FittedPolynomial.Evaluate( 0.25 ) );
}

This creates the output.

[ 0.000441173961987631 1.10644604462547 0.147104056633056
  -0.743462849129717 0.421736169818002 -0.0982995753367523
  0.00841937177923124 ]

0.276183548385269

6th degree polynomial fit using Polynomial Least Squares.

Note that the order of the coefficient vector in NMath is reversed relative that returned from MATLAB’s polyfit() function. In NMath, the constant is at index 0 and the leading coefficient is at index Coefficients.Length - 1. If you wish, you can use the Reverse() method on the coefficient vector to reverse the ordering.

-Ken

Notes
Visualization provided in partnership with Nevron