  Gravitational lensing and stellar dynamics are two independent methods, based
solely on gravity, to study the mass distributions of galaxies. Both methods
suffer from degeneracies, however, that are difficult to break. In this paper,
we present a new framework that self-consistently unifies gravitational lensing
and stellar dynamics. This approach breaks some of classical degeneracies that
have limited their individual usage, in particular in the study of
high-redshift galaxies. The methodology is based on the premise that, for any
given galaxy potential, the mapping of both the unknown lensed source
brightness distribution and the stellar phase-space distribution function on to
the photometric and kinematic observables, can be cast as a single set of
coupled linear equations. This set of linear equations is solved, maximizing
the likelihood penalty function. The evidence penalty function, as derived from
Bayesian statistics, subsequently allows the best potential-model parameters to
be found and potential-model families, or other model assumptions (e.g. PSF),
to be quantitatively ranked. We have implemented a fast algorithm that solves
for the maximum-likelihood pixelized lensed source brightness distribution and
the two-integral stellar phase-space distribution function f(E, L_z), assuming
axisymmetric potentials. To make the method practical, we have devised a new
Monte-Carlo approach to Schwarzschild's orbital superposition method, based on
the superposition of two-integral (E and L_z) toroidal components, to find the
maximum-likelihood two-integral distribution function in a matter of seconds in
any axisymmetric potential. The non-linear parameters of the potential are
subsequently found through a hybrid MCMC and Simplex optimization of the
evidence. (Abridged)
