We introduce a general phase-space representation that includes global quantum symmetries in the basis expansion. This method, called matrix phase-space, projects the basis onto a reduced Hilbert space, which can greatly reduce sampling errors of many-body quantum simulations when there is a symmetry, and unifies several previous phase-space methods. The purpose of this paper is to provide detailed proofs of basic theorems and operator identities. We also treat several different types of symmetries, corresponding to parity, number and momentum conservation laws. To illustrate these methods, we give a detailed derivation of previous applications to the topical problem of validating the outputs of Gaussian boson sampling (GBS) quantum computers (QC) with photon-number-resolving (PNR) detectors. Exact evaluation of GBS PNR output probabilities has exponential classical complexity. Fortunately, special cases of output GBS grouped probabilities have tractable solutions, providing a way to check the results. Using matrix phase-space with parity symmetry, we observe a reduction of over two orders of magnitude in sampling errors compared with earlier methods in low-loss cases. Additional nonlinear examples are also treated. These include the dynamical evolution of a nonlinear damped one-dimensional lattice, and the periodic Bose-Hubbard model. More details will be given elsewhere. As well as direct application to experiments, this scalable first-principles method may be useful for validating other QC calculations on large Hilbert spaces.