Improve computational performance for the Whittaker algorithm

[InnocenteSimone]

Hi,

I found out that the current implementation for the Whittaker algorithm is quite heavy and it takes time to compute for a good amount of spectra.
For example for 374 spectra, the current implementation took around 74 seconds to complete the preprocessing.

I found an alternative implementation https://github.com/mhvwerts/whittaker-eilers-smoother/blob/master/whittaker_smooth.py under the CeCILL-B license thus you need to cite the author.

The code will be something like this:

def _whittaker(intensity_data, spectral_axis, lam, d):
    m = len(intensity_data)
    E = sparse.eye(m, format='csc')
    D = _speyediff(m, d, format='csc')
    coefmat = E + lam * D.conj().T.dot(D)
    z = splu(coefmat).solve(intensity_data)
    return z, spectral_axis

This solution run 374 preprocessing in 0.47 seconds.

Hope it can helps, thanks!

[InnocenteSimone]

Hi, I just realized that there are missing part in the code that I provided, here the updated:

import scipy.sparse as sparse
from scipy.sparse.linalg import splu

def _speyediff(N, d, format='csc'):
    """
    (utility function)
    Construct a d-th order sparse difference matrix based on
    an initial N x N identity matrix

    Final matrix (N-d) x N
    """

    assert not (d < 0), "d must be non negative"
    shape = (N - d, N)
    diagonals = np.zeros(2 * d + 1)
    diagonals[d] = 1.
    for i in range(d):
        diff = diagonals[:-1] - diagonals[1:]
        diagonals = diff
    offsets = np.arange(d + 1)
    spmat = sparse.diags(diagonals, offsets, shape, format=format)
    return spmat

def _whittaker(intensity_data, spectral_axis, lam, d):
    m = len(intensity_data)
    E = sparse.eye(m, format='csc')
    D = _speyediff(m, d, format='csc')
    coefmat = E + lam * D.conj().T.dot(D)
    z = splu(coefmat).solve(intensity_data)
    return z, spectral_axis