Hurst Exponent calculator

I am aiming to compute the Hurst Exponent of a 1-D signal time series in Python. For now, I have one existing function hurst(sig) which returns the Hurst exponent of sig as a float. Unfortunately, the code runs very slowly even for signals with only ~7500 data points. My main concerns are:

1. Calculate the computational complexity of my algorithm given the size of sig.
2. Reduce computational complexity of the algorithm.
3. Reduce run-time (~1s).
4. Improve readability and address style convention as applies.

My code, awaiting your feedback:

import numpy as np

def hurst(sig: np.ndarray) -> float:
"""Compute the Hurst exponent of sig.

arguments
sig -- 1D signal
returns
hurst_exponent -- float
"""

n = sig.size # num timesteps
t = np.arange(1, n + 1)
y = sig.cumsum() # marginally more efficient than: np.cumsum(sig)
mean_t = y / t # running mean

s_t = np.zeros(n)
r_t = np.zeros(n)

for i in range(n):
s_t[i] = np.std(sig[:i + 1])
x_t = y - t * mean_t[i]
r_t[i] = np.ptp(x_t[:i + 1])

r_s = r_t / s_t
r_s = np.log(r_s)[1:]
n = np.log(t)[1:]
a = np.column_stack((n, np.ones(n.size)))
[hurst_exponent, c] = np.linalg.lstsq(a, r_s)
return hurst_exponent
• Does the signal have any special properties or would it be appropriate to test with random values of a given length? – AlexV Jul 17 at 18:43
• @AlexV you could test with random values of a given length. Currently, the length is anywhere between 7680 and 153600. – Matthew Anderson Jul 17 at 18:45

A few observations:

s_t = np.zeros(n) and r_t = np.zeros(n) are more than you need. Since you don't actually use the array values but solely overwrite them, you can use np.empty here.

You're doing quite a bit of redundant work in the for loop.

When calculating s_t[i], numpy basically has to repeat all the computations it has already done for s_t[i-1]. If you have a look at the definition of np.std, you can see that the mean plays a role here. Since you have already come up with a nice solution for the mean value, you could reuse it here:

s_t = np.empty(n)

for i in range(n):
s_t[i] = np.mean((sig[:i+1] - mean_t[i])**2)
...
s_t = np.sqrt(s_t)

I validated this against the original code with np.allclose.

r_t also has some extra work that you can get rid of. At the moment your code does the following:

x_t = y - t * mean_t[i]
r_t[i] = np.ptp(x_t[:i + 1])

All the values of x_t at [i+1:n] are computed, but never used, so you can do

r_t[i] = np.ptp(y[:i+1] - t[:i+1] * mean_t[i])

instead, without changing the outcome. You could also write this as two list comprehensions:

s_t = np.sqrt(
np.array([np.mean((sig[:i+1] - mean_t[i])**2) for i in range(n)])
)
r_t = np.array([np.ptp(y[:i+1] - t[:i+1] * mean_t[i]) for i in range(n)])

Or even a single one if you'd like to get a little bit creative:

rs_t = np.array([
[
np.mean((sig[:i + 1] - mean_t[i])**2),
np.ptp(y[:i + 1] - t[:i + 1] * mean_t[i])
] for i in range(n)
])
rs_t[:, 0] = np.sqrt(rs_t[:, 0])

A fully vectorized computation of x_t I have come up with, starting from

x_t = y.reshape(1, -1) - mean_t.reshape(-1, 1) @ t.reshape(1, -1)

is considerably slower than the looped version.

I'm actually not quite sure about the rest of the code from this point on at the moment. I will have to take a closer look to give a meaningful review performancewise. As a first non-functional feedback I would like to invite you to think about more meaningful variable names in the future. Especially there at the end of the code around the least-squares optimization, they are quite hard to understand.

I did some timing of different versions of this algorithm. The full benchmark code can be found in this gist. This is what I got so far (for n=15360):

base:              3.748915s
better_loop:       2.605374s
double_lc:         2.430255s
single_lc:         2.445943s
-----------------------
base_numba:        0.866583s
better_loop_numba: 0.813630s
double_lc_numba:   0.776361s
single_lc_numba*:  N/A

The lower group of results uses numba, a just-in-time compiler for Python code, that can sometimes help to speed up loops. Unfortunately, numba does not support all types of array operations in its fast nonpython mode, e.g. I could not get the single list comprehension version to properly work with numba.

The results are actually a little bit worse than I had originally expected, especially without numba, but there is still a good amount profiling work to be done here.

• I wish I could upvote this twice. By the way, the variable names match the symbols of the numbers used in the hurst exponent calculation E[R(n)/S(n)]=Cn^h – Matthew Anderson Jul 17 at 23:00
• If you refer to a mathematical formula, I find it quite helpful to link to the resource you are getting it from, or at least a short "natural language" description in a comment (where necessary). Sometimes authors use notation that slightly differs from "the usual way", though I have to admit I've never heard of the Hurst exponent before this question, so I'm maybe not the best person to judge this. – AlexV Jul 18 at 6:50