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.

        sig -- 1D signal
        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
  • \$\begingroup\$ Does the signal have any special properties or would it be appropriate to test with random values of a given length? \$\endgroup\$
    – AlexV
    Jul 17 '19 at 18:43
  • \$\begingroup\$ @AlexV you could test with random values of a given length. Currently, the length is anywhere between 7680 and 153600. \$\endgroup\$ Jul 17 '19 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.

  • \$\begingroup\$ 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 \$\endgroup\$ Jul 17 '19 at 23:00
  • 2
    \$\begingroup\$ 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. \$\endgroup\$
    – AlexV
    Jul 18 '19 at 6:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.