spm_pinv.m 1.4 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253
  1. function X = spm_pinv(A,TOL)
  2. % pseudo-inverse for sparse matrices
  3. % FORMAT X = spm_pinv(A,TOL)
  4. %
  5. % A - matrix
  6. % TOL - Tolerance to force singular value decomposition
  7. % X - generalised inverse
  8. %__________________________________________________________________________
  9. % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
  10. % Karl Friston
  11. % $Id: spm_pinv.m 5877 2014-02-11 20:03:34Z karl $
  12. % check A
  13. %--------------------------------------------------------------------------
  14. [m,n] = size(A);
  15. if isempty(A), X = sparse(n,m); return, end
  16. % try generalised inverse
  17. %--------------------------------------------------------------------------
  18. if nargin < 2
  19. sw = warning('off','MATLAB:nearlySingularMatrix');
  20. warning('off', 'MATLAB:singularMatrix');
  21. X = spm_inv(A'*A);
  22. warning(sw);
  23. % check everything is finite
  24. %----------------------------------------------------------------------
  25. if all(isfinite(X(:)))
  26. X = X*A';
  27. return
  28. end
  29. end
  30. % pseudo-inverse
  31. %--------------------------------------------------------------------------
  32. [U,S,V] = spm_svd(A,0);
  33. S = full(diag(S));
  34. if nargin < 2, TOL = max(m,n)*eps(max(S)); end
  35. % tolerance
  36. %------------------------------------------------------_-------------------
  37. r = sum(abs(S) > TOL);
  38. if ~r
  39. X = sparse(n,m);
  40. else
  41. i = 1:r;
  42. S = sparse(i,i,1./S(i),r,r);
  43. X = V(:,i)*S*U(:,i)';
  44. end