krylovsolve() error with same input as sesolve(), mesolve()

[Khamul]

Hello, 

I am trying to use krylovsolve() for state evolution. If I input the same exact Hamiltonian, initial state vector, time list, and expectation operators (with the addition of a krylov_dim parameter) as I do for sesolve() and mesolve(), I am unable to perform the evolution.

I get the following error:

ValueError: f(a) and f(b) must have different signs

Which appears to trace back to the brentq function call from the Scipy package (scipy > optimize >  _zeros_py.py line 784).

I have tried varying the krylov_dim parameter as it is the only difference, but I still cannot get it to work. I am hoping for suggestions, or if it is a bug to report the bug.

Thank you,

K

[Éric Giguère]

Could you post a code that fail and the version of qutip you are using?

[Khamul]

Hello,

I accidentally sent as a PM instead of a post, so will rewrite the message. My apologies for the duplication.

Qutip version is 5.0.1. I will address your questions below, please let me know if you need additional information.

Below is the code that fails. There is additional code beforehand to define system size, Hamiltonian, etc, but you can see here that the same inputs to sesolve(), mesolve(), and krylovsolve() (aside from c_ops and krylov_dim) result in failure for the krylovsolve() function call. I have worked extensively with sesolve() and mesolve() such that I know the Hamiltonian, initial state, e_ops, and c_ops are all 'good.'

local_coefficients = {

'omega': pt.omegas,

'delta': pt.deltas,

'r': pt.rs,

}

H_ryd_nn = Hamiltonian(local_coefficients, pt.c6, 1).rydberg()

H_ryd_ata = Hamiltonian(local_coefficients, pt.c6, test_length-1).rydberg()

# define initial and ideal final state

g = basis(2,0)

r = basis(2,1)

psi_0 = tensor([r] + [g]*(pt.length-1))

psi_f = tensor([g]*(pt.length-1) + [r])

exp_ops = [psi_f*psi_f.dag()]

c_ops = Hamiltonian.rydberg_collapse_operators(pt.length, gamma)

result_ryd_nn_se = sesolve(H_ryd_nn, psi_0, pt.times, e_ops=exp_ops)

result_ryd_ata_se = sesolve(H_ryd_ata, psi_0, pt.times, e_ops=exp_ops)

result_ryd_nn_me = mesolve(H_ryd_nn, psi_0, pt.times, c_ops=c_ops, e_ops=exp_ops)

result_ryd_ata_me = mesolve(H_ryd_ata, psi_0, pt.times, c_ops=c_ops, e_ops=exp_ops)

result_ryd_nn_krylov = krylovsolve(H_ryd_nn, psi_0, pt.times, e_ops=exp_ops, krylov_dim=10)

result_ryd_ata_krylov = krylovsolve(H_ryd_ata, psi_0, pt.times, e_ops=exp_ops, krylov_dim=10)

Additionally, below is the full error output statement:

--------------------------------------------------------------------------- ValueError Traceback (most recent call last) Cell In[7], line 7 4 result_ryd_nn_me = mesolve(H_ryd_nn, psi_0, pt.times, c_ops=c_ops, e_ops=exp_ops) 5 result_ryd_ata_me = mesolve(H_ryd_ata, psi_0, pt.times, c_ops=c_ops, e_ops=exp_ops) ----> 7 result_ryd_nn_krylov = krylovsolve(H_ryd_nn, psi_0, pt.times, e_ops=exp_ops, krylov_dim=10) 8 result_ryd_ata_krylov = krylovsolve(H_ryd_ata, psi_0, pt.times, e_ops=exp_ops, krylov_dim=10) File ~/miniconda3/envs/simulation/lib/python3.11/site-packages/qutip/solver/krylovsolve.py:101, in krylovsolve(H, psi0, tlist, krylov_dim, e_ops, args, options) 99 options["method"] = "krylov" 100 options["krylov_dim"] = krylov_dim --> 101 solver = SESolver(H, options=options) 102 return solver.run(psi0, tlist, e_ops=e_ops) File ~/miniconda3/envs/simulation/lib/python3.11/site-packages/qutip/solver/sesolve.py:150, in SESolver.__init__(self, H, options) 148 if not rhs.isoper: 149 raise ValueError("The hamiltonian must be an operator") --> 150 super().__init__(rhs, options=options) File ~/miniconda3/envs/simulation/lib/python3.11/site-packages/qutip/solver/solver_base.py:52, in Solver.__init__(self, rhs, options) 50 TypeError("The rhs must be a QobjEvo") 51 self.options = options ---> 52 self._integrator = self._get_integrator() 53 self._state_metadata = {} 54 self.stats = self._initialize_stats() File ~/miniconda3/envs/simulation/lib/python3.11/site-packages/qutip/solver/solver_base.py:226, in Solver._get_integrator(self) 224 else: 225 raise ValueError("Integrator method not supported.") --> 226 integrator_instance = integrator(self.rhs, self.options) 227 self._init_integrator_time = time() - _time_start 228 return integrator_instance File ~/miniconda3/envs/simulation/lib/python3.11/site-packages/qutip/solver/integrator/integrator.py:76, in Integrator.__init__(self, system, options) 74 self._options = self.integrator_options.copy() 75 self.options = options ---> 76 self._prepare() File ~/miniconda3/envs/simulation/lib/python3.11/site-packages/qutip/solver/integrator/krylov.py:59, in IntegratorKrylov._prepare(self) 57 self._max_step = np.inf 58 else: ---> 59 self._max_step = self._compute_max_step(krylov_tridiag, 60 krylov_basis) File ~/miniconda3/envs/simulation/lib/python3.11/site-packages/qutip/solver/integrator/krylov.py:156, in IntegratorKrylov._compute_max_step(self, krylov_tridiag, krylov_basis, krylov_state) 153 if dt > self.options["max_step"]: 154 return self.options["max_step"] --> 156 sol = root_scalar(f=krylov_error, bracket=[dt, dt * 10], 157 method="brentq", xtol=self.options['atol']) 158 if sol.converged: 159 return sol.root File ~/miniconda3/envs/simulation/lib/python3.11/site-packages/scipy/optimize/_root_scalar.py:275, in root_scalar(f, args, method, bracket, fprime, fprime2, x0, x1, xtol, rtol, maxiter, options) 272 raise ValueError('Bracket needed for %s' % method) 274 a, b = bracket[:2] --> 275 r, sol = methodc(f, a, b, args=args, **kwargs) 276 elif meth in ['secant']: 277 if x0 is None: File ~/miniconda3/envs/simulation/lib/python3.11/site-packages/scipy/optimize/_zeros_py.py:784, in brentq(f, a, b, args, xtol, rtol, maxiter, full_output, disp) 782 if rtol < _rtol: 783 raise ValueError("rtol too small (%g < %g)" % (rtol, _rtol)) --> 784 r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp) 785 return results_c(full_output, r) ValueError: f(a) and f(b) must have different signs

[Éric Giguère]

Try settings `min_step` too something smaller than the default in the options:

```

options = {"min_step": 1e-8}

qutip.krylovsolve(..., options=options)

```

When the desired tolerance cannot be reached with the minimum step, that error will be raised.

ps. I made a patch to improve the error message (https://github.com/qutip/qutip/pull/2402).

[Khamul]

This indeed worked. For my purposes, I needed to set min_step to 1e-9.

Thank you for the speedy responses!