Least-squares model optimization methods are commonly used to estimate physical media properties by fitting geophysical data with nonlinear models. I extend this formulation to joint estimation of physical properties and lithological description of the media. Incorporation of petrophysical information within the inversion scheme provides the coupling between lithology and media physics by describing the geostatistical relation between them. The resulting procedure adjusts iteratively the joint model to simultaneously fit geophysical data, the petrophysical statistical medium description, and prior information on the lithology, following equations derived for the Newton's optimization method. Although more calculations are required to incorporate the additional information and estimate the model update, the algebraic system of linearized equations can be transformed appropriately to remain within the same dimensions of the conventional inverse formulation. In the particular case when the petrophysical transform is linear (i.e., the function of lithological parameters that provides the expected values of the physical parameters), the lithological inversion equations are equivalent to the corresponding equations of a conventional inversion followed with the inverse petrophysical transform. I illustrate the methodology with synthetic examples of porosity-impedance estimation from zero offset seismic data, using Wyllie's transform to construct the statistical relationship between porosity and impedance. When the porosity range of the test models is on the nonlinear part of the petrophysical transform, the lithological inversion performs significantly better than the conventional inversion. When the porosity range is on an almost linear part of the transform, the performances are equivalent for both approaches.