We introduce a new class of quantum Monte Carlo methods, based on a Gaussian quantum operator representation of fermionic states. The methods enable first-principles dynamical or equilibrium calculations in many-body Fermi systems, and, combined with the existing Gaussian representation for bosons, provide a unified method of simulating Bose-Fermi systems. As an application, we calculate finite-temperature properties of the two dimensional Hubbard model.Calculating the quantum many-body physics of interacting Fermi systems is one of the great challenges in modern theoretical physics. These issues appear in physical problems at all energy scales, from ultra-cold atomic physics to high-energy lattice QCD. In even the simplest cases, first-principles calculations are inhibited by the complexity of the fermionic wavefunction, manifest notoriously in the Fermi sign problem. In previous quantum Monte Carlo (QMC) techniques, the sign problem appears as trajectories with negative weights, which contribute to a large sampling error [1]. QMC methods are also complicated by the calculation of large determinants.In this letter, we introduce a new QMC method for simulating many-body fermion systems, based on a Gaussian phase-space representation. As an application to condensed matter and AMO physics, we study the well-known Hubbard model. Although it is the simplest model of interacting fermions on a lattice, it is rich in physics and may even describe high-temperature superconductivity [2]. We show that for the Hubbard model the Gaussian representation leads to imaginarytime equations with no negative probabilities or weights. We demonstrate that this removes the well-known Fermi sign problem [1,2,3], by first principles numerical simulation without fixed-node[4] or variational approximations.Phase-space methods[5] provide a way to simulate quantum many-body systems both dynamically and at finite temperature, and have proved useful in bosonic cases. These methods sample the time evolution of a positive distribution on an overcomplete basis set, which is usually the set of coherent states. However, whereas coherent state representations are well-defined in the bosonic case, the only known coherent state techniques for fermions involve Grassmann algebra [6], which has an enormous computational complexity.Here we introduce a phase-space method that overcomes the problem of Grassmann complexity, using a Gaussian expansion for fermions. The operator basis is constructed from pairs of Fermi operators. Because these pairs obey commutation relations, a natural solution of the Grassmann problem is achieved. Furthermore, the resulting equations obviate the need to evaluate large determinants. The elimination of anti-commutators means that the technique is far more efficient than previous QMC and stochastic fermion methods [7]. We give examples in cases of experimental relevance involving the dynamical problem of Pauli blocking in molecular dissociation, and finite temperature correlations of fermions in an optical lattice, where the ...