Matlab 2018a: svd implementation changed and now crashes

[ftadel]

The following functions now crash with some specific noise covariance matrices:
bst_inverse_linear_2016, line 915
bst_inverse_linear_2018, line 1093

Matlab 2017b:

K>> [Un,Sn2] = svd(Cov,'econ')
% Returns valid results 

Matlab 2018a:

K>> [Un,Sn2] = svd(Cov,'econ')
Error using svd
SVD did not converge.

Example file:
Cov.zip

The only documentation I could find about changes in this functions is the following, but it is about input parameters, it doesn't say anything about the algorithms changing.
https://fr.mathworks.com/help/matlab/release-notes.html?rntext=svd&startrelease=R2017b&endrelease=R2018a&groupby=release&sortby=descending&searchHighlight=

[jcmosher]

The link you show is for svds, "sparse".

Matlab 2016b is also a valid result. The covariance is effectively singular, with the 64th deep in the bit noise. But it is not precisely "zero" (at least in 2016b), which I mention because the SVD function can be crashed, rarely, if several singular values are perfectly zero, which only happens in theoretical cases. I have not seen the SVD function fail to converge in experimental data.

How was this Covariance matrix generated?

I'll update my license to 2018a and see if I get the same problem.

[jcmosher]

Very interesting discussion here:

https://www.mathworks.com/matlabcentral/answers/285861-different-svd-results-with-r2015b-and-r2016a

version -blas

may give the different libraries used across versions and architectures.

The rank of Cov, in this example, is actually 48 out of 64. Singular values 49-63 are about 10^-25, 64 is 10^-28, relative to the first singular value at 10^-11. Relative Epsilon is 10^-16, so eigenvalues 49-63 are on the hairy edge of being in the bit noise. If the libraries changed so that they are now considered to be zero, then the SVD is challenged with finding eigenvectors in a null space of 16 zeros. Hence the "failure to converge", which I have seen happen before.

If you have time before I install 2018a, try Sn = svd(Cov) vs [Un,Sn] = svd(Cov,0) vs [Un,Sn] = svd(Cov,'econ'). I don't know if "econ" and 0 trigger different conversion routines, but the first case should still work anyway, since it only finds the singular values, and does not iterate to find the eigenvectors, which is where I believe the SVD function fails to converge.

[ftadel]

The data comes from the dataset attached to an article describing an EEGLAB+Brainstorm pipeline to process simple EEG experiments (soon to be published in the same Frontiers research topic as the new brainstorm article). The noise covariance matrix is estimated from the pre-stim baseline of the individual trials, as indicated in our guidelines: http://neuroimage.usc.edu/brainstorm/Tutorials/NoiseCovariance#Variations_on_how_to_estimate_sample_noise_covariance.
The EEG data was pre-processed in EEGLAB with ICA (at least 3 components related with subject artifacts were removed)

Here is the dataset, ready to load in Brainstorm with File > Load protocol > Load from zip file:
https://www.dropbox.com/s/d3agycnq49lljiq/TestSvd.zip?dl=0

Wee need to add ways to handle correctly these cases. Right now it just crashes abruptly with red errors in the command window. If you find a way to adjust the code (or at least to return a readable error that tells the user what should be fixed in the data).

Output in Matlab 2018a of the commands you suggested:

• Sn = svd(Cov): Works
• [Un,Sn] = svd(Cov,0): Crashes
• [Un,Sn] = svd(Cov,'econ'): Crashes

[jcmosher]

More background. Occasionally happens in Python.

http://thread.gmane.org/gmane.comp.python.numeric.general/45614

where their solution was to up the number of iterations. But in general, it has been difficult to get this convergence error. I'll keep digging.

[jcmosher]

Confirmed problem here, "SVD did not converge", which shouldn't happen.

PC: Windows 7 Service Pack 1 - X64, Intel i5-4590

Matlab 2018a on a PC:
version -blas
'Intel(R) Math Kernel Library Version 2017.0.31 Product Build 20170606 for Intel(R) 64 architecture applications, CNR branch AVX2'

Matlab 2016a on same PC:
Intel(R) Math Kernel Library Version 11.2.3 Product Build 20150413 for Intel(R) 64 architecture applications, CNR branch AVX2

Last 4 singular values (raw):

2016a
1.0e-25 *
0.402933770771953 0.361475587020278 0.344527487341072 0.000519328705814

2018a
1.0e-25 *
0.402523446650878 0.360721231051308 0.344213997389103 0.000471957411196

Scaled to first singular value:

2016a
1.0e-13 *
0.120326473509523 0.107945984677843 0.102884842586973 0.000155085019663

2018a
1.0e-13 *
0.120203940085754 0.107720714422279 0.102791226357366 0.000140938722578

and eps reports out the same as 2.220446049250313e-16

So there are subtle differences in the noise floor below eps, but that should not be the issue with the "SVD not converge here".

Testing by adding and eps to the whole matrix
S = svd(Cov); % note, converges, since singular values only

[Un,Sn] = svd(Cov); % does not converge
[Un,Sn,Vn] = svd(Cov + eps(S(1))*eye(64)); % converges

So extreme sensitivity here, have we somehow created a degenerate case for our Cov? Add a tiny amount (eps(S(1)) = 4.038967834731580e-28) to each diagonal. At least "try catch" works:

FAIL = 0;
for i = 1:64,
C = Cov;C(i,i) = Cov(i,i)+1*eps(S(1));
try
[Un,Sn]=svd(C);
catch
FAIL = FAIL+1;
end,
end,
fprintf('Failed %.0f out of 64\n',FAIL)

yields "Failed 31 out of 64"

So we have a remarkably sensitive and degenerate matrix that should nonetheless not be crashing SVD, which has been a robust algorithm for decades. I'm guessing something in the library changed. Definitely worth reporting to Mathworks ASAP.

[jcmosher]

Same problem on iMac (iOS 10.13.3 High Sierra, Intel Core i7), but it is running the same library
version -blas
'Intel(R) Math Kernel Library Version 2017.0.31 Product Build 20170606 for Intel(R) 64 architecture applications, CNR branch AVX2
'

[ftadel]

I posted a bug report on the mathworks website.

[ftadel]

Response from the Mathworks:

Thank you for sending us in this example. We are able to reproduce the issue that you are seeing. Between R2017b and R2018a, we changed which LAPACK function is called from SVD for some cases. We are filing a bug report with LAPACK and/or MKL about this, and reexamining which functions we should be calling from SVD.

[ftadel]

Additional comment from the technical support:

Unfortunately, there is no easy way to get back the R2017b behavior in R2018a. I would suggest using a small perturbation of the matrix, as mentioned in the github thread

S = svd(Cov);
[Un,Sn,Vn] = svd(Cov + eps(S(1))*eye(64));

[jcmosher]

Except that complicates the detection of rank, since we just moved all singular values up off the noise floor. However, our detection of effective rank is supposed to be more robust than relying on bit noise. And we calculate the SVD twice in this fix. We could also put “try” “catch” around these, or calculate the singular vectors separately from the singular values. Actually, that may be the way to substitute an SVD function if the user is using 2018a. Sketch: If 2018a and nargout > 1, S = svd(A); [U,~,V] = svd(A + eps(S(1))*eye(length(S))); return U, diag(S), V Expand to handle cases of vargin options (e.g. economy), nargout of 2 vs 3 for efficiencies. But can we really put in our own SVD, without taking a performance hit? We call the SVD 10s of thousands of times in a single source estimation, relying on calculations of subspaces from matrices that are often singular. For example, the MEG spherical head model is dimension three rank two, will this sometimes not converge, in our calculation of 15,000 cortical points? If not 2018a, we would want the other versions of Matlab to not have to do these additional calls sketched above. Is there a way to "pass through" the SVD arguments with very little performance hit, if not 2018a? Or, if the user is running 2018a, require them to put into their private path an alternative SVD calculation, such as above, so that only these users take the performance hit?

[ftadel]

Yes, it's easy to check for the version of Matlab.
This is already done in many places in the code, because Matlab changes quite a lot.
You can test for Matlab 2018a with:
if (bst_get('MatlabVersion') == 904)

[jcmosher]

Agreed that it's easy to test the version of Matlab. But we call the SVD in many places, so I don't see rewriting the existing code everytime we call SVD, but rather writing a local version of SVD along the lines sketched above. But how do we do that without taking a hit on all of the versions of Matlab?

In other words, I'm unsure how simple it is to write a code that just passes through the argout and argin with only a minor performance hit, except in the case that it is 2018a. Is that straightforward?

[ftadel]

We could replace the svd() calls with a bst_svd() that redirects to the appropriate function.
I think the computational cost of this would be negligible.

[jcmosher]

Would "bst_svd" be preferred over shadowing svd with a custom version (but I see Matlab generates warning messages when you try this with SVD)? I'll leave that to your expert judgement. Presumably they will fix this bad problem in svd and future versions of Matlab won't have it. So we're should only be looking at 2018a as the problem.

Francois, I can flesh out the above function, or do you have enough to go on to write one? I'm not facile at passing vargin and vargout, but SVD doesn't have a lot of options anyway. It is important to distinguish 'econ' vs 0 as an option, and to distinguish two vs three outputs, for efficiency, since we sometimes call svd with very wide matrices, for which we don't want the third output.

But in general I just want to call SVD. Is it an option to tell users to simply stay away from 2018a until this bug is fixed?

[ftadel]

Yes, it's definitely better to have a separate function bst_svd.
Shadowing Matlab functions is something extremely annoying, it makes the scripts difficult to read, the execution changes depending on the folder, etc.

If you know what to do, please write this function. I'm not so familiar with these concepts.

To test for the number of output arguments, use for instance:

if (nargout == 1)
   ...
else if (nargout == 2)
   ...
else if (nargout == 3)
   ...
end

To pass all the input arguments to a function, see for instance bst_call.m:

function varargout = bst_call( varargin )
    if (nargout)
        [varargout{1:nargout}] = feval(varargin{:});
    else
        feval(varargin{:});
    end
end

Don't worry about efficiency. Adding one or two tests and one extra level of function calls will be almost undetectable, unless you're calling this function millions of times.

[ftadel]

Bug still exists in 2018a, but was fixed in 2018b.
We simply display a warning for when there is a crash with 2018a: maybe it's because of the version of Matlab.