DyCA - Dynamical Component Analysis

[Datseris]

Dynamical Component Analysis is a new dimensionality-reduction method proposed recently in this paper: https://arxiv.org/abs/1807.10629 by Christian Uhl and coworkers. Other common dimensionality-reduction methods include PCA and DMC, dynamic mode decomposition.

DyCA is particularly suited for identifying how many linear dimensions the multivariate timeseries is composed of. The authors followed up their algorithm with an application to EEG timeseries here: https://arxiv.org/abs/1902.01777

After meeting the authors in a conference they were very kind to share an implementation of their algorithm in Python, that we can use if we want to have a Julia implementation in DynamicalSystems.jl . We thank them very much! I attach the relevant files in this issue.

---

Want to back this issue? Post a bounty on it! We accept bounties via Bountysource.

[Datseris]

One can find the notebook here: https://gist.github.com/Datseris/8121e5019071fe9e3cb9b1a0811ac26a

and the data (for using in testing the implementation) can be found here:
dyca.zip

The core algorithm is described by the following python function:

def DyCA(data, eig_threshold = 0.98):
    "Given data of the form (time, sensors) returns the DyCA projection and projected data with eigenvalue threshold eig_threshold"

    derivative_data = np.gradient(data,axis=0,edge_order=1) #get the derivative of the data
    time_length = data.shape[0] #for time averaging

    #construct the correlation matrices
    C0 = np.matmul(data.transpose(), data) / time_length 
    C1 = np.matmul(derivative_data.transpose(), data) / time_length    
    C2 = np.matmul(derivative_data.transpose(), derivative_data) / time_length   

    #solve generalized eigenproblem
    eigvalues, eigvectors = sp.linalg.eig(np.matmul(np.matmul(C1,np.linalg.inv(C0)),np.transpose(C1)), C2)
    eigvectors = eigvectors[:,np.array(eigvalues > eig_threshold) &  np.array(eigvalues <= 1)] # eigenvalues > 1 are artifacts of singularity of C0

    if eigvectors.shape[1] > 0:
        C3 = np.matmul(np.linalg.inv(C1), C2)
        proj_mat = np.concatenate((eigvectors, np.apply_along_axis(lambda x: np.matmul(C3,x), 0, eigvectors)),axis=1)
    else:
        raise ValueError('No generalized eigenvalue fulfills threshold!')        
    return proj_mat, np.matmul(data,proj_mat)

which is fortunately very compact, making this a very straight forward implementation :)

[anoojpatel]

I'm interested in implementing this in Julia!

[Datseris]

@FuturizeHandgun that is awesome! You can start with a WIP PR and ask for feedback there I think

[kimlaberinto]

@Datseris I would be very interested in implementing this DyCA algorithm if that's alright with you! I wasn't able to find any PRs that linked back to this so I'm assuming it might still be available?

I was thinking of maybe adding a new file tor the algorithm in the folder /src/dimensions , since DyCA pertains to dimensionality reduction. But if you have any thoughts on where, definitely let me know. And if you have any other suggestions or feedback definitely let me know! I'm still learning on how to best get started with doing effective PRs and contribute.

[Datseris]

Hey @kimlaberinto , that sounds like a great plan! No, there isn't any ongoing work on this so feel free to take it! Don't worry too much about where the files are, we can adjust this as a last step.

[Anantha-Rao12]

Dear @Datseris and @kimlaberinto, is this issue still open? I've had a good look at the code and the dataset and, if I'm not wrong, we only need a Julia implementation of the given pythonic DyCA. If you are fine with it, then I can open a PR.

[kimlaberinto]

Sorry! I never got around to it unfortunately. Yes certainly feel free to

[Datseris]

sure

[Anantha-Rao12]

I have implemented a major part of the proposed DyCA but have some queries that I wanted to clarify for this feature request. I will be glad if the concerned developers can advise:

1. I couldn't find any clean implementation of finding the gradient (np.gradient) in Julia. I found this post on Julia discourse but will be glad if someone can suggest how we can use this for a general data matrix? (I can do a straightforward implementation as described here on StackOverflow, but it'll be better if we use the available libraries I feel)
2. This may be a basic question, but anyway, I wanted to clarify what type of error to throw at the end. After some looking up, I found this closed Julia issue that suggests us to use BoundsError. Any help on this is also greatly appreciated.

Thanks.

[Datseris]

Seems that what gradient does is quite simple: https://numpy.org/doc/stable/reference/generated/numpy.gradient.html some central differences only. I'm sure you could utilize some Julia package for it as well, e.g. https://github.com/JuliaDiff/FiniteDifferences.jl. The stack overflow post is also super simple. In general, you should go with what takes less of your time. If you have to spend searching one hour to find some Julia package that does exactly what you need, you probably should just write it yourself instead. The Stack Overflow implementation shouldn't take more than 30 minutes.

This may be a basic question, but anyway, I wanted to clarify what type of error to throw at the end

Throw error in what context? It is probably to discuss this further via an open PR.

[Anantha-Rao12]

Thank you @Datseris. On the way to open a PR.

[Datseris]

I'm re-opening this because there is some problem with DyCA. Its test do not pass (currently commetned out).

[KalelR]

Hey. Apparently, the errors are just due to not importing all the modules. The tests run successfully if I include using DynamicalSystemsBaseand using DelayEmbeddings.

I can open a PR with these additions, and uncommenting the dyca_tests in the runtests.jl.

I could also maybe include further tests, using the EEG seizure data you used in your DyCA notebook.

[Datseris]

Yes please, sounds good to me!