[unitaryhack] Add a hybrid Gaussian/non-Gaussian compiler that merges Gaussian gates

josh146 commented 3 years ago

This issue has been tagged for a bounty during unitaryHACK

It would be great to have a compiler in Strawberry Fields that:

Applies to arbitrary programs, containing both Gaussian and non-Gaussian operations
Combines the Gaussian operations into a single symplectic operation, by multiplying together their symplectic forms.

This could significantly reduce overhead during computation.

As an example, consider the following diagram:

This circuit contains a mixture of Gaussian gates (green) and non-Gaussian gates (red). While the cross-Kerr gate in the middle of the circuit prevents us from combining all Gaussian gates, we can create two Gaussian gates, S1 and S2:

Note that this grouping is not unique! We could have potentially included the second beamsplitter within S1, rather than S2 as we did.

lilianchih commented 3 years ago

I am interested in working on this issue. Which file/folder should I start looking into?

co9olguy commented 3 years ago

Hi @lilianchih, that's great! :slightly_smiling_face:

The best place to start is to look at the existing compilers in SF, located here: https://github.com/XanaduAI/strawberryfields/tree/master/strawberryfields/compilers

The file compiler.py defines the abstract base class shared by all compilers
The most important ingredient to add in a new compiler is to replace the compile method that manages that specific case
You'll also notice that some files (e.g., fock.py) define "Compilers" but they don't actually replace the default compile method. The ones that do are in the files gaussian_unitary.py, gbs.py, xcov.py, xstrict.py, xunitary.py. These are the best references for this issue

maxtremblay commented 3 years ago

I would also be interested in giving a shot at this one. And I am open for collaboration @lilianchih.

I have a few ideas.

I guess, the main goal is the minimize the total number of combined gaussian gates at the end. A simple heuristic, would be to start from the beginning of the circuit, combine as most gates as possible and then start a new group when it is not possible to include gate anymore in the previous group.

Also, I think it is possible to use a dynamic programming approach to get the optimal solution. But, this would required a bit more thinking.

lilianchih commented 3 years ago

@maxtremblay I am happy to collaborate.

I am also thinking about the greedy algorithm that you described. I think it is optimal in the sense that it results in the minimum number of groups.

nquesada commented 3 years ago

Neat ideas @lilianchih and @maxtremblay ! You could also represent the program as a directed acyclic graph and use some of the tools from networkx to start reducing nodes until you are only left with Gaussian transformation interleaved with non Gaussian ones.

nquesada commented 3 years ago

A very basic test for your_compiler would go as follows:

# Import and preliminaries
import strawberryfields as sf
from strawberryfields.ops import Ket, BSgate, Interferometer, Kgate, CKgate,Rgate
import numpy as np

modes = 2
cutoff_dim = 3  # (1+ total number of photons)

initial_state = np.zeros([cutoff_dim] * modes, dtype=np.complex)
# The ket below corresponds to a single horizontal photon in each of the modes
initial_state[1, 1] = 1

# Here is the main program
# We create the input state and then send it through a network of beamsplitters, rotations and Kerr gates
# Since all the gates preserve the number of photons we need not to worry about cutoff issues
prog = sf.Program(2)
with prog.context as q:
    Ket(initial_state) | q  # Initial state preparation
    # Gaussian layer
    BSgate(0.5,0.9)|(q[0],q[1])
    Rgate(0.6)|q[0]
    BSgate(1.9,0.7)|(q[0],q[1])
    Rgate(0.5)|q[1]
    BSgate(1.5,2.9)|(q[0],q[1])
    BSgate(1.8,0.9)|(q[0],q[1])
    BSgate(0.5,0.9)|(q[0],q[1])
    Rgate(0.6)|q[1]
    BSgate(1.9,0.7)|(q[0],q[1])
    Rgate(0.5)|q[1]
    BSgate(1.5,2.9)|(q[0],q[1])
    BSgate(1.8,0.9)|(q[0],q[1])
    Rgate(0.5)|q[1]
    # Non Gaussian layer
    Kgate(0.5)|q[0]
    Kgate(0.4)|q[1]
    # Gaussian layer
    BSgate(0.5,0.9)|(q[0],q[1])
    Rgate(0.6)|q[0]
    BSgate(1.9,0.7)|(q[0],q[1])
    Rgate(0.5)|q[1]
    BSgate(1.5,2.9)|(q[0],q[1])
    BSgate(1.8,0.9)|(q[0],q[1])
    Rgate(0.5)|q[1]
    BSgate(0.5,0.9)|(q[0],q[1])
    Rgate(0.6)|q[1]
    BSgate(1.9,0.7)|(q[0],q[1])
    Rgate(0.5)|q[1]
    BSgate(1.5,2.9)|(q[0],q[1])
    BSgate(1.8,0.9)|(q[0],q[1])
    Rgate(0.5)|q[1]
    # Non Gaussian layer
    CKgate(0.5)|(q[0],q[1])
# We run the simulation
eng = sf.Engine("fock", backend_options={"cutoff_dim": cutoff_dim})
results_norm = eng.run(prog)
prog_merged = prog.compile(compiler="your_compiler")
results_merged = eng.run(prog_merged)
ket = results_norm.state.ket()
ket_merged = results_merged.state.ket()
np.allclose(np.abs(np.sum(np.conj(ket) * ket_merged)), 1)

federico0112 commented 3 years ago

Heres an edit to that test to make it more complex, it uses three qumodes:

# Import and preliminaries
import strawberryfields as sf
from strawberryfields.ops import Ket, BSgate, Interferometer, Kgate, CKgate,Rgate
import numpy as np

modes = 3
cutoff_dim = 4  # (1+ total number of photons)

initial_state = np.zeros([cutoff_dim] * modes, dtype=np.complex)
# The ket below corresponds to a single horizontal photon in each of the modes
initial_state[1, 1, 1] = 1

# Here is the main program
# We create the input state and then send it through a network of beamsplitters, rotations and Kerr gates
# Since all the gates preserve the number of photons we need not to worry about cutoff issues
prog = sf.Program(modes)
with prog.context as q:
    Ket(initial_state) | q  # Initial state preparation
    # Gaussian layer
    BSgate(0.5,0.9)|(q[0],q[1])
    Rgate(0.6)|q[0]
    Rgate(0.2)|q[2]
    BSgate(0.5,0.9)|(q[1],q[2])
    BSgate(1.9,0.7)|(q[0],q[1])
    Rgate(0.5)|q[1]
    BSgate(1.5,2.9)|(q[0],q[1])
    BSgate(1.8,0.9)|(q[0],q[1])
    BSgate(1.9,0.7)|(q[1],q[2])
    BSgate(0.5,0.9)|(q[0],q[1])
    Rgate(0.6)|q[1]
    Rgate(0.9)|q[2]
    BSgate(1.9,0.7)|(q[0],q[1])
    BSgate(1.8,0.9)|(q[1],q[2])
    Rgate(0.5)|q[1]
    Rgate(0.3)|q[2]
    BSgate(1.5,2.9)|(q[0],q[1])
    BSgate(1.8,0.9)|(q[0],q[1])
    BSgate(0.5,0.9)|(q[1],q[2])
    Rgate(0.5)|q[1]
    # Non Gaussian layer
    Kgate(0.5)|q[0]
    Kgate(0.4)|q[1]
    Kgate(0.6)|q[2]
    # Gaussian layer
    BSgate(0.5,0.9)|(q[0],q[1])
    Rgate(0.6)|q[0]
    Rgate(0.2)|q[2]
    BSgate(0.5,0.9)|(q[1],q[2])
    BSgate(1.9,0.7)|(q[0],q[1])
    Rgate(0.5)|q[1]
    BSgate(1.5,2.9)|(q[0],q[1])
    BSgate(1.8,0.9)|(q[0],q[1])
    BSgate(1.9,0.7)|(q[1],q[2])
    BSgate(0.5,0.9)|(q[0],q[1])
    Rgate(0.6)|q[1]
    Rgate(0.9)|q[2]
    BSgate(1.9,0.7)|(q[0],q[1])
    BSgate(1.8,0.9)|(q[1],q[2])
    Rgate(0.5)|q[1]
    Rgate(0.3)|q[2]
    BSgate(1.5,2.9)|(q[0],q[1])
    BSgate(1.8,0.9)|(q[0],q[1])
    BSgate(0.5,0.9)|(q[1],q[2])
    Rgate(0.5)|q[1]
    # Non Gaussian layer
    CKgate(0.5)|(q[0],q[1])
    CKgate(0.3)|(q[1],q[2])

# We run the simulation
eng = sf.Engine("fock", backend_options={"cutoff_dim": cutoff_dim})
results_norm = eng.run(prog)
prog_merged = prog.compile(compiler="gaussian_merge")
results_merged = eng.run(prog_merged)
ket = results_norm.state.ket()
ket_merged = results_merged.state.ket()
if np.allclose(np.abs(np.sum(np.conj(ket) * ket_merged)), 1):
    print("Wow its working!!!")

federico0112 commented 3 years ago

Here is my commit with a working implementation of Gaussian Merge: https://github.com/federico0112/StrawberyFields_GaussianMerge/commit/633052d292a9092b29cfbce5436ba367be9a773e

nquesada commented 3 years ago

Implemented in https://github.com/XanaduAI/strawberryfields/pull/591

XanaduAI / strawberryfields

[unitaryhack] Add a hybrid Gaussian/non-Gaussian compiler that merges Gaussian gates #574