Implementation of a network of interconnected beamsplitters for coherent state evolution (strawberryfields)

Hello, good morning dear colleagues. I’m trying to evolve a set of coherent states through a network of 50/50 symmetric divisors. My problem is related to the amount of RAM required by my code. Are there any ways around this problem? What would they be?

import strawberryfields as sf
from strawberryfields.ops import *

import numpy as np
from numpy import pi, sqrt

# set the random seed

prog = sf.Program(6)

alpha = 1+0.5j
r = np.abs(alpha)
phi = np.angle(alpha)

with prog.context as q:
    # prepare initial states
    Coherent(r, phi) | q[0]
    Coherent(r, phi) | q[1]
    Coherent(r, phi) | q[2]
    Coherent(r, phi) | q[3]
    Coherent(r, phi) | q[4]
    Coherent(r, phi) | q[5]
    # apply gates
    # Primeira bateria de BS
    BSgate() | (q[0], q[1])
    BSgate() | (q[2], q[3])
    BSgate() | (q[4], q[5])

    # Segunda bateria de BS
    BSgate() | (q[2], q[4])
    BSgate() | (q[0], q[5])
    BSgate() | (q[1], q[3])

eng = sf.Engine('fock', backend_options={"cutoff_dim": 6})

result =, shots=1, modes=None, compile_options={})

eng = sf.Engine('fock', backend_options={"cutoff_dim": 6})

result =, shots=1, modes=None, compile_options={})

If you want help with diagnosing an error, please put the full error message below:

MemoryError                               Traceback (most recent call last)
Cell In[9], line 1
----> 1 result =, shots=1, modes=None, compile_options={})

File ~\anaconda3\lib\site-packages\strawberryfields\, in, program, args, compile_options, **kwargs)
    565         if c.op.measurement_deps and eng_run_options["shots"] > 1:
    566             raise NotImplementedError(
    567                 "Feed-forwarding of measurements cannot be used together with multiple shots."
    568             )
--> 570 return super()._run(
    571     program_lst, args=args, compile_options=compile_options, **eng_run_options
    572 )

File ~\anaconda3\lib\site-packages\strawberryfields\, in BaseEngine._run(self, program, args, compile_options, **kwargs)
    303 p.bind_params(args)
    304 p.lock()
--> 306 _, self.samples, self.samples_dict = self._run_program(p, **kwargs)
    307 self.run_progs.append(p)
    309 if isinstance(p, TDMProgram) and received_rolled:

File ~\anaconda3\lib\site-packages\strawberryfields\, in LocalEngine._run_program(self, prog, **kwargs)
    427 for cmd in prog.circuit:
    428     try:
    429         # try to apply it to the backend and, if op is a measurement, store it in values
--> 430         val = cmd.op.apply(cmd.reg, self.backend, **kwargs)
    431         if val is not None:
    432             for i, r in enumerate(cmd.reg):

File ~\anaconda3\lib\site-packages\strawberryfields\, in Operation.apply(self, reg, backend, **kwargs)
    226 temp = [rr.ind for rr in reg]
    227 # call the child class specialized _apply method
--> 228 return self._apply(temp, backend, **kwargs)

File ~\anaconda3\lib\site-packages\strawberryfields\, in Coherent._apply(self, reg, backend, **kwargs)
    627 phi = par_evaluate(self.p[1])
    629 self._check_for_complex_args([r, phi], "Coherent(r, phi)")
--> 631 backend.prepare_coherent_state(r, phi, *reg)

File ~\anaconda3\lib\site-packages\strawberryfields\backends\fockbackend\, in FockBackend.prepare_coherent_state(self, r, phi, mode)
    152 def prepare_coherent_state(self, r, phi, mode):
--> 153     self.circuit.prepare_mode_coherent(r, phi, self._remap_modes(mode))

File ~\anaconda3\lib\site-packages\strawberryfields\backends\fockbackend\, in Circuit.prepare_mode_coherent(self, r, phi, mode)
    501 """
    502 Prepares a mode in a coherent state.
    503 """
    504 if self._pure:
--> 505     self.prepare(ops.coherentState(r, phi, self._trunc), mode)
    506 else:
    507     st = ops.coherentState(r, phi, self._trunc)

File ~\anaconda3\lib\site-packages\strawberryfields\backends\fockbackend\, in Circuit.prepare(self, state, mode)
    485 if isinstance(mode, int):
    486     mode = [mode]
--> 487 self.prepare_multimode(state, mode)

File ~\anaconda3\lib\site-packages\strawberryfields\backends\fockbackend\, in Circuit.prepare_multimode(self, state, modes)
    455     reduced_state = ops.partial_trace(self._state, self._num_modes, modes)
    457     # Insert state at the end (I know there is also tensor() from ops but it has extra aguments wich only confuse here)
--> 458     self._state = np.tensordot(reduced_state, state, axes=0)
    460     # unless the preparation was meant to go into the last modes in the standard order, we need to swap indices around
    461 if modes != list(range(self._num_modes - len(modes), self._num_modes)):

File <__array_function__ internals>:180, in tensordot(*args, **kwargs)

File ~\anaconda3\lib\site-packages\numpy\core\, in tensordot(a, b, axes)
   1136 at = a.transpose(newaxes_a).reshape(newshape_a)
   1137 bt = b.transpose(newaxes_b).reshape(newshape_b)
-> 1138 res = dot(at, bt)
   1139 return res.reshape(olda + oldb)

File <__array_function__ internals>:180, in dot(*args, **kwargs)

MemoryError: Unable to allocate 32.4 GiB for an array with shape (60466176, 36) and data type complex128

conda 23.3.1
strawberryfields 0.23.0
numpy 1.23.5

Hey @Antonio_Aguiar! You can reduce the cutoff_dim, but that isn’t ideal :confused:. Do you need the fock backend? The gaussian backend might be better here — it’s exact and less computationally intensive than an equivalent fock calculation. You would just need to do the following:

eng = sf.Engine('gaussian')

Although this isn’t a problem with your example, the gaussian backend cannot represent non-gaussian states :smile:. You can find out more information on the various backends available here: Circuits — Strawberry Fields 0.23.0 documentation

1 Like

Thank you for your attention. The multimode solution suits me best. About the last code, when you defined the states, there are six coherent states that will then pass through the beam splitters, correct? My next step is to build a function that switches the inputs between |alpha> and |-alpha> so that I evaluate all 64 possible inputs and evolution of my system. Do you know a way to do this analysis?

I don’t speak spanish :sweat_smile:, but I assume your question is similar to your follow-up question here. This is possible in SF! You just need to construct a proper for-loop and correctly change the sign of r :slight_smile:

1 Like

Thank you for your patience, I must have got confused in the translation of my message. I’m glad you understand and thanks for the help.

1 Like

My pleasure! Glad this helps :smiley:

1 Like

Hi good morning, before I start to evaluate all the possible entries and their evolution of entries I would like to evaluate an entry and an exit. How do I rate my outings? How can I perform my measurements to find out if |alpha>, |-alpha> or some different coherent state came out?

I suppose this might be a little redundant, but you could just apply a Displacement operator at the end of your circuit with the positive value of alpha then do a fock measurement:

prog = sf.Program(1)

alpha = 1+0.5j
r = np.abs(alpha)
phi = np.angle(alpha)

eng = sf.Engine('gaussian')

for _ in range(2):
    results = []
    for rval in [-r, r]:
        prog = sf.Program(1)

        with prog.context as q:
            # prepare initial states
            Coherent(rval, phi) | q[0]
            Dgate(r, phi) | q[0]
            MeasureFock() | q[0]

        result =, shots=1).samples

[0, 4]
[0, 2]

If the value of \alpha is negative, then the result will be 0. If the value of \alpha is positive, then the result will be non-zero.

1 Like

There are six input modes and each one receives a state composed of a superposition of coherent states that will be used for binary coding where |𝛼> is 1 and |-𝛼> corresponds to 0, considering that any result different from these is discarded and will contribute to decreasing the efficiency of the circuit. The system starts with the tensor product of the six input states, three auxiliaries and three for encoding the output states, producing a matrix with six columns by 64 rows where each column indicates which input mode at that position. While each line indicates a combination of input states out of 64 possible. Then that input state goes to the first battery of (3) beam splitters (BS).


How would I do this with six modes and six beam splitters?

In this case, I would know that a positive α came out, but I wouldn’t know if it was exactly the same α that entered one of the inputs? It is interesting that I know when α, -α exited, which will be used for my encoding of zeros and ones, and I also know when none of the two states exited.

You can find some guidance here on conventions for the BSGate operator here under “Details and Conventions”.

And to try and make sure that the circuit you create in Strawberryfields is what you have in mind, you can use the draw_circuit method :smile:

Hello good morning, at the moment my biggest difficulty is to loop my circuit. So that I have all the outputs depending on my inputs |α> and |-α>.

prog = sf.Program(6)

alpha = 1+0.5j
r = np.abs(alpha)
phi = np.angle(alpha)

eng = sf.Engine(‘gaussian’)

for _ in range(2):
results =
for rval in [-r, r]:
prog = sf.Program(6)

    with prog.context as q:
        # prepare initial states
        Coherent(rval, phi) | q[0]
        Dgate(r, phi) | q[0]
        MeasureFock() | q[0]
        Coherent(rval, phi) | q[1]
        Dgate(r, phi) | q[1]
        MeasureFock() | q[1]
        Coherent(rval, phi) | q[2]
        Dgate(r, phi) | q[2]
        MeasureFock() | q[2]
        Coherent(rval, phi) | q[3]
        Dgate(r, phi) | q[3]
        MeasureFock() | q[3]
        Coherent(rval, phi) | q[4]
        Dgate(r, phi) | q[4]
        MeasureFock() | q[4]
        Coherent(rval, phi) | q[5]
        Dgate(r, phi) | q[5]
        MeasureFock() | q[5]
        # BS first battery
        BSgate() | (q[0], q[1])
        BSgate() | (q[2], q[3])
        BSgate() | (q[4], q[5])

        # Second battery of BS
        BSgate() | (q[2], q[4])
        BSgate() | (q[0], q[5])
        BSgate() | (q[1], q[3])

    result =, shots=1).samples

This code generates the following error:


And trying to build the circuit is new to me.


The .item() error is happening because result[0] is already a scalar :slight_smile:. You shouldn’t need to call .item().

For the draw_circuit method, it’s a method of the prog you created :slight_smile:. So, you have to call it like so: prob.draw_circuit(). Note that it will generate LaTex code that can then be rendered in a .tex document with the Qcircuit package. If you’re not too familiar with Latex, then I would recommend not to worry about it :slight_smile:

I would like to understand this result better.

As for the attempt to build the circuit still unsuccessful.

UnsupportedGateException Traceback (most recent call last)
Cell In[7], line 1
----> 1 prog.draw_circuit().tex

File ~\anaconda3\lib\site-packages\strawberryfields\, in Program.draw_circuit(self, tex_dir, write_to_file)
810 r""“Draw the circuit using the Qcircuit :math:\LaTeX package.
812 This will generate the LaTeX code required to draw the quantum circuit
838 list[str]: filename of the written tex document and the written tex content
839 “””
840 drawer = sfcd.Circuit(wires=self.init_num_subsystems)
→ 841 self.print(drawer.parse_op)
842 tex = drawer.dump_to_document()
844 document = None

File ~\anaconda3\lib\site-packages\strawberryfields\, in Program.print(self, print_fn)
273 “”“Print the program contents using Blackbird syntax.
275 Example:
299 print_fn (function): optional custom function to use for string printing
300 “””
301 for k in self.circuit:
→ 302 print_fn(k)

File ~\anaconda3\lib\site-packages\strawberryfields\, in Circuit.parse_op(self, op)
212 wires = list(map(lambda register: register.ind, op.reg))
214 if method is None:
→ 215 raise UnsupportedGateException(
216 “Unsupported operation {0} not printable by circuit builder!”.format(str(op))
217 )
218 if mode == len(wires):
219 method(*wires)

UnsupportedGateException: Unsupported operation Coherent(-1.118, 0.4636) | (q[0]) not printable by circuit builder!

Woops! It’s not a scalar. result[0] is an array, so calling .item() doesn’t make any sense.

For the drawer, it will look like this: prog.draw_circuit("circuit_tex"). However, I just realized that it doesn’t support drawing Coherent operations. So, you can forget about using draw_circuit here :slight_smile: