Python scipy.linalg.fractional_matrix_power() Examples

def test_larger_abs_fractional_matrix_powers(self):
        np.random.seed(1234)
        for n in (2, 3, 5):
            for i in range(10):
                M = np.random.randn(n, n) + 1j * np.random.randn(n, n)
                M_one_fifth = fractional_matrix_power(M, 0.2)
                # Test the round trip.
                M_round_trip = np.linalg.matrix_power(M_one_fifth, 5)
                assert_allclose(M, M_round_trip)
                # Test a large abs fractional power.
                X = fractional_matrix_power(M, -5.4)
                Y = np.linalg.matrix_power(M_one_fifth, -27)
                assert_allclose(X, Y)
                # Test another large abs fractional power.
                X = fractional_matrix_power(M, 3.8)
                Y = np.linalg.matrix_power(M_one_fifth, 19)
                assert_allclose(X, Y) 

def test_random_matrices_and_powers(self):
        # Each independent iteration of this fuzz test picks random parameters.
        # It tries to hit some edge cases.
        np.random.seed(1234)
        nsamples = 20
        for i in range(nsamples):
            # Sample a matrix size and a random real power.
            n = random.randrange(1, 5)
            p = np.random.randn()

            # Sample a random real or complex matrix.
            matrix_scale = np.exp(random.randrange(-4, 5))
            A = np.random.randn(n, n)
            if random.choice((True, False)):
                A = A + 1j * np.random.randn(n, n)
            A = A * matrix_scale

            # Check a couple of analytically equivalent ways
            # to compute the fractional matrix power.
            # These can be compared because they both use the principal branch.
            A_power = fractional_matrix_power(A, p)
            A_logm, info = logm(A, disp=False)
            A_power_expm_logm = expm(A_logm * p)
            assert_allclose(A_power, A_power_expm_logm) 

def test_al_mohy_higham_2012_experiment_1(self):
        # Fractional powers of a tricky upper triangular matrix.
        A = _get_al_mohy_higham_2012_experiment_1()

        # Test remainder matrix power.
        A_funm_sqrt, info = funm(A, np.sqrt, disp=False)
        A_sqrtm, info = sqrtm(A, disp=False)
        A_rem_power = _matfuncs_inv_ssq._remainder_matrix_power(A, 0.5)
        A_power = fractional_matrix_power(A, 0.5)
        assert_array_equal(A_rem_power, A_power)
        assert_allclose(A_sqrtm, A_power)
        assert_allclose(A_sqrtm, A_funm_sqrt)

        # Test more fractional powers.
        for p in (1/2, 5/3):
            A_power = fractional_matrix_power(A, p)
            A_round_trip = fractional_matrix_power(A_power, 1/p)
            assert_allclose(A_round_trip, A, rtol=1e-2)
            assert_allclose(np.tril(A_round_trip, 1), np.tril(A, 1)) 

def test_type_conversion_mixed_sign_or_complex_spectrum(self):
        complex_dtype_chars = ('F', 'D', 'G')
        for matrix_as_list in (
                [[1, 0], [0, -1]],
                [[0, 1], [1, 0]],
                [[0, 1, 0], [0, 0, 1], [1, 0, 0]]):

            # check that the spectrum has the expected properties
            W = scipy.linalg.eigvals(matrix_as_list)
            assert_(any(w.imag or w.real < 0 for w in W))

            # Check various positive and negative powers
            # with absolute values bigger and smaller than 1.
            for p in (-2.4, -0.9, 0.2, 3.3):

                # check complex->complex
                A = np.array(matrix_as_list, dtype=complex)
                A_power = fractional_matrix_power(A, p)
                assert_(A_power.dtype.char in complex_dtype_chars)

                # check float->complex
                A = np.array(matrix_as_list, dtype=float)
                A_power = fractional_matrix_power(A, p)
                assert_(A_power.dtype.char in complex_dtype_chars) 

def test_singular(self):
        # Negative fractional powers do not work with singular matrices.
        for matrix_as_list in (
                [[0, 0], [0, 0]],
                [[1, 1], [1, 1]],
                [[1, 2], [3, 6]],
                [[0, 0, 0], [0, 1, 1], [0, -1, 1]]):

            # Check fractional powers both for float and for complex types.
            for newtype in (float, complex):
                A = np.array(matrix_as_list, dtype=newtype)
                for p in (-0.7, -0.9, -2.4, -1.3):
                    A_power = fractional_matrix_power(A, p)
                    assert_(np.isnan(A_power).all())
                for p in (0.2, 1.43):
                    A_power = fractional_matrix_power(A, p)
                    A_round_trip = fractional_matrix_power(A_power, 1/p)
                    assert_allclose(A_round_trip, A) 

def test_round_trip_random_complex(self):
        np.random.seed(1234)
        for p in range(1, 5):
            for n in range(1, 5):
                M_unscaled = np.random.randn(n, n) + 1j * np.random.randn(n, n)
                for scale in np.logspace(-4, 4, 9):
                    M = M_unscaled * scale
                    M_root = fractional_matrix_power(M, 1/p)
                    M_round_trip = np.linalg.matrix_power(M_root, p)
                    assert_allclose(M_round_trip, M) 

def test_round_trip_random_float(self):
        # This test is more annoying because it can hit the branch cut;
        # this happens when the matrix has an eigenvalue
        # with no imaginary component and with a real negative component,
        # and it means that the principal branch does not exist.
        np.random.seed(1234)
        for p in range(1, 5):
            for n in range(1, 5):
                M_unscaled = np.random.randn(n, n)
                for scale in np.logspace(-4, 4, 9):
                    M = M_unscaled * scale
                    M_root = fractional_matrix_power(M, 1/p)
                    M_round_trip = np.linalg.matrix_power(M_root, p)
                    assert_allclose(M_round_trip, M) 

def __pow__(self, t: float) -> 'Gate':
        """Return this gate raised to the given power."""
        # Note: This operation cannot be performed within the tensorflow or
        # torch backends in general. Subclasses of Gate may override
        # for special cases.
        N = self.qubit_nb
        matrix = asarray(self.vec.flatten())
        matrix = matpow(matrix, t)
        matrix = np.reshape(matrix, ([2]*(2*N)))
        return Gate(matrix, self.qubits)

    # TODO: Refactor functionality into QubitVector 

def quantum_chernoff_bound(rho: np.ndarray,
                           sigma: np.ndarray,
                           tol: float = 1000) -> Tuple[float, float]:
    r"""
    Computes the quantum Chernoff bound between rho and sigma.

    It is defined as

    .. math::

        ξ_{QCB}(\rho, \sigma) = - \log[ \min_{0\le s\le 1} tr(\rho^s \sigma^{1-s}) ]

    It is also common to study the non-logarithmic variety of the quantum Chernoff bound

    .. math::

        Q_{QCB}(\rho, \sigma) = \min_{0\le s\le 1} tr(\rho^s \sigma^{1-s})

    The quantum Chernoff bound has many nice properties, see [QCB]_. Importantly it is
    operationally important in the following context. Given n copies of rho or sigma the minimum
    error probability for discriminating rho from sigma is :math:`P_{e,min,n} ~ exp[-n ξ_{QCB}]`.

    .. [QCB] The Quantum Chernoff Bound.
          Audenaert et al.
          Phys. Rev. Lett. 98, 160501 (2007).
          https://dx.doi.org/10.1103/PhysRevLett.98.160501
          https://arxiv.org/abs/quant-ph/0610027

    :param rho: Is a dim by dim positive matrix with unit trace.
    :param sigma: Is a dim by dim positive matrix with unit trace.
    :param tol: Tolerance in machine epsilons for np.real_if_close.
    :return: The non-logarithmic quantum Chernoff bound and the s achieving the minimum.
    """

    def f(s):
        s = np.real_if_close(s)
        return np.trace(
            np.matmul(fractional_matrix_power(rho, s), fractional_matrix_power(sigma, 1 - s)))

    f_min = minimize_scalar(f, bounds=(0, 1), method='bounded')
    s_opt = np.real_if_close(f_min.x, tol)
    qcb = np.real_if_close(f_min.fun, tol)
    return qcb, s_opt 

def test_diamond_norm_distance():
    if int(os.getenv('SKIP_SCS', 0)) == 1:
        return pytest.skip('Having issues with SCS, skipping for now')

    # Test cases borrowed from qutip,
    # https://github.com/qutip/qutip/blob/master/qutip/tests/test_metrics.py
    # which were in turn generated using QuantumUtils for MATLAB
    # (https://goo.gl/oWXhO9) by Christopher Granade
    choi0 = kraus2choi(I_MAT)
    choi1 = kraus2choi(X_MAT)
    dnorm = dm.diamond_norm_distance(choi0, choi1)
    assert np.isclose(2.0, dnorm, rtol=0.01)

    turns_dnorm = [[1.000000e-03, 3.141591e-03],
                   [3.100000e-03, 9.738899e-03],
                   [1.000000e-02, 3.141463e-02],
                   [3.100000e-02, 9.735089e-02],
                   [1.000000e-01, 3.128689e-01],
                   [3.100000e-01, 9.358596e-01]]

    for turns, target in turns_dnorm:
        choi0 = kraus2choi(X_MAT)
        choi1 = kraus2choi(matpow(X_MAT, 1 + turns))
        dnorm = dm.diamond_norm_distance(choi0, choi1)
        assert np.isclose(target, dnorm, rtol=0.01)

    hadamard_mixtures = [[1.000000e-03, 2.000000e-03],
                         [3.100000e-03, 6.200000e-03],
                         [1.000000e-02, 2.000000e-02],
                         [3.100000e-02, 6.200000e-02],
                         [1.000000e-01, 2.000000e-01],
                         [3.100000e-01, 6.200000e-01]]

    for p, target in hadamard_mixtures:
        chan0 = kraus2superop(I_MAT) * (1 - p) + kraus2superop(H_MAT) * p
        chan1 = kraus2superop(I_MAT)

        choi0 = superop2choi(chan0)
        choi1 = superop2choi(chan1)
        dnorm = dm.diamond_norm_distance(choi0, choi1)
        assert np.isclose(dnorm, target, rtol=0.01)

    choi0 = kraus2choi(I_MAT)
    choi1 = kraus2choi(matpow(Y_MAT, 0.5))
    dnorm = dm.diamond_norm_distance(choi0, choi1)
    assert np.isclose(dnorm, np.sqrt(2), rtol=0.01) 

def test_watrous_bounds():
    # Test cases borrowed from qutip,
    # https://github.com/qutip/qutip/blob/master/qutip/tests/test_metrics.py
    # which were in turn generated using QuantumUtils for MATLAB
    # (https://goo.gl/oWXhO9) by Christopher Granade

    choi0 = kraus2choi(I_MAT)
    choi1 = kraus2choi(X_MAT)
    wbounds = dm.watrous_bounds(choi0-choi1)
    assert wbounds[0]/2 <= 2.0 or np.isclose(wbounds[0]/2, 2.0, rtol=1e-2)
    assert wbounds[1]/2 >= 2.0 or np.isclose(wbounds[1]/2, 2.0, rtol=1e-2)

    turns_dnorm = [[1.000000e-03, 3.141591e-03],
                   [3.100000e-03, 9.738899e-03],
                   [1.000000e-02, 3.141463e-02],
                   [3.100000e-02, 9.735089e-02],
                   [1.000000e-01, 3.128689e-01],
                   [3.100000e-01, 9.358596e-01]]

    for turns, target in turns_dnorm:
        choi0 = kraus2choi(X_MAT)
        choi1 = kraus2choi(matpow(X_MAT, 1 + turns))
        wbounds = dm.watrous_bounds(choi0-choi1)
        assert wbounds[0]/2 <= target or np.isclose(wbounds[0]/2, target, rtol=1e-2)
        assert wbounds[1]/2 >= target or np.isclose(wbounds[1]/2, target, rtol=1e-2)
                          
    hadamard_mixtures = [[1.000000e-03, 2.000000e-03],
                         [3.100000e-03, 6.200000e-03],
                         [1.000000e-02, 2.000000e-02],
                         [3.100000e-02, 6.200000e-02],
                         [1.000000e-01, 2.000000e-01],
                         [3.100000e-01, 6.200000e-01]]

    for p, target in hadamard_mixtures:
        chan0 = kraus2superop(I_MAT) * (1 - p) + kraus2superop(H_MAT) * p
        chan1 = kraus2superop(I_MAT)

        choi0 = superop2choi(chan0)
        choi1 = superop2choi(chan1)
        wbounds = dm.watrous_bounds(choi0-choi1)
        assert wbounds[0]/2 <= target or np.isclose(wbounds[0]/2, target, rtol=1e-2)
        assert wbounds[1]/2 >= target or np.isclose(wbounds[1]/2, target, rtol=1e-2)

    choi0 = kraus2choi(I_MAT)
    choi1 = kraus2choi(matpow(Y_MAT, 0.5))
    wbounds = dm.watrous_bounds(choi0-choi1)
    assert wbounds[0]/2 <= np.sqrt(2) or np.isclose(wbounds[0]/2, np.sqrt(2), rtol=1e-2)
    assert wbounds[1]/2 >= np.sqrt(2) or np.isclose(wbounds[0]/2, np.sqrt(2), rtol=1e-2) 

def one_electron_integrals(self):
        """Routine to set up and compute one-electron integrals"""
        N = self.nbasis
        # core integrals
        self.S = np.zeros((N,N)) 
        self.V = np.zeros((N,N)) 
        self.T = np.zeros((N,N)) 
        # dipole integrals
        self.M = np.zeros((3,N,N)) 
        self.mu = np.zeros(3,dtype='complex') 
 
        # angular momentum
        self.L = np.zeros((3,N,N)) 

        self.nuc_energy = 0.0
        # Get one electron integrals
        #print "One-electron integrals"

        for i in (range(N)):
            for j in range(i+1):
                self.S[i,j] = self.S[j,i] \
                    = S(self.bfs[i],self.bfs[j])
                self.T[i,j] = self.T[j,i] \
                    = T(self.bfs[i],self.bfs[j])
                self.M[0,i,j] = self.M[0,j,i] \
                    = Mu(self.bfs[i],self.bfs[j],self.center_of_charge,'x')
                self.M[1,i,j] = self.M[1,j,i] \
                    = Mu(self.bfs[i],self.bfs[j],self.center_of_charge,'y')
                self.M[2,i,j] = self.M[2,j,i] \
                    = Mu(self.bfs[i],self.bfs[j],self.center_of_charge,'z')
                for atom in self.atoms:
                    self.V[i,j] += -atom.charge*V(self.bfs[i],self.bfs[j],atom.origin)
                self.V[j,i] = self.V[i,j]

                # RxDel is antisymmetric
                self.L[0,i,j] \
                    = RxDel(self.bfs[i],self.bfs[j],self.center_of_charge,'x')
                self.L[1,i,j] \
                    = RxDel(self.bfs[i],self.bfs[j],self.center_of_charge,'y')
                self.L[2,i,j] \
                    = RxDel(self.bfs[i],self.bfs[j],self.center_of_charge,'z')
                self.L[:,j,i] = -1*self.L[:,i,j] 

        # Compute nuclear repulsion energy 
        for pair in itertools.combinations(self.atoms,2):
            self.nuc_energy += pair[0].charge*pair[1].charge \
                              / np.linalg.norm(pair[0].origin - pair[1].origin)
           
        # Preparing for SCF
        self.Core       = self.T + self.V
        self.X          = mat_pow(self.S,-0.5)
        self.U          = mat_pow(self.S,0.5)