Can not build circuit unitary for a three-qubit gate in Cirq

Question

I am trying to implement a three-qubit gate in an eight-qubit circuit. The method I use is the same by which I create my two-qubit gates with no issues. I produce the unitary of the gate and then introduce it to the _unitary_(self) method. I did the same for the three-qubit gate:

import cirq
import numpy as np

class CsAgate(cirq.ThreeQubitGate):
    def __init__(self, theta):
        self.theta = theta

    def _unitary_(self):
        st = np.sin(self.theta)
        ct = np.cos(self.theta)

        return np.array([
        [1., 0, 0, 0,   0, 0, 0, 0,     0, 0, 0, 0,     0, 0, 0, 0],
        [0, ct, st,0,   0, 0, 0, 0,     0, 0, 0, 0,     0, 0, 0, 0],
        [0, st,-ct,0,   0, 0, 0, 0,     0, 0, 0, 0,     0, 0, 0, 0],
        [0, 0, 0, 1.,   0, 0, 0, 0,     0, 0, 0, 0,     0, 0, 0, 0],

        [0, 0, 0, 0,    1., 0, 0, 0,    0, 0, 0, 0,     0, 0, 0, 0],
        [0, 0, 0, 0,    0, ct,st, 0,    0, 0, 0, 0,     0, 0, 0, 0],
        [0, 0, 0, 0,    0, st,-ct,0,    0, 0, 0, 0,     0, 0, 0, 0],
        [0, 0, 0, 0,    0, 0, 0, 1.,    0, 0, 0, 0,     0, 0, 0, 0],

        [0, 0, 0, 0,    0, 0, 0, 0,     0, 0, 0, 0,     0, 0, 0, 0],
        [0, 0, 0, 0,    0, 0, 0, 0,     0, 0, 0, 0,     0, 0, 0, 0],
        [0, 0, 0, 0,    0, 0, 0, 0,     0, 0, 0, 0,     0, 0, 0, 0],
        [0, 0, 0, 0,    0, 0, 0, 0,     0, 0, 0, 0,     0, 0, 0, 0],

        [0, 0, 0, 0,    0, 0, 0, 0,     0, 0, 0, 0,     0, 0, 0, 0],
        [0, 0, 0, 0,    0, 0, 0, 0,     0, 0, 0, 0,     0, 0, 0, 0],
        [0, 0, 0, 0,    0, 0, 0, 0,     0, 0, 0, 0,     0, 0, 0, 0],
        [0, 0, 0, 0,    0, 0, 0, 0,     0, 0, 0, 0,     0, 0, 0, 0],
        ]
        )


    def _circuit_diagram_info_(self, args):
        return 'o','theta({})'.format(round(self.theta, 2)), '0'

    def __str__(self):
        return 'test'

n_qubits = 8
qubits = cirq.LineQubit.range(n_qubits)
my_circuit=cirq.Circuit([CsAgate(2.)(qubits[0], qubits[1], qubits[2])])
print(my_circuit)
print(my_circuit.unitary())

The circuit is generated successfully but Cirq throws an error when calculating the circuit unitary. Any idea how to solve this issue?

Traceback (most recent call last):
  File "/Users/mjahanpo/Desktop/opt_api/test.py", line 46, in <module>
    print(my_circuit.unitary())
  File "/Users/mjahanpo/opt/anaconda3/envs/opt_api/lib/python3.8/site-packages/cirq/circuits/circuit.py", line 1418, in unitary
    result = _apply_unitary_circuit(self, state, qs, dtype)
  File "/Users/mjahanpo/opt/anaconda3/envs/opt_api/lib/python3.8/site-packages/cirq/circuits/circuit.py", line 1960, in _apply_unitary_circuit
    return protocols.apply_unitaries(
  File "/Users/mjahanpo/opt/anaconda3/envs/opt_api/lib/python3.8/site-packages/cirq/protocols/apply_unitary_protocol.py", line 492, in apply_unitaries
    result = apply_unitary(unitary_value=op,
  File "/Users/mjahanpo/opt/anaconda3/envs/opt_api/lib/python3.8/site-packages/cirq/protocols/apply_unitary_protocol.py", line 347, in apply_unitary
    result = strat(unitary_value, args)
  File "/Users/mjahanpo/opt/anaconda3/envs/opt_api/lib/python3.8/site-packages/cirq/protocols/apply_unitary_protocol.py", line 384, in _strat_apply_unitary_from_apply_unitary
    sub_result = func(sub_args)
  File "/Users/mjahanpo/opt/anaconda3/envs/opt_api/lib/python3.8/site-packages/cirq/ops/gate_operation.py", line 113, in _apply_unitary_
    return protocols.apply_unitary(self.gate, args, default=None)
  File "/Users/mjahanpo/opt/anaconda3/envs/opt_api/lib/python3.8/site-packages/cirq/protocols/apply_unitary_protocol.py", line 347, in apply_unitary
    result = strat(unitary_value, args)
  File "/Users/mjahanpo/opt/anaconda3/envs/opt_api/lib/python3.8/site-packages/cirq/protocols/apply_unitary_protocol.py", line 419, in _strat_apply_unitary_from_unitary
    matrix.reshape(val_qid_shape * 2),
ValueError: cannot reshape array of size 256 into shape (2,2,2,2,2,2)

Answer 1

You specified a matrix whose size is consistent with a four qubit gate, not a three qubit gate. Also, the matrix isn't unitary (it has columns containing nothing but zeros). The matrix' size must be consistent with the number of qubits that the gate applies to.

Do this:

    def _unitary_(self):
        st = np.sin(self.theta)
        ct = np.cos(self.theta)

        return np.array([
            [1., 0, 0, 0,   0, 0, 0, 0],
            [0, ct, st,0,   0, 0, 0, 0],
            [0, st,-ct,0,   0, 0, 0, 0],
            [0, 0, 0, 1.,   0, 0, 0, 0],

            [0, 0, 0, 0,    1., 0, 0, 0],
            [0, 0, 0, 0,    0, ct,st, 0],
            [0, 0, 0, 0,    0, st,-ct,0],
            [0, 0, 0, 0,    0, 0, 0, 1.],
        ])

Also, based on that matrix, it looks like this is actually a two qubit gate tensored with an identity gate on a third qubit.