sorting_analyzer.compute("amplitude_scalings") returned numpy.linalg.LinAlgError: Matrix is singular.

[chiyu1203]

Dear Spikeinterface community, I used sorting_analyzer.compute("amplitude_scalings") on Spikeinterface 0.101.0 (installing from source).
numpy version is 1.26.4 Kilosort 4.0.4 Window11

I first did spike sorting on Kilosort4 and manually curated on phy. Then I created the sorting_analyser with binary_folder format, and sparse =True, and I used the default setting to compute this "ampliteude_scalings". However, it returned this error.

concurrent.futures.process._RemoteTraceback: """ Traceback (most recent call last): File "c:\Users\neuroPC\anaconda3\envs\spike_interface\lib\concurrent\futures\process.py", line 246, in _process_worker r = call_item.fn(*call_item.args, call_item.kwargs) File "c:\Users\neuroPC\anaconda3\envs\spike_interface\lib\concurrent\futures\process.py", line 205, in _process_chunk return [fn(args) for args in chunk] File "c:\Users\neuroPC\anaconda3\envs\spike_interface\lib\concurrent\futures\process.py", line 205, in return [fn(args) for args in chunk] File "C:\Users\neuroPC\Documents\GitHub\spikeinterface\src\spikeinterface\core\job_tools.py", line 466, in function_wrapper return _func(segment_index, start_frame, end_frame, _worker_ctx) File "C:\Users\neuroPC\Documents\GitHub\spikeinterface\src\spikeinterface\core\node_pipeline.py", line 556, in _compute_peak_pipeline_chunk node_output = node.compute(traces_chunk, node_input_args) File "C:\Users\neuroPC\Documents\GitHub\spikeinterface\src\spikeinterface\postprocessing\amplitude_scalings.py", line 349, in compute scaled_amps = fit_collision( File "C:\Users\neuroPC\Documents\GitHub\spikeinterface\src\spikeinterface\postprocessing\amplitude_scalings.py", line 530, in fit_collision reg = LinearRegression(fit_intercept=True, positive=True).fit(X, y) File "c:\Users\neuroPC\anaconda3\envs\spike_interface\lib\site-packages\sklearn\base.py", line 1474, in wrapper return fit_method(estimator, args, kwargs) File "c:\Users\neuroPC\anaconda3\envs\spike_interface\lib\site-packages\sklearn\linearmodel_base.py", line 611, in fit self.coef = optimize.nnls(X, y)[0] File "c:\Users\neuroPC\anaconda3\envs\spike_interface\lib\site-packages\scipy\optimize_nnls.py", line 91, in nnls x, rnorm, mode = _nnls(A, b, maxiter, tol=atol) File "c:\Users\neuroPC\anaconda3\envs\spike_interface\lib\site-packages\scipy\optimize_nnls.py", line 135, in _nnls s[P] = solve(AtA[P_ind[:, None], P_ind[None, :]], Atb[P], File "c:\Users\neuroPC\anaconda3\envs\spike_interface\lib\site-packages\scipy\linalg_basic.py", line 247, in solve _solve_check(n, info) File "c:\Users\neuroPC\anaconda3\envs\spike_interface\lib\site-packages\scipy\linalg_basic.py", line 44, in _solve_check raise LinAlgError('Matrix is singular.') numpy.linalg.LinAlgError: Matrix is singular. """

The above exception was the direct cause of the following exception:

Traceback (most recent call last): File "", line 1, in File "C:\Users\neuroPC\Documents\GitHub\spikeinterface\src\spikeinterface\core\sortinganalyzer.py", line 884, in compute self.compute_several_extensions(extensions=extensions, save=save, **job_kwargs) File "C:\Users\neuroPC\Documents\GitHub\spikeinterface\src\spikeinterface\core\sortinganalyzer.py", line 1027, in compute_several_extensions results = run_node_pipeline( File "C:\Users\neuroPC\Documents\GitHub\spikeinterface\src\spikeinterface\core\node_pipeline.py", line 505, in run_node_pipeline processor.run() File "C:\Users\neuroPC\Documents\GitHub\spikeinterface\src\spikeinterface\core\job_tools.py", line 428, in run for res in results: File "c:\Users\neuroPC\anaconda3\envs\spike_interface\lib\site-packages\tqdm\std.py", line 1181, in iter for obj in iterable: File "c:\Users\neuroPC\anaconda3\envs\spike_interface\lib\concurrent\futures\process.py", line 575, in _chain_from_iterable_of_lists for element in iterable: File "c:\Users\neuroPC\anaconda3\envs\spike_interface\lib\concurrent\futures_base.py", line 621, in result_iterator yield _result_or_cancel(fs.pop()) File "c:\Users\neuroPC\anaconda3\envs\spike_interface\lib\concurrent\futures_base.py", line 319, in _result_or_cancel return fut.result(timeout) File "c:\Users\neuroPC\anaconda3\envs\spike_interface\lib\concurrent\futures_base.py", line 458, in result return self.get_result() File "c:\Users\neuroPC\anaconda3\envs\spike_interface\lib\concurrent\futures_base.py", line 403, in get_result raise self._exception numpy.linalg.LinAlgError: Matrix is singular.

Could anyone suggest me what might be the issue and how to fix it? There is no much information about calculating this one so I do not know how to fix it myself. All the best,

[alejoe91]

Hi @chiyu1203

I think that this simply means that the LinearRegression to fit a specific collision. I think we should protect against this and either use NaN for the failed spikes or revert to the non-collision option when this happens.

[JoeZiminski]

Hey @chiyu1203 thanks for this bug report. The safest option like Alessio said is to make sure that execution doesn't crash if fitting fails.

It would be nice to get a bit more information on the input variables X, y that are causing the failure to fit. However, this is not necessary so don't worry if you don't have much time to do the below, it would just add a little background information. Are you familiar with debugging or entering breakpoints? This can be useful to get more information on errors. To break into the code you can:

1) pip show spikeinterface just to confirm where spikeinterface is installed. From above, it should be installed at: "C:\Users\neuroPC\Documents\GitHub\spikeinterface\src\spikeinterface". 2) Next open the script at "C:\Users\neuroPC\Documents\GitHub\spikeinterface\src\spikeinterface\postprocessing\amplitude_scalings.py" 3) at line 530 can you replace the line:

reg = LinearRegression(fit_intercept=True, positive=True).fit(X, y) with

try:
    reg = LinearRegression(fit_intercept=True, positive=True).fit(X, y)
except:
    breakpoint()

When you run the script, execution will stop at the crash point and you will be able to interact with the variables. You can print them, but they will probably be very large. However, even printing them may show that something strange is occurring with these input variables (e.g. they have nan or are empty). Alternatively you can save them to disk np.save("/path/to/save/variable_X.npy", X). If you do this for X and y and they are not too large you could upload them to dropbox or google drive to share.

[chiyu1203]

Dear @JoeZiminski and @alejoe91 Thank you for your quick response and clear instruction!! I have followed your step and save X and y when the execution stopped at the crash point. variable_X.zip I have attached them here. Thanks!

[JoeZiminski]

Thanks @chiyu1203! Indeed the matrix X has two columns that are both identical. X^TX is not invertible, which is required for computation of OLS.

The template (y) and traces are shown below. I'm trying to look a bit deeper into the code to understand why this might occur. Wrapping the LinearRegression will definitely handle this but I wonder if it is indicative of any other problems. @alejoe91 does the problem seem to be something that could be expected to occur randomly from time to time and handled at this level? Or might it be an edge case that could be caught higher-up the processing chain?

[JoeZiminski]

@alejoe91 @samuelgarcia I understand a bit more what the code is doing, this must mean that somehow there are two spikes, from the sample template, occurring at exactly the same time (?), that were assigned as a collision in find_collisions? I guess this could happen if a single spike was somehow split in two, and registered to the same unit?

Maybe the easiest thing to do, is when building X, either check collisions contains no duplicates, or if the template to be added to X is not already in X.

[alejoe91]

@JoeZiminski thank you so much for looking into this!

I agree that the problem is that there are to spikes from the same unit exactly at the same time.

@chiyu1203 to fix this, you can try to run this curation step before creating the analyzer:

import spikeinterface.curation as scur

sorting_deduplicated = scur.remove_duplicated_spikes(sorting)

Let me know if it works!

[chiyu1203]

Thank you for your reply and solution! This method works with me. I still do not know how these duplicated spikes generated though. After I ran spike sorting with kilosort, I manually curated the spikes on phy and merged clusters from the same channel (and same template). Maybe there is a bigger issue about how I ran kilosort?!

@JoeZiminski thank you so much for looking into this!

I agree that the problem is that there are to spikes from the same unit exactly at the same time.

@chiyu1203 to fix this, you can try to run this curation step before creating the analyzer:

import spikeinterface.curation as scur

sorting_deduplicated = scur.remove_duplicated_spikes(sorting)

Let me know if it works!

[JoeZiminski]

That's great it's working @chiyu1203! Could you please actually re-open this as it would be useful to perform a check when filling the X matrix that there are no duplicate templates. If there are it will cause strange results in the fitting.

It may also be worth looking into the duplicate spike issue, but my guess is that if there is already a function for it is a sorter-side problem that SpikeInterface has to compensate for with the remove_duplicated_spikes function. So I doubt it is anything to do with how you performing the sorting / curation, @alejoe91 will be able to confirm though.