Detailed walkthrough

In this tutorial, we will try to cover any basic code to get our hands on Perceval.

I. Introduction

1. Perceval installation

[1]:
import perceval as pcvl
pcvl.__version__
[1]:
'0.10.0'
[2]:
from perceval.components.unitary_components import PS, BS, PERM
import numpy as np

## Use the symbolic skin for display
from perceval.rendering.circuit import DisplayConfig, SymbSkin
DisplayConfig.select_skin(SymbSkin)

2. BasicStates

In Linear Optical Circuits, photons can have many discrete degrees of freedom, called modes. It can be the frequency, the polarisation, the position, or all of them.

We represent these degrees of freedom with Fock states. If we have \(n\) photons over \(m\) modes, the Fock state \(|s_1,s_2,...,s_m\rangle\) means we have \(s_i\) photons in the \(i^{th}\) mode. Note that \(\sum_{i=1}^m s_i =n\).

In Perceval, we will use the module pcvl.BasicState

[3]:
## Syntax of different BasicState (list, string, etc)
bs1 = pcvl.BasicState('|0,2,0,1>')
bs2 = pcvl.BasicState([0, 2, 0, 1])

if bs1==bs2:
    print("Those are the same states")

## You can iterate on modes
for i, photon_count in enumerate(bs2):
    print(f"There is {photon_count} photon in mode {i}")

Those are the same states
There is 0 photon in mode 0
There is 2 photon in mode 1
There is 0 photon in mode 2
There is 1 photon in mode 3

3. LO-Components

The linear optical components are the elementary blocks which will act on our Fock states.

It’s important to know all the possible components that can be found in Perceval and understand their effects.

[4]:
## Permutation

perm=PERM([2,0,1])

print(perm.name)
print(perm.describe())
pcvl.pdisplay(perm.definition())
pcvl.pdisplay(perm)
PERM
PERM([2, 0, 1])
$\left[\begin{matrix}0 & 1 & 0\\0 & 0 & 1\\1 & 0 & 0\end{matrix}\right]$
[4]:
../_images/notebooks_Tutorial_7_2.svg
[5]:
## Phase shifter
ps = PS(phi=np.pi)

print(ps.name)
print(ps.describe())
pcvl.pdisplay(ps.definition())
pcvl.pdisplay(ps)  # A pdisplay call on a circuit/processor needs to be the last line of a cell

PS
PS(phi=pi)
$\left[\begin{matrix}e^{i \phi}\end{matrix}\right]$
[5]:
../_images/notebooks_Tutorial_8_2.svg
[6]:
## Beam splitters

bs_rx = BS.Rx()  # By default a beam splitter follows the Rx gate convention, so bs=BS() has the same matrix

# But other conventions exist too:
bs_h = BS.H()
bs_ry = BS.Ry()

## Check the difference in the unitary definition:
print("BS.Rx() unitary matrix")
pcvl.pdisplay(bs_rx.definition())
print("BS.H() unitary matrix")
pcvl.pdisplay(bs_h.definition())
print("BS.Ry() unitary matrix")
pcvl.pdisplay(bs_ry.definition())
print("BS displays its convention as a small label")
pcvl.pdisplay(bs_ry)
BS.Rx() unitary matrix
$\left[\begin{matrix}e^{i \left(\phi_{tl} + \phi_{tr}\right)} \cos{\left(\frac{\theta}{2} \right)} & i e^{i \left(\phi_{bl} + \phi_{tr}\right)} \sin{\left(\frac{\theta}{2} \right)}\\i e^{i \left(\phi_{br} + \phi_{tl}\right)} \sin{\left(\frac{\theta}{2} \right)} & e^{i \left(\phi_{bl} + \phi_{br}\right)} \cos{\left(\frac{\theta}{2} \right)}\end{matrix}\right]$
BS.H() unitary matrix
$\left[\begin{matrix}e^{i \left(\phi_{tl} + \phi_{tr}\right)} \cos{\left(\frac{\theta}{2} \right)} & e^{i \left(\phi_{bl} + \phi_{tr}\right)} \sin{\left(\frac{\theta}{2} \right)}\\e^{i \left(\phi_{br} + \phi_{tl}\right)} \sin{\left(\frac{\theta}{2} \right)} & - e^{i \left(\phi_{bl} + \phi_{br}\right)} \cos{\left(\frac{\theta}{2} \right)}\end{matrix}\right]$
BS.Ry() unitary matrix
$\left[\begin{matrix}e^{i \left(\phi_{tl} + \phi_{tr}\right)} \cos{\left(\frac{\theta}{2} \right)} & - e^{i \left(\phi_{bl} + \phi_{tr}\right)} \sin{\left(\frac{\theta}{2} \right)}\\e^{i \left(\phi_{br} + \phi_{tl}\right)} \sin{\left(\frac{\theta}{2} \right)} & e^{i \left(\phi_{bl} + \phi_{br}\right)} \cos{\left(\frac{\theta}{2} \right)}\end{matrix}\right]$
BS displays its convention as a small label
[6]:
../_images/notebooks_Tutorial_9_7.svg
[7]:
# You can ask for the symbolic matrix value of your component with the attribute U
my_ps = PS(phi=np.pi/8)
pcvl.pdisplay(my_ps.U)
# And for the numerical value with the method compute_unitary
pcvl.pdisplay(my_ps.compute_unitary())
print("")

# If you do it for a Beam-Splitter, you can see that by default theta=pi/2, and the phi's are 0
print("A default beam-splitter:")
pcvl.pdisplay(BS().compute_unitary())  #this is a balanced Beamsplitter
print("")

# To control the value of the parameters of a component, several choices are possible:
#  - by setting a numerical value during the creation of the component
print("A Beam-Splitter with a numerical value for theta:")
bs_rx = BS.Rx(theta=10)
pcvl.pdisplay(bs_rx.U)
pcvl.pdisplay(bs_rx.compute_unitary())
print("")

#  - by using the syntax pcvl.P to create a symbolic variable
#    (note that you cannot compute the numerical value of your component anymore)
print("A Phase Shifter with a symbolic value for phi:")
ps = PS(phi=pcvl.P('\psi'))
pcvl.pdisplay(ps.U)
print("")

#  - you can still modify the value of a symbolic variable after its creation
#    This is not true for a numerical variable!
print("A beam-splitter with a symbolic variable...")
bs_rx = BS(theta=pcvl.P('toto'))
pcvl.pdisplay(bs_rx.U)
bs_rx.assign({'toto':5})
bs_rx.assign({'toto':10})
print("... set to a numerical value")
pcvl.pdisplay(bs_rx.compute_unitary())
print("")

# To check which parameters can be modified, you can call the method get_parameters
# You can also directly change the output of get_parameters to change the values of the parameters
bs_rx = BS(theta=pcvl.P('toto'), phi_tl = pcvl.P('tata'), phi_tr = -1)
parameters = bs_rx.get_parameters()
parameters[0].set_value(np.pi)
print("Modified parameters...")
for param in parameters:
    print("    ", param)
print("... and successfully modified Beam-Splitter:")
pcvl.pdisplay(bs_rx.U)
$\left[\begin{matrix}e^{0.392699081698724 i}\end{matrix}\right]$
$\left[\begin{matrix}0.92388 + 0.382683 i\end{matrix}\right]$

A default beam-splitter:
$\left[\begin{matrix}\frac{\sqrt{2}}{2} & \frac{\sqrt{2} i}{2}\\\frac{\sqrt{2} i}{2} & \frac{\sqrt{2}}{2}\end{matrix}\right]$

A Beam-Splitter with a numerical value for theta:
$\left[\begin{matrix}\cos{\left(5 \right)} & i \sin{\left(5 \right)}\\i \sin{\left(5 \right)} & \cos{\left(5 \right)}\end{matrix}\right]$
$\left[\begin{matrix}0.283662 & - 0.958924 i\\- 0.958924 i & 0.283662\end{matrix}\right]$

A Phase Shifter with a symbolic value for phi:
$\left[\begin{matrix}e^{i \psi}\end{matrix}\right]$

A beam-splitter with a symbolic variable...
$\left[\begin{matrix}\cos{\left(\frac{toto}{2} \right)} & i \sin{\left(\frac{toto}{2} \right)}\\i \sin{\left(\frac{toto}{2} \right)} & \cos{\left(\frac{toto}{2} \right)}\end{matrix}\right]$
... set to a numerical value
$\left[\begin{matrix}0.283662 & - 0.958924 i\\- 0.958924 i & 0.283662\end{matrix}\right]$

Modified parameters...
     Parameter(name='toto', value=3.141592653589793, min_v=0.0, max_v=12.566370614359172)
     Parameter(name='tata', value=None, min_v=0.0, max_v=6.283185307179586)
... and successfully modified Beam-Splitter:
$\left[\begin{matrix}6.12323399573677 \cdot 10^{-17} e^{i \left(tata + 5.28318530717959\right)} & 1.0 i e^{5.28318530717959 i}\\1.0 i e^{i tata} & 6.12323399573677 \cdot 10^{-17}\end{matrix}\right]$
[8]:
## to understand the conventions, you can note that a BS.Rx with the 4 phases phi (top left/right and bottom left/right) can be represented like that

bs_rx_circuit=pcvl.Circuit(2) // (0,PS(phi=pcvl.P("phi_tl"))) // (1,PS(phi=pcvl.P("phi_bl"))) // BS(theta=pcvl.P('theta')) // (0,PS(phi=pcvl.P("phi_tr"))) // (1,PS(phi=pcvl.P("phi_br")))


pcvl.pdisplay(bs_rx_circuit.U)

# we can check it's the same as bs_rx.definition()
pcvl.pdisplay(bs_rx_circuit)

## For this cell, we needed the syntax to builds circuits... Good transition !
$\left[\begin{matrix}e^{i \left(\phi_{tl} + \phi_{tr}\right)} \cos{\left(\frac{\theta}{2} \right)} & i e^{i \left(\phi_{bl} + \phi_{tr}\right)} \sin{\left(\frac{\theta}{2} \right)}\\i e^{i \left(\phi_{br} + \phi_{tl}\right)} \sin{\left(\frac{\theta}{2} \right)} & e^{i \left(\phi_{bl} + \phi_{br}\right)} \cos{\left(\frac{\theta}{2} \right)}\end{matrix}\right]$
[8]:
../_images/notebooks_Tutorial_11_1.svg

II. LO-Circuits

From the LO-components, we can build a LO-circuit, i.e. a sequence of those components acting on our different modes.

1. Syntax

[9]:
circuit = pcvl.Circuit(3)  # Create a 3 mode circuit


circuit.add(0, BS())
circuit.add(0, PS(phi=np.pi/2)).add(1, PS(phi=pcvl.P('phi'))).add(1, BS())

# Equivalent syntax:
# circuit // BS() // PS(phi=np.pi/2) // (1, PS(phi=pcvl.P('phi'))) // (1, BS())


pcvl.pdisplay(circuit.U)
pcvl.pdisplay(circuit)
$\left[\begin{matrix}\frac{\sqrt{2} e^{1.5707963267949 i}}{2} & \frac{\sqrt{2} i e^{1.5707963267949 i}}{2} & 0\\\frac{i e^{i \phi}}{2} & \frac{e^{i \phi}}{2} & \frac{\sqrt{2} i}{2}\\- \frac{e^{i \phi}}{2} & \frac{i e^{i \phi}}{2} & \frac{\sqrt{2}}{2}\end{matrix}\right]$
[9]:
../_images/notebooks_Tutorial_13_1.svg

The syntax pcvl.P('phi') allows you to use parameters in the circuit, where you can assign a value or not. The behavior of the parameters of a circuit is similar to the case of the components.

For instance, you can use :

[10]:
params=circuit.get_parameters()
print(params) #list of the parameters

# the value is None, but we can change that with :

params[0].set_value(np.pi)
pcvl.pdisplay(circuit)

[Parameter(name='phi', value=None, min_v=0.0, max_v=6.283185307179586)]
[10]:
../_images/notebooks_Tutorial_15_1.svg

2. Mach-Zehnder Interferometers

The beamsplitter’s angle \(\theta\) can also be defined as a parameter.

However, as the reflexivity depends on the mirror, it’s hard to have adaptibility on the angle. Therefore, in practice, we use a Mach-Zehnder Interferometer.

The beamsplitter with a parameterised \(\theta\) is therefore implemented with a parameterised phase shifter \(\phi\) between two fixed beamsplitters.

[11]:
## TO-DO: build a circuit implementing the mzi

mzi = pcvl.Circuit(2) // BS() // (1,PS(phi=pcvl.P("phi"))) // BS()


## TO-DO: Check that the parameterised phase allows you to change the reflexivity of your MZI

import matplotlib.pyplot as plt

## We create a list of all different values for theta
X = np.linspace(0, 2*np.pi, 1000)
Y = []
for theta in X:
    phase = mzi.get_parameters()[0]
    phase.set_value(theta)
    Y.append(abs(mzi.compute_unitary()[0,0])**2)  #compute_unitary is numerical, so far faster that mzi.U, which uses symbolic expressions.

plt.plot(X, Y)
plt.xlabel("phi")
plt.ylabel("R")
plt.show()

pcvl.pdisplay(mzi)
## Note: If you need to create a BS directly from the reflexivity value, please use:
## BS(BS.r_to_theta(reflectivity_value))
## However, be aware that only theta value is stored inside the BS object
../_images/notebooks_Tutorial_17_0.png
[11]:
../_images/notebooks_Tutorial_17_1.svg

3. Universal Circuits

An operation on the modes of our circuit can also be expressed as a unitary.

For three modes, the unitary \(U=\begin{pmatrix} a_{1,1} & a_{1,2} & a_{1,3}\\ a_{2,1} & a_{2,2} & a_{2,3} \\ a_{3,1} & a_{3,2} & a_{3,3} \end{pmatrix}\) performs the following operation on the Fock state basis:

\[\begin{split}\begin{array}{rcl} |1,0,0\rangle & \mapsto& a_{1,1}|1,0,0\rangle + a_{1,2}|0,1,0\rangle + a_{1,3}|0,0,1\rangle\\ |0,1,0\rangle & \mapsto& a_{2,1}|1,0,0\rangle + a_{2,2}|0,1,0\rangle + a_{2,3}|0,0,1\rangle\\ |0,0,1\rangle & \mapsto& a_{3,1}|1,0,0\rangle + a_{3,2}|0,1,0\rangle + a_{3,3}|0,0,1\rangle \end{array}\end{split}\]

Since 1994, we know that any \(U\) on the modes can be implemented as an LO-circuit Reck’s et al.

This decomposition can be done easily in Perceval using beamsplitters and phase-shifters as follows.

[12]:
## From any unitary
n = 3
U = pcvl.Matrix.random_unitary(n)

circuit_u = pcvl.Circuit.decomposition(U, BS(theta=pcvl.P('theta'),phi_tr=pcvl.P('phi')), phase_shifter_fn=PS)

pcvl.pdisplay(circuit_u)
[12]:
../_images/notebooks_Tutorial_20_0.svg
[13]:
print("The error between the two unitaries is", np.linalg.norm(U-circuit_u.compute_unitary()))
The error between the two unitaries is 8.546891021709984e-09
[14]:
## TO-DO: decompose your unitary with only phase shifters and balanced beamsplitters.
mzi=pcvl.Circuit(2) // BS() // PS(pcvl.P("phi1")) //  BS()  //  PS(pcvl.P("phi2"))

circuit_u=pcvl.Circuit.decomposition(U,mzi,phase_shifter_fn=PS)

## Note: you can use a MZI. Be careful to put the phase on the right, as the full layer of phase_shifter_fn is on the left of the circuit

[15]:
## TO-DO: check the norm of the difference to be sure it has worked well
print("The error between the two unitaries is",np.linalg.norm(U-circuit_u.compute_unitary()))
The error between the two unitaries is 3.4586946619580422e-09

4. Black Box

To improve readibility, the circuit can be constructed in multiple steps which are then combined as black boxes. This will also help when we’ll need generic operations.

[16]:
pre_MZI = (pcvl.Circuit(4, name="Bell State Prepar.")
           .add(0, BS())
           .add(2, BS())
           .add(1, PERM([1, 0])))

upper_MZI = (pcvl.Circuit(2, name="upper MZI")
             .add(0, PS(phi=pcvl.P('phi_0')))
             .add(0, BS())
             .add(0, PS(phi=pcvl.P('phi_2')))
             .add(0, BS()))

lower_MZI = (pcvl.Circuit(2, name="lower MZI")
             .add(0, PS(phi=pcvl.P('phi_1')))
             .add(0, BS())
             .add(0, PS(phi=pcvl.P('phi_3')))
             .add(0, BS()))

chip = (pcvl.Circuit(4)
              .add(0, pre_MZI)
              .add(0, upper_MZI, merge=False)
              .add(2, lower_MZI, merge=False))

pcvl.pdisplay(chip)
[16]:
../_images/notebooks_Tutorial_25_0.svg
[17]:
## You can still display the inside of black boxes with:
pcvl.pdisplay(chip, recursive=True)
[17]:
../_images/notebooks_Tutorial_26_0.svg

III. Simulation

Up to this point, we have focused on creating circuits. It’s time to learn how to sample from them or describe their output distribution, on many different inputs.

1. Computing probabilities

For this part, we will take the Hong-Ou-Mandel experience as an example.

It’s one of the simplest experiments and yet it is very useful.

Making two indistinguishable photons, one in each mode, enter one balanced beamsplitter \(BS=\frac{1}{\sqrt{2}} \left[\begin{matrix}1 & 1\\1& -1\end{matrix}\right]\), we expect the outcome to be:

\[|1,1\rangle \mapsto \frac{|2,0\rangle - |0,2\rangle}{\sqrt{2}}\]

We will show how to verify this in the next steps using the Naive backend to recover the full probability distribution.

[18]:
## TO-DO: build the circuit with the convention above


circuit = pcvl.Circuit(2) // BS.H()



[19]:
# Syntax to compute the amplitudes
backend = pcvl.BackendFactory.get_backend("Naive")
backend.set_circuit(circuit)
backend.set_input_state(pcvl.BasicState([1,1]))
print(backend.prob_amplitude(pcvl.BasicState([2,0])))  #note that it's the amplitude !
print(backend.prob_amplitude(pcvl.BasicState([0,2])))
print(backend.probability(pcvl.BasicState([0,2])))


## We can also use the Analyser module to compute a table of probabilities
## The Analyser uses a Processor to work with. A Processor aims at simulating a photonic source plugged into a circuit
## with a given backend.
## The main syntax is :
## >>> p = pcvl.Processor(backend_name, circuit, source)
p = pcvl.Processor("Naive", BS())
analyzer = pcvl.algorithm.Analyzer(p, [pcvl.BasicState([1,1])], '*')
pcvl.pdisplay(analyzer)

(0.7071067811865477+0j)
(-0.7071067811865477+0j)
0.5000000000000002
|1,1>|0,2> |2,0>
|1,1> 01/2 1/2
[20]:
## TO-DO:  Choose a random unitary 3x3 U, and output the table probablities when the input |1,1,0> passes through the LO-Circuit of unitary U.

randU=pcvl.Unitary(pcvl.Matrix.random_unitary(3))
input=pcvl.BasicState([1,1,0])

p=pcvl.Processor("SLOS", randU)  #We can put in the processor a pcvl.Unitary instead of a circuit ! We don't need to use pcvl.decomposition
analyzer= pcvl.algorithm.Analyzer(p, [input], '*')
pcvl.pdisplay(analyzer)
|0,2,0> |0,1,1> |2,0,0> |1,1,0> |1,0,1> |0,0,2>
|1,1,0> 0.043518 0.302556 0.072713 0.54482 0.014496 0.021897

2. Sampling

Although it’s crucial to compute the output distribution, it’s not what we can expect from a photonic chip. Indeed, realistically, we only can obtain a single sample from the distribution each time we run the circuit. This can be done using the backend SLOS.

[21]:
p = pcvl.Processor("SLOS", BS())
p.with_input(pcvl.BasicState([1,1]))

# The sampler holds 'probs', 'sample_count' and 'samples' calls. You can use the one that fits your needs!
sampler = pcvl.algorithm.Sampler(p)

# A sampler call will return a Python dictionary containing sampling results, and two performance scores
# sample_count = sampler.sample_count(1000)
# sample_count contains {'results': <actual count>, 'physical_perf': float [0.0 - 1.0], 'logical_perf': float [0.0 - 1.0]}
sample_count = sampler.sample_count(1000)
print(sample_count['results'])

{
  |2,0>: 499
  |0,2>: 501
}
[22]:
## TO-DO: implement the code to sample from the 3x3 Unitary of earlier

p=pcvl.Processor("CliffordClifford2017",randU)
p.with_input(pcvl.BasicState([1,1,0]))
p.min_detected_photons_filter(0)  # Do not filter out any output state

sampler = pcvl.algorithm.Sampler(p)
sample_count = sampler.sample_count(1000)
print(sample_count['results'])

## Question: how many states do we have for 3 modes and 2 photons?
## There are 6 different states


## Question : how many states do we have for m modes and n photons?
## There are m+n-1 choose n different states. Cf Bar and Star problems.
{
  |1,1,0>: 552
  |2,0,0>: 73
  |0,1,1>: 301
  |0,0,2>: 20
  |0,2,0>: 44
  |1,0,1>: 10
}

Note : to approximate with decent precision a distribution over \(M\) different states, we would need \(M^2\) samples. This can be shown by Hoeffding’s inequality.

3. Performance and output state filtering

Perceval Processors have a built-in way of computing performance scores.

There are two different performance scores: * Physical performance * Logical performance

These performance scores help measure the real duration of a data acquisition on a real QPU.

a. Physical performance

This score is related to the number of detections (on a QPU: number of clicks). It drops output states where photons have been lost, or finish in the same mode.

For instance, an imperfect source makes this score drop.

However, you can choose not to filter any output state by lowering the expected clicks with: > proc.min_detected_photons_filter(0)

Processor.min_detected_photons_filter method

Perceval aims at being an interface for the QPU and as such, proc.min_detected_photons_filter(int k) post selects on having at least k photons detected (for threshold detection: photons on at least k different modes). By default, this value is set to n where n is the expected number of input photons. This is useful for retrieving a logical interpretation, making sure that no photon has been lost due to noise and coherent with the use of threshold detectors. However, for various applications (for instance machine learning where we use the full Fock space and resolve the number of photons, you will have to set it to 0 (and you may introduce you own post selection scheme if needed).

[23]:
# Create an empty circuit (each input mode is directly connected to a detector without interacting with any other)
empty_circuit = pcvl.Circuit(4)

perfect_proc = pcvl.Processor("Naive", empty_circuit)
imperfect_proc = pcvl.Processor("Naive", empty_circuit, pcvl.Source(emission_probability=0.3))

# Set the same input in both processors
input_state = pcvl.BasicState([1,0,1,0])
perfect_proc.with_input(input_state)
imperfect_proc.with_input(input_state)

perfect_sampler = pcvl.algorithm.Sampler(perfect_proc)
perfect_probs = perfect_sampler.probs()
imperfect_sampler = pcvl.algorithm.Sampler(imperfect_proc)
imperfect_probs = imperfect_sampler.probs()

print('Physical perf of perfect processor =', perfect_probs['physical_perf'])
print('Physical perf of imperfect processor =', imperfect_probs['physical_perf'])  # source emission probability**2

# You can still disable output state filtering
imperfect_proc.min_detected_photons_filter(0)
imperfect_probs = imperfect_sampler.probs()
print('Physical perf of imperfect processor (without selection) =', imperfect_probs['physical_perf'])
Physical perf of perfect processor = 1
Physical perf of imperfect processor = 0.08999999999999997
Physical perf of imperfect processor (without selection) = 1

b. Logical performance

This performance computation is set up by heralded modes and/or post-selection function set in a processor.

Depending on the circuit used, on the post-selection function, you may observe that physical and logical performance score interact. So, if you’re interested on a theoretical gate performance, you should disable physical post-selection with: > proc.min_detected_photons_filter(0)

Here is a quick example of the heralding / post-selection syntax in Perceval. You will see the result later on in this notebook.

[24]:
circuit = pcvl.Circuit(3) // BS() // (1, BS()) // BS()
p = pcvl.Processor("Naive", circuit)
p.add_herald(2,0)  # Third mode is heralded (0 photon in, 0 photon expected out)

# After a mode is heralded, you must not take it into account when setting an input to the processor
p.with_input(pcvl.BasicState([1, 0]))
sampler = pcvl.algorithm.Sampler(p)
probs = sampler.probs()
print("With herald only")
print("Logical perf =", probs['logical_perf'])
print(probs['results'])

# A post-selection function can be created like this:
postselect_func = pcvl.PostSelect("[1] == 1")  # meaning we required 1 photon detection in mode #1

p.set_postselection(postselect_func)  # Add post-selection
probs = sampler.probs()
print("With herald + post-selection function")
print("Logical perf =", probs['logical_perf'])
print(probs['results'])
With herald only
Logical perf = 0.75
{
  |1,0>: 0.02859547920896832
  |0,1>: 0.9714045207910317
}
With herald + post-selection function
Logical perf = 0.7285533905932737
{
  |0,1>: 1.0
}

4. Variational algorithm

In variational algorithms, the samples from a quantum circuit allow us to approximate an expectation value, which is then used to determine the value of a loss function. This loss function is chosen such that minimising it yields a solution to a given problem. By changing the values of the parameters in our quantum circuit, we can search for this minimum.

We won’t go into the details of variational algorithms. However, it may be useful to see how to perform an optimisation with Perceval.

We will use the library scipy.optimise.

The following code solves the problem of finding an LO-Circuit which, given a Fock State \(|1,1,1,1\rangle\), maximises the probability of outputting \(|4,0,0,0\rangle\). The solution below works for an arbitrary \(n\).

[3]:
import random
from scipy import optimize

# Data
n = 4
input = pcvl.BasicState([1]*n)
output_to_max = pcvl.BasicState([n]+[0]*(n-1))
backend = pcvl.BackendFactory.get_backend("SLOS")

# TO-DO: implement a generic circuit of size n with parameters. Code the loss function to maximise the good output. Launch the optimisation procedure. Output the probability and circuit obtained

# We take a universal circuit
circuit = pcvl.GenericInterferometer(n,
    lambda i: BS(theta=pcvl.P(f"theta{i}"),
    phi_tr=pcvl.P(f"phi_tr{i}")),
    phase_shifter_fun_gen=lambda i: PS(phi=pcvl.P(f"phi{i}")))
param_circuit = circuit.get_parameters()
params_init = [random.random()*np.pi for _ in param_circuit]


def loss_function(params):
    for i, value in enumerate(params):
        param_circuit[i].set_value(value)
    backend.set_circuit(circuit)
    backend.set_input_state(input)
    return -backend.probability(output_to_max)  # we want to maximise the prob, so we want to minimise the -prob


# We run the otpimisation
o = optimize.minimize(loss_function, params_init, method="Powell")

print(f"The maximum probability is {-loss_function(o.x)}")

# For n=4, the probability should be 3/32
# The maximum can also be obtained with the Hadamard matrix :

H4 = (1/2)*np.array([[1,1,1,1], [1,-1,1,-1], [1,1,-1,-1], [1,-1,-1,1]])
backend.set_circuit(pcvl.Unitary(pcvl.Matrix(H4)))
backend.set_input_state(input)
backend.probability(output_to_max)
The maximum probability is 0.09374999914144536
[3]:
0.09374999999999999

5. To go further : connect to a chip

Perceval is also connected to real/physical chips. Here’s the syntax to sample directly from them !

doc to connect to a chip

IV. Encoding Qubits

1. Path encoding

To perform quantum computations using photons, we need an encoding: a correspondance between our Fock states and our qubit states.

We therefore want to associate each qubit state with one of our Fock states.

One natural way to encode qubits is the path encoding. A qubit is a two-level quantum state, so we will use two spatial modes to encode it: this is the dual-rail or path encoding.

The logical qubit state \(|0\rangle_L\) will correspond to a photon in the upper mode, as in the Fock state \(|1,0\rangle\), while \(|1\rangle_L\) will be encoded as \(|0,1\rangle\).

We can extend this to multiple qubits by having twice as many modes as there are qubits. For example the \(3\)-qubit state \(\frac{1}{\sqrt{2}}(|000\rangle_L+|111\rangle_L)\) can be encoded with \(3\) photons and \(3\times 2=6\) modes : \(\frac{1}{\sqrt{2}}(|1,0,1,0,1,0\rangle+|0,1,0,1,0,1\rangle)\)

2. Single-qubit gates

Using the dual-rail enconding, single-qubit gates only deal with one photon and are straightforward. Can you give the LO-circuits for the gates below?

\[\begin{split}X=\left[\begin{matrix}0 & 1\\1& 0\end{matrix}\right]\end{split}\]
\[\begin{split}Y=\left[\begin{matrix}0 & -i\\i& 0\end{matrix}\right]\end{split}\]
\[\begin{split}Z=\left[\begin{matrix}1 & 0\\0& -1\end{matrix}\right]\end{split}\]
\[\begin{split}H=\frac{1}{\sqrt{2}} \left[\begin{matrix}1 & 1\\1& -1\end{matrix}\right]\end{split}\]
\[\begin{split}R_X=\left[\begin{matrix}\cos{\left(\frac{\theta}{2} \right)} & -i \sin{\left(\frac{\theta}{2} \right)}\\-i \sin{\left(\frac{\theta}{2} \right)} & \cos{\left(\frac{\theta}{2} \right)}\end{matrix}\right]\end{split}\]
\[\begin{split}R_Y=\left[\begin{matrix}\cos{\left(\frac{\theta}{2} \right)} & - \sin{\left(\frac{\theta}{2} \right)}\\ \sin{\left(\frac{\theta}{2} \right)} & \cos{\left(\frac{\theta}{2} \right)}\end{matrix}\right]\end{split}\]
\[\begin{split}R_Z=\left[\begin{matrix}e^{-i\frac{\theta}{2}} & 0 \\ 0 & e^{i\frac{\theta}{2}}\end{matrix}\right]\end{split}\]
[26]:
## TO-DO: find the LO-circuits for each gate

circuit_x = PERM([1,0])  #it's not the only way
circuit_y = PERM([1,0]) // (0,PS(-np.pi/2)) // (1,PS(np.pi/2))
circuit_z = pcvl.Circuit(2) // (1,PS(np.pi))
circuit_h = BS.H()

circuit_rx = pcvl.Circuit(2) // (0, PS(np.pi)) // BS.Rx(theta=pcvl.P("theta")) // (0, PS(np.pi)) #Be careful for the minus ! We use a convention
circuit_ry = BS.Ry(theta=pcvl.P("theta"))
circuit_rz = BS.H() // circuit_rx // BS.H()  # Indeed, Rz = H Rx H. Of course, we would like to be able to have many parameters with the same name. It will come soon :)

3. Two-qubit gates

On the other hand, in dual-rail encoding, it can be shown that two-qubit gates can’t be deterministic, and have a probability to fail.

There are two ways to detect that failure:

  • We can use additional photons called ancillas, which we can measure independently from the main circuit photons. Depending on the state obtained on the ancilla, we know whether the gate has succeeded or not on the main qubits. Those gates will be called heralded.

  • We can also directly measure the main circuit qubits, and depending on the result, assess whether the gate has succeeded or not. Those gates will be called postselected.

The CNOT gate acts on two qubits, a control and a target, and flips the value of the target if the control qubit is in state \(|1\rangle_L\). In the following two exercices, we will see the two types of CNOT gates: - the postselected CNOT of Ralph et al. - the heralded CNOT of the KLM protocol

[27]:
## We introduce the component catalog. It contains both CNOT gates.
from perceval.components import catalog
print(catalog.list())
['klm cnot', 'postprocessed cnot', 'heralded cz', 'generic 2 mode circuit', 'mzi phase first', 'mzi phase last']
[28]:
## Ralph's et al. CNot

print(catalog['postprocessed cnot'].doc)
ralph_cnot = catalog['postprocessed cnot'].build_processor()
## You can set its input state with a LogicalState
ralph_cnot.with_input(pcvl.LogicalState([1, 0]))

pcvl.pdisplay(ralph_cnot, recursive=True, render_size=1.25)
POSTPROCESSED CNOT DOCUMENTATION
---------------------------------

CNOT gate with 2 heralded modes and a post-selection function

Scientific article reference: https://journals.aps.org/pra/abstract/10.1103/PhysRevA.65.062324

Schema:
                      ╭─────╮
ctrl (dual rail) ─────┤     ├───── ctrl (dual rail)
                 ─────┤     ├─────
                      │     │
data (dual rail) ─────┤     ├───── data (dual rail)
                 ─────┤     ├─────
                      ╰─────╯
[28]:
../_images/notebooks_Tutorial_50_1.svg
[29]:
## TO-DO: Check/Convince yourself that the circuit above is performing a CNOT in the dual rail encoding
[30]:
## You can sample some output states
ralph_cnot.min_detected_photons_filter(0)
cnot_sampler = pcvl.algorithm.Sampler(ralph_cnot)
samples = cnot_sampler.probs()
print(samples['results'])
print("Some output states were not selected because of heralds and post-processing => you can check the logical performance")
print("Logical performance = ", samples['logical_perf'])
{
  |0,1,0,1>: 1.0
}
Some output states were not selected because of heralds and post-processing => you can check the logical performance
Logical performance =  0.11111111111111131
[31]:
## You can sample some output states
h_cnot = catalog['klm cnot'].build_processor()
cnot_sampler = pcvl.algorithm.Sampler(h_cnot)
h_cnot.with_input(pcvl.LogicalState([0, 0]))

samples = cnot_sampler.samples(10)
print(samples['results'])
print("Some output states were not selected because of heralds and post-processing => you can check the logical performance")
print("Logical performance =", samples['logical_perf'])
[|1,0,1,0>, |1,0,1,0>, |1,0,1,0>, |1,0,1,0>, |1,0,1,0>, |1,0,1,0>, |1,0,1,0>, |1,0,1,0>, |1,0,1,0>, |1,0,1,0>]
Some output states were not selected because of heralds and post-processing => you can check the logical performance
Logical performance = 0.05156442163200023
[32]:
## TO-DO: Check it perfoms a CNOT, and explicit the difference between the two types of CNOT

Exercise

The next circuit comes from the following paper.

image0

[33]:
## KLM's CNOT

print(catalog['klm cnot'].doc)
knill_cnot = catalog['klm cnot'].build_processor()

pcvl.pdisplay(knill_cnot, recursive=True)
HERALDED CNOT DOCUMENTATION
----------------------------

CNOT gate with 4 heralded modes

Scientific article reference: https://doi.org/10.1073/pnas.1018839108

Schema:
                      ╭─────╮
ctrl (dual rail) ─────┤     ├───── ctrl (dual rail)
                 ─────┤     ├─────
                      │     │
data (dual rail) ─────┤     ├───── data (dual rail)
                 ─────┤     ├─────
                      ╰─────╯
[33]:
../_images/notebooks_Tutorial_57_1.svg
[34]:
# TO-DO : reproduce it in the encoding seen above

# Let's try to implement that circuit properly.
# First, the quantum gates, as coded above :

Rx = lambda i: pcvl.Circuit(2) // (0,PS(np.pi)) // BS.Rx(theta=pcvl.P(f"theta{i}")) // (0,PS(np.pi)) #Be careful for the minus ! We use a convention
Ry = lambda i: pcvl.Circuit(2,name=f"Ry{i}") // BS.Ry(theta=pcvl.P(f"theta{i}"))
Rz = lambda i: pcvl.Circuit(2,name=f"Rz{i}") // BS.H() // circuit_rx // BS.H()
cnot = catalog['klm cnot'].build_processor()

# Our qubits in the dual rail encoding
q1, q2, q3 = [0,1], [2,3], [4,5]

p = pcvl.Processor("SLOS",6)

for i in range(3):
    p.add(2*i,Ry(i+1)).add(2*i,Rz(i+4))
p.add(q1+q2, cnot)
p.add(q1+q3, cnot)
p.add(q2+q3, cnot)

for i in range(3):
    p.add(2*i, Ry(i+7)).add(2*i, Rz(i+10))

pcvl.pdisplay(p,recursive=False)
[34]:
../_images/notebooks_Tutorial_58_0.svg