Forcepho infers the fluxes and shapes of galaxies from astronomical images. It models the appearance of multiple sources in multiple bands simultaneously and compares to observed data via a likelihood function. Gradients of this likelihood allow for efficient maximization of the posterior probability or sampling of the posterior probability distribution via Hamiltonian Monte Carlo. The model intrinsic galaxy shapes and positions are shared across the different bands, but the fluxes are fit separately for each band. Forcepho does not perform detection; initial locations and (very rough) parameter estimates must be supplied by the user.