qig.duhamel
Tools for evaluating the Duhamel (Wilcox) formula for quantum exponential-family derivatives,
where \(\rho = e^{H}\) with \(H = \sum_a \theta_a F_a - \psi(\theta) I\).
Conceptually there are three evaluation strategies:
Quadrature-based Duhamel (
duhamel_derivative,duhamel_derivative_simpson):Treat the integral over \(s \in [0,1]\) literally and approximate it with trapezoid or Simpson rules.
Works for any Hermitian \(H\), with accuracy controlled by the number of quadrature points.
Used as a high-precision reference in tests and for legacy
method='duhamel'paths.
Spectral/BCH Duhamel (
duhamel_derivative_spectral):Use the eigen-decomposition \(H = U \operatorname{diag}(\lambda) U^\dagger\) and evaluate the Fréchet derivative of the exponential in that basis.
In the eigenbasis of \(H\) the derivative has closed form
\[\begin{split}\bigl(D\exp_H[X]\bigr)_{ij} = \begin{cases} e^{\lambda_i} X_{ii}, & i = j,\\[4pt] \displaystyle X_{ij}\, \frac{e^{\lambda_i}-e^{\lambda_j}}{\lambda_i-\lambda_j}, & i \neq j, \end{cases}\end{split}\]which can be written as a matrix function \(f(\mathrm{ad}_H)\) applied to \(X = F_a - \langle F_a \rangle I\), with \(f(z) = (e^z - 1)/z\).
This realises the Lie-closure/BCH observation from the paper in a concrete numerical form: when the operator basis closes under commutation, the Duhamel integral becomes a finite-dimensional linear map on the Lie algebra, and can be evaluated analytically (up to diagonalisation error) without explicit \(s\)-quadrature.
Block-matrix Duhamel (
duhamel_derivative_block):Uses Higham’s block-matrix identity to compute the Fréchet derivative via a single \(2n \times 2n\) matrix exponential:
\[\begin{split}\exp\begin{pmatrix} H & F_a - \langle F_a \rangle I \\ 0 & H \end{pmatrix} = \begin{pmatrix} e^H & D\exp_H[F_a - \langle F_a \rangle I] \\ 0 & e^H \end{pmatrix}\end{split}\]The (1,2) block of the result is the Duhamel integral.
Avoids both explicit quadrature and eigendecomposition, instead “compiling away” the integral into the exponential itself.
More robust than spectral method for ill-conditioned \(H\) where eigendecomposition may be unstable.
Cost: one call to a highly-optimized
expmroutine (Padé approximation with scaling and squaring) on a \(2n \times 2n\) matrix.Best for small to medium systems (\(n \leq 100\)) when numerical robustness is important.
Method Selection
Method |
Cost |
Accuracy |
Best For |
|---|---|---|---|
Quadrature |
50 |
~10⁻⁵ |
Validation |
Spectral |
1 |
Machine precision |
Well-conditioned H |
Block |
1 |
Machine precision |
Ill-conditioned H |
SLD |
2 evaluations |
~10⁻³ |
Fast approximation |
The quadrature, spectral, and block implementations are numerically cross-validated in
the test suite (see TestRhoDerivativeNumerical in
tests/test_pair_exponential_family.py and TestBlockFrechet in
tests/test_block_frechet.py).
For small finite-dimensional examples, both the spectral and block methods achieve machine precision. Choose based on your needs:
Spectral: Faster when eigendecomposition is cheap and well-conditioned; aligns with the adjoint/BCH structure in theory
Block: More robust for ill-conditioned problems; leverages highly-optimized
expmwithout eigendecomposition
See CIP-000A for detailed comparison and implementation notes.
References
Higham, N. J. (2008). Functions of Matrices: Theory and Computation. SIAM. Chapter 10: The Fréchet Derivative.
Al-Mohy, A. H., & Higham, N. J. (2009). Computing the Fréchet derivative of the matrix exponential, with an application to condition number estimation. SIAM J. Matrix Anal. Appl., 30(4), 1639–1657.
Duhamel formula for quantum exponential family derivatives.
For ρ = exp(H) where H = ∑ θ_a F_a - ψ(θ)I, the exact derivative is:
∂ρ/∂θ_a = ∫₀¹ exp(sH) (F_a - ⟨F_a⟩I) exp((1-s)H) ds
This is the Duhamel (or Dalecki-Krein exponential) formula. The SLD is just the trapezoid rule approximation (average of s=0 and s=1 endpoints).
QUANTUM DERIVATIVE PRINCIPLE: This is the EXACT formula that respects operator ordering and preserves Hermiticity.
- qig.duhamel.duhamel_derivative(rho, H, F_centered, n_points=50)[source]
Compute ∂ρ/∂θ using the Duhamel formula with numerical integration.
∂ρ/∂θ = ∫₀¹ exp(sH) F_centered exp((1-s)H) ds
where F_centered = F - ⟨F⟩I.
- Parameters:
rho (
ndarray,shape (D,D)) – Density matrixH (
ndarray,shape (D,D)) – Hamiltonian H = ∑ θ_a F_a - ψ(θ)IF_centered (
ndarray,shape (D,D)) – Centered operator F - ⟨F⟩In_points (
int) – Number of quadrature points for integration Default: 50 (gives ~5e-05 error for typical cases)
- Returns:
drho – Derivative ∂ρ/∂θ (exactly Hermitian)
- Return type:
ndarray,shape (D,D)
Notes
Uses trapezoid rule for integration
For n_points=2, recovers the SLD formula
Convergence (typical case): n=50→6e-05, n=100→1.5e-05, n=200→3.6e-06
For high precision, use n_points ≥ 200 or use theta_only method instead
- qig.duhamel.duhamel_derivative_block(rho, H, F_centered)[source]
Compute ∂ρ/∂θ using Higham’s block-matrix identity (Fréchet derivative of exp).
This evaluates the same Duhamel integral as duhamel_derivative_spectral:
D exp_H[F_centered] = ∫₀¹ e^{(1-s)H} F_centered e^{sH} ds
but without eigendecomposition or explicit quadrature. Instead, form the 2n×2n block matrix:
- B = [[H, F_centered],
[0, H]]
Then:
- exp(B) = [[exp(H), D exp_H[F_centered]],
[0, exp(H)]]
so the (1,2) block is the desired derivative.
Notes
Cost is one expm on a (2D)×(2D) matrix.
For Hermitian H and Hermitian F_centered, the exact result is Hermitian; we symmetrise to guard against tiny numerical asymmetry.
References
Higham, N. J. (2008). Functions of Matrices: Theory and Computation. SIAM. Chapter 10.
Al-Mohy, A. H., & Higham, N. J. (2009). Computing the Fréchet derivative of the matrix exponential, with an application to condition number estimation. SIAM J. Matrix Anal. Appl., 30(4), 1639–1657.
- qig.duhamel.expm_frechet_block_1(A, E)[source]
First Fréchet derivative of exp at A in direction E via 2×2 block trick.
Returns D exp_A[E].
- qig.duhamel.expm_frechet_block_2(A, E, F)[source]
Second Fréchet derivative of exp at A in directions (E, F) via 3×3 block trick.
Constructs the 3n×3n block matrix:
- [[A, E, 0],
[0, A, F], [0, 0, A]]
Then (1,3) block of exp is D² exp_A[E, F].
- qig.duhamel.expm_frechet_block_3(A, E, F, G)[source]
Third Fréchet derivative of exp at A in directions (E, F, G) via 4×4 block trick.
Constructs the 4n×4n block matrix:
- [[A, E, 0, 0],
[0, A, F, 0], [0, 0, A, G], [0, 0, 0, A]]
Then (1,4) block of exp is D³ exp_A[E, F, G].
- qig.duhamel.duhamel_derivative_simpson(rho, H, F_centered, n_points=51)[source]
Compute ∂ρ/∂θ using Duhamel formula with Simpson’s rule.
More accurate than trapezoid rule for smooth integrands.
- Parameters:
rho (
ndarray,shape (D,D)) – Density matrixH (
ndarray,shape (D,D)) – HamiltonianF_centered (
ndarray,shape (D,D)) – Centered operatorn_points (
int,odd) – Number of quadrature points (must be odd)
- Returns:
drho – Derivative ∂ρ/∂θ
- Return type:
ndarray,shape (D,D)
- qig.duhamel.duhamel_derivative_spectral(rho, H, F_centered)[source]
Compute ∂ρ/∂θ using the Duhamel formula via the spectral representation of H.
- This evaluates the Fréchet derivative of the matrix exponential
D exp_H[F_centered] = ∫₀¹ e^{(1-s)H} F_centered e^{sH} ds
exactly (up to diagonalisation error), by working in the eigenbasis of H.
- In that basis we have, for H = U diag(λ) U†,
- (D exp_H[F_centered])_{ij} =
- { e^{λ_i} F_{ij} if i = j,
F_{ij} (e^{λ_i} - e^{λ_j}) / (λ_i - λ_j) if i ≠ j }.
- Parameters:
rho (
ndarray,shape (D,D)) – Density matrix (unused, included for API symmetry with duhamel_derivative)H (
ndarray,shape (D,D)) – Hamiltonian H = ∑ θ_a F_a - ψ(θ)I (Hermitian)F_centered (
ndarray,shape (D,D)) – Centered operator F - ⟨F⟩I
- Returns:
drho – Derivative ∂ρ/∂θ (exact up to eigen-decomposition accuracy)
- Return type:
ndarray,shape (D,D)