Doubling algorithm not working for solving discrete Lyapunov equation

[shizejin]

For a specific example borrowed from RMT5 (example 3 in section 2.4), the doubling algorithm is unable to find the correct solution to the discrete Lyapunov equation while Bartels-Stewart algorithm (the scipy routine) is able to find one.

The demonstration of the failure is here.

[oyamad]

I have no knowledge about the theory, but the docstring says that "We assume in order for convergence that the eigenvalues of A have moduli bounded by unity". Your A matrix has some eigenvalues larger than one in magnitude:

w, v = np.linalg.eig(A)
np.abs(w)

array([1.14605125e+00, 1.14605125e+00, 8.50498522e-01, 8.50498522e-01,
       9.50203544e-01, 9.50203544e-01, 8.16036422e-01, 1.93952061e-17])

It might be an implicit premise that it is the user's responsibility to switch between the methods depending on the eigenvalues of A.

Maybe rewrite the docstring to say something like, "the method "doubling" assumes that the eigenvalues of A ...", and "if this is not the case, use the method "bartels-stewart""?

[shizejin]

Thanks @oyamad for your clear explanation. I guess the point is whether we intend to make this solve_discrete_lyapunov function only work for finding a "unique" solution, because having moduli bounded by one is not a necessary condition for the existence of the solution.

If we aim to make this general to use, then it might be burdensome to let users to check the eigenvalues each time before they call solve_discrete_lyapunov. Could we set the default method to be None and select one by checking the magnitudes of the eigenvalues inside the function?

But again, I am not sure about in what shape do we want method to be. For the purpose of economics research, having it work only with stable linear state models seems to be enough.

[rht]

Alternatively, if the method selected is doubling, there could be an assertion on the eigenvalues condition before the solving algorithm starts.