2
votes

I'm using MATLAB to do compute conditional covariance and mean using Gaussian mixture model, which always relates to Schur complement. It is suggested in Wiki that if the matrix C is singular, generalized inverse of C can be used to compute the Schur complement.

In MATLAB, pinv is for this aim. As my matrix is very large (more than 1000 columns) and will result in covariance matrix with size > 1000*1000, it can be much faster to use eig instead of svd to compute the pinv. However, this can lose remarkable precision as it truncates eigenvectors corresponding to small eigenvalues under set threshold.

Another way is to use rmdivide function to compute BC^(-1) as B/C since the inverse of a matrix can be considered as a least square problem. In my problem, this can obtain much higher precision and run much faster than using B*pinv(C). Further, rmdivide can deal with some singular matrix, therefore, this method is preferable. But in some cases the warning Matrix is singular to working precision can arise, which result in NaNs if rmdivide is used. So, is there a way to identify when this warning will arise, thus I can use pinv instead?

Update

Additional to @Dohyun's answer, what I'm now doing is to check the obtained results, basing on the fact that NaN can be obtained in the results if the matrix is singular.

warning('off','MATLAB:singularMatrix')
x = b/C; % in my codes, vector is obtained, I think matrix can also be checked in this way
if isnan(sum(x))
    x = b*pinv(C);
end
1

1 Answers

3
votes

If your matrix C is 'usually' nonsingular and mrdvide is much faster than pinv, then you can try mrdivide and then catch warning to switch to pinv.

However, in MATLAB, we cannot catch the warning by using try. Fortunately, there is undocumented solution to catch the warning.

The basic idea is turn warning to error for specific warning ID (in your case, MATLAB:singularMatrix) you will have, and then use try and catch.

myWarn = warning('error','MATLAB:singularMatrix'); % turn singular matrix warning to error
try
    u = mrdivide(B,C);
catch
    u = B*pinv(C);
end
warning(myWarn); % return to warning