Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Robust Shared Response Model #302

Merged
merged 33 commits into from Dec 12, 2017
Merged

Robust Shared Response Model #302

merged 33 commits into from Dec 12, 2017

Conversation

TuKo
Copy link
Contributor

@TuKo TuKo commented Nov 23, 2017

The code introduces the Robust SRM, tests and an example with synthetic data.

@TuKo
Copy link
Contributor Author

TuKo commented Nov 23, 2017

@CameronTEllis and @hejiaz, please do you mind reviewing this PR?

Thank you and happy thanksgiving!

Copy link
Contributor

@CameronTEllis CameronTEllis left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Overall it is looking good, my comments are relatively minor and mainly just focused on fleshing details out

"""Robust Shared Response Model (RSRM)

The implementation is based on the following publications:

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It would be good to provide a paragraph explaining what RSRM offers and how it compares to SRM, SSRM etc. You could copy and paste one of the paragraphs from manuscript. Also should mention the influence of RPCA

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added an extended description of the method. However, I don't think that the relation to RPCA is relevant to the code/implementation. The idea is that people can understand the code and relate to the paper if there is something unclear.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fair enough


gamma : float, default: 1.0
Regularization parameter for the sparseness of the individual
components.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would be good to add direction to this description (Is high more or less sparse?)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done

The shared response.

s_ : list of array, element i has shape=[voxels_i, samples]
The individual components with outlying data for each subject.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You use outlier here and elsewhere but sometimes refer to this as the individual component. It is important to be consistent so check that the same word is used everywehere

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

----

The number of voxels may be different between subjects. However, the
number of samples for the alignment data must be the same across
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

'samples' is not a typical way of describing timepoints in fMRI data. Consider changing this to say that 'timepoints must be aligned across participants'

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done (everywhere)

X : list of 2D arrays, element i has shape=[voxels_i, samples]
Each element in the list contains the fMRI data of one subject.

y : not used
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is this for then?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed.

Parameters
----------

X : array, shape=[voxels, samples]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe use a different letter, or lower case since this is no longer a list like it is in other functions. same with S below

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

"""
A = X.dot(R.T)
A -= S.dot(R.T)
U, _, V = np.linalg.svd(A, full_matrices=False)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is probably obvious but I am unsure why you are doing the SVD. A comment would help

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

" \n",
"S = [None] * subjects\n",
"for s in range(subjects):\n",
" S[s] = (np.random.rand(data[s].shape[0], samples) < p) * ((np.random.rand(data[s].shape[0], samples) * amplitude - amplitude/2) )\n",
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why use uniform noise? I would have thought gaussian noise would be more reasonable?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

S is not the noise but the sparse individual component. If we use a Normal distribution, S will mix with the shared signal because sums of Normal distributions are also Normal distributions. In that case you will have a difficult time to identify what comes from where (Shared vs individual). This is also a common assumption for the RPCA.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fine

"features = 3\n",
"snr = 20 # in dB\n",
"amplitude = 8.0 # approximate amplitude used for the values in Si\n",
"p = 0.05 # probability of being non-zero in Si"
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How does this relate to the gamma value used later? Why are they not the same?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

p is the probability of an element to be non-zero in the original Si (drawn from the model).
gamma is selected to approximately have the same number of non-zeros in the estimated Si.
However, there is no reason for them to have the same values. Ones is the probability and the other is a balancing value between regularization terms.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I might still be misunderstanding but shouldn't the non-zero elements of the original Si be thought of as the individual component that RSRM is aiming to extract? Why would there be a different number in the input vs the output?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The math is not translated 1:1. Remember that this is probability, and when you draw a signal there are not exactly 5% of non-zeros. There could be anywhere near the 0.05. Added to that this will be differently with all subjects. So every subject will have approximately that, but not exactly.
Therefore, even if you can obtain a 1:1 translation from p (range 0 to 1) to gamma (range 0 to +inf), still the signal itself will not necessarily behave in that way.

Also, this is a particular case of drawing Si's, and you can change it by other distributions as well. In such case, you will have parameters to generate the signals and still gamma to recover it. Dividing between gamma and p, it is a good idea.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Okay fair enough I understand

" n = noise_level * np.random.randn(voxels, samples)\n",
" data[s] += n\n",
" \n",
"S = [None] * subjects\n",
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Might be good to make another example where there is no S variability and show that RSRM does as well as SRM

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The idea of the example is to show how to use it. I chose the synthetic example as it is easy to replicate and visually show the improvements. If you set gamma to "infinity", you will be solving exactly the (deterministic) SRM.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good

@TuKo
Copy link
Contributor Author

TuKo commented Dec 6, 2017

@CameronTEllis do you have any suggestion or comment for the newest version?

@TuKo
Copy link
Contributor Author

TuKo commented Dec 11, 2017

@CameronTEllis do you think there are more changes to be done? If we are done, then don't forget to approve the review. Otherwise let me know.
Thanks!

@mihaic mihaic merged commit f708dfb into brainiak:master Dec 12, 2017
@TuKo TuKo deleted the rsrm branch December 12, 2017 22:41
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants