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
def two_qubit_kak ( unitary ) : 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.
Decompose a two - qubit gate over SU ( 2 ) + CNOT using the KAK decomposition .