The full polarization property of volume holographic grating diffraction is investigated theoretically. With a simple volume grating model, the diffracted fields and Mueller matrices are first derived from Maxwell's equations by using the Green's function algorithms. The formalism is derived for the general case that the diffraction beam and the grating wave vector are not in the plane of incidence, where s waves and p waves are not decoupled. The derived photon-momentum relations determine the Bragg angle selectivity. The parameters of diffraction strength related to the hologram-writing process and material are defined and are not necessarily small in general. The diffracted-beam profiles are analytically calculated by using the known grating shape function. This theory has provided a fundamental understanding of the polarization phenomena of a real holographic diffraction grating device. The derived algorithm would provide a simulation-analysis tool for the engineering design of real holographic beam combiner/splitter devices.