Commit 3facea88 authored by kerry-he's avatar kerry-he
Browse files

Added tests for Renyi entropies

parent 5da67fad
Loading
Loading
Loading
Loading
Loading
+2 −2
Original line number Diff line number Diff line
@@ -186,7 +186,7 @@ class BaseTrRenyiEntrEpiConstraint(Constraint):
        assert isinstance(divergence.X, required_type)
        assert isinstance(divergence.Y, required_type)

        super(BaseRenyiEntrConstraint, self).__init__(divergence._typeStr)
        super(BaseTrRenyiEntrEpiConstraint, self).__init__(divergence._typeStr)

    def _required_type(self):
        from ..expressions import AffineExpression
@@ -277,7 +277,7 @@ class TrSandRenyiEntrEpiConstraint(BaseTrRenyiEntrEpiConstraint):
        return 1 <= alpha and alpha <= 2


class ComplexTrSandRenyiEntrEpiConstraint(TrRenyiEntrEpiConstraint):
class ComplexTrSandRenyiEntrEpiConstraint(TrSandRenyiEntrEpiConstraint):
    """Upper bound of complex convex trace function used to define sandwiched  
    Renyi entropies.
    """

tests/ptest_renyi.py

0 → 100644
+138 −0
Original line number Diff line number Diff line
# ------------------------------------------------------------------------------
# Copyright (C) 2024 Kerry He
#
# This file is part of PICOS Testbench.
#
# PICOS Testbench is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# PICOS Testbench is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along with
# this program.  If not, see <http://www.gnu.org/licenses/>.
# ------------------------------------------------------------------------------

"""Test quantum relative entropy programs."""

import cvxopt
import math
import numpy as np

import picos

from .ptest import ProductionTestCase


def mpower(A, p):
    D, U = np.linalg.eigh(A)
    return U @ np.diag(np.power(D, p)) @ U.conj().T

class RMI_REAL(ProductionTestCase):
    """Renyi mutual information."""

    def setUp(self):  # noqa
        np.random.seed(42)

        # Primal problem.
        self.P = picos.Problem()
        self.X = picos.SymmetricVariable("X", 4)

        A = np.random.randn(16, 16)
        A = A @ A.T
        self.A = A / np.trace(A)
        self.tr2_A = picos.partial_trace(A, 1, (4, 4))

    def _opt_renyi_mutual_information(self, alpha):
        A, tr2_A = self.A, self.tr2_A
        temp = mpower(tr2_A @ picos.I(4), 1 - alpha) @ mpower(A, alpha)
        temp = mpower(picos.partial_trace(temp, 0, (4, 4)), 1 / alpha)
        return temp / np.trace(temp)
    
    def _opt_sand_renyi_mutual_information(self, Xstar, alpha):
        A, tr2_A = self.A, self.tr2_A
        temp = mpower(tr2_A @ Xstar, (1 - alpha) / (2 * alpha))
        temp = mpower(temp @ A @ temp, alpha)
        temp = picos.partial_trace(temp, 0, (4, 4))
        return temp / np.trace(temp)

    def testRenyi(self):
        P, X, A, tr2_A = self.P, self.X, self.A, self.tr2_A

        alpha = 0.5

        P.set_objective("min", picos.renyientr(A, tr2_A @ X, alpha))
        P.add_constraint(picos.trace(X) == 1)

        Xstar = self._opt_renyi_mutual_information(alpha)

        self.primalSolve(self.P)
        self.expectVariable(self.X, cvxopt.matrix(Xstar))

    def _test_trrenyi(self, alpha, direction):
        P, X, A, tr2_A = self.P, self.X, self.A, self.tr2_A

        P.set_objective(direction, picos.trrenyientr(A, tr2_A @ X, alpha))
        P.add_constraint(picos.trace(X) == 1)

        Xstar = self._opt_renyi_mutual_information(alpha)

        self.primalSolve(self.P)
        self.expectVariable(self.X, cvxopt.matrix(Xstar))

    def testTraceRenyi1(self):
        self._test_trrenyi(-0.5, "min")

    def testTraceRenyi2(self):
        self._test_trrenyi(0.5, "max")

    def testTraceRenyi3(self):
        self._test_trrenyi(1.5, "min")

    def testSandwichedRenyi(self):
        P, X, A, tr2_A = self.P, self.X, self.A, self.tr2_A

        alpha = 0.5

        P.set_objective("min", picos.sandrenyientr(A, tr2_A @ X, alpha))
        P.add_constraint(picos.trace(X) == 1)

        self.primalSolve(self.P)

        RHS = self._opt_sand_renyi_mutual_information(self.X, alpha)
        self.expectVariable(self.X, RHS.value)

    def _test_trsandrenyi(self, alpha, direction):
        P, X, A, tr2_A = self.P, self.X, self.A, self.tr2_A

        P.set_objective(direction, picos.trsandrenyientr(A, tr2_A @ X, alpha))
        P.add_constraint(picos.trace(X) == 1)

        self.primalSolve(self.P)

        RHS = self._opt_sand_renyi_mutual_information(self.X, alpha)
        self.expectVariable(self.X, RHS.value)

    def testTraceSandwichedRenyi1(self):
        self._test_trsandrenyi(0.75, "max")

    def testTraceSandwichedRenyi2(self):
        self._test_trsandrenyi(1.5, "min")

class RMI_COMPLEX(RMI_REAL):
    """Renyi mutual information."""

    def setUp(self):  # noqa
        np.random.seed(42)

        # Primal problem.
        self.P = picos.Problem()
        self.X = picos.HermitianVariable("X", 4)

        A = np.random.randn(16, 16) + np.random.randn(16, 16) * 1j
        A = A @ A.conj().T
        self.A = A / np.trace(A)
        self.tr2_A = picos.partial_trace(A, 1, (4, 4))