regparamrange crashes due to gsvd

[luisfabib]

The regparamrange employs the deerlab.utils.gsvd function which under certain circumstances leads to a crash. Until now crashes have been caused:

• Numerical errors internally of Numpy:

~\AppData\Local\Programs\Python\Python37-32\lib\site-packages\numpy\core\fromnumeric.py in amax(a, axis, out, keepdims, initial, where)
   2704     """
   2705     return _wrapreduction(a, np.maximum, 'max', axis, None, out,
-> 2706                           keepdims=keepdims, initial=initial, where=where)
   2707 
   2708 

~\AppData\Local\Programs\Python\Python37-32\lib\site-packages\numpy\core\fromnumeric.py in _wrapreduction(obj, ufunc, method, axis, dtype, out, **kwargs)
     85                 return reduction(axis=axis, out=out, **passkwargs)
     86 
---> 87     return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
     88 
     89 

ValueError: zero-size array to reduction operation maximum which has no identity

• Invalid values being sent to specific memory adresses in Scipy's Fortran subroutines:

/opt/hostedtoolcache/Python/3.8.6/x64/lib/python3.8/site-packages/scipy/linalg/decomp_qr.py:140: in qr
    qr, tau = safecall(geqrf, "geqrf", a1, lwork=lwork,
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 

f = <fortran object>, name = 'geqrf'
args = (array([], shape=(100, 0), dtype=float64),)
kwargs = {'lwork': 0, 'overwrite_a': False}, lwork = None
ret = (array([], shape=(100, 0), dtype=float64), array([], dtype=float64), array([0.]), -7)

    def safecall(f, name, *args, **kwargs):
        """Call a LAPACK routine, determining lwork automatically and handling
        error return values"""
        lwork = kwargs.get("lwork", None)
        if lwork in (None, -1):
            kwargs['lwork'] = -1
            ret = f(*args, **kwargs)
            kwargs['lwork'] = ret[-2][0].real.astype(numpy.int_)
        ret = f(*args, **kwargs)
        if ret[-1] < 0:
>           raise ValueError("illegal value in %dth argument of internal %s"
                             % (-ret[-1], name))
E           ValueError: illegal value in 7th argument of internal geqrf

Thus far these errors seem to depend on the OS, the precision of the hardware and the BLAS linked to Numpy.

Since neither numpy nor scipy have a working gsvd algorithm (that I could find), there are two options:

• Invest time into debugging and fixing deerlab.utils.gsvd
• Deprecate regparamrange. The function has lost some utility now that the regularization parameter is optimized via Golden-search algorithms rather than grid searches.

For now the best approach is to just patch regparamrange and catch the exception when running gsvd(K,L). Once caught we can just run svd(K) as an approximating alternative.

[luisfabib]

Here is a MWE to reproduce the first error described above:

import numpy as np 
import deerlab as dl
t = np.linspace(0,5,100)
r = np.linspace(2,6,150)
P = dl.dd_gauss(r,[4.5, 0.2])
K = dl.dipolarkernel(t,r)
V = K@P
fit = dl.fitsignal(V,t,r,'P',None,None,uqanalysis=False)

resulting in

deerlab\fitsignal.py:425: in fitsignal
    fit = dl.fitregmodel(Vexp,Ks,r,regtype,regparam, weights=weights,uqanalysis=uqanalysis)
deerlab\fitregmodel.py:131: in fitregmodel
    alpha = dl.selregparam(V,K,r,regtype,regparam,regorder=regorder,weights=weights, nonnegativity=nonnegativity,huberparam=huberparam)
deerlab\selregparam.py:111: in selregparam
    alphaCandidates = dl.regparamrange(K,L)
deerlab\regparamrange.py:41: in regparamrange
    singularValues = gsvd(K,L)
deerlab\utils\utils.py:245: in gsvd
    U,_,_,C,S = csd(Q1,Q2)
deerlab\utils\utils.py:278: in csd
    V,U,Z,S,C = csd(Q2,Q1)
deerlab\utils\utils.py:311: in csd
    k = max(0, np.max(np.where(np.diag(C) <= 1/math.sqrt(2))))
<__array_function__ internals>:6: in amax
    ???
C:\Users\luis\AppData\Local\Programs\Python\Python36\lib\site-packages\numpy\core\fromnumeric.py:2706: in amax
    keepdims=keepdims, initial=initial, where=where)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _  

obj = (array([], dtype=int64),), ufunc = <ufunc 'maximum'>, method = 'max', axis = None, dtype = None, out = None, kwargs = {'initial': <no value>, 'keepdims': <no value>, 'where': <no value>}, passkwargs = {}

    def _wrapreduction(obj, ufunc, method, axis, dtype, out, **kwargs):
        passkwargs = {k: v for k, v in kwargs.items()
                      if v is not np._NoValue}

        if type(obj) is not mu.ndarray:
            try:
                reduction = getattr(obj, method)
            except AttributeError:
                pass
            else:
                # This branch is needed for reductions like any which don't
                # support a dtype.
                if dtype is not None:
                    return reduction(axis=axis, dtype=dtype, out=out, **passkwargs)
                else:
                    return reduction(axis=axis, out=out, **passkwargs)

>       return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
E       ValueError: zero-size array to reduction operation maximum which has no identity

C:\Users\luis\AppData\Local\Programs\Python\Python36\lib\site-packages\numpy\core\fromnumeric.py:87: ValueError

[luisfabib]

After a short benchmark to get a statistical sense of these errors, it seems to happen too frequently:

import numpy as np 
import deerlab as dl
import matplotlib.pyplot as plt


alphamin = []
alphamax = []
Npass = 0
Nfails = 0
N = 500

for _ in range(N):

    tmax = np.random.uniform(2,8)
    rmax = np.random.uniform(2,8)
    nt = np.random.randint(80,900)
    nr = np.random.randint(80,900)

    t = np.linspace(-0.3,tmax,nt)
    r = np.linspace(1,rmax,nr)
    try:
        K = dl.dipolarkernel(t,r)
        L = dl.regoperator(r,2)
        alphas = dl.regparamrange(K,L)
        Npass += 1
    except: 
        Nfails += 1

[luisfabib]

Using the same benchmark we can see that the edges of the alpha-ranges are statistically quite localized. So the best solution would be to just deprecate the regparamrange function and use a fixed range for the golden-search in selregparam.