# Parametric spectrum estimation using aryule: difference between Python and Matlab

I have some code illustrating my problem below. As you can see from the results- Matlab is significantly faster than python. It seems that the function which performs the levinson durbin recursion in Matlab is a Mex file. Could this be why Matlab is significantly faster? How would I make the Python implementation on par with Matlab, or does anyone know if such a version already exists?

On a separate issue, the python aryule function does not seem to be able to handle white noise that is not unit variance. For instance if I were to use

X = 0.2*randn(1, 256000)


the python estimates for the AR coefficients would be less accurate. Matlab doesn't seem to have this problem. Is there a way to implement this in Python?

## Matlab Implementation

A=[1 -2.7607 3.8106 -2.6535 0.9238];
% AR(4) coefficients
y=filter(1,A,0.2*randn(256000,1));
%filter a white noise input to create AR(4) process
t = cputime;
[ar_coeffs v]=aryule(y,4);
e = cputime-t;
disp(['Elapsed Time: ' num2str(e)]);
%compare the results in ar_coeffs to the vector A.
ar_coeffs


The result:

Elapsed Time: 0.14


## Python Implementation

from pylab import *
import scipy.signal
from spectrum import *
import time

# Create an AR model
a = [1, -2.7607, 3.8106, -2.6535, 0.9238]
# create some data based on these AR parameters
X = randn(1, 256000)
y = scipy.signal.lfilter([1], a, X)
# now, let us try to estimate the original AR parameters
t0 = time.time()
AR, P, k = aryule(y[0], 4)
print "Time elapsed: ", time.time()-t0
%compare the results in ar_coeffs to the vector A.
print(AR)


The result:

Time elapsed:  2.48965501785

• As to your second question: This is hard to answer without knowing the algorithm and its exact implementation of the algorithm both in Spectrum (which is open source) and Matlab (which isn't). In any case, this is really an independent question -- could you ask this as a separate question? It would help if you could give some brief example codes for this as well, like you did for the timing (+1 for that, by the way). It always helps to have an explicit comparison between what is expected and what is observed. – Christian Clason Aug 13 '14 at 20:15

The reason that Python code is slower than the equivalent Matlab code is usually that in Matlab, more computations are carried out by libraries that are

1. written in a lower-level language, and
2. heavily optimized.

Luckily, you can do the same in Python.

For example, Matlab is linked against (and bundles) Intel's hand-tuned MKL library, while Python (that is, NumPy/SciPy) usually links against the generic system BLAS/LAPACK. With a bit of wrangling (and after paying Intel for a license), you can build NumPy/SciPy yourself so that it calls MKL instead; after that, you get roughly the same performance for anything that relies mostly on BLAS/LAPACK calls. You can do the same with OpenBLAS, which has similar performance as MKL but is open source and easier to build NumPy with. With this setup, I get 0.66s for your Python code on my machine (which runs your Matlab code in 0.09s).

The remaining performance difference is probably due, as you point out, to the fact that Matlab's version is implemented as a MEX file (i.e., written in C and pre-compiled). There are several options to do the same thing in Python, e.g., Cython. As awesomebytes so eloquently suggested, you might want to avoid Spectrum (which is no longer actively developed) and try the correlation function in scipy.signal.correlate if you don't want to get your hands dirty with C. (The LD recursion is implemented in scipy.linalg.solve_toeplitz). You might also want to take a look at statsmodel and nitime.

How to make that process faster? Well, firstly if you feel like having the time you could check the implementation, which I'm doing just right now. (Warning: I need to truncate links as stackexchange does not let me have more than 2 links)

Code is here: assembla . com /code/PySpectrum/subversion/nodes/37/trunk/src/spectrum

Giving a look at the package it already contains some .c code to fasten up some operations. That may ease the implementation in C of some part to later on import it via Python if needed.

I did svn checkout https://subversion.assembla.com/svn/PySpectrum/

And I got: ~/spectrumlib/PySpectrum/trunk/src/cpp\$ ls init.py mydpss.c

mydpss.c implements some subroutines. Something could be added there maybe.

But let's take a look at the code that is actually being executed: assembla . com /code/PySpectrum/subversion/nodes/37/trunk/src/spectrum/yulewalker.py

from correlation import CORRELATION
from levinson import LEVINSON

r = CORRELATION(X, maxlags=order, norm=norm)
A, P, k = LEVINSON(r, allow_singularity=allow_singularity)
return A, P, k


This is what is being executed. Let's check what is the slow part here.

So I modified yulewalker.py like:

import time
print "Correlation calculation"
t0 = time.time()
r = CORRELATION(X, maxlags=order, norm=norm)
print "Time elapsed: ", time.time()-t0
print "Levinson calculation"
t1 = time.time()
A, P, k = LEVINSON(r, allow_singularity=allow_singularity)
print "Time elapsed: ", time.time()-t1
return A, P, k


I ran your piece of code of stackexchange and I got:

./aryulespeed.py
Correlation calculation
Time elapsed:  1.50509214401
Levinson calculation
Time elapsed:  0.000110864639282
Time elapsed:  1.50538802147
[-2.7616739   3.81292317 -2.65578205  0.9248034 ]


So the correlation calculation is the slow part! aha! This guy: assembla . com /code/PySpectrum/subversion/nodes/37/trunk/src/spectrum/correlation.py seems to be the fat guy.

Oh, look, a funny comment in the code:

Provides two correlation functions. :func:CORRELATION is slower than
:func:xcorr. However, the output is as expected by some other functions.
Ultimately, it should be replaced by :func:xcorr.

For real data, the behaviour of the 2 functions is identical. However, for
complex data, xcorr returns a 2-sides correlation.


Maybe we should try to interchange this CORRELATION for xcorr!

So I tried:

#! /usr/bin/env python
from pylab import *
import scipy.signal
from spectrum import *
import time

# The original imports
from spectrum.correlation import CORRELATION
from spectrum.levinson import LEVINSON
# The faster corr
from spectrum.correlation import xcorr

def superawesomearyule(X, order, norm='biased', allow_singularity=True):

assert norm in ['biased', 'unbiased']
import time
print "Correlation calculation in superawesomearyule"
t0 = time.time()
#r = CORRELATION(X, maxlags=order, norm=norm) # Commenting you, slowwww     code!
r = xcorr(X, maxlags=order, norm=norm)
print "Time elapsed: ", time.time()-t0
print "Levinson calculation"
t1 = time.time()
A, P, k = LEVINSON(r, allow_singularity=allow_singularity)
print "Time elapsed: ", time.time()-t1
return A, P, k

# Create an AR model
a = [1, -2.7607, 3.8106, -2.6535, 0.9238]
# create some data based on these AR parameters
X = randn(1, 256000)
y = scipy.signal.lfilter([1], a, X)
# now, let us try to estimate the original AR parameters
t0 = time.time()
#AR, P, k = aryule(y[0], 4)
AR, P, k = superawesomearyule(y[0], 4)
print "Time elapsed: ", time.time()-t0
#%compare the results in ar_coeffs to the vector A.
print(AR)


And executed:

./aryulespeed.py
Correlation calculation in superawesomearyule
Time elapsed:  69.9396679401
Levinson calculation
Traceback (most recent call last):
File "./aryulespeed.py", line 38, in <module>
AR, P, k = superawesomearyule(y[0], 4)
File "./aryulespeed.py", line 24, in superawesomearyule
A, P, k = LEVINSON(r, allow_singularity=allow_singularity)
File "/usr/local/lib/python2.7/dist-packages/spectrum/levinson.py", line     112, in LEVINSON
if P <= 0 and allow_singularity==False:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()


And...... oh noes :( it took 69 seconds!! and also crashed after that on levinson.

I think giving it a bit of love something may be done here. Maybe sending a mail to it's developer Thomas Cokelaer would help.

And I'm sorry but I can't invest more time in this issue right now. I hope this helps you in some way :)

I also found this github repo: https://github.com/RhysU/ar

Which has an ar implementation with some python bindings. It may also help you.

Good luck!

• So TL;DR: Yes, Matlab is faster because the recursion is implemented in C as opposed to pure Python, but you can speed the Python version up by computing the correlation in C as well (e.g., via Cython)? – Christian Clason Aug 13 '14 at 13:44
• @ChristianClason That's what I got out of it. I'd convert your comment to an answer. As much as I appreciate awesomebytes' effort, brevity counts. – Geoff Oxberry Aug 13 '14 at 18:54
• This was a really good answer and very useful for my purposes, but as @GeoffOxberry pointed out, I think other people will find ChristianClason answer to be more useful. – Dipole Aug 13 '14 at 23:31