Error in running of phiD Wilson coefficient

[peterstangl]

When trying to run Wilson coefficients containing a non-zero value of phiD, I receive the following error I am not sure where it is coming from:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-9-effb4070f91a> in <module>()
----> 1 mywilson.match_run(1e3, 'SMEFT', 'Warsaw')

~/.local/lib/python3.6/site-packages/wilson/classes.py in match_run(self, scale, eft, basis, sectors)
     90         if self.wc.eft == 'SMEFT':
     91             if eft == 'SMEFT':
---> 92                 smeft = SMEFT(self.wc)
     93                 # if input and output EFT ist SMEFT, just run.
     94                 wc_out = smeft.run(scale)

~/.local/lib/python3.6/site-packages/wilson/run/smeft/classes.py in __init__(self, wc, get_smpar)
     34         self.C_in = None
     35         if wc is not None:
---> 36             self._set_initial_wcxf(wc, get_smpar=get_smpar)
     37 
     38     def _set_initial(self, C_in, scale_in):

~/.local/lib/python3.6/site-packages/wilson/run/smeft/classes.py in _set_initial_wcxf(self, wc, get_smpar)
     78             self.C_in.update(C)
     79         if get_smpar:
---> 80             self.C_in.update(self._get_sm_scale_in())
     81 
     82 

~/.local/lib/python3.6/site-packages/wilson/run/smeft/classes.py in _get_sm_scale_in(self, scale_sm)
    177         _smeft = SMEFT(self.wc, get_smpar=False)
    178         # Step 1: run the SM up, using the WCs at scale_input as (constant) estimate
--> 179         _smeft.C_in.update(self._run_sm_scale_in(self.C_in, scale_sm=scale_sm))
    180         # Step 2: run the WCs down in LL approximation
    181         C_out = _smeft._rgevolve_leadinglog(scale_sm)

~/.local/lib/python3.6/site-packages/wilson/run/smeft/classes.py in _run_sm_scale_in(self, C_out, scale_sm)
    153         C_in_sm = beta.C_array2dict(np.zeros(9999))
    154         # set the SM parameters to the values obtained from smpar.smeftpar
--> 155         C_SM = smpar.smeftpar(scale_sm, C_out, basis='Warsaw')
    156         SM_keys = set(smeftutil.SM_keys)  # to speed up lookup
    157         C_SM = {k: v for k, v in C_SM.items() if k in SM_keys}

~/.local/lib/python3.6/site-packages/wilson/run/smeft/smpar.py in smeftpar(scale, C, basis)
     98     vb = sqrt(1 / sqrt(2) / GF)
     99     v = vb  # TODO
--> 100     _d = vMh2_to_m2Lambda(v=v, Mh2=Mh**2, C=C)
    101     m2 = _d['m2'].real
    102     Lambda = _d['Lambda'].real

~/.local/lib/python3.6/site-packages/wilson/run/smeft/smpar.py in vMh2_to_m2Lambda(v, Mh2, C)
     66         x0 = np.array([dSM['m2'], dSM['Lambda']])
     67         try:
---> 68             xres = scipy.optimize.newton_krylov(f0, x0)
     69         except scipy.optimize.nonlin.NoConvergence:
     70             raise ValueError("No solution for m^2 and Lambda found")

~/.local/lib/python3.6/site-packages/scipy/optimize/nonlin.py in newton_krylov(F, xin, iter, rdiff, method, inner_maxiter, inner_M, outer_k, verbose, maxiter, f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search, callback, **kw)

~/.local/lib/python3.6/site-packages/scipy/optimize/nonlin.py in nonlin_solve(F, x0, jacobian, iter, verbose, maxiter, f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search, callback, full_output, raise_exception)
    310 
    311         if norm(dx) == 0:
--> 312             raise ValueError("Jacobian inversion yielded zero vector. "
    313                              "This indicates a bug in the Jacobian "
    314                              "approximation.")

ValueError: Jacobian inversion yielded zero vector. This indicates a bug in the Jacobian approximation.

This error is generated by the following code:

from wilson import Wilson
mywilson = Wilson({'phiD': 1}, 1e3, 'SMEFT', 'Warsaw')
mywilson.match_run(1e2, 'SMEFT', 'Warsaw')

[DavidMStraub]

Similarly to #7, this is due essentially to a crazily large WC ... in this case, the determination of the Higgs potential parameters fails. In this case there is even a try ... except block in lines 67-70, but the error thrown here is not caught; this is why the error message is not very enlightening.

[peterstangl]

Ok, so my comment in #7 also applies here.

I will keep the issue open as a reminder to perhaps catch the error and add a corresponding error message.

[peterstangl]

The NoConvergence error actually also appears with far less crazily large WCs.

E.g.

mywilson = Wilson({'phiD': 9.678e-06}, 1e3, 'SMEFT', 'Warsaw')

and

mywilson = Wilson({'phiD': 9.679e-06}, 1e3, 'SMEFT', 'Warsaw')

yield an error when trying to run the WCs.

However, for slightly smaller or larger values

mywilson = Wilson({'phiD': 9.677e-06}, 1e3, 'SMEFT', 'Warsaw')

and

mywilson = Wilson({'phiD': 9.680e-06}, 1e3, 'SMEFT', 'Warsaw')

the optimization algorithm does converge.

Changing the algorithm in https://github.com/wilson-eft/wilson/blob/a4bb7946218e2f569470ff080e7ec28805d661c7/wilson/run/smeft/smpar.py#L68 from the standard method='lgmres' to method='gmres' leads to convergence in all of the above cases. This is kind of unexpected as according to https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.lgmres.html,

The LGMRES algorithm [BJM] [BPh] is designed to avoid some problems in the convergence in restarted GMRES, and often converges in fewer iterations.

An opposite behaviour is observed here.

Only changing the method from method='lgmres' to method='gmres' leads to some tests failing (test_Cphi0 with AssertionError: 246.00000012925454 != 246 within 7 places and test_smpar_roundtrip with AssertionError: 0.9999936694445539 != 1 within 1e-06 delta : Failed for m_e). However, this can be resolved by also reducing the tolerance in the optimization by adding f_tol=1e-7, which does not significantly increase the computing time for the four above cases.

I've tried the other available algorithms ‘bicgstab’, ‘cgs’, ‘minres’ but they lead to non-convergence or other errors when running the tests. So only method='lgmres' and method='gmres', f_tol=1e-7 pass all tests, while the latter also leads to convergence in all of the four above cases. I would thus suggest to either replace the former version by the latter one or to modify the try ... except block such that in case of non-convergence of LGMRES, the GMRES algorithm is tried before throwing an error.

[DavidMStraub]

Thanks, very interesting. I would avoid modifying the default chosen by scipy unless necessary, so indeed I think the best way to solve this would be to nest another try ... except into the original except that tries with gmres. Your examples could be added as unit tests.

[peterstangl]

I have solved this issue in PR #10.

The nested try/except construction may lead to a doubling of computing time. To see how often this actually happens, a warning is issued if the standard algorithm fails to converge. If this should happen frequently, the computational speed could be increased by always using the LGMRES algorithm.

[DavidMStraub]

Perfect!