# Boson Sampling with MPS

In this notebook, we explain how to use the MPS (Matrix Product State) backend to simulate a linear circuit. MPS simulation is based on a type of tensor network simulation, which gives an approximation of the output states [1] [2]. It does the computation on each component of the circuits one-by-one, and not on the whole unitary. The states are represented by tensors, which are then updated at each component. These tensors can be seen as a big set of matrices,
and the approximation is done by chosing the dimension of these matrices, called the *bond dimension*. For this example, we simulate a simple boson sampling problem, with 6 modes and 3 photons [3].

```
[1]:
```

```
import perceval as pcvl
import matplotlib.pyplot as plt
```

## Definition of the circuit

Just like in the Boson Sampling notebook, we generate a Haar-random unitary and its decomposition in a circuit :

```
[2]:
```

```
m = 6
unitary = pcvl.Matrix.random_unitary(m)
mzi = (pcvl.BS() // (0, pcvl.PS(phi=pcvl.Parameter("φ_a")))
// pcvl.BS() // (1, pcvl.PS(phi=pcvl.Parameter("φ_b"))))
linear_circuit = pcvl.Circuit.decomposition(unitary, mzi,
phase_shifter_fn=pcvl.PS,
shape="triangle")
print(m, " modes triangular Boson Sampler :")
pcvl.pdisplay(linear_circuit, compact = True)
```

```
6 modes triangular Boson Sampler :
```

```
[2]:
```

## MPS simulation

Let us now define the MPS simulation, using the MPS backend in Perceval.

```
[3]:
```

```
mps = pcvl.MPSBackend()
mps.set_circuit(linear_circuit)
```

We now choose the size of the matrices (the bond dimension) for our simulation.

```
[4]:
```

```
chi = 8
mps.set_cutoff(chi)
```

And finally, we get the output probability distribution from a given input state with 3 photons.

```
[5]:
```

```
n = 3
input_state = pcvl.BasicState([1]*n + [0]*(m-n))
mps.set_input_state(input_state)
probs = mps.prob_distribution()
pcvl.pdisplay(probs, max_v=20)
```

state | probability |
---|---|

|0,0,1,0,2,0> | 0.089692 |

|1,1,1,0,0,0> | 0.074157 |

|0,1,0,0,0,2> | 0.066353 |

|0,0,1,0,1,1> | 0.057677 |

|0,1,1,0,0,1> | 0.057048 |

|0,0,0,0,1,2> | 0.052347 |

|0,1,1,0,1,0> | 0.048586 |

|0,1,2,0,0,0> | 0.039534 |

|1,1,0,0,0,1> | 0.034403 |

|2,0,1,0,0,0> | 0.034118 |

|1,0,2,0,0,0> | 0.029711 |

|0,0,0,1,2,0> | 0.029622 |

|0,0,0,0,3,0> | 0.027699 |

|0,0,2,0,1,0> | 0.026586 |

|0,0,0,1,1,1> | 0.022117 |

|2,0,0,0,0,1> | 0.021964 |

|0,1,0,1,0,1> | 0.020708 |

|1,1,0,1,0,0> | 0.020591 |

|1,0,0,0,0,2> | 0.018927 |

|0,0,2,0,0,1> | 0.01776 |

## Certification of the MPS method :

As we make an approximation by chosing the bond dimension, we have to check when does this approximation becomes good enough. Unfortunately, there is no formula giving the minimal size for a given approximation error. What we can do though is to compute the *Total Variance Distance* (TVD) between an ideal simulation of Boson Sampling, and an approximated one. To compute the ideal one, we can for instance use the *SLOS* backend on Perceval :

```
[6]:
```

```
slos = pcvl.SLOSBackend()
slos.set_circuit(pcvl.Unitary(unitary))
slos.set_input_state(input_state)
probs_slos = slos.prob_distribution()
pcvl.pdisplay(probs_slos, max_v=20)
```

state | probability |
---|---|

|0,0,1,0,2,0> | 0.089692 |

|1,1,1,0,0,0> | 0.074157 |

|0,1,0,0,0,2> | 0.066353 |

|0,0,1,0,1,1> | 0.057677 |

|0,1,1,0,0,1> | 0.057048 |

|0,0,0,0,1,2> | 0.052347 |

|0,1,1,0,1,0> | 0.048586 |

|0,1,2,0,0,0> | 0.039534 |

|1,1,0,0,0,1> | 0.034403 |

|2,0,1,0,0,0> | 0.034118 |

|1,0,2,0,0,0> | 0.029711 |

|0,0,0,1,2,0> | 0.029622 |

|0,0,0,0,3,0> | 0.027699 |

|0,0,2,0,1,0> | 0.026586 |

|0,0,0,1,1,1> | 0.022117 |

|2,0,0,0,0,1> | 0.021964 |

|0,1,0,1,0,1> | 0.020708 |

|1,1,0,1,0,0> | 0.020591 |

|1,0,0,0,0,2> | 0.018927 |

|0,0,2,0,0,1> | 0.01776 |

We also have to define the TVD function.

```
[7]:
```

```
def tvd(probs1, probs2):
tvd = 0
for state, prob in probs1.items():
tvd += abs(prob - probs2[state])
return tvd
```

Now we compute the TVD between the two simulations for different bond dimensions, going from 1 to 10.

```
[8]:
```

```
TVD = []
for chi in range(1,11):
mps.set_cutoff(chi)
mps.set_input_state(input_state)
probs_mps = mps.prob_distribution()
TVD.append(tvd(probs_mps, probs_slos))
plt.plot(TVD)
plt.xlabel('Bond dimension')
plt.ylabel('TVD')
plt.title('Evolution of the TVD between SLOS and MPS probability distributions');
```

We can see that the TVD decreases as the size of the matrices increases, untill reaching 0 for a bond dimension of 7.

## References

[1] Ulrich Schollwöck. The density-matrix renormalization group in the age of matrix product states. Annals of Physics, 326(1):96–192, jan 2011.

[2] Changhun Oh, Kyungjoo Noh, Bill Fefferman, and Liang Jiang. Classical simulation of lossy boson sampling using matrix product operators. Physical Review A, 104(2), aug 2021.

[3] Hui Wang, et al. Boson Sampling with 20 Input Photons and a 60-Mode Interferometer in a \(10^{14}\)-Dimensional Hilbert Space. Physical Review Letters, 123(25):250503, December 2019. Publisher: American Physical Society.