qiskit/quantum_info/synthesis/two_qubit_kak.py
def two_qubit_kak(unitary):
    """Decompose a two-qubit gate over SU(2)+CNOT using the KAK decomposition.

    Args:
        unitary (Operator): a 4x4 unitary operator to decompose.

    Returns:
        QuantumCircuit: a circuit implementing the unitary over SU(2)+CNOT

    Raises:
        QiskitError: input not a unitary, or error in KAK decomposition.
    """
    if hasattr(unitary, 'to_operator'):
        # If input is a BaseOperator subclass this attempts to convert
        # the object to an Operator so that we can extract the underlying
        # numpy matrix from `Operator.data`.
        unitary = unitary.to_operator().data
    if hasattr(unitary, 'to_matrix'):
        # If input is Gate subclass or some other class object that has
        # a to_matrix method this will call that method.
        unitary = unitary.to_matrix()
    # Convert to numpy array incase not already an array
    unitary_matrix = np.array(unitary, dtype=complex)
    # Check input is a 2-qubit unitary
    if unitary_matrix.shape != (4, 4):
        raise QiskitError("two_qubit_kak: Expected 4x4 matrix")
    if not is_unitary_matrix(unitary_matrix):
        raise QiskitError("Input matrix is not unitary.")
    phase = la.det(unitary_matrix)**(-1.0/4.0)
    # Make it in SU(4), correct phase at the end
    U = phase * unitary_matrix
    # B changes to the Bell basis
    B = (1.0/math.sqrt(2)) * np.array([[1, 1j, 0, 0],
                                       [0, 0, 1j, 1],
                                       [0, 0, 1j, -1],
                                       [1, -1j, 0, 0]], dtype=complex)

    # We also need B.conj().T below
    Bdag = B.conj().T
    # U' = Bdag . U . B
    Uprime = Bdag.dot(U.dot(B))
    # M^2 = trans(U') . U'
    M2 = Uprime.T.dot(Uprime)

    # Diagonalize M2
    # Must use diagonalization routine which finds a real orthogonal matrix P
    # when M2 is real.
    D, P = la.eig(M2)
    D = np.diag(D)
    # If det(P) == -1 then in O(4), apply a swap to make P in SO(4)
    if abs(la.det(P)+1) < 1e-5:
        swap = np.array([[1, 0, 0, 0],
                         [0, 0, 1, 0],
                         [0, 1, 0, 0],
                         [0, 0, 0, 1]], dtype=complex)
        P = P.dot(swap)
        D = swap.dot(D.dot(swap))

    Q = np.sqrt(D)  # array from elementwise sqrt
    # Want to take square root so that Q has determinant 1
    if abs(la.det(Q)+1) < 1e-5:
        Q[0, 0] = -Q[0, 0]

    # Q^-1*P.T = P' -> QP' = P.T (solve for P' using Ax=b)
    Pprime = la.solve(Q, P.T)
    # K' now just U' * P * P'
    Kprime = Uprime.dot(P.dot(Pprime))

    K1 = B.dot(Kprime.dot(P.dot(Bdag)))
    A = B.dot(Q.dot(Bdag))
    K2 = B.dot(P.T.dot(Bdag))
    # KAK = K1 * A * K2
    KAK = K1.dot(A.dot(K2))

    # Verify decomp matches input unitary.
    if la.norm(KAK - U) > 1e-6:
        raise QiskitError("two_qubit_kak: KAK decomposition " +
                          "does not return input unitary.")

    # Compute parameters alpha, beta, gamma so that
    # A = exp(i * (alpha * XX + beta * YY + gamma * ZZ))
    xx = np.array([[0, 0, 0, 1],
                   [0, 0, 1, 0],
                   [0, 1, 0, 0],
                   [1, 0, 0, 0]], dtype=complex)

    yy = np.array([[0, 0, 0, -1],
                   [0, 0, 1, 0],
                   [0, 1, 0, 0],
                   [-1, 0, 0, 0]], dtype=complex)

    zz = np.array([[1, 0, 0, 0],
                   [0, -1, 0, 0],
                   [0, 0, -1, 0],
                   [0, 0, 0, 1]], dtype=complex)

    A_real_tr = A.real.trace()
    alpha = math.atan2(A.dot(xx).imag.trace(), A_real_tr)
    beta = math.atan2(A.dot(yy).imag.trace(), A_real_tr)
    gamma = math.atan2(A.dot(zz).imag.trace(), A_real_tr)

    # K1 = kron(U1, U2) and K2 = kron(V1, V2)
    # Find the matrices U1, U2, V1, V2

    # Find a block in K1 where U1_ij * [U2] is not zero
    L = K1[0:2, 0:2]
    if la.norm(L) < 1e-9:
        L = K1[0:2, 2:4]
        if la.norm(L) < 1e-9:
            L = K1[2:4, 2:4]
    # Remove the U1_ij prefactor
    Q = L.dot(L.conj().T)
    U2 = L / math.sqrt(Q[0, 0].real)

    # Now grab U1 given we know U2
    R = K1.dot(np.kron(np.identity(2), U2.conj().T))
    U1 = np.zeros((2, 2), dtype=complex)
    U1[0, 0] = R[0, 0]
    U1[0, 1] = R[0, 2]
    U1[1, 0] = R[2, 0]
    U1[1, 1] = R[2, 2]

    # Repeat K1 routine for K2
    L = K2[0:2, 0:2]
    if la.norm(L) < 1e-9:
        L = K2[0:2, 2:4]
        if la.norm(L) < 1e-9:
            L = K2[2:4, 2:4]
    Q = np.dot(L, np.transpose(L.conjugate()))
    V2 = L / np.sqrt(Q[0, 0])
    R = np.dot(K2, np.kron(np.identity(2), np.transpose(V2.conjugate())))

    V1 = np.zeros_like(U1)
    V1[0, 0] = R[0, 0]
    V1[0, 1] = R[0, 2]
    V1[1, 0] = R[2, 0]
    V1[1, 1] = R[2, 2]

    if la.norm(np.kron(U1, U2) - K1) > 1e-4:
        raise QiskitError("two_qubit_kak: K1 != U1 x U2")
    if la.norm(np.kron(V1, V2) - K2) > 1e-4:
        raise QiskitError("two_qubit_kak: K2 != V1 x V2")

    test = la.expm(1j*(alpha * xx + beta * yy + gamma * zz))
    if la.norm(A - test) > 1e-4:
        raise QiskitError("two_qubit_kak: " +
                          "Matrix A does not match xx,yy,zz decomposition.")

    # Circuit that implements K1 * A * K2 (up to phase), using
    # Vatan and Williams Fig. 6 of quant-ph/0308006v3
    # Include prefix and suffix single-qubit gates into U2, V1 respectively.

    V2 = np.array([[np.exp(1j*np.pi/4), 0],
                   [0, np.exp(-1j*np.pi/4)]], dtype=complex).dot(V2)
    U1 = U1.dot(np.array([[np.exp(-1j*np.pi/4), 0],
                          [0, np.exp(1j*np.pi/4)]], dtype=complex))

    # Corrects global phase: exp(ipi/4)*phase'
    U1 = U1.dot(np.array([[np.exp(1j*np.pi/4), 0],
                          [0, np.exp(1j*np.pi/4)]], dtype=complex))
    U1 = phase.conjugate() * U1

    # Test
    g1 = np.kron(V1, V2)
    g2 = np.array([[1, 0, 0, 0],
                   [0, 0, 0, 1],
                   [0, 0, 1, 0],
                   [0, 1, 0, 0]], dtype=complex)

    theta = 2*gamma - np.pi/2

    Ztheta = np.array([[np.exp(1j*theta/2), 0],
                       [0, np.exp(-1j*theta/2)]], dtype=complex)

    kappa = np.pi/2 - 2*alpha
    Ykappa = np.array([[math.cos(kappa/2), math.sin(kappa/2)],
                       [-math.sin(kappa/2), math.cos(kappa/2)]], dtype=complex)
    g3 = np.kron(Ztheta, Ykappa)
    g4 = np.array([[1, 0, 0, 0],
                   [0, 1, 0, 0],
                   [0, 0, 0, 1],
                   [0, 0, 1, 0]], dtype=complex)

    zeta = 2*beta - np.pi/2
    Yzeta = np.array([[math.cos(zeta/2), math.sin(zeta/2)],
                      [-math.sin(zeta/2), math.cos(zeta/2)]], dtype=complex)
    g5 = np.kron(np.identity(2), Yzeta)
    g6 = g2
    g7 = np.kron(U1, U2)

    V = g2.dot(g1)
    V = g3.dot(V)
    V = g4.dot(V)
    V = g5.dot(V)
    V = g6.dot(V)
    V = g7.dot(V)

    if la.norm(V - U*phase.conjugate()) > 1e-6:
        raise QiskitError("two_qubit_kak: " +
                          "sequence incorrect, unknown error")

    v1_param = euler_angles_1q(V1)
    v2_param = euler_angles_1q(V2)
    u1_param = euler_angles_1q(U1)
    u2_param = euler_angles_1q(U2)

    v1_gate = U3Gate(v1_param[0], v1_param[1], v1_param[2])
    v2_gate = U3Gate(v2_param[0], v2_param[1], v2_param[2])
    u1_gate = U3Gate(u1_param[0], u1_param[1], u1_param[2])
    u2_gate = U3Gate(u2_param[0], u2_param[1], u2_param[2])

    q = QuantumRegister(2)
    return_circuit = QuantumCircuit(q)

    return_circuit.append(v1_gate, [q[1]])

    return_circuit.append(v2_gate, [q[0]])

    return_circuit.append(CnotGate(), [q[0], q[1]])

    gate = U3Gate(0.0, 0.0, -2.0*gamma + np.pi/2.0)
    return_circuit.append(gate, [q[1]])

    gate = U3Gate(-np.pi/2.0 + 2.0*alpha, 0.0, 0.0)
    return_circuit.append(gate, [q[0]])

    return_circuit.append(CnotGate(), [q[1], q[0]])

    gate = U3Gate(-2.0*beta + np.pi/2.0, 0.0, 0.0)
    return_circuit.append(gate, [q[0]])

    return_circuit.append(CnotGate(), [q[0], q[1]])

    return_circuit.append(u1_gate, [q[1]])

    return_circuit.append(u2_gate, [q[0]])

    return return_circuit
def two_qubit_kak(unitary):
    """Decompose a two-qubit gate over SU(2)+CNOT using the KAK decomposition.

    Args:
        unitary (Operator): a 4x4 unitary operator to decompose.

    Returns:
        QuantumCircuit: a circuit implementing the unitary over SU(2)+CNOT

    Raises:
        QiskitError: input not a unitary, or error in KAK decomposition.
    """
    if hasattr(unitary, 'to_operator'):
        # If input is a BaseOperator subclass this attempts to convert
        # the object to an Operator so that we can extract the underlying
        # numpy matrix from `Operator.data`.
        unitary = unitary.to_operator().data
    if hasattr(unitary, 'to_matrix'):
        # If input is Gate subclass or some other class object that has
        # a to_matrix method this will call that method.
        unitary = unitary.to_matrix()
    # Convert to numpy array incase not already an array
    unitary_matrix = np.array(unitary, dtype=complex)
    # Check input is a 2-qubit unitary
    if unitary_matrix.shape != (4, 4):
        raise QiskitError("two_qubit_kak: Expected 4x4 matrix")
    if not is_unitary_matrix(unitary_matrix):
        raise QiskitError("Input matrix is not unitary.")
    phase = la.det(unitary_matrix)**(-1.0/4.0)
    # Make it in SU(4), correct phase at the end
    U = phase * unitary_matrix
    # B changes to the Bell basis
    B = (1.0/math.sqrt(2)) * np.array([[1, 1j, 0, 0],
                                       [0, 0, 1j, 1],
                                       [0, 0, 1j, -1],
                                       [1, -1j, 0, 0]], dtype=complex)

    # We also need B.conj().T below
    Bdag = B.conj().T
    # U' = Bdag . U . B
    Uprime = Bdag.dot(U.dot(B))
    # M^2 = trans(U') . U'
    M2 = Uprime.T.dot(Uprime)

    # Diagonalize M2
    # Must use diagonalization routine which finds a real orthogonal matrix P
    # when M2 is real.
    D, P = la.eig(M2)
    D = np.diag(D)
    # If det(P) == -1 then in O(4), apply a swap to make P in SO(4)
    if abs(la.det(P)+1) < 1e-5:
        swap = np.array([[1, 0, 0, 0],
                         [0, 0, 1, 0],
                         [0, 1, 0, 0],
                         [0, 0, 0, 1]], dtype=complex)
        P = P.dot(swap)
        D = swap.dot(D.dot(swap))

    Q = np.sqrt(D)  # array from elementwise sqrt
    # Want to take square root so that Q has determinant 1
    if abs(la.det(Q)+1) < 1e-5:
        Q[0, 0] = -Q[0, 0]

    # Q^-1*P.T = P' -> QP' = P.T (solve for P' using Ax=b)
    Pprime = la.solve(Q, P.T)
    # K' now just U' * P * P'
    Kprime = Uprime.dot(P.dot(Pprime))

    K1 = B.dot(Kprime.dot(P.dot(Bdag)))
    A = B.dot(Q.dot(Bdag))
    K2 = B.dot(P.T.dot(Bdag))
    # KAK = K1 * A * K2
    KAK = K1.dot(A.dot(K2))

    # Verify decomp matches input unitary.
    if la.norm(KAK - U) > 1e-6:
        raise QiskitError("two_qubit_kak: KAK decomposition " +
                          "does not return input unitary.")

    # Compute parameters alpha, beta, gamma so that
    # A = exp(i * (alpha * XX + beta * YY + gamma * ZZ))
    xx = np.array([[0, 0, 0, 1],
                   [0, 0, 1, 0],
                   [0, 1, 0, 0],
                   [1, 0, 0, 0]], dtype=complex)

    yy = np.array([[0, 0, 0, -1],
                   [0, 0, 1, 0],
                   [0, 1, 0, 0],
                   [-1, 0, 0, 0]], dtype=complex)

    zz = np.array([[1, 0, 0, 0],
                   [0, -1, 0, 0],
                   [0, 0, -1, 0],
                   [0, 0, 0, 1]], dtype=complex)

    A_real_tr = A.real.trace()
    alpha = math.atan2(A.dot(xx).imag.trace(), A_real_tr)
    beta = math.atan2(A.dot(yy).imag.trace(), A_real_tr)
    gamma = math.atan2(A.dot(zz).imag.trace(), A_real_tr)

    # K1 = kron(U1, U2) and K2 = kron(V1, V2)
    # Find the matrices U1, U2, V1, V2

    # Find a block in K1 where U1_ij * [U2] is not zero
    L = K1[0:2, 0:2]
    if la.norm(L) < 1e-9:
        L = K1[0:2, 2:4]
        if la.norm(L) < 1e-9:
            L = K1[2:4, 2:4]
    # Remove the U1_ij prefactor
    Q = L.dot(L.conj().T)
    U2 = L / math.sqrt(Q[0, 0].real)

    # Now grab U1 given we know U2
    R = K1.dot(np.kron(np.identity(2), U2.conj().T))
    U1 = np.zeros((2, 2), dtype=complex)
    U1[0, 0] = R[0, 0]
    U1[0, 1] = R[0, 2]
    U1[1, 0] = R[2, 0]
    U1[1, 1] = R[2, 2]

    # Repeat K1 routine for K2
    L = K2[0:2, 0:2]
    if la.norm(L) < 1e-9:
        L = K2[0:2, 2:4]
        if la.norm(L) < 1e-9:
            L = K2[2:4, 2:4]
    Q = np.dot(L, np.transpose(L.conjugate()))
    V2 = L / np.sqrt(Q[0, 0])
    R = np.dot(K2, np.kron(np.identity(2), np.transpose(V2.conjugate())))

    V1 = np.zeros_like(U1)
    V1[0, 0] = R[0, 0]
    V1[0, 1] = R[0, 2]
    V1[1, 0] = R[2, 0]
    V1[1, 1] = R[2, 2]

    if la.norm(np.kron(U1, U2) - K1) > 1e-4:
        raise QiskitError("two_qubit_kak: K1 != U1 x U2")
    if la.norm(np.kron(V1, V2) - K2) > 1e-4:
        raise QiskitError("two_qubit_kak: K2 != V1 x V2")

    test = la.expm(1j*(alpha * xx + beta * yy + gamma * zz))
    if la.norm(A - test) > 1e-4:
        raise QiskitError("two_qubit_kak: " +
                          "Matrix A does not match xx,yy,zz decomposition.")

    # Circuit that implements K1 * A * K2 (up to phase), using
    # Vatan and Williams Fig. 6 of quant-ph/0308006v3
    # Include prefix and suffix single-qubit gates into U2, V1 respectively.

    V2 = np.array([[np.exp(1j*np.pi/4), 0],
                   [0, np.exp(-1j*np.pi/4)]], dtype=complex).dot(V2)
    U1 = U1.dot(np.array([[np.exp(-1j*np.pi/4), 0],
                          [0, np.exp(1j*np.pi/4)]], dtype=complex))

    # Corrects global phase: exp(ipi/4)*phase'
    U1 = U1.dot(np.array([[np.exp(1j*np.pi/4), 0],
                          [0, np.exp(1j*np.pi/4)]], dtype=complex))
    U1 = phase.conjugate() * U1

    # Test
    g1 = np.kron(V1, V2)
    g2 = np.array([[1, 0, 0, 0],
                   [0, 0, 0, 1],
                   [0, 0, 1, 0],
                   [0, 1, 0, 0]], dtype=complex)

    theta = 2*gamma - np.pi/2

    Ztheta = np.array([[np.exp(1j*theta/2), 0],
                       [0, np.exp(-1j*theta/2)]], dtype=complex)

    kappa = np.pi/2 - 2*alpha
    Ykappa = np.array([[math.cos(kappa/2), math.sin(kappa/2)],
                       [-math.sin(kappa/2), math.cos(kappa/2)]], dtype=complex)
    g3 = np.kron(Ztheta, Ykappa)
    g4 = np.array([[1, 0, 0, 0],
                   [0, 1, 0, 0],
                   [0, 0, 0, 1],
                   [0, 0, 1, 0]], dtype=complex)

    zeta = 2*beta - np.pi/2
    Yzeta = np.array([[math.cos(zeta/2), math.sin(zeta/2)],
                      [-math.sin(zeta/2), math.cos(zeta/2)]], dtype=complex)
    g5 = np.kron(np.identity(2), Yzeta)
    g6 = g2
    g7 = np.kron(U1, U2)

    V = g2.dot(g1)
    V = g3.dot(V)
    V = g4.dot(V)
    V = g5.dot(V)
    V = g6.dot(V)
    V = g7.dot(V)

    if la.norm(V - U*phase.conjugate()) > 1e-6:
        raise QiskitError("two_qubit_kak: " +
                          "sequence incorrect, unknown error")

    v1_param = euler_angles_1q(V1)
    v2_param = euler_angles_1q(V2)
    u1_param = euler_angles_1q(U1)
    u2_param = euler_angles_1q(U2)

    v1_gate = U3Gate(v1_param[0], v1_param[1], v1_param[2])
    v2_gate = U3Gate(v2_param[0], v2_param[1], v2_param[2])
    u1_gate = U3Gate(u1_param[0], u1_param[1], u1_param[2])
    u2_gate = U3Gate(u2_param[0], u2_param[1], u2_param[2])

    q = QuantumRegister(2)
    return_circuit = QuantumCircuit(q)

    return_circuit.append(v1_gate, [q[1]])

    return_circuit.append(v2_gate, [q[0]])

    return_circuit.append(CnotGate(), [q[0], q[1]])

    gate = U3Gate(0.0, 0.0, -2.0*gamma + np.pi/2.0)
    return_circuit.append(gate, [q[1]])

    gate = U3Gate(-np.pi/2.0 + 2.0*alpha, 0.0, 0.0)
    return_circuit.append(gate, [q[0]])

    return_circuit.append(CnotGate(), [q[1], q[0]])

    gate = U3Gate(-2.0*beta + np.pi/2.0, 0.0, 0.0)
    return_circuit.append(gate, [q[0]])

    return_circuit.append(CnotGate(), [q[0], q[1]])

    return_circuit.append(u1_gate, [q[1]])

    return_circuit.append(u2_gate, [q[0]])

    return return_circuit
Decompose a two-qubit gate over SU(2)+CNOT using the KAK decomposition.

    Args:
        unitary (Operator): a 4x4 unitary operator to decompose.

    Returns:
        QuantumCircuit: a circuit implementing the unitary over SU(2)+CNOT

    Raises:
        QiskitError: input not a unitary, or error in KAK decomposition.