From daa5b21ef6c5363c668d16b2144dab999a78165b Mon Sep 17 00:00:00 2001 From: Victor Yu Date: Thu, 21 May 2026 10:46:30 -0500 Subject: [PATCH 1/3] Use eig + inv to compute exponential for non-Hermitian L Faster than expm. --- pycce/run/mastereq.py | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/pycce/run/mastereq.py b/pycce/run/mastereq.py index 1262300..256f975 100644 --- a/pycce/run/mastereq.py +++ b/pycce/run/mastereq.py @@ -1,5 +1,4 @@ import numpy as np -import scipy from pycce.bath.map import process_key_operator, process_key_dissipator from pycce.constants import PI2 from pycce.h import total_hamiltonian, projected_addition @@ -22,16 +21,13 @@ def simple_incoherent_propagator(timespace, lindbladian): ndarray with shape (n, N, N): Master equation propagators, evaluated at each timepoint. Use with vector form of density matrix. """ - return scipy.linalg.expm(timespace[:, np.newaxis, np.newaxis] * lindbladian[np.newaxis, :] * PI2) - - # This is how I did it for coherent operator, doesn't work with non-Hermitian L - # evalues, evec = np.linalg.eig(lindbladian * PI2) - # - # eigexp = np.exp(np.outer(timespace, evalues), - # dtype=np.complex128) - # - # return np.matmul(np.einsum('...ij,...j->...ij', evec, eigexp, dtype=np.complex128), - # evec.conj().T) + # Using eig is faster than expm for non-Hermitian L + # return scipy.linalg.expm(timespace[:, np.newaxis, np.newaxis] * lindbladian[np.newaxis, :] * PI2) + + evals, V = np.linalg.eig(lindbladian * PI2) + Vinv = np.linalg.inv(V) + eigexp = np.exp(np.outer(timespace, evals)) + return (V[None] * eigexp[:, None, :]) @ Vinv def simple_dissipator(left, right): From c4cb96d40e0487e543e2873fd6352461a04ae870 Mon Sep 17 00:00:00 2001 From: Victor Yu Date: Thu, 21 May 2026 17:53:19 -0500 Subject: [PATCH 2/3] Fall back to expm if eig + inv fails --- pycce/run/mastereq.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/pycce/run/mastereq.py b/pycce/run/mastereq.py index 256f975..26974c8 100644 --- a/pycce/run/mastereq.py +++ b/pycce/run/mastereq.py @@ -1,4 +1,5 @@ import numpy as np +import scipy from pycce.bath.map import process_key_operator, process_key_dissipator from pycce.constants import PI2 from pycce.h import total_hamiltonian, projected_addition @@ -21,13 +22,15 @@ def simple_incoherent_propagator(timespace, lindbladian): ndarray with shape (n, N, N): Master equation propagators, evaluated at each timepoint. Use with vector form of density matrix. """ - # Using eig is faster than expm for non-Hermitian L - # return scipy.linalg.expm(timespace[:, np.newaxis, np.newaxis] * lindbladian[np.newaxis, :] * PI2) - - evals, V = np.linalg.eig(lindbladian * PI2) - Vinv = np.linalg.inv(V) - eigexp = np.exp(np.outer(timespace, evals)) - return (V[None] * eigexp[:, None, :]) @ Vinv + # Non-Hermitian L, eig + inv faster than expm + try: + evals, V = np.linalg.eig(lindbladian * PI2) + Vinv = np.linalg.inv(V) + eigexp = np.exp(np.outer(timespace, evals)) + return (V[None] * eigexp[:, None, :]) @ Vinv + except np.linalg.LinAlgError as e: + print(f"eig + inv failed: {e}, falling back to expm") + return scipy.linalg.expm(timespace[:, np.newaxis, np.newaxis] * lindbladian[np.newaxis, :] * PI2) def simple_dissipator(left, right): From 231da121ddf15e9e74dfed82e09831bdccec3cd6 Mon Sep 17 00:00:00 2001 From: Victor Yu Date: Fri, 22 May 2026 10:04:28 -0500 Subject: [PATCH 3/3] Check condition number of eigenvector of L --- pycce/run/mastereq.py | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/pycce/run/mastereq.py b/pycce/run/mastereq.py index 26974c8..dbd3872 100644 --- a/pycce/run/mastereq.py +++ b/pycce/run/mastereq.py @@ -22,15 +22,26 @@ def simple_incoherent_propagator(timespace, lindbladian): ndarray with shape (n, N, N): Master equation propagators, evaluated at each timepoint. Use with vector form of density matrix. """ + L = lindbladian * PI2 + # Non-Hermitian L, eig + inv faster than expm try: - evals, V = np.linalg.eig(lindbladian * PI2) - Vinv = np.linalg.inv(V) + evals, V = np.linalg.eig(L) + # Reject ill-conditioned V + condV = np.linalg.cond(V) + if not np.isfinite(condV) or condV > 1e6: + raise np.linalg.LinAlgError(f"Eigenvector of L ill-conditioned (cond={condV:.2e})") + + Vinv = np.linalg.solve(V, np.eye(V.shape[0], dtype=V.dtype)) eigexp = np.exp(np.outer(timespace, evals)) + + if not np.all(np.isfinite(eigexp)): + raise np.linalg.LinAlgError("Non-finite values in exp(eigvals * t)") + return (V[None] * eigexp[:, None, :]) @ Vinv except np.linalg.LinAlgError as e: print(f"eig + inv failed: {e}, falling back to expm") - return scipy.linalg.expm(timespace[:, np.newaxis, np.newaxis] * lindbladian[np.newaxis, :] * PI2) + return scipy.linalg.expm(timespace[:, None, None] * L[None, :]) def simple_dissipator(left, right):