Technical fixes to roi_extract.

CHANGELOG.md

@@ -9,6 +9,16 @@ neat.
 ### Changed
 - Changed takesnap inf view angle to fix shadow from @judyseinkim.
 - Changed scripts to record version number from @clane9.
+- Changed method for computing beta values in roi_extract. Previously, betas were
+  computed the naive way, by explicitly inverting `X' * X`. Although correct,
+  this is slow and numerically unstable. It's better to use the SVD, as
+  discussed [on wikipedia](https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)#Orthogonal_decomposition_methods)
+  (from @clane9).
+
+### Fixed
+- Redundant spike covariates are now removed before fitting the GLM. This
+  prevents the "Design matrix is rank-deficient" error from coming up in this
+  case (from @clane9).
 
 ## [1.0] - 2016-08-09
 

bin/roi_extract

@@ -632,14 +632,25 @@ def read_run_cfds(run_cfds, run_len):
 
     # Determine type of confound
     # Normalize continuous covariates (in particular, center)
+    # Also, remove any redundant spike covariates
     spike_mask = np.ones(run_cfd.shape[1], dtype=int)
+    keep_mask = np.ones(run_cfd.shape[1], dtype=int)
+    spike_inds = [];
     for i in range(run_cfd.shape[1]):
       if np.sum(run_cfd[:, i] == np.round(run_cfd[:, i])) < run_len:
         run_cfd[:, i] = normalize(run_cfd[:, i])
         spike_mask[i] = 0
       elif np.sum(run_cfd[:, i]) != 1.:
         raise NPDLError(('Only continuous covariates and spike regressors' +
                         'accepted as confounds.'))
+      else:
+        # Check if spike covariate is redundant.
+        spike_ind = np.argmax(run_cfd[:, i])
+        if spike_ind in spike_inds:
+          keep_mask[i] = 0
+        else:
+          spike_inds.append(spike_ind)
+    run_cfd = run_cfd[:, keep_mask==1]
   return run_cfd, spike_mask
 
 class Events(object):

lib/npdl_utils.py

@@ -413,7 +413,7 @@ def prep_fig(leadbuff, backbuff, bottombuff, topbuff, figw, figh):
 
 def fit_glm(TS, DM, cfd_mat=None, intercept=True, outdir=None,
             out_prefix=None, logger=None):
-  """Calculate GLM beta weights. Return baseline and PSC betas.
+  """Calculate GLM beta weights using SVD. Return baseline and PSC betas.
   """
   # prepend intercept
   if intercept:
@@ -426,25 +426,25 @@ def fit_glm(TS, DM, cfd_mat=None, intercept=True, outdir=None,
     np.savetxt('{}/{}_design.txt'.format(outdir, out_prefix), DM, delimiter=' ')
   # calculate singular values
   try:
-    S = np.linalg.svd(DM, compute_uv=False)
+    U, S, VT = np.linalg.svd(DM, full_matrices=False, compute_uv=True)
   except LinAlgError:
     raise NPDLError('Bad design matrix: SVD did not converge.')
   S_str = ' '.join(['{0:.0f}'.format(s) for s in S])
   if logger is not None and out_prefix is not None:
     logger.log('Singular values for {} design mat: {}.'.format(out_prefix, S_str))
   # determine if rank deficient as in matrix_rank function
   eps = np.finfo('f8').eps
-  thresh = S.max() * np.max(DM.shape) * eps
+  thresh = S[0] * np.max(DM.shape) * eps
   rank = np.sum(S > thresh)
   if logger is not None and out_prefix is not None:
     logger.log('Design matrix width for {}: {}; rank: {}.'.format(out_prefix, DM.shape[1], rank))
   if rank < S.size:
     raise NPDLError('Design matrix is rank deficient.')
-  # calculate betas
-  # y = X*beta
-  # beta = (X'X)^-1 * X' * y
+  # calculate betas using numerically stable SVD method.
+  # See: https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)#Orthogonal_decomposition_methods
   TS = TS.reshape((-1, 1))
-  beta = np.dot(np.linalg.inv(np.dot(DM.T, DM)), np.dot(DM.T, TS))
+  Sinv = (S ** -1).reshape(-1, 1)
+  beta = np.dot(VT.T, Sinv * np.dot(U.T, TS))
   if outdir is not None and out_prefix is not None:
     np.savetxt('{}/{}_beta.txt'.format(outdir, out_prefix), beta)
   beta = beta.reshape(-1)