staging > main: refactored numerical experiment as module

src/experiments/__init__.py

@@ -0,0 +1,6 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+"""
+Module containing experiments
+"""

src/experiments/choi_cholesky/__init__.py

@@ -0,0 +1,20 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+"""
+Module containing experiments for Choi--Cholesky decomposition
+"""
+
+# ----------------------------------------------------------------
+# IMPORTS
+# ----------------------------------------------------------------
+
+from .verify import *
+
+# ----------------------------------------------------------------
+# EXPORTS
+# ----------------------------------------------------------------
+
+__all__ = [
+    "verify_choi_cholesky",
+]

src/experiments/choi_cholesky/verify.py

@@ -0,0 +1,134 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+"""
+- author:      Raj Dahya
+- institute:   Universität Leipzig
+- department:  Fakultät für Mathematik und Informatik
+- created:     2025-12-27
+- updated:     2026-03-29
+- description:
+    Verifies the Choi-Cholesky decomposition of a CP/CPTP map numerically,
+    as presented in the paper _The Choi-Cholesky algorithm for completely positive maps_
+    (see <https://arxiv.org/abs/2603.19444>).
+
+    Verification occurs by running the experiment multiple times without failure.
+
+    Failure criteria:
+
+    - [critical] algorithm fails to produce positive entries for the diagonal operator (within tolerance)
+    - in final result L·L^* is not close to Choi matrix
+
+    NOTE: Due to the discontinuous nature of pseudo-inverses (cf Remark 3.3 in §3),
+    the algorithm can get close to a critical error under low numerical precision (e.g. Float32).
+    For this reason we use double precision (Float64) numbers.
+
+    NOTE: The option `map_kind = "CP"` establishes a proper confirmation.
+    For completeness we have included the option `"CPTP"`,
+    but this provides a circular test, as the generation of CPTP-maps here assumes results of paper in advance.
+
+"""
+
+# ----------------------------------------------------------------
+# IMPORTS
+# ----------------------------------------------------------------
+
+import logging
+
+import torch
+from tabulate import tabulate
+
+from ..._core.utils.linalg import *
+from ...algorithms.choi_cholesky import *
+from ...algorithms.constants import *
+from ...algorithms.generation import *
+
+# ----------------------------------------------------------------
+# SETTINGS, CONSTANTS
+# ----------------------------------------------------------------
+
+TOL = 1e-6
+# PRECISION = torch.float32 # works but to lower precision
+PRECISION = torch.float64
+# PRECISION = torch.double # same as f64
+
+# ----------------------------------------------------------------
+# EXECUTION
+# ----------------------------------------------------------------
+
+
+def verify_choi_cholesky(
+    *,
+    d1: int,
+    d2: int,
+    num_experiments: int,
+    map_kind: MAP_CATEGORY,
+    algorithm_choice: ALGORITHM_CHOICE,
+):
+    """
+    Experimental verification of Choi--Cholesky decomposition of CP-/CPTP-maps
+    """
+    # execute main method
+    tol = TOL
+    diffs: list[float] = []
+    for index in range(1, 1 + num_experiments):
+        logging.info(f"run experiment {index}")
+
+        """
+        Create a random CP/CPTP-map
+        """
+
+        match map_kind:
+            case "CPTP":
+                logging.info("generate Choi matrix C of a random CPTP map Φ")
+                C = random_cptp(d1, d2, precision=PRECISION)
+                diff = torch.norm(trace_2(C) - torch.eye(d1)).item()
+                assert diff < TOL, \
+                    f"generated map failed to be trace preserving within tolerance {tol:.4g}"  # fmt: skip
+                logging.info(f"generated map is trace preserving within tolerance: ‖tr₂(C) – I‖ = {diff:.4g} < {tol:.4g}")  # fmt: skip
+
+            case _:
+                logging.info("generate Choi matrix C of a random CP map Φ")
+                C = random_cp(d1, d2, precision=PRECISION)
+
+        """
+        Perform Choi-Cholesky decomposition
+        """
+
+        logging.info("compute Choi-Cholesky decomposition for Φ")
+        match algorithm_choice:
+            case "CHOI-CHOLESKY-MODIFIED":
+                L = algorithm_choi_cholesky_modified(C, tol=TOL)
+
+            # case "CHOI-CHOLESKY":
+            case _:
+                L = algorithm_choi_cholesky(C, tol=TOL)
+
+        """
+        Verify validity of Choi-Cholesky decomposition
+        """
+
+        C_recovered = tensor_multiply(L, tensor_conj(L))
+        diff = tensor_diff(C, C_recovered, rel=True)
+        diffs.append(diff)
+
+    """
+    Perform statistical evaluation of success rate
+    """
+
+    diffs = torch.asarray(diffs)
+    levels = [0.25, 0.5, 0.75, 0.975]
+    quantiles = torch.quantile(diffs, q=torch.asarray(levels)).tolist()
+    table = tabulate(
+        zip(levels, quantiles),
+        headers=("quantile", "Diff"),
+        tablefmt="rst",
+        colalign=("right", "right"),
+        floatfmt=(".2%", ".3e"),
+    )
+    logging.info(f"quantiles for rel. differences ‖L·L^* – C(Φ)‖/‖C(Φ)‖:\n{table}")
+    diff = diffs[-1]  # get 97.5% quantile of rel differences
+    assert diff < 0.5*TOL, \
+        f"Algorithm failed to produce decomposition that recovers original Choi matrix within tolerance {TOL:.4g}."  # fmt: skip
+    logging.info(f"Algorithm succeeded in producing decomposition that recovers original Choi matrix within tolerance {TOL:.4g}.")  # fmt: skip
+    return