Parametric spectrum estimation using aryule: difference between Python and Matlab

Question

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

Answer 1

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.

Answer 2

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!