  We propose a framework to design and optimize generic photonic quantum
circuits composed of Gaussian objects (pure and mixed Gaussian states, Gaussian
unitaries, Gaussian channels, Gaussian measurements) as well as non-Gaussian
effects such as photon-number-resolving measurements. In this framework, we
parametrize a phase space representation of Gaussian objects using elements of
the symplectic group (or the unitary or orthogonal group in special cases), and
then we transform it into the Fock representation using a single linear
recurrence relation that computes the Fock amplitudes of any Gaussian object
recursively. We also compute the gradient of the Fock amplitudes with respect
to phase space parameters by differentiating through the recurrence relation.
We can then use Riemannian optimization on the symplectic group to optimize
M-mode Gaussian objects, avoiding the need to commit to particular realizations
in terms of fundamental gates. This allows us to "mod out" all the different
gate-level implementations of the same circuit, which now can be chosen after
the optimization has completed. This can be especially useful when looking to
answer general questions, such as bounding the value of a property over a class
of states or transformations, or when one would like to worry about hardware
constraints separately from the circuit optimization step. Finally, we make our
framework extendable to non-Gaussian objects that can be written as linear
combinations of Gaussian ones, by explicitly computing the change in global
phase when states undergo Gaussian transformations. We implemented all of these
methods in the freely available open-source library MrMustard, which we use in
three examples to optimize the 216-mode interferometer in Borealis, and 2- and
3-modes circuits (with Fock measurements) to produce cat states and cubic phase
states.
