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:


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


My code, awaiting your feedback:

import numpy as np

def hurst(sig: np.ndarray) -&gt; 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)[0]
    return hurst_exponent
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.