Optimizer that decomposes two qubit operations into sqrt(iSWAP)

[Dripto]

Is your feature request related to a use case or problem? Please describe.
As it stands, cirq.two_qubit_matrix_to_operations() outputs a decomposition into CZ (or partial CZ) gates. If you want a decomposition into a different gate set, the best (automated) way of doing it is to first decompose into CZ's, and then use cirq.google.ConvertToSqrtIswapGates().optimize_circuit(), however this leads to inefficient circuit decompositions.

In general, any two qubit unitary should be doable with just 3 sqrt(iSWAP) gates, but with the above workflow you end up with circuits that use 6.

Describe the solution you'd like
Bill's proposed solution (#4059) would encompass this one if it recognizes the optimal method for a given gate set. Alternatively, an additional gate_set parameter to input into cirq.two_qubit_matrix_to_operations() which changes the method of decomposition it uses would also suffice,

[optional] Describe alternatives/workarounds you've considered
The current workaround is described above, using cirq.two_qubit_matrix_to_operations() and cirq.google.ConvertToSqrtIswapGates().optimize_circuit() in succession.

What is the urgency from your perspective for this issue? Is it blocking important work?

P1/P2, not vital but also shouldn't be difficult.

Thanks!

[mpharrigan]

the decomposition into CZs is analytic and very fast. Using the gate tabulation could be slow or inaccurate (depending on your grid spacing).

Has anyone tried to pencil out a decomposition a la CZ?

[balopat]

The gate_set parameter can be done in cirq.two_qubit_matrix_to_operations, we will have to work out the sqrt(iSWAP) version of the _non_local_part function:

Cirq/cirq-core/cirq/optimizers/two_qubit_decompositions.py

Lines 242 to 266 in a41df7f

	def _non_local_part(
	q0: 'cirq.Qid',
	q1: 'cirq.Qid',
	interaction_coefficients: Tuple[float, float, float],
	allow_partial_czs: bool,
	atol: float = 1e-8,
	):
	"""Yields non-local operation of KAK decomposition."""
	x, y, z = interaction_coefficients
	if allow_partial_czs or all(_is_trivial_angle(e, atol) for e in [x, y, z]):
	return [
	_parity_interaction(q0, q1, x, atol, ops.Y ** -0.5),
	_parity_interaction(q0, q1, y, atol, ops.X ** 0.5),
	_parity_interaction(q0, q1, z, atol),
	]
	if abs(z) >= atol:
	return _xx_yy_zz_interaction_via_full_czs(q0, q1, x, y, z)
	if y >= atol:
	return _xx_yy_interaction_via_full_czs(q0, q1, x, y)
	return _xx_interaction_via_full_czs(q0, q1, x)

It should give an equivalent (up to global phase) unitary to e^{i * (k_1 * XX + k_2 * YY + k_3 * ZZ)}:

xx,yy,zz = cirq.unitary(cirq.XX), cirq.unitary(cirq.YY), cirq.unitary(cirq.ZZ)
u1 = cirq.map_eigenvalues( sum([term[0] * term[1] for term in zip([xx,yy,zz], interaction_coefficients)]), lambda v: np.exp(1j * v))
u2 = cirq.Circuit(_non_local_part(a,b, interaction_coefficients)).unitary()
cirq.testing.assert_allclose_up_to_global_phase(u1,u2, atol=1e-8)

[balopat]

I went down the rabbit-hole of "how hard this can be". Here are my learnings so far:

---

Sqrt(ISWAP) has a (pi/8, pi/8, 0) KAK vector (the coordinate in the a+ Weyl chamber). I couldn't find yet a paper that describes an optimized decomposition of arbitrary 2 qubit unitaries using sqrt(ISWAP) as the base entangling gate.

The current 6 sqrt(ISWAP) decomposition seems like an upper bound for a generic entangling gate (as long as it's not equivalent to local unitaries or the SWAP gate): https://arxiv.org/pdf/quant-ph/0212109.pdf

There is work by that explores optimal decompositions with base entangling gates that have the (pi/4, alpha, 0) KAK vector (super-controlled gates) - 3 applications should be sufficient to simulate an aribtrary 2 qubit unitary. https://arxiv.org/pdf/quant-ph/0407108.pdf - but that is still not the (pi/8, pi/8, 0) class!

Other frameworks:

• Qiskit has a TwoQubitBasisDecomposer, but even that only works with super-controlled base entangling gates. They also use some interesting approximation tricks to reduce the number of two qubit gates described in Appendix B of https://arxiv.org/abs/1811.12926.
• trueq seems to have something similar, but unclear what they support, also not open source

Not sure how useful this is, but also here's the "path" the unitary takes in the a+ Weyl chamber with our 6 sqrt(ISWAP) to get from identity (light green) to a random unitary(dark blue):

Each step is a sqrt(ISWAP) or its inverse.

---

To me it sounds like someone will need to sit down and deduce a more optimal decomposition, but it doesn't seem like an easy task.

[Dripto]

I'll have to think about this, if I figure out the math I'll make another comment on this issue.

In a (maybe) related note, I've written up this little method for what I needed, which prepares any arbitrary 2q state using a single sqrt(iSWAP). I don't know if this is useful enough for us to put it into cirq, but I'll leave it here in case someone finds this issue and can use it:

def prep_2q_state(q: list,
                  state: np.ndarray):

  if not np.allclose(sum([abs(x)**2 for x in state]), 1):
    print('state not normalized: ', sum([abs(x)**2 for x in state]))

  state_matrix = np.array([[state[0], state[1]], [state[2], state[3]]])

  u, s, vh = np.linalg.svd(state_matrix)

  if np.allclose(s[0], 1):
    alpha = 0
  else:
    alpha = np.arccos(np.sqrt(1-s[0]*(np.sqrt(4 - 4*s[0]**2, dtype='complex64'))))
  
  iSWAP_state_matrix = np.array([[np.cos(alpha), -1j*np.sin(alpha)/np.sqrt(2)],[np.sin(alpha)/np.sqrt(2), 0]])

  u_iSWAP, s_iSWAP, vh_iSWAP = np.linalg.svd(iSWAP_state_matrix)

  prep_circ = cirq.Circuit()

  prep_circ.append(cirq.ry(2*alpha).on(q[0]))
  prep_circ.append(cirq.ISwapPowGate(exponent=-0.5).on(*q))

  gate_0 = np.dot(u, np.linalg.inv(u_iSWAP))
  gate_1 = np.dot(vh.T, np.linalg.inv(vh_iSWAP.T))

  prep_circ.append(cirq.single_qubit_matrix_to_phxz(gate_0).on(q[0]))
  prep_circ.append(cirq.single_qubit_matrix_to_phxz(gate_1).on(q[1]))

  return prep_circ 

[mpharrigan]

@Strilanc also pointed out that we have decompose_two_qubit_interaction_into_four_fsim_gates which takes four FSimGate's that are close enough to ISWAP (unfortunately sqrt-iswap is too far away).

[cduck]

This new preprint (page 14) appears to give a complete algorithm for KAK synthesis with sqrt-iSWAP.

[balopat]

Wow, that's awesome, thanks @cduck, very timely.
Also, pretty involved, I guess there are multiple compilation paths based on where the unitary is in the weyl chamber, significantly more complicated than the CZ compilation - but this is definitely a great source, it looks implementable.

[mrwojtek]

Just mention that there is an approximate decomposition of FSim(0, phi) into sqrt(iSWAP)-like gate in ReCirq that we have used for Fermi-Hubbard. It's analytical and fast but fails if the underlying sqrt(iSWAP)-like gate has a large c-phase angle.

[cduck]

I'm looking into implementing the decomposition in that paper. It looks reasonable. I'll post here or make a PR if I make any progress.

[balopat]

That's awesome - thank you, looking forward to it!

[cduck]

For anyone curious, here are the decompositions for some common gates:

import cirq

for gate in [
            cirq.IdentityGate(2),
            cirq.SQRT_ISWAP,
            cirq.SQRT_ISWAP_INV,
            cirq.ISWAP,  # The extra S gates could be canceled
            cirq.SWAP,
            cirq.CZ,
            cirq.CZ**0.5,
            cirq.optimizers.two_qubit_to_fsim._BGate(),
            cirq.H.controlled(),
            cirq.XX,
            cirq.XX**0.5,
        ]:
    u_target = cirq.unitary(gate)
    q0, q1 = cirq.LineQubit.range(2)
    c = cirq.Circuit(
        cirq.two_qubit_matrix_to_sqrt_iswap_operations(q0, q1, u_target)
    )
    print(f'Decompose({gate}) =')
    print(c)
    print()

Decompose(I(2)) =


Decompose(ISWAP**0.5) =
0: ───iSwap───────
      │
1: ───iSwap^0.5───

Decompose(ISWAP**-0.5) =
0: ───Z───iSwap───────Z───
          │
1: ───────iSwap^0.5───────

Decompose(ISWAP) =
0: ───S───iSwap───────iSwap──────────────
          │           │
1: ───────iSwap^0.5───iSwap^0.5───S^-1───

Decompose(SWAP) =
0: ───S───────Y^0.5───iSwap───────Z^-0.17───X^0.5───iSwap───────Y^0.5───iSwap───────Y^0.83───S^-1──────
                      │                             │                   │
1: ───Y^0.5───────────iSwap^0.5───Z^-0.17───X^0.5───iSwap^0.5───Y^0.5───iSwap^0.5───Z────────Y^-0.83───

Decompose(CZ) =
0: ───Z^-0.75───X^0.5───iSwap───────X───Z───iSwap───────X^-0.5───T─────────
                        │                   │
1: ───Z^-0.75───X^0.5───iSwap^0.5───────────iSwap^0.5───X^-0.5───Z^-0.75───

Decompose(CZ**0.5) =
0: ───Z^0.367───X^0.5─────iSwap───────Z───────────iSwap───────X^0.5─────Z^0.883───
                          │                       │
1: ───Z^-0.83───X^0.136───iSwap^0.5───X^(-4/11)───iSwap^0.5───X^0.136───Z^-0.92───

Decompose(<cirq.optimizers.two_qubit_to_fsim._BGate object at 0x1460d8ca0>) =
0: ───────────────iSwap───────Z───────────────iSwap───────Z───────────────
                  │                           │
1: ───X^(-4/11)───iSwap^0.5───Z───X^(-7/11)───iSwap^0.5───Z───X^(-4/11)───

Decompose(CH) =
0: ───Z^0.75─────X^0.5─────Z^0.196───iSwap───────Y───Z^0.392───iSwap───────Z^0.804───X^0.5─────T^-1───────
                                     │                         │
1: ───Z^-0.196───X^(1/3)─────────────iSwap^0.5─────────────────iSwap^0.5───Z^0.608───X^(1/3)───Z^-0.196───

Decompose(XX) =
0: ───X───

1: ───X───

Decompose(XX**0.5) =
0: ───X^0.5─────Z^-0.192───iSwap───────X───Z^(-5/13)───iSwap───────Z^0.192────X^0.5──────────────
                           │                           │
1: ───Z^0.808───X^0────────iSwap^0.5───────────────────iSwap^0.5───Z^-0.784───X^0─────Z^-0.024───

[Dripto]

looks great, thanks Casey!

This will make it a lot easier for users to build things for our sqrt(iSWAP) grids.