svd crash during HCCM Jacobian inversion

[ritsch-i]

Hi, while fitting a series of data, one dataset somehow seemed to be problematic, and I got the error report that one sub-routine in fitsignal did not converge (see below):

1. What can cause this? Is there anything I can do to 'help' the fit algorithm?
2. Would be nice that even when this error is caught somehow, such that the fit of the series continues anyways with the remaining data in the set (otherwise I have to figure out where it crashed, remove and restart manually).

python errormessage

/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/numpy/linalg/linalg.py in _raise_linalgerror_svd_nonconvergence(err, flag)
     95 
     96 def _raise_linalgerror_svd_nonconvergence(err, flag):
---> 97     raise LinAlgError("SVD did not converge")
     98 
     99 def _raise_linalgerror_lstsq(err, flag):

LinAlgError: SVD did not converge

[luisfabib]

Duplicate of #42

This is related to multiple errors we have already encountered in regparamrange. What kind of OS do you use? MacOS?

This will be solved as soon as #42 is closed. There is not much you can do at the moment except to set a fix regularization parameter manually to avoid entering regparamrange.

[luisfabib]

Can you please provide the full error stack (the full error message) to see where exactly in the call tree this is crashing?

[ritsch-i]

Here you go:

---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
<ipython-input-479-5763123bbc87> in <module>
      6     r[i] = np.linspace(1.5,6,len(trid[i]))# distance axis, ns
      7     # standard DEER kernel
----> 8     ridmefit[i] = dl.fitsignal(Vrid[i],trid[i],r[i],'P',dl.bg_strexp,dl.ex_4pdeer)
      9     lambda0 = ridmefit[i].exparam[0]
     10     # extend to (static) RIDME kernel defined in RIDME exp. pathways

/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/deerlab/fitsignal.py in fitsignal(Vexp, t, r, dd_model, bg_model, ex_model, par0, lb, ub, verbose, weights, uqanalysis, regparam, regtype)
    469 
    470         # Separable non-linear least squares (SNNLS)
--> 471         fit = dl.snlls(Vexp,lambda par: multiPathwayModel(par)[0],par0,lb,ub,lbl, reg=True,
    472                             regparam=regparam, lin_maxiter=[], uqanalysis=uqanalysis, lin_tol=[], weights=weights)
    473         parfit_ = fit.nonlin

/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/deerlab/snlls.py in snlls(y, Amodel, par0, lb, ub, lbl, ubl, nnlsSolver, reg, weights, regtype, regparam, multistart, regorder, alphareopt, nonlin_tol, nonlin_maxiter, lin_tol, lin_maxiter, huberparam, uqanalysis)
    374 
    375         # Calculate the heteroscedasticity consistent covariance matrix
--> 376         covmatrix = hccm(J, res, 'HC1')
    377 
    378         # Get combined parameter sets and boundaries

/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/deerlab/utils/utils.py in hccm(J, *args)
    171 
    172     # Hat matrix
--> 173     H = J@np.linalg.pinv(J.T@J)@J.T
    174     # Get leverage
    175     h = np.diag(H)

<__array_function__ internals> in pinv(*args, **kwargs)

/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/numpy/linalg/linalg.py in pinv(a, rcond, hermitian)
   2001         return wrap(res)
   2002     a = a.conjugate()
-> 2003     u, s, vt = svd(a, full_matrices=False, hermitian=hermitian)
   2004 
   2005     # discard small singular values

<__array_function__ internals> in svd(*args, **kwargs)

/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/numpy/linalg/linalg.py in svd(a, full_matrices, compute_uv, hermitian)
   1659 
   1660         signature = 'D->DdD' if isComplexType(t) else 'd->ddd'
-> 1661         u, s, vh = gufunc(a, signature=signature, extobj=extobj)
   1662         u = u.astype(result_t, copy=False)
   1663         s = s.astype(_realType(result_t), copy=False)

/Library/Frameworks/Python.framework/Versions/3.8/lib/python3.8/site-packages/numpy/linalg/linalg.py in _raise_linalgerror_svd_nonconvergence(err, flag)
     95 
     96 def _raise_linalgerror_svd_nonconvergence(err, flag):
---> 97     raise LinAlgError("SVD did not converge")
     98 
     99 def _raise_linalgerror_lstsq(err, flag):

LinAlgError: SVD did not converge

[luisfabib]

Thank you! It seems the error does not arise from regparamrange as I thought, but rather from the covariance uncertainty quantification.

This is a new bug, a not a duplicate from #42

For now, you can get rid of this error by disabling the uncertainty quantification via the keyword: uqanalysis=False

[ritsch-i]

Tried it, when using: ridmefit[i] = dl.fitsignal(Vrid[i],trid[i],r[i],'P',dl.bg_strexp,dl.ex_4pdeer,uqanalysis=False), I get a new error:

AttributeError                            Traceback (most recent call last)
<ipython-input-480-3f535c4cc6fc> in <module>
     28     plt.plot(r[i],ridmefit[i].P,color='k')
     29     Puq = ridmefit[i].Puncert# Uncertainty quantification of the distance distribution
---> 30     Puq_95=Puq.ci(95)
     31     plt.fill_between(r[i],Puq_95[:,0],Puq_95[:,1],color='r',alpha=0.5)
     32     #plot RIDME-kernel result

AttributeError: 'list' object has no attribute 'ci'

[luisfabib]

The error LinAlgError: SVD did not converge in numpy.linalg.svd occurs most commonly when one of the terms passed to the function contains invalid values (nan,inf,...). Since the HCCM function just takes the Jacobians, these values must have appeared while computing the Jacobian which might happen when values are evaluated outside of the boundaries of a model.

This might have been solved now that the Jacobian engine has been modified in #55. Still this it is difficult to say without a working example, so more analysis is required.