I was wondering about the complexity of cvxeda.

[Momellouky]

The cvxeda method is included in the neurokit2 python package. I wanted to use it to decompose the phasic and tonic signals from the raw EDA signal. However, as the sampling rate was quite high (700 Hz, and there is no possibility of downsampling), The algorithm takes minutes (I estimate it takes hours in my case, as I have a total of 60 hours of EDA recordings). So I'd like to know the complexity (the big O notation) of this algorithm. This information will help me justify the long execution of the algorithm on my data.

Yours sincerely

[lciti]

Assessing the complexity in big O notation can be difficult and even misleading in this case, for a number of reasons:

• Since the optimisation uses an iterative method (in the case of the python version used in neurokit2, an interior-point method is used), the overall complexity scales with the number of iterations times the complexity of each iteration. However the number of iterations is hard to determine. Empirically, the interior point method tends to converge in a roughly constant number of iterations (in the order of 50) regardless of the problem size but there's no guarantee that this will be true for every problem and every set of data.
• The complexity of each iteration depends on the sparsity pattern of the quadratic form matrix. For a banded matrix, a linear system can be solved in $O(b^2 n)$ where $b$ is the bandwidth and $n$ is the length. In our case $n$ is the length of the signal (length of the phasic component), plus 2 (length of the linear drift component), plus the number of spline knots for the tonic component (one every 10s). Determining $b$ is much harder. The simplest thing to do would be to copy the cvxEDA function and insert print(len(H), band(H)) after the line H = cv.sparse(...) where band is defined somewhere as follows.

  def band(M):
  permutation = cvxopt.amd.order(M, uplo='L')
  inv_permutation = np.empty(len(permutation), dtype=int)
  inv_permutation[np.reshape(permutation, -1)] = np.arange(len(inv_permutation))
  ImJ = inv_permutation[M.I] - inv_permutation[M.J]
  return 1 + np.max(ImJ) + np.max(-ImJ)

I have put this together from bit and pieces that I had around and I have not tried it myself. You could have a go and see if it works.

In your case, the easiest gains would be gained by downsampling. A 700 sps sampling frequency is unnecessarily high for a slow signal like EDA. You say you can't downsample it but provide no reason for that. If at all possible, I would reconsider that.