  Linear optical quantum circuits with photon number resolving (PNR) detectors
are used for both Gaussian Boson Sampling (GBS) and for the preparation of
non-Gaussian states such as Gottesman-Kitaev-Preskill (GKP), cat and NOON
states. They are crucial in many schemes of quantum computing and quantum
metrology. Classically optimizing circuits with PNR detectors is challenging
due to their exponentially large Hilbert space, and quadratically more
challenging in the presence of decoherence as state vectors are replaced by
density matrices. To tackle this problem, we introduce a family of algorithms
that calculate detection probabilities, conditional states (as well as their
gradients with respect to circuit parametrizations) with a complexity that is
comparable to the noiseless case. As a consequence we can simulate and optimize
circuits with twice the number of modes as we could before, using the same
resources. More precisely, for an $M$-mode noisy circuit with detected modes
$D$ and undetected modes $U$, the complexity of our algorithm is $O(M^2
\prod_{i\in U} C_i^2 \prod_{i\in D} C_i)$, rather than $O(M^2 \prod_{i \in
D\cup U} C_i^2)$, where $C_i$ is the Fock cutoff of mode $i$. As a particular
case, our approach offers a full quadratic speedup for calculating detection
probabilities, as in that case all modes are detected. Finally, these
algorithms are implemented and ready to use in the open-source photonic
optimization library MrMustard.
