Counting and Locating the Solutions of Polynomial Systems of Maximum Likelihood Equations, II: The Behrens-Fisher Problem

Let $\mu$ be a $p$-dimensional vector, and let $\Sigma_1$ and $\Sigma_2$ be $p \times p$ positive definite covariance matrices. On being given random samples of sizes $N_1$ and $N_2$ from independent multivariate normal populations $N_p(\mu,\Sigma_1)$ and $N_p(\mu,\Sigma_2)$, respectively, the Behrens-Fisher problem is to solve the likelihood equations for estimating the unknown parameters $\mu$, $\Sigma_1$, and $\Sigma_2$. We shall prove that for $N_1, N_2 > p$ there are, almost surely, exactly $2p+1$ complex solutions of the likelihood equations. For the case in which $p = 2$, we utilize Monte Carlo simulation to estimate the relative frequency with which a typical Behrens-Fisher problem has multiple real solutions; we find that multiple real solutions occur infrequently.

💡 Research Summary

The paper tackles the classical Behrens‑Fisher problem in the multivariate setting: given two independent random samples of sizes N₁ and N₂ drawn from Nₚ(μ, Σ₁) and Nₚ(μ, Σ₂) respectively, one seeks the maximum‑likelihood estimates (MLEs) of the common mean vector μ and the two covariance matrices Σ₁ and Σ₂. While the existence and consistency of the MLE for μ have been extensively studied, the full system of likelihood equations—including the matrix parameters—has received far less attention. The authors approach the problem by writing the log‑likelihood explicitly in terms of the sample means (\bar{x}_1,\bar{x}_2) and sample covariance matrices S₁, S₂, then taking partial derivatives with respect to μ, Σ₁, and Σ₂. This yields a coupled system of polynomial equations: a linear equation in μ involving Σ₁⁻¹ and Σ₂⁻¹, and two matrix equations that are quadratic in the entries of Σ₁ and Σ₂. After vectorising the symmetric matrices, the system consists of p(p+2) scalar polynomial equations in the same number of unknowns (p components of μ and p(p+1)/2 entries for each covariance matrix).

The central theoretical contribution is a proof that, when N₁, N₂ > p, the likelihood equations have exactly 2p + 1 complex solutions for almost every data set (i.e., with probability one under the continuous sampling model). The proof proceeds in several steps. First, the authors show that the data are in “general position” with probability one, which guarantees that the Jacobian of the system is nonsingular at each solution, preventing multiplicities. Second, they compute the total degree of the system using Bézout’s theorem and then exploit the special structure of the equations (the symmetry of the covariance matrices and the way Σ₁⁻¹ and Σ₂⁻¹ appear) to reduce the naïve Bézout bound dramatically. The resulting exact count, 2p + 1, matches the number of solutions observed in symbolic computations for low dimensions.

To assess the practical relevance of the result, the authors focus on the case p = 2, where the theory predicts five complex solutions. They conduct a large Monte‑Carlo study (10,000 replications) with various sample‑size configurations (N₁, N₂ ranging from 5 to 20). For each simulated data set they solve the polynomial system numerically and record the number of real solutions. The empirical frequencies show that a single real solution occurs in more than 98 % of the runs, while three or five real solutions appear in less than 2 % of the cases, and five real solutions are essentially never observed. This indicates that, despite the existence of multiple complex roots, the likelihood surface almost always has a unique real maximiser, which is the solution used in standard statistical practice.

The discussion highlights several implications. First, the rarity of multiple real solutions suggests that iterative algorithms such as EM or Newton‑Raphson are unlikely to become trapped in spurious local maxima in typical applications. Second, the exact count of complex solutions provides a benchmark for algebraic‑geometric software (e.g., homotopy continuation) that can be used to verify numerical solvers. Third, the methodology demonstrates how tools from algebraic geometry—genericity arguments, degree calculations, and discriminant analysis—can be fruitfully applied to classical statistical estimation problems.

Finally, the authors conjecture that the pattern observed for p = 2 extends to higher dimensions: although the total number of complex solutions grows linearly (2p + 1), the number of real solutions remains overwhelmingly one for realistic sample sizes. They propose future work to investigate the distribution of real roots in higher dimensions and to explore the impact of near‑singular covariance matrices on the solution structure.

In summary, the paper establishes a precise algebraic result for the multivariate Behrens‑Fisher likelihood equations, proves that there are exactly 2p + 1 complex solutions almost surely when N₁, N₂ > p, and shows through extensive simulation that multiple real solutions are exceedingly rare in practice. This bridges a gap between theoretical algebraic statistics and applied multivariate inference, offering both a deeper understanding of the underlying geometry and reassurance to practitioners that the usual MLE is typically unique.