# Blog

## Optimal Portfolio Allocation

Wednesday, January 22nd, 2014

The problem of optimal portfolio allocation, in its simplest form, asks the question of how to fully allocate a given amount of wealth across a fixed universe of investments to achieve a minimum-risk goal-expected return. The known quantities are the potential field of investments, their performance history, and the goal rate of return; The unknown is the wealth allocation across the investments. In this standard formulation, the phrase ‘minimum-risk’ is understood to be equivalent to ‘minimum-variance’ of the portfolio return. Said another way, in this model the investor will expect to get a ‘risk-premium’ (higher return) for investing in assets which have a higher historic return variance. This understanding of risk is at the heart of Harry Markowitz’s portfolio theory [5].

In the following brief development of optimal portfolio allocation theory we follow the paper by P. J. Atzberger [2]. If we have n investments we represent a fully invested portfolio as,

where the fraction of wealth invested in the ith asset is wi. We use this normalization to avoid working with the absolute magnitudes of the assets and resultant portfolios.

The rate of return of the entire portfolio, ρp at a given time t is,

or simply the weighted sum of investment returns.

And because expectation is a linear operator, the expected rate of return of a given portfolio is,

Where μi is expected rate of return of the ith investment option.

And finally from the expected return of return, the portfolio variance will be given by,

where w is the vector of normalized fractional investments defined above. (See eq. (13) in [2] for the full derivation). This is the minimum theoretical framework needed to setup the optimization problem.

### Optimization with Markowitz Theory

Using our established problem framework we can cast the optimal portfolio allocation problem as a constrained quadratic optimization problem by minimizing the portfoilo’s variance (risk) subject to two linear constraints.

subject to,

and,

The minimized objective function is the variance of the portfolio – which is thought of as risk equivalent in this model. The first constraint insures that the built portfolio achieves a goal return of μp and the second constraint enforces the full allocation of wealth. There are many variants of this basic setup such as including a zero-risk investment (cash) to the pool of possible investments or adding constraints to allocate various percentages of wealth into different risk-classified investments [2][4].

### Solving the Quadratic Programming Problem

This constrained optimization problem can be solved using NMath‘s quadratic programming class, ActiveSetQPSolver(). The basic steps of setting up the QP solver involve (1) Reading and parsing the investment data, (2) Setting up the problem constraints, (3) Executing the solve. The data for this example come from a data repository for exploring various Operations Research (OR) problems, run by J.E. Beasley. The specific data we study here has a universe of 225 potential assets, is found here, and is the largest test data set used in the paper “Heuristics for cardinality constrained portfolio optimisation” [3] (A brief description of the data set can be found on the same server here). The data set provides the basic information we need for optimizing the portfolio – the mean return and standard deviation of return of each investment and the correlation between all possible pairs of assets. The code to parse this file is appended to the end of this article.

Before generating the QP’s constraints we need to set up a QuadraticProgrammingProblem instance.

 // Set up QR problem using port5.txt data FloatMatrix covariance = correlation * ( NMathFunctions.Product( new FloatMatrix( stdDevReturn ), ( new FloatMatrix( stdDevReturn ).Transpose() ) ) ); QuadraticProgrammingProblem problem = new QuadraticProgrammingProblem(); problem.H = covariance; problem.c = new DoubleVector(n, 0);

First, the covariance matrix is computed from the investment standard deviations and correlations. Next we create and initialize a QuadraticProgrammingProblem instance. Since this QP solver class solves programs of the form: Minimize x'Hx + x'c, the c variable is set to zero.

The following code snippet sets up the QP constraints.

 ... QuadraticProgrammingProblem problem = new QuadraticProgrammingProblem(); ... // Set up constraints LinearConstraint sumToOne = new LinearConstraint( new DoubleVector( n, 1 ), 1, ConstraintType.EqualTo ); LinearConstraint returnRate = new LinearConstraint( meanReturn, targetReturn, ConstraintType.GreaterThanOrEqualTo ); problem.AddConstraint( returnRate ); problem.AddConstraint( sumToOne );   // Add the bound on each var, w: w > 0 for ( int i = 0; i < n; i++ ) problem.AddLowerBound( i, 0.0 );

The first two constraints, sumToOne and returnRate were outlined above in the brief Markowitz theory section. The third constraint has been added to insure that we have a positive ownership in all assets; in this way we have disallowed shorting assets. This constraint can be deleted to include shorted asset in the optimized portfolio, if desired. All of these constraints are added to the QuadraticProgrammingProblem instance.

We can now put all of this together and solve the constrained 225 variable QP using the ActiveSetQPSolver class [0].

Asset Allocation 8, 0.0795 39, 0.0866 42, 0.0812 59, 0.1201 61, 0.2567 96, 0.0593 128, 0.0741 170, 0.0573 195, 0.0980 214, 0.0688 224, 0.0183 (All other assets have an allocation of zero.)   Total allocation: 1.0000 Total return-rate: 0.2000%, Target return-rate 0.2000% Total solver time required: 0.1830 seconds

The optimal allocation includes 11 of our 225 investments and closely adheres to all constraints. The complete code listing is provided at the end of this article.

Happy Computing,

Paul Shirkey

#### References, Resources, Notes

[0] Due to the size of this QP problem and the number of constraints (227) the current release of NMath‘s ActiveSetQPSolver() fails to correctly find the optimal portfolio. However, the next release of NMath due in May (NMath 6.0) includes a major upgrade to many of our LP, QP, and NLP classes. Many of these classes in NMath 6.0 are backed by the Microsoft Solver Foundation, providing the ability to solve much larger problems correctly and efficiently.

[1] The port5 data set can be found at http://people.brunel.ac.uk/~mastjjb/jeb/orlib/files/port5.txt . This data set is part of the OR-Library, a set of standard data sets maintained by J. E. Beasley at http://people.brunel.ac.uk/~mastjjb/jeb/info.html.

[2] This blog note follows the notation of this clearly written paper, by Paul J. Atzberger, on optimal portfolio allocation, variants, and the efficient frontier: http://www.math.ucsb.edu/~atzberg/finance/portfolioTheory.pdf.

[3] Chang, T.-J., Meade, N., Beasley, J.E. and Sharaiha, Y.M., “Heuristics for cardinality constrained portfolio optimisation” Comp. & Opns. Res. 27 (2000) 1271-1302.

[4] MATLAB uses this same data set in their optimization documentation.

[5] Markowitz, H.M. (March 1952). “Portfolio Selection”. The Journal of Finance 7 (1): 77–91. [paper in pdf]

### .NET Code Listings

Full QP solution code. Note that NMath 6.0 (due May 2014) is required.

using System; using CenterSpace.NMath.Core; using CenterSpace.NMath.Analysis; using CenterSpace.NMath.Charting.Microsoft; using System.Windows.Forms.DataVisualization.Charting;   public void OptimalPortAllocation() { // Get port5 data set : http://people.brunel.ac.uk/~mastjjb/jeb/orlib/files/port5.txt // File Format of port5: // Line 1: Text description of data // Line 2: number of assets, N // Line 2 - 2+N : list of mean return _ standard return // Remaining 1/2 N*N lines: correlation matrix between asset i and j FloatVector meanReturn = null ; FloatVector stdDevReturn = null; FloatMatrix correlation = null; int n = loadData( ref meanReturn, ref stdDevReturn, ref correlation );   double targetReturn = 0.002;   // Set up QR problem using port5.txt data FloatMatrix covariance = correlation * ( NMathFunctions.Product( new FloatMatrix( stdDevReturn ), ( new FloatMatrix( stdDevReturn ).Transpose() ) ) ); QuadraticProgrammingProblem problem = new QuadraticProgrammingProblem(); problem.H = covariance; problem.c = new DoubleVector(n, 0);   // Set up constraints LinearConstraint sumToOne = new LinearConstraint( new DoubleVector( n, 1 ), 1, ConstraintType.EqualTo ); LinearConstraint returnRate = new LinearConstraint( meanReturn, targetReturn, ConstraintType.GreaterThanOrEqualTo ); problem.AddConstraint( returnRate ); problem.AddConstraint( sumToOne );   // Add the bound on each var, w: w > 0 for ( int i = 0; i < n; i++ ) problem.AddLowerBound( i, 0.0 );   // Solve constrainted QP var solver = new ActiveSetQPSolver(); solver.StepSizeEpsilon = 1e-10;   var startEqualAllocation = new DoubleVector( n, 1.0 / n ); var optAllocation = new DoubleVector();   long before = DateTime.Now.Ticks; if ( solver.Solve( problem, startEqualAllocation ) ) { optAllocation = solver.OptimalX; } else { throw new Exception( "Solve failed to find a solution." ); } long after = DateTime.Now.Ticks; long ticks = after - before; double millis = ( (double) ticks ) / 10000;   Console.WriteLine("Asset # Allocation"); double wealth = 0.0; double measuredReturnRate = 0.0; for ( int i = 0; i < optAllocation.Length; i++ ) { if ( optAllocation[i] > 0.0000001 ) { Console.WriteLine( String.Format( "{0,4}, {1,7:0.0000}",i, optAllocation[i])); wealth += optAllocation[i]; measuredReturnRate += optAllocation[i] * meanReturn[i]; }   } Console.WriteLine(String.Format("Total allocation: {0,7:0.0000}", wealth) ); Console.WriteLine( String.Format( "Total return rate: {0,7:0.0000}%, Target return rate {1,7:0.0000}%", 100 * measuredReturnRate, 100 * targetReturn ) ); Console.WriteLine( String.Format( "Total solver time required: {0,7:0.0000} seconds", millis/1000 ) );     // Chart optimal portfolio allocation DoubleVector x = new DoubleVector(n, 0, 1); DoubleVector w = optAllocation; var chart = NMathChart.ToChart( x, w ); chart.Series[0].MarkerStyle = MarkerStyle.Square; chart.Series[0].ChartType = SeriesChartType.Bar; NMathChart.Show( chart );   }

Below is the full parsing code for the port5.txt data file.

 private int loadData( ref FloatVector meanReturn, ref FloatVector stdDevReturn, ref FloatMatrix correlation ) { int n;   // Fix up the path here to point to your data file copy. using ( StreamReader file = new StreamReader( "./data/PortfolioOptimzation_port5.data.txt" ) ) { file.ReadLine(); // skip first line of explanatory text n = Int32.Parse( file.ReadLine() );   // Now get the next n mean return and standard return values. meanReturn = new FloatVector( n ); stdDevReturn = new FloatVector( n ); string[] values = new string[3]; for ( int k = 0; k < n; k++ ) { values = file.ReadLine().Split( ' ' ); meanReturn[k] = float.Parse( values[1] ); stdDevReturn[k] = float.Parse( values[2] ); }   // Now get the correlation matrix. correlation = new FloatMatrix( n, n ); int cnt = 0; for ( int i = cnt; i < n; i++ ) { for ( int j = cnt; j < n; j++ ) { values = file.ReadLine().Split( ' ' ); correlation[i, j] = float.Parse( values[3] ); correlation[j, i] = float.Parse( values[3] ); }   cnt++; } } return n; }

## Chebyshev Filters in C#

Wednesday, December 11th, 2013

There are three classes of widely used IIR (recursive) filters in signal processing: Butterworth, Chebyshev, and elliptical. In this article I will give a short introduction to the Chebyshev filter, present a code implementation, and end with a usage example. The Butterworth filter was discussed in a previous blog article.

Chebyshev filters come in two flavors defined by either allowing ripple in the pass-band (type 1) or ripple the stop-band (type 2). We will only discuss the type 1 filters here, as the type 2 filters are rarely used. Type 1 Chebyshev filters trade-off steeper roll-off with ripple in the pass band and in the limit of zero ripple, mimic Butterworth filters. Simply put, the more ripple allowed in the passband, the steep roll-off the filter achieves.
Chebyshev filters derive their name from Chebyshev polynomials [1]. Finding Chebyshev polynomials for a given order represent a computational challenge, and when computed often exhibit various numerical instabilities (often overflow) if not done properly. For an excellent review and comparison of computational techniques for computing Chebyshev polynomials see this paper by Wolfram Koepf. Here, we will side step the entire computational issue, as legions of engineers have done before us, and simply look up the Chebyshev polynomials, by order, in a switch statement. Because Chebyshev filters with an order n above 10 or 12 would be unusual, this table look-up approach can be considered, practically speaking, complete.

Wikipedia has a excellent overview article discussing the Chebyshev filter, so I’ll refer you there for a background treatment of Chebyshev filtering theory.

## Chebyshev Filters

The gain expression for Chebyshev filters has a very similar structure to Butterworth filters. As with the Butterworth filter, there are the three familiar design parameters, DC Gain, cutoff frequency f0, and filter order n, with an additional parameter, epsilon, for specifying the desired pass-band ripple.

Chebyshev gain equation

The Tn signifies a Chebyshev polynomial of order n. So for example when n = 2, it happens that Tn = 2*x*x - 1, where when substituted into the gain equation we replace x with f/f0. The ripple factor can be equated to pass-band ripple in Db with the equation:

Passband ripple in Db

Typical values for epsilon would range between 0.1 and 1.

## The Chebyshev Filter in Code

We take the identical approach to implementing the Chebyshev filter in code as we did with the Butterworth filter. We will first compute the input signal’s FFT, then multiply that by the above filter gain, and then take the inverse FFT of that product resulting in our filtered signal. This is a O( n*log(n)) operation.

But before we do that, we need a way to find the Chebyshev polynomials of a given order n. To that end I’ve build a simple method that returns a Chebyshev polynomial given the order; if the requested order is larger than the look-up table, a NotImplementedException is thrown. This look-up approach can easily be extended as far as would be needed for most any practicable Chebyshev filter.

### Chebyshev Polynomial Lookup

 private Polynomial ChebyshevPolynomial( int order) { switch( order ) { case 0: return new Polynomial( new DoubleVector( "1" ) ); case 1: return new Polynomial( new DoubleVector( "0 1" ) ); case 2: return new Polynomial( new DoubleVector( "-1 0 2" ) ); case 3: return new Polynomial( new DoubleVector( "0 -3 0 4" ) ); case 4: return new Polynomial( new DoubleVector( "1 0 -8 0 8" ) ); case 5: return new Polynomial( new DoubleVector( "0 5 0 -20 0 16" ) ); case 6: return new Polynomial( new DoubleVector( "-1 0 18 0 -48 0 32" ) ); case 7: return new Polynomial( new DoubleVector( "0 7 0 56 0 -112 0 64" ) ); case 8: return new Polynomial( new DoubleVector( "1 0 -32 0 160 0 -256 0 128" ) ); default: throw new NotImplementedException(); } }

An excellent on-line resource for Chebyshev polynomials, or any other integer sequence can be found at OEIS, the Online Encyclopaedia of Integer Sequences. If you wish to delve into computing the coefficients of Chebyshev polynomials, a good review of techniques can be found here.

### Chebyshev Filter

The following code implements the Chebyshev filter, following the three steps mentioned above: [1] FFT(signal), [2] FFT(signal)*G(f), [3] IFFT.

private DoubleComplexVector ChebeshevFilter( DoubleComplexVector signal, double sampleFrequency, int order, double f0, double ripple, double DCGain ) { // Get the Chebyshev polynomial Polynomial T = ChebyshevPolynomial( order );   int N = signal.Length;   // Apply forward FFT in place DoubleComplexForward1DFFT fft = new DoubleComplexForward1DFFT( N ); DoubleComplexVector signalFFT = fft.FFT( signal );   // Remove DC offset signalFFT[0] = 0;   if ( f0 > 0 ) {   double numBins = N / 2; // Half the length of the FFT by symmetry double binWidth = sampleFrequency / numBins; // Hz   // Filter System.Threading.Tasks.Parallel.For( 1, N / 2, i => { double binFreq = binWidth * i; double gain = DCGain / ( Math.Sqrt( ( 1 + Math.Pow( T.Evaluate( binFreq / f0) , 2.0 ) ) ) ); signalFFT[i] *= gain; signalFFT[N - i] *= gain; } );   }   // Reverse filtered signal DoubleComplexBackward1DFFT ifft = new DoubleComplexBackward1DFFT( N ); ifft.SetScaleFactorByLength(); // Needed to get the correct amplitude signal = ifft.FFT( signalFFT );   return signal; }

## Filtering Example

To test the Chebyshev filtering code above, I’ve generated a test signal which is the sum to two sinusoids of frequences 1 Hz and 8 Hz. For comparison, this is the same test signal I used to test the Butterworth filter in my previous post. The test signal is sampled at 25 Hz and then filtered using a third order Chebyshev filter with a frequency cut-off of 1.5 Hz and ripple factor (epsilon) of 0.2.

public void ChebyshevFilterTest() { DoubleComplexVector signal = new DoubleComplexVector( 100 );   double Fs = 25; // Sample rate, 25Hz double f1 = 1; // 1 Hz Signal double f2 = 8; // 8 Hz Signal for ( int t = 0; t < 100; t++ ) { signal[t] = new DoubleComplex( Math.Sin( 2 * Math.PI * f1 / Fs * t ) + Math.Sin( 2 * Math.PI * f2 / Fs * t ), 0.0 ); }   // Filter DoubleComplexVector filteredSignal = this.ChebeshevFilter( signal, 25, 3, 1.5, 0.2, 1.0 );   // Display results Chart chart = NMathChart.ToChart( new DoubleVector( signal.ToRealDataTable() ) ); chart.Titles[0].Text = "Input signal, sampled at 25 Hz"; chart.ChartAreas[0].AxisY.Title = "Input signal magnitude"; chart.ChartAreas[0].AxisX.Minimum = 0; Chart chartfft = NMathChart.ToChart( new DoubleVector( filteredSignal.ToRealDataTable() ) ); chartfft.Titles[0].Text = "Filtered signal"; chartfft.ChartAreas[0].AxisY.Title = "Filtered signal magnitude"; chartfft.ChartAreas[0].AxisX.Minimum = 0;   List<Chart> charts = new List<Chart>() { chart, chartfft };   Chart compositeChart = NMathChart.Compose( charts, 1, 2, NMathChart.AreaLayoutOrder.ColumnMajor ); compositeChart.Width = 650; compositeChart.Height = 300;   NMathChart.Show( compositeChart );   }

This example test code, using NMath charts, generates the following chart of the unfiltered and filtered signal side by side.

-Happy Computing,

Paul

[1] Chebyshev is a transliteration of a Russian name and in the past has been variously spelled as Chebychev, Tschebyscheff, Tchebysheff and Tchebichef.

## NMath Tutorial videos

Thursday, October 31st, 2013

We are proud to announce a series of NMath tutorial videos on how to use CenterSpace’s math library products. We are starting, naturally, with Getting Started with NMath.

Please let us know which topics you want us to cover. Email support@centerspace.net

Cheers,
Trevor

## Butterworth Filters in C#

Wednesday, October 30th, 2013

There are three classes of widely used IIR (recursive) filters in signal processing: Butterworth, Chebyshev, and elliptical. In this article I will discuss the Butterworth filter and provide example code implementing and using the filter. The Chebyshev and elliptical filters will be discussed in follow up articles.

Butterworth filters are desirable for their ease of implementation, good phase response, and their smooth monotonic frequency response in both the pass-band and the stop-band. However they don’t achieve the steep roll-off of Chebyshev filters, but Chebyshev filters achieve this steeper roll-off by allowing either some ripple in the pass-band (type 1) or the stop-band (type 2). The type 2 Chebyshev filters are rarely used so only type 1 Chebyshev filters will be discussed in follow up articles. The elliptical filters strike a balance between the type 1 and type 2 Chebyshev filters by allowing ripple in both the stop and pass-bands while achieving the steepest roll-off gain of all three filter types. Wikipedia has excellent overview articles discussing the Butterworth, the Chebyshev, and the elliptical filters, so I’ll refer you to those articles for a thorough background treatment of each filter type instead of repeating that information here.

## Butterworth Filters

Butterworth filters exhibited a ripple free frequency response with a -20*n Db/decade roll-off at the cutoff frequency, where n is the order of the filter. There are only three design parameters for a Butterworth filter, the order n, the cut-off frequency , and the DC gain, , or the gain at zero frequency. So the gain of any Butterworth filter can be written in terms of these three parameters.

Butterworth filter gain

Engineers typically think in terms of Hertz, so we can relate the linear frequency cutoff we are accustom to, to the angular frequency we have above, with the familiar equation, . Now rewriting the above equation in terms of frequency we have:

Butterworth filter gain in terms of f

## The Butterworth Filter in C#

To digitally implement a Butterworth filter in code, we work in the frequency domain (for performance reasons), following the three steps outlined in the pseudo code below.

function ButterworthFilter { fftSignal = FFT( signal );   // Scale FFT by butterworth gain for( int w = 0; w < N/2; w++) filteredSignal[w] = fftSignal[w] * butterworthGain[w];   return inverseFFT( filteredSignal ); }

Given the dual relationship between convolution and multiplication in the time and frequency domains, we could equally well implement this filter by convolving in the time domain the inverse fourier of the Butterworth’s frequency response with the signal, but this would require O(n*n) time. By using the FFT, and working in the frequency domain, we achieve the same effect in O(n* log(n)) time, where n is the length of the signal to be filtered.

DoubleComplexVector ButterworthFilter( DoubleComplexVector signal, double sampleFrequency, int order, double f0, double DCGain ) {   var N = signal.Length;   // Apply forward FFT var fft = new DoubleComplexForward1DFFT( N ); var signalFFT = fft.FFT( signal );   if ( f0 > 0 ) {   var numBins = N / 2; // Half the length of the FFT by symmetry var binWidth = sampleFrequency / numBins; // Hz   // Filter System.Threading.Tasks.Parallel.For( 1, N / 2, i => { var binFreq = binWidth * i; var gain = DCGain / ( Math.Sqrt( ( 1 + Math.Pow( binFreq / f0, 2.0 * order ) ) ) ); signalFFT[i] *= gain; signalFFT[N - i] *= gain; } );   }   // Reverse filtered signal var ifft = new DoubleComplexBackward1DFFT( N ); ifft.SetScaleFactorByLength(); // Needed to get the correct amplitude signal = ifft.FFT( signalFFT );   return signal;

This general Butterworth filtering code begs a bit of explanation. The FFT provides the Fourier transform for N/2 positive frequencies, and symmetrically (via a complex conjugation) for N/2 negative frequencies. Because of this we can loop just N/2 times scaling the symmetric halves simultaneously by our Butterworth gain in our Parallel.For loop for improved performance. When scaling the FFT by the gain of the Butterworth filter, we must know the angular frequency of each bin we scale, which is captured by the binFreq variable; recalling that the frequency resolution of a length N FFT is ( Fs / N / 2 ) Hz, where Fs is the sampling frequency.

Although this code will execute very efficiently it can be significantly optimized if it is specialized to filtering only signals of a fixed length and fixed sampling frequency – a very common situation in many systems. In this case the Butterworth gains can all be computed ahead, avoiding the expensive Math.Sqrt() calls, also the NMath FFT classes would only need to be created once outside the filtering routine. Finally we could do all of the computation in-place and eliminate the copy to a new result vector, with:

 ... // Apply forward FFT in place var fft = new DoubleComplexForward1DFFT( N ); fft.FFTInPlace( signal ); ...

## Filtering Example

To test our new Butterworth filter, I generated a test signal which is simply the sum of two sinusoids of frequencies 1 Hz and 8 Hz. This continuous signal is then sampled at 25 Hz and filtered with a third order Butterworth filter using a cut-off frequency of 1.5 Hz and a DC gain of one. Finally, a chart is generated of the input signal and of the filtered signal to visualize the action of the filter in the time domain.

 public void FilterTest() { var signal = new DoubleComplexVector( 100 );   var Fs = 25; // Sample rate, 25Hz var f1 = 1; // 1 Hz Signal var f2 = 8; // 8 Hz Signal for ( int t = 0; t < 100; t++ ) { signal[t] = new DoubleComplex( Math.Sin( 2 * Math.PI * f1 / Fs * t ) + Math.Sin( 2 * Math.PI * f2 / Fs * t ), 0.0 ); }   // Filter var filteredSignal = ButterworthFilter( signal, 25, 3, 1.5, 1.0 );   // Display results var chart = NMathChart.ToChart( new DoubleVector( signal.ToRealDataTable() ) ); chart.Titles[0].Text = "Input signal, sampled at 25 Hz"; chart.ChartAreas[0].AxisY.Title = "Input signal magnitude"; var chartfft = NMathChart.ToChart( new DoubleVector( filteredSignal.ToRealDataTable() ) ); chartfft.Titles[0].Text = "Filtered signal"; chartfft.ChartAreas[0].AxisY.Title = "Filtered signal magnitude";   List<Chart> charts = new List<Chart>() { chart, chartfft };   var compositeChart = NMathChart.Compose( charts, 1, 2, NMathChart.AreaLayoutOrder.ColumnMajor ); compositeChart.Width = 650; compositeChart.Height = 300; NMathChart.Show( compositeChart );   }

This example code generates the following composite chart showing on the left the unfiltered signal, and on the right the filtered signal.

Butterworth filtering example

As expected, almost all of the high frequency sinusoid at 8Hz has been eliminated from the original signal. In this case where the ‘noisy’ 8 Hz signal is well separated from our 1Hz signal, the Butterworth filter does an excellent job. The next signal processing blog entry will cover the Chebyshev type 1 filter.

-Happy Computing,

Paul Shirkey

## Smoothing Cubic Splines

Tuesday, October 15th, 2013

Cubic smoothing splines embody a curve fitting technique which blends the ideas of cubic splines and curvature minimization to create an effective data modeling tool for noisy data. Traditional cubic splines represent the tabulated data as a piece-wise continuous curve which passes through each value in the data table. The curve spanning each data interval is represented by a cubic polynomial, with the requirement that the endpoints of adjacent cubic polynomials match in location and in their first and second derivatives. If we need to fit a spline to some non-uniformly sampled noisy data, the results can be rather unsatisfying because every data point is visited by the spline (see chart below).

Cubic spline fit to noisy data

If a smoothed cubic spline is fitted to the same noisy dataset, the general sense of the data is captured without the spline mimicking the noise.

Smoothed cubic spline with p = 0.98

## Understanding Smoothing Cubic Splines

Both cubic splines and cubic smoothing splines are constrained to have contiguous first and second derivatives, however only (interpolating) cubic splines are required to pass though all data points. Freed of these N data-point constraints, new additional constraints are now needed to have a fully determined solution. In the case of cubic smoothing splines we have the following two functional constraints: one which minimizes the squared error between the data and spline, and a second which minimizes the curvature.

To minimize the squared error we minimize the following function, where the yi‘s are the data’s ordinates, and the g(x)‘s are cubic polynomials.

Squared error

To minimize the curvature we minimize:

Curvature

Since we are not concerned with the convexity or concavity of the curvature, only the amount of curvature, we square the curvature eliminating the sign.

For most data sets, these two constraint functions act in opposition to one another. The first constraint pushes the spline as close as possible to the data points while the second attempts to keep the spline as free of curvature as possible. In the final formulation, these two countervailing constraints are joining by a weighting parameter, p, which balances these two opposing constraints and results in a smoothing cubic spline.

Cubic Smoothing Spline Constraints

The two degenerate cases are worth pondering. If p = 1 the curvature constraint is nullified and the smoothing cubic spline then passes through all data points resulting in an interpolation spline. In other words, with p = 1 we have reverted back to a natural cubic spline solution. If p = 0 only the curvature of the spline is minimized. Therefore we are minimizing the mean squared error over linear functions resulting in a linear (straight-line) least squares fit. A few smoothing cubic splines are shown below demonstrating the effect of various choices of p.

## Smoothing Cubic Splines C# code example.

The NMath library currently supports two spline classes for natural and clamped splines. The class SmoothCubicSpline is not currently released but will be available in the spring 2014 release of NMath 6.0. If you have a need for this class before our next release, please contact us.

NMath Class Name
Description
Natural cubic spline
Clamped cubic spline
SmoothCubicSpline
Smoothing cubic spline

The following example code creates a noisy sinusoid signal then generates a smoothing cubic spline fit and finally charts the results.

 // Build a noisy signal to fit with the smooth cubic spline. var numPoints = 64; x = new DoubleVector( numPoints, 0, 2.0 * Math.PI / ( numPoints - 1.0 ) ); yNoisy = new DoubleVector( numPoints ); RandGenUniform rng = new RandGenUniform( 0x162 ); for ( int i = 0; i < x.Length; i++ ) { yNoisy[i] = ( 2.3 * Math.Cos( 3.0 * x[i] ) + 1.2 * Math.Sin( 4.5 * x[i] ) + Math.Cos( 1.92 * x[i] ) ) + rng.Next(); }   // Compute the cubic smoothing spline. double p = 0.75; SmoothCubicSpline spline = new SmoothCubicSpline( x, yNoisy, p );   // Chart the noisy signal and smoothing spline. var numInterpolatedValues = 512; var chart = NMathChart.ToChart( spline, numInterpolatedValues ); chart.Titles.Add( String.Format( "Smoothed Cubic Spline with P = {0}", p.ToString() ) ); chart.Series[0].Name = "Noisy Data"; chart.Series[1].Name = "Spline";   NMathChart.Show( chart );

The above code was used to generate the following four charts which compare cubic spline fits of the same dataset with different weighting factors. The four weighting factors used were p = 0.05, 0.5, 0.75, 0.98.

It’s worth noting that if you need to filter uniformly sampled data Savitzky-Golay filters work similarly to a cubic smoothing spline by fitting piecewise continuous polynomials to the data. To more learn about Savitzky-Golay filtering have a look at our blog article.

-Happy Computing,

Paul Shirkey

References
P. Craven and G. Wahba, “Smoothing Noisy Data with Spline Functions”, Numerische Mathematik, vol. 31, pp. 377-403, 1979.