Quantum dynamic programming model

examples/transpiler/KAK_decomposition.py

@@ -0,0 +1,236 @@
+import numpy as np
+
+
+def orth_decomp_of_unitary(X):
+    """Decomposes unitary X = Ql @ exp(i Theta) @ Qr'
+    into orthogonal matrices Ql, Qr and angles Theta.
+
+    """
+
+    # one can also directly decompose XX^T and infer Ql from Qr. This might save about 3 eigh / qrs
+    # I will here try to implement something where I understand the proof, why it works.
+    # The benefit here is that we reliably work with real orthogonal matrices all the time.
+    Xr = X.real
+    Xi = X.imag
+    _, Ql = np.linalg.eigh(Xr @ Xi.transpose())
+    _, Qr = np.linalg.eigh(Xr.transpose() @ Xi)
+
+    Dr = Ql.transpose() @ Xr @ Qr
+    # if Xi injective this is diagonal, if not there is an arbitrary SO(dim ker(xi)) too much
+    # fixing the kernels, I don't know if there is a smarted way to do this!
+    if not np.allclose(Dr, np.diag(np.diag(Dr))):
+        Q, R = np.linalg.qr(Dr)
+        Dr = np.diag(R).copy()
+        Ql = Ql @ Q
+    else:
+        Dr = np.diag(Dr).copy()
+
+    Di = Ql.transpose() @ Xi @ Qr
+    # if Xr injective this is diagonal, if not there is an arbitrary SO(dim ker(Xr)) too much
+    if not np.allclose(Di, np.diag(np.diag(Di))):
+        Q, R = np.linalg.qr(Di.T)
+        Di = np.diag(R).copy()
+        Qr = Qr @ Q
+
+    else:
+        Di = np.diag(Di).copy()
+
+    # ensure Ql, Qr are in SO
+    if np.linalg.det(Ql) < 0:
+        Ql[:, 0] = -Ql[:, 0]
+        Dr[0] = -Dr[0]
+        Di[0] = -Di[0]
+
+    if np.linalg.det(Qr) < 0:
+        Qr[:, 0] = -Qr[:, 0]
+        Dr[0] = -Dr[0]
+        Di[0] = -Di[0]
+
+    Theta = np.angle(Dr + 1j * Di)
+
+    return Theta, Qr, Ql
+
+
+def unit_kronecker_rank_approx(X, verbose=True):
+    """Approximates a n^2 x m^2 matrix X  as A otimes B"""
+    # better check for square dimensions
+
+    n, m = np.sqrt(X.shape).astype(int)
+
+    Y = X.reshape(n, n, m, m).transpose((0, 2, 1, 3)).reshape(n * m, n * m)
+    U, S, Vh = np.linalg.svd(Y)
+
+    if verbose:
+        print("Truncating the spectrum", S)
+
+    U = U @ np.sqrt(np.diag(S))
+    Vh = np.sqrt(np.diag(S)) @ Vh
+    A = U[:, 0].reshape(n, m)
+    B = Vh[0, :].reshape(n, m)
+    return A, B
+
+
+MAGIC_BASIS = np.array(
+    [[1, 0, 0, 1j], [0, 1j, 1, 0], [0, 1j, -1, 0], [1, 0, 0, -1j]]
+) / np.sqrt(2)
+
+HADAMARD = np.array([[1, 1, -1, 1], [1, 1, 1, -1], [1, -1, -1, -1], [1, -1, 1, 1]]) / 2
+
+
+def KAK_decomp(U):
+    """Calculates the KAK decomposition of U in U(4).
+    U = np.kron(A0, A1) @ heisenberg_unitary(K) @ np.kron(B0.conj().T, B1.conj().T)
+    """
+
+    Theta, Qr, Ql = orth_decomp_of_unitary(MAGIC_BASIS.conj().T @ U @ MAGIC_BASIS)
+    A0, A1 = unit_kronecker_rank_approx(MAGIC_BASIS @ Ql @ MAGIC_BASIS.conj().T)
+    B0, B1 = unit_kronecker_rank_approx(MAGIC_BASIS @ Qr @ MAGIC_BASIS.conj().T)
+    K = HADAMARD.T @ Theta / 2
+
+    print("Theta", Theta)
+    print("K", K)
+
+    # Make A0, A1, B0, B1 special
+    def make_special(X):
+        return np.sqrt(np.linalg.det(X).conj()) * X
+
+    A0, A1, B0, B1 = map(make_special, (A0, A1, B0, B1))
+
+    return A0, A1, K, B0, B1
+
+
+# ----
+# Tests
+# ----
+
+
+# some convenient functions and constants
+def haar_rand_unitary(dim, real=False):
+    """Returns a random unitary dim x dim matrix drawn from the Haar measure.
+
+    Args:
+        dim (int): the dimension of the unitary
+        real (boolean): If true returns an orthogonal matrix.
+
+    Returns:
+        numpy.array: random unitary matrix
+    """
+
+    # Draw a random Gaussian
+    gaussianMatrix = np.random.randn(dim, dim)
+    if not real:  # Add a random imaginary part
+        gaussianMatrix = gaussianMatrix + 1j * np.random.randn(dim, dim)
+
+    # project to Haar random unitary
+    Q, R = np.linalg.qr(gaussianMatrix)
+    D = np.diag(R)
+    D = D / np.abs(D)
+    R = np.diag(D)
+    Q = Q.dot(R)
+
+    return Q
+
+
+def is_unitary(X, special=False, eps=1e-8):
+    is_special = np.abs(np.linalg.det(X) - 1) < 1e-8 if special else True
+    return is_special and np.linalg.norm(X @ X.conj().T - np.eye(X.shape[0])) < eps
+
+
+PAULI_X = np.array([[0, 1], [1, 0]])
+PAULI_Y = np.array([[0, -1j], [1j, 0]])
+PAULI_Z = np.array([[1, 0], [0, -1]])
+
+PXX = np.kron(PAULI_X, PAULI_X)
+PYY = np.kron(PAULI_Y, PAULI_Y)
+PZZ = np.kron(PAULI_Z, PAULI_Z)
+
+
+def heisenberg_unitary(k):
+    from scipy.linalg import expm
+
+    H = k[0] * np.eye(4) + k[1] * PXX + k[2] * PYY + k[3] * PZZ
+    return expm(1j * H)
+
+
+# test low_kronecker_rank_approx
+
+n, m = 2, 3
+
+A = np.random.randn(n, m)
+B = np.random.randn(n, m)
+
+C = np.kron(A, B)
+
+Arec, Brec = unit_kronecker_rank_approx(C)
+
+assert np.linalg.norm(C - np.kron(Arec, Brec)) < 1e-8
+
+
+# test orth_decomp_of_unitary
+dim = 10
+
+X = haar_rand_unitary(dim)
+
+
+def test_orth_decomp_of_unitary(X):
+
+    Theta, Qr, Ql = orth_decomp_of_unitary(X)
+    assert np.linalg.norm(Ql @ np.diag(np.exp(1j * Theta)) @ Qr.T - X) < 1e-8
+
+
+test_orth_decomp_of_unitary(X)
+
+# test KAK_decomp(U)
+
+
+def test_KAK_decomposition(U):
+    A0, A1, K, B0, B1 = KAK_decomp(U)
+
+    # Correctly decomposed
+    assert (
+        np.linalg.norm(
+            np.kron(A0, A1) @ heisenberg_unitary(K) @ np.kron(B0.conj().T, B1.conj().T)
+            - U
+        )
+        < 1e-8
+    )
+
+    # A0, A1, B0, B1 are in SU(2)
+    assert is_unitary(A0, special=True)
+    assert is_unitary(A1, special=True)
+    assert is_unitary(B0, special=True)
+    assert is_unitary(B1, special=True)
+
+
+# test KAK_decomp of random unitary
+
+test_KAK_decomposition(haar_rand_unitary(4))
+
+# test on Marek's example
+from scipy.linalg import expm
+
+H = 0.1 * (
+    PXX + PYY + 0.5 * (np.kron(PAULI_Z, np.eye(2)) + np.kron(np.eye(2), PAULI_Z))
+)
+U = expm(1j * H)
+UM = MAGIC_BASIS.conj().T @ U @ MAGIC_BASIS
+assert is_unitary(U)
+assert is_unitary(UM)
+
+test_orth_decomp_of_unitary(UM)
+
+test_KAK_decomposition(U)
+
+H2 = 0.1 * (
+    PXX + PZZ + PYY + (np.kron(PAULI_Z, np.eye(2)) + np.kron(np.eye(2), PAULI_Z))
+)
+U2 = expm(1j * H2)
+UM2 = MAGIC_BASIS.conj().T @ U2 @ MAGIC_BASIS
+test_orth_decomp_of_unitary(UM2)
+
+test_KAK_decomposition(expm(1j * H2))
+
+# test Heisenberg type Hamiltonian
+
+K = [0.1, 0.3, 0.6, 0.2]
+test_KAK_decomposition(heisenberg_unitary(K))

src/qibo/models/qdp/memory_usage_query.py

@@ -0,0 +1,122 @@
+import numpy as np
+import scipy
+
+from qibo import gates
+from qibo.models.qdp.quantum_dynamic_programming import (
+    MeasurementEmulation,
+    MeasurementReset,
+    SequentialInstruction,
+)
+from qibo.transpiler.unitary_decompositions import two_qubit_decomposition
+
+
+class DensityMatrixExponentiation(SequentialInstruction):
+    """
+    Subclass of AbstractQuantumDynamicProgramming for density matrix exponentiation,
+    where we attempt to instruct the work qubit to do an X rotation, using SWAP gate.
+
+    Args:
+        theta (float): Overall rotation angle.
+        N (int): Number of steps.
+        num_work_qubits (int): Number of work qubits.
+        num_instruction_qubits (int): Number of instruction qubits.
+        number_muq_per_call (int): Number of memory units per call.
+
+    Example:
+        import numpy as np
+        from qibo.models.qdp.dynamic_programming import DensityMatrixExponentiation
+        my_protocol = DensityMatrixExponentiation(theta=np.pi,N=3,num_work_qubits=1,num_instruction_qubits=3,number_muq_per_call=1)
+        my_protocol.memory_call_circuit(num_instruction_qubits_per_query=3)
+        print('DME, q0 is target qubit, q1,q2 and q3 are instruction qubit')
+        print(my_protocol.c.draw())
+        my_protocol.c.execute(nshots=1000).frequencies()
+    """
+
+    def __init__(
+        self, theta, N, num_work_qubits, num_instruction_qubits, number_muq_per_call
+    ):
+        super().__init__(
+            num_work_qubits, num_instruction_qubits, number_muq_per_call, circuit=None
+        )
+        self.theta = theta  # overall rotation angle
+        self.N = N  # number of steps
+        self.delta = theta / N  # small rotation angle
+        self.id_current_work_reg = self.list_id_work_reg[0]
+
+    def memory_usage_query_circuit(self):
+        """Defines the memory usage query circuit."""
+        delta_swap = scipy.linalg.expm(
+            -1j
+            * gates.SWAP(
+                self.id_current_work_reg, self.id_current_instruction_reg
+            ).matrix()
+            * self.delta
+        )
+        for decomposed_gate in two_qubit_decomposition(
+            self.id_current_work_reg,
+            self.id_current_instruction_reg,
+            unitary=delta_swap,
+        ):
+            self.c.add(decomposed_gate)
+
+    def instruction_qubits_initialization(self):
+        """Initializes the instruction qubits."""
+        for instruction_qubit in self.list_id_current_instruction_reg:
+            self.c.add(gates.X(instruction_qubit))
+
+
+class DME_reset(MeasurementReset):
+    """
+    Warning: Functional, but without a way to actually do reset.
+    DME using reset method.
+    """
+
+    def __init__(
+        self, theta, N, num_work_qubits, num_instruction_qubits, number_muq_per_call
+    ):
+        super().__init__(
+            num_work_qubits, num_instruction_qubits, number_muq_per_call, circuit=None
+        )
+        self.theta = theta  # overall rotation angle
+        self.N = N  # number of steps
+        self.delta = theta / N  # small rotation angle
+        self.id_current_work_reg = self.list_id_work_reg[0]
+
+    def current_register_reset(self):
+        """
+        Resets a single register.
+
+        Args:
+            register (int): The register index.
+            _c = self.c.copy()
+        """
+        # todo: find a way to do reset
+        # result = _c.execute(nshots=1).samples(binary=True)[0][0]
+        result = 1
+        if result == 1:
+            self.c.add(gates.X(self.id_current_instruction_reg))
+        elif result == 0:
+            pass
+        else:
+            print("Warning: qubit wasn't reset")
+
+    def memory_usage_query_circuit(self):
+        """Defines the memory usage query circuit."""
+        delta_SWAP = scipy.linalg.expm(
+            -1j
+            * gates.SWAP(
+                self.id_current_work_reg, self.id_current_instruction_reg
+            ).matrix()
+            * self.delta
+        )
+        for decomposed_gate in two_qubit_decomposition(
+            self.id_current_work_reg,
+            self.id_current_instruction_reg,
+            unitary=delta_SWAP,
+        ):
+            self.c.add(decomposed_gate)
+
+    def instruction_qubits_initialization(self):
+        """Initializes the instruction qubits."""
+        for instruction_qubit in self.list_id_current_instruction_reg:
+            self.c.add(gates.X(instruction_qubit))