feat(cones): composable cones for tsallis, conditional entropy, and quadrature/integral relative entropy (#1868)

Add composable CVXPY cone builders for the remaining value functions: tsallis_entropy (hypo), tsallis_relative_entropy (epi), quantum_conditional_entropy (hypo), relative_entropy_quadrature (epi), and integral_relative_entropy (lower + upper). The tsallis and conditional-entropy builders compose existing cones; the quadrature/integral builders factor the SDP constraints out of the former internal solve, reproducing it exactly. cones stays a strict leaf (the conditional-entropy builder uses cvxpy.partial_trace, no matrix_ops import).

BREAKING CHANGE: relative_entropy_quadrature no longer accepts the m and k parameters; they moved to relative_entropy_quadrature_epi_cone.

Closes #1846.
Closes #1867.

Author: simsaidan

toqito/channel_metrics/tests/test_channel_relative_entropy.py

@@ -38,7 +38,7 @@ def failing_eigh(*args, **kwargs):
         raise LinAlgError("singular pencil")
 
     monkeypatch.setattr(
-        "toqito.state_props.integral_relative_entropy._generalized_eigenvalues",
+        "toqito.cones._integral_relative_entropy_helpers._generalized_eigenvalues",
         failing_eigh,
     )
 

toqito/cones/__init__.py

@@ -7,11 +7,27 @@
 
 from toqito.cones.geometric_mean_epi_cone import geometric_mean_epi_cone
 from toqito.cones.geometric_mean_hypo_cone import geometric_mean_hypo_cone
+from toqito.cones.integral_relative_entropy_lower_cone import (
+    integral_relative_entropy_lower_cone,
+)
+from toqito.cones.integral_relative_entropy_upper_cone import (
+    integral_relative_entropy_upper_cone,
+)
 from toqito.cones.lieb_ando_epi_cone import lieb_ando_epi_cone
 from toqito.cones.lieb_ando_hypo_cone import lieb_ando_hypo_cone
 from toqito.cones.ln_quantum_entropy_hypo_cone import ln_quantum_entropy_hypo_cone
 from toqito.cones.operator_relative_entropy_epi_cone import operator_relative_entropy_epi_cone
+from toqito.cones.quantum_conditional_entropy_hypo_cone import (
+    quantum_conditional_entropy_hypo_cone,
+)
 from toqito.cones.quantum_relative_entropy_epi_cone import quantum_relative_entropy_epi_cone
+from toqito.cones.relative_entropy_quadrature_epi_cone import (
+    relative_entropy_quadrature_epi_cone,
+)
 from toqito.cones.trace_matrix_log_hypo_cone import trace_matrix_log_hypo_cone
 from toqito.cones.trace_matrix_power_epi_cone import trace_matrix_power_epi_cone
 from toqito.cones.trace_matrix_power_hypo_cone import trace_matrix_power_hypo_cone
+from toqito.cones.tsallis_entropy_hypo_cone import tsallis_entropy_hypo_cone
+from toqito.cones.tsallis_relative_entropy_epi_cone import (
+    tsallis_relative_entropy_epi_cone,
+)

toqito/cones/_integral_relative_entropy_helpers.py

@@ -0,0 +1,95 @@
+"""Numeric helpers for integral relative-entropy sandwich grids."""
+
+import cvxpy
+import numpy as np
+from scipy.linalg import LinAlgError, eig, eigh
+
+
+def _generalized_eigenvalues(a: np.ndarray, b: np.ndarray) -> np.ndarray:
+    """Return real generalized eigenvalues for the pencil ``(a, b)``."""
+    try:
+        return np.real(eigh(a, b, check_finite=False)[0])
+    except LinAlgError:
+        return np.real(eig(a, b, left=False, right=False)[0])
+
+
+def _sandwich_parameters(rho: np.ndarray, sigma: np.ndarray) -> tuple[float, float]:
+    r"""Return sandwich bounds \(\mu\) and \(\lambda\) for PSD matrices \(X\) and \(Y\)."""
+    try:
+        w_xy = _generalized_eigenvalues(rho, sigma)
+        w_yx = _generalized_eigenvalues(sigma, rho)
+    except LinAlgError as exc:
+        raise ValueError("Failed to compute sandwich parameters from generalized eigenvalues.") from exc
+    finite_xy = w_xy[np.isfinite(w_xy)]
+    finite_yx = w_yx[np.isfinite(w_yx)]
+    if finite_xy.size == 0 or finite_yx.size == 0:
+        raise ValueError("Failed to compute sandwich parameters from generalized eigenvalues.")
+    lam = float(np.max(finite_xy))
+    mu = float(np.min(finite_yx))
+    return mu, lam
+
+
+def _make_grid(mu: float, lam: float, epsilon: float) -> np.ndarray:
+    r"""Make a grid of points for the integral representation of the relative entropy.
+
+    The first point of the grid is set to \(\mu\). For the k-th point \(t_k\), where
+    \(k \gt 1\), \(t_k = t_{k-1} + \sqrt{8 \epsilon t_{k-1}}\).
+
+    This formula yields \(O(\sqrt{\lambda/\epsilon})\) points in the grid.
+    """
+    grid = [mu]
+    curr = mu + np.sqrt(epsilon * mu * 8)
+    while curr < lam:
+        grid.append(curr)
+        curr = curr + np.sqrt(epsilon * 8 * curr)
+    return np.array(grid + [lam])
+
+
+def _make_delta(t: np.ndarray) -> np.ndarray:
+    r"""Make the delta coefficients for the integral representation."""
+    delta = np.zeros(len(t))
+    delta[0] = t[0] * ((1 + t[0] / (t[1] - t[0])) * np.log(t[1] / t[0]) - 1)
+    delta[-1] = t[-1] * (1 - (np.log(t[-1] / t[-2]) * t[-2] / (t[-1] - t[-2])))
+    for i in range(1, len(t) - 1):
+        delta[i] = t[i] * (
+            (1 + t[i] / (t[i + 1] - t[i])) * np.log(t[i + 1] / t[i])
+            - t[i - 1] * np.log(t[i] / t[i - 1]) / (t[i] - t[i - 1])
+        )
+    return delta
+
+
+def _make_gamma(t: np.ndarray) -> np.ndarray:
+    r"""Make the gamma coefficients for the integral representation."""
+    gamma = np.zeros(len(t))
+    gamma[0] = -1 * ((1 + t[0] / (t[1] - t[0])) * np.log(t[1] / t[0]) - 1)
+    gamma[-1] = -1 * (1 - t[-2] * np.log(t[-1] / t[-2]) / (t[-1] - t[-2]))
+    for i in range(1, len(t) - 1):
+        gamma[i] = -1 * (
+            (1 + t[i] / (t[i + 1] - t[i])) * np.log(t[i + 1] / t[i])
+            - np.log(t[i] / t[i - 1]) * t[i - 1] / (t[i] - t[i - 1])
+        )
+    return gamma
+
+
+def _integral_correction(lam: float) -> float:
+    return float(np.log(lam) + 1 - lam)
+
+
+def _require_valid_sandwich(mu: float, lam: float) -> None:
+    if mu <= 0 or lam <= mu:
+        raise ValueError(
+            "The integral representation requires 0 < mu < lambda. "
+            "This typically means the matrices are too close for the bound "
+            "(support of X may not be contained in support of Y)."
+        )
+
+
+def _numeric_pair_for_sandwich(
+    mat_x: cvxpy.Expression,
+    mat_y: cvxpy.Expression,
+) -> tuple[np.ndarray, np.ndarray]:
+    if mat_x.value is None or mat_y.value is None:
+        raise ValueError(
+            "Sandwich parameters require numeric `.value` on mat_x and mat_y, or pass mu and lam explicitly."
+        )
+    return np.asarray(mat_x.value), np.asarray(mat_y.value)

toqito/cones/integral_relative_entropy_lower_cone.py

@@ -0,0 +1,84 @@
+"""CVXPY constraints for the integral lower approximation of quantum relative entropy."""
+
+import cvxpy
+import numpy as np
+
+from toqito.cones._integral_relative_entropy_helpers import (
+    _integral_correction,
+    _make_grid,
+    _numeric_pair_for_sandwich,
+    _require_valid_sandwich,
+    _sandwich_parameters,
+)
+from toqito.cones._utils import _require_square_2d, _symmetric_like_variable
+
+
+def integral_relative_entropy_lower_cone(
+    mat_x: cvxpy.Expression,
+    mat_y: cvxpy.Expression,
+    t: cvxpy.Expression,
+    *,
+    epsilon_dec: float = 1e-2,
+    mu: float | None = None,
+    lam: float | None = None,
+    hermitian: bool = False,
+) -> list[cvxpy.Constraint]:
+    r"""Return CVXPY constraints for the integral lower approximation of \(D(X\|Y)\).
+
+    Discretizes the integral representation of quantum relative entropy
+    [@kossmann2024optimisingrelativeentropy] and enforces
+
+    \[
+        t \geqslant L_{\varepsilon}(X, Y),
+    \]
+
+    where \(L_{\varepsilon}\) is the lower SDP bound (with \(n \times n\) auxiliaries,
+    not the \(n^2\) Kronecker lift of ``quantum_relative_entropy_epi_cone``).
+
+    Sandwich endpoints \(\mu, \lambda\) are taken from ``mu`` / ``lam`` when provided;
+    otherwise they are computed from ``mat_x.value`` and ``mat_y.value``.
+
+    Args:
+        mat_x: A CVXPY expression for an ``n x n`` PSD matrix \(X\).
+        mat_y: A CVXPY expression for an ``n x n`` PSD matrix \(Y\).
+        t: A CVXPY scalar (or ``1 x 1``) epigraph variable.
+        epsilon_dec: Grid refinement parameter \(\varepsilon\).
+        mu: Optional sandwich lower endpoint.
+        lam: Optional sandwich upper endpoint.
+        hermitian: Whether the matrices are Hermitian or symmetric.
+
+    Raises:
+        ValueError: If ``mat_x`` or ``mat_y`` is not square 2D.
+        ValueError: If shapes differ, sandwich is degenerate, or numeric values
+            are missing when ``mu`` / ``lam`` are omitted.
+
+    Returns:
+        A list of CVXPY constraints.
+
+    """
+    _require_square_2d(mat_x, "mat_x")
+    _require_square_2d(mat_y, "mat_y")
+    if mat_x.shape != mat_y.shape:
+        raise ValueError("mat_x and mat_y must have the same shape")
+
+    if mu is None or lam is None:
+        if mu is not None or lam is not None:
+            raise ValueError("mu and lam must both be provided or both omitted")
+        x_val, y_val = _numeric_pair_for_sandwich(mat_x, mat_y)
+        mu, lam = _sandwich_parameters(x_val, y_val)
+    _require_valid_sandwich(mu, lam)
+
+    grid = _make_grid(mu, lam, epsilon_dec)
+    r = len(grid)
+    alpha = [float(np.log(grid[k] / grid[k + 1])) for k in range(r - 1)]
+    beta = [float(grid[k + 1] - grid[k]) for k in range(r - 1)]
+
+    n = int(mat_x.shape[0])
+    mu_vars = [_symmetric_like_variable(n, hermitian=hermitian) for _ in range(r - 1)]
+    constraints: list[cvxpy.Constraint] = [mu_vars[k] >> 0 for k in range(r - 1)]
+    constraints.extend(mu_vars[k] - alpha[k] * mat_x - beta[k] * mat_y >> 0 for k in range(r - 1))
+    bound = cvxpy.sum([cvxpy.trace(mu_vars[k]) for k in range(r - 1)]) + _integral_correction(lam)
+    if hermitian:
+        bound = cvxpy.real(bound)
+    constraints.append(t >= bound)
+    return constraints

toqito/cones/integral_relative_entropy_upper_cone.py

@@ -0,0 +1,85 @@
+"""CVXPY constraints for the integral upper approximation of quantum relative entropy."""
+
+import cvxpy
+
+from toqito.cones._integral_relative_entropy_helpers import (
+    _integral_correction,
+    _make_delta,
+    _make_gamma,
+    _make_grid,
+    _numeric_pair_for_sandwich,
+    _require_valid_sandwich,
+    _sandwich_parameters,
+)
+from toqito.cones._utils import _require_square_2d, _symmetric_like_variable
+
+
+def integral_relative_entropy_upper_cone(
+    mat_x: cvxpy.Expression,
+    mat_y: cvxpy.Expression,
+    t: cvxpy.Expression,
+    *,
+    epsilon_dec: float = 1e-2,
+    mu: float | None = None,
+    lam: float | None = None,
+    hermitian: bool = False,
+) -> list[cvxpy.Constraint]:
+    r"""Return CVXPY constraints for the integral upper approximation of \(D(X\|Y)\).
+
+    Discretizes the integral representation of quantum relative entropy
+    [@kossmann2024optimisingrelativeentropy] and enforces
+
+    \[
+        t \geqslant U_{\varepsilon}(X, Y),
+    \]
+
+    where \(U_{\varepsilon}\) is the upper SDP bound (with \(n \times n\) auxiliaries,
+    not the \(n^2\) Kronecker lift of ``quantum_relative_entropy_epi_cone``).
+
+    Sandwich endpoints \(\mu, \lambda\) are taken from ``mu`` / ``lam`` when provided;
+    otherwise they are computed from ``mat_x.value`` and ``mat_y.value``.
+
+    Args:
+        mat_x: A CVXPY expression for an ``n x n`` PSD matrix \(X\).
+        mat_y: A CVXPY expression for an ``n x n`` PSD matrix \(Y\).
+        t: A CVXPY scalar (or ``1 x 1``) epigraph variable.
+        epsilon_dec: Grid refinement parameter \(\varepsilon\).
+        mu: Optional sandwich lower endpoint.
+        lam: Optional sandwich upper endpoint.
+        hermitian: Whether the matrices are Hermitian or symmetric.
+
+    Raises:
+        ValueError: If ``mat_x`` or ``mat_y`` is not square 2D.
+        ValueError: If shapes differ, sandwich is degenerate, or numeric values
+            are missing when ``mu`` / ``lam`` are omitted.
+
+    Returns:
+        A list of CVXPY constraints.
+
+    """
+    _require_square_2d(mat_x, "mat_x")
+    _require_square_2d(mat_y, "mat_y")
+    if mat_x.shape != mat_y.shape:
+        raise ValueError("mat_x and mat_y must have the same shape")
+
+    if mu is None or lam is None:
+        if mu is not None or lam is not None:
+            raise ValueError("mu and lam must both be provided or both omitted")
+        x_val, y_val = _numeric_pair_for_sandwich(mat_x, mat_y)
+        mu, lam = _sandwich_parameters(x_val, y_val)
+    _require_valid_sandwich(mu, lam)
+
+    grid = _make_grid(mu, lam, epsilon_dec)
+    gamma = _make_gamma(grid)
+    delta = _make_delta(grid)
+    r = len(grid)
+
+    n = int(mat_x.shape[0])
+    nu_vars = [_symmetric_like_variable(n, hermitian=hermitian) for _ in range(r)]
+    constraints: list[cvxpy.Constraint] = [nu_vars[k] >> 0 for k in range(r)]
+    constraints.extend(nu_vars[k] - float(gamma[k]) * mat_x - float(delta[k]) * mat_y >> 0 for k in range(r))
+    bound = cvxpy.sum([cvxpy.trace(nu_vars[k]) for k in range(r)]) + _integral_correction(lam)
+    if hermitian:
+        bound = cvxpy.real(bound)
+    constraints.append(t >= bound)
+    return constraints