def canonical_decomposition(gate: Gate) -> Circuit:
    """Decompose a 2-qubit gate by removing local 1-qubit gates to leave
    the non-local canonical two-qubit gate. [1]_ [2]_ [3]_ [4]_

    Returns: A Circuit of 5 gates: two initial 1-qubit gates; a CANONICAL
    gate, with coordinates in the Weyl chamber; two final 1-qubit gates

    The canonical coordinates can be found in circ.elements[2].params

    More or less follows the algorithm outlined in [2]_.

    .. [1] A geometric theory of non-local two-qubit operations, J. Zhang,
        J. Vala, K. B. Whaley, S. Sastry quant-ph/0291120
    .. [2] An analytical decomposition protocol for optimal implementation of
        two-qubit entangling gates. M. Blaauboer, R.L. de Visser,
        cond-mat/0609750
    .. [3] Metric structure of two-qubit gates, perfect entangles and quantum
        control, P. Watts, M. O'Conner, J. Vala, Entropy (2013)
    .. [4] Constructive Quantum Shannon Decomposition from Cartan Involutions
        B. Drury, P. Love, arXiv:0806.4015
    """

    # Implementation note: The canonical decomposition is easy. Constraining
    # canonical coordinates to the Weyl chamber is easy. But doing the
    # canonical decomposition with the canonical gate in the Weyl chamber
    # proved to be surprisingly tricky.

    # Unitary transform to Magic Basis of Bell states
    Q = np.asarray([[1, 0, 0, 1j],
                    [0, 1j, 1, 0],
                    [0, 1j, -1, 0],
                    [1, 0, 0, -1j]]) / np.sqrt(2)
    Q_H = Q.conj().T

    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)  # U is in SU(4) so det U = 1

    U_mb = Q_H @ U @ Q   # Transform gate to Magic Basis [1, (eq. 17, 18)]
    M = U_mb.transpose() @ U_mb         # Construct M matrix [1, (eq. 22)]

    # Diagonalize symmetric complex matrix
    eigvals, eigvecs = _eig_complex_symmetric(M)

    lambdas = np.sqrt(eigvals)          # Eigenvalues of F
    # Lambdas only fixed up to a sign. So make sure det F = 1 as it should
    det_F = np.prod(lambdas)
    if det_F.real < 0:
        lambdas[0] *= -1

    coords, signs, perm = _constrain_to_weyl(lambdas)

    # Construct local and canonical gates in magic basis
    lambdas = (lambdas*signs)[perm]
    O2 = (np.diag(signs) @ eigvecs.transpose())[perm]
    F = np.diag(lambdas)
    O1 = U_mb @ O2.transpose() @ F.conj()

    # Sanity check: Make sure O1 and O2 are orthogonal
    assert np.allclose(np.eye(4), O2.transpose() @ O2)  # Sanity check
    assert np.allclose(np.eye(4), O1.transpose() @ O1)  # Sanity check

    # Sometimes O1 & O2 end up with det = -1, instead of +1 as they should.
    # We can commute a diagonal matrix through F to fix this up.
    neg = np.diag([-1, 1, 1, 1])
    if np.linalg.det(O2).real < 0:
        O2 = neg @ O2
        O1 = O1 @ neg

    # Transform gates back from magic basis
    K1 = Q @ O1 @ Q_H
    A = Q @ F @ Q_H
    K2 = Q @ O2 @ Q_H

    assert gates_close(Gate(U), Gate(K1 @ A @ K2))  # Sanity check
    canon = CANONICAL(coords[0], coords[1], coords[2], 0, 1)

    # Sanity check
    assert gates_close(Gate(A, qubits=gate.qubits), canon, tolerance=1e-4)

    # Decompose local gates into the two component 1-qubit gates
    gateK1 = Gate(K1, qubits=gate.qubits)
    circK1 = kronecker_decomposition(gateK1)
    assert gates_close(gateK1, circK1.asgate())  # Sanity check

    gateK2 = Gate(K2, qubits=gate.qubits)
    circK2 = kronecker_decomposition(gateK2)
    assert gates_close(gateK2, circK2.asgate())  # Sanity check

    # Build and return circuit
    circ = Circuit()
    circ += circK2
    circ += canon
    circ += circK1

    return circ