Pure Inverse Problems: Maximum Entropy and Cross-Entropy

Ganbaatar Jambal

2026-08-30

1 The pure inverse problem

Suppose we face the following inverse problem

\[\begin{equation} \mathbf{y} = \mathbf{X}\mathbf{p} \tag{1.1} \end{equation}\]

where \(\mathbf{y} = (y_1, \dots, y_T)'\) is a \(T\)-dimensional vector of observations (moments), \(\mathbf{X}\) is a known \((T \times K)\) matrix, and \(\mathbf{p} = (p_1, \dots, p_K)'\) is a \(K\)-dimensional vector of unknown probabilities. The goal is to recover \(\mathbf{p}\), but in many cases this task is not directly feasible — for example, when \(\mathbf{X}\) has no inverse. A classical example is the under-determined case in which the number of observations is smaller than the number of unknowns, \(T < K\). There are then infinitely many \(\mathbf{p}\) consistent with (1.1), so we need a criterion to select among them.

Throughout, inverse_ce() solves exactly this problem through its concentrated (dual) formulation.

library(infometrics)

2 Maximum entropy

Following Shannon (1948), (Jaynes, 1957a, 1957b), and Levine (1980), one natural criterion is to choose the \(\mathbf{p}\) that maximises the Shannon entropy

\[\begin{equation} H(\mathbf{p}) = -\sum_{k=1}^K p_k \log(p_k), \tag{2.1} \end{equation}\]

with the convention \(p_k \log(p_k) = 0\) whenever \(p_k = 0\). Entropy (2.1) is maximised when \(\mathbf{p}\) is uniform, i.e. \(p_k = 1/K\) for all \(k\), so maximising it selects the least committal distribution consistent with the data. We also impose the additivity (normalization) constraint

\[\begin{equation} \sum_{k=1}^K p_k = 1. \tag{2.2} \end{equation}\]

We therefore maximise (2.1) subject to (1.1) and (2.2). Writing the moment residual as \(\big(\sum_k p_k x_{tk} - y_t\big)\), the Lagrangian is

\[\begin{equation} \mathcal{L} = -\sum_{k=1}^K p_k \log(p_k) + \sum_{t=1}^T \lambda_t\!\left(\sum_{k=1}^K p_k x_{tk} - y_t\right) + \mu\!\left(1 - \sum_{k=1}^K p_k\right), \tag{2.3} \end{equation}\]

where the \(\lambda_t\) and \(\mu\) are Lagrange multipliers. The multiplier \(\lambda_t\) measures the information content of moment \(t\): if \(\lambda_t = 0\), the \(t\)-th moment carries no information and dropping it would not change the solution or the level of remaining uncertainty. The first-order conditions are

\[\begin{equation} \begin{aligned} \frac{\partial \mathcal{L}}{\partial p_k} &= -\log(p_k) - 1 + \sum_{t=1}^T \lambda_t x_{tk} - \mu = 0, & k &= 1,\dots,K,\\ \frac{\partial \mathcal{L}}{\partial \lambda_t} &= \sum_{k=1}^K p_k x_{tk} - y_t = 0, & t &= 1,\dots,T,\\ \frac{\partial \mathcal{L}}{\partial \mu} &= 1 - \sum_{k=1}^K p_k = 0. \end{aligned} \tag{2.4} \end{equation}\]

Solving the first condition and imposing (2.2) gives the exponential-family solution

\[\begin{equation} \hat{p}_k = \frac{\exp\!\big(\sum_t \hat\lambda_t x_{tk}\big)} {\sum_{k'} \exp\!\big(\sum_t \hat\lambda_t x_{tk'}\big)} = \frac{\exp\!\big(\sum_t \hat\lambda_t x_{tk}\big)}{\Omega(\hat{\boldsymbol\lambda})}, \tag{2.5} \end{equation}\]

where \(\Omega(\boldsymbol\lambda) = \sum_k \exp\!\big(\sum_t \lambda_t x_{tk}\big)\) is the normalization factor, usually called the partition function. Substituting (2.5) back into (2.3) — the normalization is already satisfied, so the \(\mu\) term drops — yields the concentrated (dual) model

\[\begin{equation} M(\boldsymbol\lambda) = -\sum_t \lambda_t y_t + \log \Omega(\boldsymbol\lambda), \tag{2.6} \end{equation}\]

which is minimised over the \(T\)-dimensional vector \(\boldsymbol\lambda\). This is exactly the objective inverse_ce() solves, so the reported lambda_hat and objective correspond one-to-one with \(\hat{\boldsymbol\lambda}\) and \(M(\hat{\boldsymbol\lambda})\) here. Working in \(\boldsymbol\lambda\)-space (dimension \(T\)) rather than \(\mathbf{p}\)-space (dimension \(K \gg T\)) is the key computational advantage of the dual (Golan, 2008).

2.1 Solving with inverse_ce()

Consider a six-sided die for which we know only that the mean roll is \(4.5\) (rather than the fair \(3.5\)). Which distribution over the faces \(\{1,\dots,6\}\) is most consistent with that single moment while assuming as little as possible? In inverse_ce()’s formula interface the response is the moment vector and each right-hand-side term is a state (a face); - 1 drops the intercept.

dice <- data.frame(y = 4.5,
                   s1 = 1, s2 = 2, s3 = 3, s4 = 4, s5 = 5, s6 = 6)
fit_me <- inverse_ce(y ~ s1 + s2 + s3 + s4 + s5 + s6 - 1, data = dice)
round(coef(fit_me), 4)          # p_hat over the six faces
#>     s1     s2     s3     s4     s5     s6 
#> 0.0544 0.0788 0.1142 0.1654 0.2398 0.3475

The moment condition (1.1) is satisfied exactly, and the recovered distribution is the flattest one compatible with a mean of \(4.5\):

sum(1:6 * coef(fit_me))          # reproduces the mean 4.5
#> [1] 4.5
c(H_phat = shannon_entropy(coef(fit_me)), H_max = log(6))
#>   H_phat    H_max 
#> 1.613581 1.791759

3 Cross-entropy: incorporating a prior

If we hold prior beliefs about the data-generating process, summarised by a distribution \(\mathbf{p}^0 = (p_1^0, \dots, p_K^0)'\), we can incorporate them through the cross-entropy (Kullback–Leibler) divergence (Csiszar, 1991; Kullback, 1974; Levine, 1980; Shore & Johnson, 1980)

\[\begin{equation} D(\mathbf{p} \,\Vert\, \mathbf{p}^0) = \sum_{k=1}^K p_k \log(p_k / p_k^0), \tag{3.1} \end{equation}\]

which measures the divergence of \(\mathbf{p}\) from the prior \(\mathbf{p}^0\). We now minimise (3.1) subject to the same constraints (1.1) and (2.2). The analogous derivation gives

\[\begin{equation} \hat{p}_k = \frac{p_k^0 \exp\!\big(\sum_t \hat\lambda_t x_{tk}\big)} {\sum_{k'} p_{k'}^0 \exp\!\big(\sum_t \hat\lambda_t x_{tk'}\big)} = \frac{p_k^0 \exp\!\big(\sum_t \hat\lambda_t x_{tk}\big)}{\Omega(\hat{\boldsymbol\lambda})}, \tag{3.2} \end{equation}\]

with \(\Omega(\boldsymbol\lambda) = \sum_k p_k^0 \exp\!\big(\sum_t \lambda_t x_{tk}\big)\). Maximum entropy is the special case of (3.2) with a uniform prior \(p_k^0 = 1/K\): the constant \(1/K\) cancels between numerator and denominator and (3.2) collapses to (2.5). The concentrated dual has the same form as (2.6), only the partition function changes, so inverse_ce() needs a single implementation for both: uniform p0 gives ME, a non-uniform p0 gives CE (hence the name).

Since the divergence (3.1) is non-negative, a strictly positive value means the data carried information beyond the prior. Following (Zellner, 1988, 1991), (3.2) can be read as an efficient information-processing rule that combines prior and data into a posterior.

3.1 The ME special case, and a genuine prior

First, the uniform-prior identity — passing an explicit uniform p0 reproduces the maximum-entropy fit:

fit_unif <- inverse_ce(y ~ s1 + s2 + s3 + s4 + s5 + s6 - 1,
                       data = dice, p0 = rep(1 / 6, 6))
max(abs(coef(fit_unif) - coef(fit_me)))     # ~ 0: ME == CE(uniform prior)
#> [1] 0

Now a genuine prior that leans toward high faces. Cross-entropy pulls the fit toward p0 while still matching the observed mean of \(4.5\):

p0_load <- c(.05, .05, .10, .15, .25, .40)  # prior beliefs favouring high faces
fit_ce  <- inverse_ce(y ~ s1 + s2 + s3 + s4 + s5 + s6 - 1,
                      data = dice, p0 = p0_load)
round(rbind(ME = coef(fit_me), CE = coef(fit_ce)), 4)
#>        s1     s2     s3     s4     s5     s6
#> ME 0.0544 0.0788 0.1142 0.1654 0.2398 0.3475
#> CE 0.0689 0.0630 0.1153 0.1582 0.2413 0.3533
sum(1:6 * coef(fit_ce))                      # still satisfies the mean moment
#> [1] 4.5

4 The information matrix and standard errors

The gradient of the dual (2.6) returns the moment condition (1.1), and its Hessian is the Fisher information matrix for \(\boldsymbol\lambda\),

\[\begin{equation} \big[\mathbf{I}(\boldsymbol\lambda)\big]_{st} = \frac{\partial^2 M}{\partial \lambda_s\, \partial \lambda_t} = \sum_k x_{sk} x_{tk}\, \hat p_k - \Big(\sum_k x_{sk}\hat p_k\Big)\Big(\sum_k x_{tk}\hat p_k\Big) = \mathrm{Cov}_{\hat p}(\mathbf{x}_s, \mathbf{x}_t), \tag{4.1} \end{equation}\]

i.e. the covariance of the moment functions under the fitted distribution \(\hat{\mathbf p}\) (Golan Eq. 4.7). The diagonal elements are the variances \(\mathrm{Var}_{\hat p}(\mathbf{x}_t)\); the off-diagonal elements (\(s \neq t\)) are the covariances \(\mathrm{Cov}_{\hat p}(\mathbf{x}_s, \mathbf{x}_t)\). As shown in Chapter 17 of Cover & Thomas (2006) and Golan (2008) (p. 59), the same matrix can be recovered from the local curvature of the cross-entropy. The variance–covariance matrix of \(\hat{\boldsymbol\lambda}\) is its inverse,

\[\begin{equation} \mathrm{Var}(\hat{\boldsymbol\lambda}) = \mathbf{I}^{-1}(\boldsymbol\lambda), \tag{4.2} \end{equation}\]

and a delta-method step propagates it to \(\hat{\mathbf p}\). These are curvature / identification quantities, not sampling standard errors: a pure inverse problem is deterministic, so they describe how sharply the moments pin down \(\boldsymbol\lambda\), not sampling variability. For sampling inference, use the stochastic-moment sibling inverse_noise().

4.1 Curvature-based standard errors from inverse_ce()

To obtain a full \(\mathbf{I}(\boldsymbol\lambda)\) we use two moments: the mean and the second moment of the faces. We generate consistent moments from a reference distribution (the recovered \(\hat{\mathbf p}\) will still be the maximum-entropy distribution matching them, not that reference).

faces  <- 1:6
p_ref  <- c(.10, .12, .15, .18, .20, .25)
Xm     <- rbind(faces, faces^2)                 # 2 moments x 6 states
ym     <- as.numeric(Xm %*% p_ref)              # feasible (mean, 2nd moment)
d2 <- data.frame(y  = ym,
                 s1 = Xm[, 1], s2 = Xm[, 2], s3 = Xm[, 3],
                 s4 = Xm[, 4], s5 = Xm[, 5], s6 = Xm[, 6])
fit2 <- inverse_ce(y ~ s1 + s2 + s3 + s4 + s5 + s6 - 1, data = d2)
summary(fit2)
#> 
#> Call:
#> inverse_ce(formula = y ~ s1 + s2 + s3 + s4 + s5 + s6 - 1, data = d2)
#> 
#> Moment-fitting residuals (y - Xp):
#> [1]  -3.888e-07  -3.081e-06
#> 
#> Coefficients (estimated probabilities p; delta-method SE from I^-1):
#>     Estimate  Std. Error
#> s1    0.1002      0.2583
#> s2    0.1218      0.1482
#> s3    0.1470      0.1684
#> s4    0.1759      0.1985
#> s5    0.2089      0.1581
#> s6    0.2462      0.3898
#> 
#> Lagrange multipliers (lambda; info-matrix SE, Var(lambda) = I^-1):
#>        Estimate  Std. Error  z value
#> faces    0.2068      3.1054   0.0666
#>         -0.0039      0.4124  -0.0094
#> 
#> Normalized entropy S = 0.9749   pseudo-R2 = 0.02515
#> Fano (sec 7.5): modal error pe = 0.7538 >= bound 0.588
#> Convergence code: 0 (converged) 
#> 
#> Note: the SEs above are information-matrix (curvature) quantities,
#> Var(lambda) = I^-1 (Golan 2008 sec 4.2), measuring how well the moments
#> identify lambda -- not frequentist sampling SEs (a pure inverse problem
#> is deterministic). For sampling-based inference use the
#> stochastic-moment inverse_noise().

The stored vcov() is \(\mathbf{I}^{-1}\); se_lambda and se_p are the square roots of its (propagated) diagonals:

vcov(fit2)                       # I^{-1}, the T x T dual covariance
#>           faces           
#> faces  9.643578 -1.2560837
#>       -1.256084  0.1700659
fit2$se_lambda                   # sqrt(diag(vcov)) for lambda
#>     faces           
#> 3.1054111 0.4123905
fit2$se_p                        # delta-method curvature SEs for p_hat
#>        s1        s2        s3        s4        s5        s6 
#> 0.2583316 0.1482227 0.1683946 0.1984702 0.1581484 0.3898053

By (4.1), a moment whose function has a larger variance under \(\hat{\mathbf p}\) contributes more curvature and so receives a smaller se_lambda. When a moment row is constant, collinear, or the system reaches \(T \ge K\), \(\mathbf{I}(\boldsymbol\lambda)\) is singular; inverse_ce() then rank-checks the inverse and returns NA for the affected standard errors rather than a misleadingly finite pseudo-inverse value.

5 Normalized entropy and Fano bounds

The normalized entropy summarises how much uncertainty remains after conditioning on the data. Relative to a uniform reference it is

\[\begin{equation} S(\hat{\mathbf p}) = \frac{-\sum_k \hat p_k \log(\hat p_k)}{\log(K)} \in [0, 1], \tag{5.1} \end{equation}\]

equal to \(1\) under complete ignorance (uniform \(\hat{\mathbf p}\)) and \(0\) under perfect certainty (a point mass). With a genuine prior it generalises to the prior-relative form \(S = H(\hat{\mathbf p}) / H(\mathbf{p}^0)\),

\[\begin{equation} S(\hat{\mathbf p}) = \frac{-\sum_k \hat p_k \log(\hat p_k)} {-\sum_k p_k^0 \log(p_k^0)}, \tag{5.2} \end{equation}\]

which is what inverse_ce() reports as fit$S:

c(ME = fit_me$S, CE = fit_ce$S)
#>        ME        CE 
#> 0.9005568 1.0539417

The recovered \(\hat{\mathbf p}\) is a single distribution over \(K\) states, so Golan’s Fano inequality applies directly: the modal-classification error \(p_e = 1 - \max_k \hat p_k\) is bounded below by \(S_u(\hat{\mathbf p}) - \log(2)/\log(K)\), where \(S_u = H(\hat{\mathbf p})/\log(K)\) is the uniform-reference normalized entropy. fano_bounds() returns these quantities:

fano_bounds(fit2)
#>       p_max        pe        H         S  pe_lower
#> 1 0.2461965 0.7538035 1.746698 0.9748505 0.5879977

6 From exact to noisy moments

Everything above treats the moments in (1.1) as holding exactly. When the moments are instead measured with noise — \(\mathbf{y} = \mathbf{X}\mathbf{p} + \boldsymbol\varepsilon\) — the problem becomes a generalized maximum-entropy / cross-entropy problem, and the multipliers acquire genuine sampling standard errors. That stochastic-moment sibling is inverse_noise(), which shares inverse_ce()’s formula interface and adds a noise support and signal/noise weight. As the noise support shrinks to zero, its estimates converge to the inverse_ce() solution developed here.

References

Cover, T. M., & Thomas, J. (2006). Elements of information theory (2nd Edition). John Wiley & Sons.
Csiszar, I. (1991). Why least squares and maximum entropy? An axiomatic approach to inference for linear inverse problems. The Annals of Statistics, 19(4), 2032–2066.
Golan, A. (2008). Information and entropy econometrics—a review and synthesis. Foundations and Trends in Econometrics, 2(1-2), 1–145.
Jaynes, E. T. (1957a). Information theory and statistical mechanics. Physical Review, 106(4), 620.
Jaynes, E. T. (1957b). Information theory and statistical mechanics. II. Physical Review, 108(2), 171.
Kullback, S. (1974). The information in contingency tables.
Levine, R. D. (1980). An information theoretical approach to inversion problems. Journal of Physics A: Mathematical and General, 13(1), 91–108.
Shannon, C. E. (1948). A mathematical theory of communication. The Bell System Technical Journal, 27(3), 379–423.
Shore, J., & Johnson, R. (1980). Axiomatic derivation of the principle of maximum entropy and the principle of minimum cross-entropy. IEEE Transactions on Information Theory, 26(1), 26–37.
Zellner, A. (1988). Optimal information processing and bayes’s theorem. The American Statistician, 42(4), 278–280.
Zellner, A. (1991). Bayesian methods and entropy in economics and econometrics. In Maximum entropy and bayesian methods: Laramie, wyoming, 1990 (pp. 17–31). Springer.