Decompose a 2-qubit unitary composed of two 1-qubit local gates.

    Uses the "Nearest Kronecker Product" algorithm. Will give erratic
    results if the gate is not the direct product of two 1-qubit gates.

def kronecker_decomposition(gate: Gate) -> Circuit:
    """
    Decompose a 2-qubit unitary composed of two 1-qubit local gates.

    Uses the "Nearest Kronecker Product" algorithm. Will give erratic
    results if the gate is not the direct product of two 1-qubit gates.
    """
    # An alternative approach would be to take partial traces, but
    # this approach appears to be more robust.

    if gate.qubit_nb != 2:
        raise ValueError('Expected 2-qubit gate')

    U = asarray(gate.asoperator())
    rank = 2**gate.qubit_nb
    U /= np.linalg.det(U) ** (1/rank)

    R = np.stack([U[0:2, 0:2].reshape(4),
                  U[0:2, 2:4].reshape(4),
                  U[2:4, 0:2].reshape(4),
                  U[2:4, 2:4].reshape(4)])
    u, s, vh = np.linalg.svd(R)
    v = vh.transpose()
    A = (np.sqrt(s[0]) * u[:, 0]).reshape(2, 2)
    B = (np.sqrt(s[0]) * v[:, 0]).reshape(2, 2)

    q0, q1 = gate.qubits
    g0 = Gate(A, qubits=[q0])
    g1 = Gate(B, qubits=[q1])

    if not gates_close(gate, Circuit([g0, g1]).asgate()):
        raise ValueError("Gate cannot be decomposed into two 1-qubit gates")

    circ = Circuit()
    circ += zyz_decomposition(g0)
    circ += zyz_decomposition(g1)

    assert gates_close(gate, circ.asgate())  # Sanity check

    return circ