torch_vn_algebra.divergences¶
Source: torch_vn_algebra/divergences.py. Signatures and docstrings as in the source.
Module documentation¶
Quantum Renyi divergences, Petz recovery maps and noncommutative L^p norms on M = (+)_c M_{k_c}.
Conventions
-----------
States are normal positive functionals on M, represented (as everywhere in the library, see states.py) by their
densities with respect to the blunt trace Tr = sum_c Tr_c: rho = (+)_c rho_c >= 0, omega(A) = Tr(rho A). Plain
positive Operators are accepted as well (they need not be normalised); all functions are batched over the leading
dimension (batch 1 broadcasts against batch B) and return tensors of shape (batch,) in the real dtype of the
algebra. Logarithms are natural (divergences in nats). Functions of positive operators are computed by the spectral
calculus of the library, one batched eigendecomposition per distinct block size, with gradients by the
Daleckii-Krein formula (functional.spectral_apply); there are no host synchronisations in the computational paths.
Definitions (rho, sigma >= 0, rho != 0, alpha in (0, 1) u (1, inf); for states Tr rho = 1 and the normalisation
below is void)
Petz [Petz 1986]: D_alpha(rho||sigma) = (alpha-1)^{-1} log( Tr rho^alpha sigma^{1-alpha} / Tr rho )
sandwiched [MDSFT 2013, WWY 2014]:
~D_alpha(rho||sigma) = (alpha-1)^{-1} log( Tr (sigma^g rho sigma^g)^alpha / Tr rho ),
g = (1-alpha)/(2 alpha)
Umegaki (alpha -> 1 of both families):
D(rho||sigma) = Tr rho (log rho - log sigma) / Tr rho
max-relative entropy (alpha -> inf of the sandwiched family):
D_max(rho||sigma) = log min{lambda : rho <= lambda sigma}
= log lambda_max(sigma^{-1/2} rho sigma^{-1/2})
alpha -> 0 of the Petz family (Datta's "min-relative entropy"; not the sandwiched alpha = 1/2 quantity that
Tomamichel's book calls D_min):
D_0(rho||sigma) = -log( Tr(Pi_rho sigma) / Tr rho ), Pi_rho = support of rho
fidelity: ~D_{1/2}(rho||sigma) = -log F(rho, sigma) with F = ||sqrt(rho) sqrt(sigma)||_1^2 (states.fidelity).
The quantities Q_alpha = Tr rho^alpha sigma^{1-alpha} and ~Q_alpha = Tr (sigma^g rho sigma^g)^alpha are returned per
sector by sector_contributions.
Support conventions
-------------------
Powers (including negative ones) and logarithms of a positive operator are taken on its numerical support: the
eigenvalues above rtol * lambda_max (lambda_max over all sectors of the sample; default rtol = N eps with
N = sum_c k_c, the cut-off of numpy.linalg.matrix_rank) are kept, the others count as exact zeros, and f(0) := 0 for every
power p != 0 as well as for the logarithm (x^0 is the support projection). Hence Q_alpha and ~Q_alpha only involve
the compression of rho to supp sigma. The support condition supp rho <= supp sigma is tested through the weight of
rho outside supp sigma, leak = Tr rho (1 - Pi_sigma); it counts as violated when leak > supp_tol * Tr rho (default
supp_tol = 100 N eps). Then
alpha > 1 (both families), D (Umegaki) and D_max: the divergence is +inf;
alpha < 1: the divergence is finite unless Q_alpha = 0 (rho and sigma with orthogonal supports): +inf.
The dual test decides orthogonality: a sector in which the weight of rho inside supp sigma, Tr rho_c Pi_sigma_c
(computed from the eigenpairs, accurate also for orthogonal supports), is at most supp_tol * Tr rho_c contributes
Q_c = 0 (otherwise rounding-level overlaps of orthogonal eigenvectors, ~ eps^2, would give finite values ~ -log eps^2
/ (1 - alpha)); for D_0 a sector with Tr(Pi_rho_c sigma_c) <= supp_tol * Tr sigma_c contributes 0.
These are the standard conventions (Tomamichel 2016, Ch. 4); they are the limits eps -> 0 of the divergences of rho
and sigma + eps 1. The thresholds are relative, so all results are invariant under rescaling of the inputs.
Numerical evaluation
--------------------
The Renyi quantities are evaluated per sector as logarithms log Q_c, with rescaled factors (every factor <= 1, scales
per block), so that large alpha (D_max as a limit) neither overflows nor underflows, also in single precision; the
divergence is (alpha-1)^{-1} (logsumexp_c log Q_c - log Tr rho). The Petz value is the exact pair sum
Q_c = sum_ij lambda_i^alpha mu_j^{1-alpha} |<v_i|u_j>|^2 over the eigenpairs of rho_c and sigma_c (a logsumexp; the
matrix form Tr rho^alpha sigma^{1-alpha} loses all digits when its dominant pair is not (lambda_max, mu_min), as for
nearly commuting states at large alpha); its gradient is that of the matrix form (Daleckii-Krein), which is dropped
where the rescaled matrix form is below sqrt(eps). The spectrum of the sandwich sigma^g rho sigma^g is obtained as
the squared singular values of sigma^g rho^{1/2}: the SVD resolves it down to (eps ||sigma^g|| ||rho||^{1/2})^2, not
only down to eps ||sigma^g rho sigma^g|| as the eigenvalues of the product would (this matters for alpha < 1, where
eigenvalues lambda contribute lambda^alpha); singular values below rtol * max(s_max, ||sigma^g|| ||rho||^{1/2}) (per
block) count as zeros. Limit of accuracy: for alpha < 1 the value depends on all eigenvalues of the sandwich, whose
condition number can reach cond(sigma)^{2g} cond(rho) with g = (1-alpha)/(2 alpha) -> inf as alpha -> 0; eigenvalues
below the resolution above cannot be represented, so for small alpha and ill-conditioned states the computed
sandwiched divergence is not accurate (see scripts/renyi_threshold.py, which therefore restricts alpha >= 0.2 and
bounds the condition of the states).
Trace conventions
-----------------
With a trace tau_w = sum_c w_c Tr_c (Tr_blunt: w = 1, Tr_norm: w_c = 1/k_c, tau_vN: w_c = 1/(C k_c)), the density of
the state omega with respect to tau_w is d = (+)_c rho_c / w_c. Every divergence above, computed with tau_w from the
densities d (powers of d and traces tau_w), equals the one computed with Tr_blunt from rho: the central factors
w_c^{-alpha} w_c^{alpha-1} w_c (Petz), w_c^{-alpha(2g+1)} w_c = 1 (sandwiched) etc. cancel. Divergences are therefore
functions of the states alone. The argument trace='norm' / 'tau_vN' declares plain Operator inputs to be densities
with respect to that trace (d = DensityMatrix.density(trace)). The L^p norms ||A||_{p,w} do depend on the trace; the
state-weighted (Kosaki) norms ||x||_{p,rho} again do not.
Petz recovery map
-----------------
For a completely positive map Phi: M_in -> M_out (Channel or InterSectorChannel) and sigma >= 0 on M_in,
R_{sigma,Phi}(X) = sigma^{1/2} Phi^*( Phi(sigma)^{-1/2} X Phi(sigma)^{-1/2} ) sigma^{1/2},
with the adjoint Phi^* for Tr_blunt (the trace in which states are densities) and the inverse square root on
supp Phi(sigma). R is completely positive with Tr R(X) = Tr(Pi_{Phi(sigma)} X), and R(Phi(sigma)) = sigma when Phi is
trace preserving. The rotated maps R^t = U_{sigma,t} o R o U_{Phi(sigma),-t}, U_{s,t}(X) = s^{it} X s^{-it}
[JRSWW 2018], and the universal recovery map R_univ = int beta_0(t) R^{t/2} dt, beta_0(t) = pi/2 (cosh(pi t) + 1)^{-1},
satisfy D(rho||sigma) - D(Phi rho||Phi sigma) >= -log F(rho, R_univ(Phi(rho))) [JRSWW 2018, Thm. 2.1].
L^p norms
---------
||A||_{p,w} = (sum_c w_c Tr|A_c|^p)^{1/p} for p in [1, inf), ||A||_inf = max_c ||A_c|| (sectors with w_c > 0), and
the state-weighted norms of Kosaki [Kosaki 1984] ||x||_{p,rho,eta} = ||rho^{eta/p} x rho^{(1-eta)/p}||_p (eta = 1/2:
the symmetric embedding used by Beigi 2013, eta = 1 / 0: left / right), i.e. the norm in the Haagerup L^p space of
the image of x under the embedding x -> d^{eta/p} x d^{(1-eta)/p} (d the density of the state; in finite dimensions
the Haagerup L^p space is M with the Schatten norm of the trace). For example
~D_alpha(rho||sigma) = alpha/(alpha-1) log ||sigma^{-1/2} rho sigma^{-1/2}||_{alpha,sigma} [Beigi 2013].
Data processing
---------------
For every completely positive trace-preserving Phi (theorems; tested on random instances in
tests/test_divergences.py and explored by scripts/renyi_threshold.py):
~D_alpha(Phi rho||Phi sigma) <= ~D_alpha(rho||sigma) for alpha >= 1/2 [Frank-Lieb 2013; Beigi 2013 (alpha > 1);
MDSFT 2013, WWY 2014 (1 < alpha <= 2)],
D_alpha(Phi rho||Phi sigma) <= D_alpha(rho||sigma) for 0 < alpha <= 2 [Petz 1986],
and both fail in general outside these ranges. Further properties used in the tests: ~D_alpha <= D_alpha
(Araki-Lieb-Thirring), both families nondecreasing in alpha (Petz on (0, inf), sandwiched on [1/2, inf)), and
positive for states (sandwiched for alpha >= 1/2) with equality iff rho = sigma.
References
----------
D. Petz, Quasi-entropies for finite quantum systems, Rep. Math. Phys. 23, 57 (1986).
M. Mueller-Lennert, F. Dupuis, O. Szehr, S. Fehr, M. Tomamichel, On quantum Renyi entropies: a new generalization
and some properties, J. Math. Phys. 54, 122203 (2013). [MDSFT 2013]
M. M. Wilde, A. Winter, D. Yang, Strong converse for the classical capacity of entanglement-breaking and Hadamard
channels via a sandwiched Renyi relative entropy, Commun. Math. Phys. 331, 593 (2014). [WWY 2014]
R. L. Frank, E. H. Lieb, Monotonicity of a relative Renyi entropy, J. Math. Phys. 54, 122201 (2013).
S. Beigi, Sandwiched Renyi divergence satisfies data processing inequality, J. Math. Phys. 54, 122202 (2013).
M. Junge, R. Renner, D. Sutter, M. M. Wilde, A. Winter, Universal recovery maps and approximate sufficiency of
quantum relative entropy, Ann. Henri Poincare 19, 2955 (2018). [JRSWW 2018]
H. Kosaki, Applications of the complex interpolation method to a von Neumann algebra: non-commutative L^p-spaces,
J. Funct. Anal. 56, 29 (1984).
N. Datta, Min- and max-relative entropies and a new entanglement monotone, IEEE Trans. Inf. Theory 55, 2816 (2009).
M. Tomamichel, Quantum Information Processing with Finite Resources, Springer (2016).
petz_renyi_divergence¶
def petz_renyi_divergence(rho: Operator, sigma: Operator, alpha: Alpha, trace: str='blunt', rtol: Optional[float]=None, supp_tol: Optional[float]=None) -> torch.Tensor
Petz-Renyi divergence D_alpha(rho||sigma) = (alpha-1)^{-1} log(Tr rho^alpha sigma^{1-alpha} / Tr rho), shape
(batch,).
rho, sigma: states (DensityMatrix) or positive Operators of the same algebra, batch 1 or B.
alpha: a number in (0, 1) u (1, inf) (checked), or a tensor (batch,) of positive values, one per sample (not
checked; entries alpha = 1 give the Umegaki relative entropy, computed in the same pass).
trace: the trace for which plain Operator inputs are densities (module docstring; the value does not depend on it).
rtol, supp_tol: support thresholds (module docstring). +inf for alpha > 1 if supp rho is not in supp sigma, and
for alpha < 1 if the supports are orthogonal.
Differentiable in rho, sigma (and alpha) at points where the supports are stable.
sandwiched_renyi_divergence¶
def sandwiched_renyi_divergence(rho: Operator, sigma: Operator, alpha: Alpha, trace: str='blunt', rtol: Optional[float]=None, supp_tol: Optional[float]=None) -> torch.Tensor
Sandwiched Renyi divergence ~D_alpha(rho||sigma) = (alpha-1)^{-1} log(Tr (sigma^g rho sigma^g)^alpha / Tr rho),
g = (1-alpha)/(2 alpha), shape (batch,). Arguments as for petz_renyi_divergence; ~D_{1/2} = -log F (fidelity).
umegaki_relative_entropy¶
def umegaki_relative_entropy(rho: Operator, sigma: Operator, trace: str='blunt', rtol: Optional[float]=None, supp_tol: Optional[float]=None) -> torch.Tensor
Umegaki relative entropy D(rho||sigma) = Tr rho (log rho - log sigma) / Tr rho, shape (batch,), with the exact
support convention (+inf if supp rho is not in supp sigma; 0 log 0 = 0). For states with supp rho in supp sigma
it agrees with states.relative_entropy, which instead clamps the spectrum of sigma (finite values ~ -log eps).
max_relative_entropy¶
def max_relative_entropy(rho: Operator, sigma: Operator, trace: str='blunt', rtol: Optional[float]=None, supp_tol: Optional[float]=None) -> torch.Tensor
D_max(rho||sigma) = log min{lambda: rho <= lambda sigma} = log lambda_max(sigma^{-1/2} rho sigma^{-1/2}) (the
largest over all sectors), +inf if supp rho is not in supp sigma; the limit alpha -> inf of the sandwiched
divergences. Shape (batch,).
min_relative_entropy¶
def min_relative_entropy(rho: Operator, sigma: Operator, trace: str='blunt', rtol: Optional[float]=None, supp_tol: Optional[float]=None) -> torch.Tensor
D_0(rho||sigma) = -log(Tr(Pi_rho sigma) / Tr rho), the limit alpha -> 0 of the Petz divergences (Datta's
min-relative entropy), +inf if the supports are orthogonal (a sector with Tr(Pi_rho_c sigma_c) <= supp_tol *
Tr sigma_c, default 100 N eps, contributes 0). Shape (batch,). Differentiable in sigma only (the support
projection Pi_rho is piecewise constant in rho).
renyi_divergence¶
def renyi_divergence(rho: Operator, sigma: Operator, alpha: Alpha, kind: str='sandwiched', **kwargs) -> torch.Tensor
D_alpha (kind='petz') or ~D_alpha (kind='sandwiched'), including the limits for a number alpha: alpha = 1 gives
the Umegaki relative entropy (both kinds), alpha = inf gives D_max (sandwiched), alpha = 0 gives D_0 (Petz).
Keyword arguments (trace, rtol, supp_tol) as for petz_renyi_divergence.
SectorContributions¶
class SectorContributions()
Per-sector decomposition of a divergence of states on M = (+)_c M_{k_c}; all fields (batch, C) except alpha.
p, q: sector weights Tr rho_c, Tr sigma_c.
terms: Renyi kinds: Q_c = Tr rho_c^alpha sigma_c^{1-alpha} (Petz) or Tr (sigma_c^g rho_c sigma_c^g)^alpha
(sandwiched), so that Tr rho^alpha sigma^{1-alpha} = sum_c Q_c (and likewise ~Q); these may
overflow or underflow for very large alpha, the conditional divergences are computed from log Q_c
and do not;
Umegaki: Tr rho_c (log rho_c - log sigma_c), summing to Tr rho D(rho||sigma).
conditional: D(rho_c / p_c || sigma_c / q_c), the divergence of the conditional states of sector c: NaN where
p_c = 0 (rounding level: p_c <= rtol Tr rho), +inf where q_c = 0 < p_c or where the support
conventions of the module give +inf for the sector; for alpha given per sample, entries alpha = 1
give the conditional Umegaki relative entropies (terms are then Tr rho_c Pi_sigma_c).
alpha: (batch,).
Identities (sectors with p_c, q_c > 0): Q_c = p_c^alpha q_c^{1-alpha} exp((alpha-1) conditional_c), hence
D_alpha(rho||sigma) = (alpha-1)^{-1} log sum_c p_c^alpha q_c^{1-alpha} exp((alpha-1) D_alpha(rho_c/p_c||sigma_c/q_c))
for states, and for Umegaki terms_c = p_c log(p_c/q_c) + p_c conditional_c, i.e. D(rho||sigma) = D(p||q) +
sum_c p_c D(rho_c/p_c||sigma_c/q_c).
sector_contributions¶
def sector_contributions(rho: Operator, sigma: Operator, alpha: Optional[Alpha]=None, kind: str='petz', trace: str='blunt', rtol: Optional[float]=None, supp_tol: Optional[float]=None) -> SectorContributions
Per-sector quantities of D_alpha (kind='petz' / 'sandwiched') or of the Umegaki relative entropy
(kind='umegaki', alpha ignored); see SectorContributions.
PetzRecovery¶
class PetzRecovery(sigma: Operator, channel, t: float=0.0, rtol: Optional[float]=None)
Petz recovery map of sigma and Phi, optionally rotated [JRSWW 2018]:
R^t(X) = sigma^{1/2 + it} Phi^*( Phi(sigma)^{-1/2 - it} X Phi(sigma)^{-1/2 + it} ) sigma^{1/2 - it},
i.e. R^t = U_{sigma,t} o R o U_{Phi(sigma),-t} with U_{s,t}(X) = s^{it} X s^{-it}; t = 0 is the Petz map.
Phi: Channel (M -> M) or InterSectorChannel (M_in -> M_out), batch 1 or B; Phi^* its adjoint for Tr_blunt.
sigma: positive operator (e.g. a DensityMatrix) on M_in, batch 1 or B. The inverse powers act on the numerical
support of Phi(sigma) (relative threshold rtol, module docstring). t != 0 needs a complex algebra.
R^t is completely positive with Tr R^t(X) = Tr(Pi_{Phi(sigma)} X) (= Tr X if supp X is in supp Phi(sigma)); for
trace-preserving Phi, R^t(Phi(sigma)) = sigma. Calling the map on an Operator X of M_out returns
a lazy Operator of M_in (flagged positive if X is; wrap it as a DensityMatrix when X is a state supported in
supp Phi(sigma)), differentiable in X, sigma and the Kraus operators. sigma and Phi(sigma) are decomposed once
per map, also for several rotation parameters (method rotated).
PetzRecovery.factors¶
def factors(self, t: Optional[float]=None) -> Tuple[torch.Tensor, torch.Tensor]
(sigma^{1/2 + it}, Phi(sigma)^{-1/2 - it}) as padded tensors (B, C, k, k), computed once per t.
PetzRecovery.rotated¶
def rotated(self, X: Operator, t: float) -> Operator
R^t(X) for another rotation parameter t (same decompositions).
petz_recovery¶
def petz_recovery(sigma: Operator, channel, X: Operator, t: float=0.0, rtol: Optional[float]=None) -> Operator
R^t_{sigma,Phi}(X) (PetzRecovery; t = 0: the Petz recovery map).
universal_recovery_nodes¶
def universal_recovery_nodes(step: float=0.125, cutoff: float=12.0) -> Tuple[List[float], List[float]]
Nodes t_j and weights of the trapezoidal rule for int beta_0(t) g(t) dt on [-cutoff, cutoff] (weights
renormalised to sum 1; beta_0(t) = pi/2 (cosh(pi t) + 1)^{-1}).
universal_recovery¶
def universal_recovery(sigma: Operator, channel, X: Operator, step: float=0.125, cutoff: float=12.0, rtol: Optional[float]=None) -> Operator
Universal recovery map of [JRSWW 2018], R(X) = int dt beta_0(t) R^{t/2}(X), beta_0(t) = pi/2 (cosh(pi t) + 1)^{-1}
(a probability density), with R^t the rotated Petz maps of PetzRecovery. The integral is computed by the
trapezoidal rule with the given step on [-cutoff, cutoff] (universal_recovery_nodes): the integrand is analytic in
the strip |Im t| < 1, so the discretisation error decays like exp(-2 pi / step) times a factor that grows with the
condition of sigma and Phi(sigma), and the truncation error is about 2 exp(-pi cutoff) (~1e-16 with the
defaults). sigma and Phi(sigma) are decomposed once. Complex algebras only. Returns a lazy Operator of M_in
(positive if X is).
lp_norm¶
def lp_norm(A: Operator, p: float, trace: str='blunt', weights: Optional[torch.Tensor]=None) -> torch.Tensor
Noncommutative L^p norm ||A||_{p,w} = (sum_c w_c Tr|A_c|^p)^{1/p} for p in [1, inf), ||A||_inf = max_c ||A_c||_op
over the sectors with w_c > 0; shape (batch,). The weights are those of `trace` ('blunt': 1, 'norm': 1/k_c,
'tau_vN': 1/(C k_c)) or explicit non-negative `weights` of shape (C,) or (batch, C). Singular values from the
Hermitian eigensolver if A is flagged self-adjoint, otherwise from the batched SVD; differentiable (at
non-degenerate singular values).
kosaki_norm¶
def kosaki_norm(x: Operator, rho: Operator, p: float, eta: float=0.5, trace: str='blunt', rtol: Optional[float]=None) -> torch.Tensor
State-weighted (Kosaki) L^p norm ||x||_{p,rho,eta} = ||d^{eta/p} x d^{(1-eta)/p}||_{p,trace} of x in M, for the
density d of the state with respect to `trace` (rho itself for trace='blunt'); eta = 1/2: symmetric embedding
(Beigi 2013), eta = 1: left, eta = 0: right. Factors with exponent 0 are omitted (identity). The value does not
depend on the trace: the central weights cancel (module docstring). For a faithful state, ||x||_{p} <= ||x||_inf
and the norms increase with p. Shape (batch,); differentiable in x and rho.