BUG: "lbfgs" optimizer in MixedLM: fixed effects silently returned as `0.0` when `Var(RE)` falls below `tol` in `get_fe_params` but returing convergence "true" when using l

## Code Sample 1 — minimal, no optimizer involved

import numpy as np
import pandas as pd
import statsmodels
import statsmodels.formula.api as smf

print("statsmodels", statsmodels.__version__, "| numpy", np.__version__)

# 30 subjects x 3 conditions. `grp` is constant within a subject
# (between-group), `cond` varies within a subject.
rng  = np.random.default_rng(0)
subj = np.repeat(np.arange(30), 3)
cond = np.tile(["c1", "c2", "c3"], 30)
grp  = np.where(subj % 2 == 0, "A", "B")
y    = 10.0 + 5.0 * (grp == "B") + 2.0 * (cond == "c2") + rng.standard_normal(90)
df   = pd.DataFrame({"subj": subj, "grp": grp, "cond": cond, "y": y})

ols   = smf.ols("y ~ grp + cond", df).fit()
model = smf.mixedlm("y ~ grp + cond", df, groups=df["subj"])

print("terms:", model.exog_names)
print("OLS                        ", np.round(ols.params.values, 6))
for v in [1e-9, 1e-10, 9.9e-11, 1e-12, 0.0]:
    fe, singular = model.get_fe_params(np.array([[v]]), np.array([]))
    print(f"MixedLM Var(RE) = {v:8.2e} ", np.round(fe, 6), " singular =", singular)

### sample 2

## Code Sample 2 — the same failure through the public `fit()` API

#Here the residuals carry a small *negative* within-subject correlation, so the ML
#estimate of the random-effect variance lies on the boundary at zero and `lbfgs` descends
#past `tol` on its own. (`standard_normal` plus an explicit linear transform is used
#instead of `multivariate_normal`, because the latter goes through an SVD and is not
#bit-reproducible across BLAS builds; the printed checksum was identical on every machine
#I tried.)

import numpy as np
import pandas as pd
import statsmodels.formula.api as smf

# 74 subjects x 3 repeated measures, one between-subject factor `grp`
# and one within-subject factor `cond`.
C, SD, N = 0.25, 28.0, 74
rng = np.random.default_rng(12345)
Z = rng.standard_normal((N, 3))
E = SD * (Z - C * Z.mean(axis=1, keepdims=True))   # small negative within-subject rho

rows = []
for s in range(N):
    g = s % 2                                      # between-subject group
    for k, cond in enumerate(["c1", "c2", "c3"]):
        mu = 10.0 + 5.0 * g + (20.0 + 25.0 * g) * (k == 1) + (0.0 + 8.0 * g) * (k == 2)
        rows.append({"subj": s, "grp": "B" if g else "A", "cond": cond, "y": mu + E[s, k]})
df = pd.DataFrame(rows)
print("checksum sum(y) =", round(df["y"].sum(), 6), " (expected 5423.484816)")

terms = ["Intercept", "grp[T.B]", "cond[T.c2]", "grp[T.B]:cond[T.c2]"]
ols = smf.ols("y ~ grp * cond", df).fit()
out = {"OLS (reference)": ols.params[terms].to_dict() | {"Var(RE)": np.nan, "converged": True}}
for method in ["lbfgs", "bfgs", "powell", "cg", "nm"]:
    f = smf.mixedlm("y ~ grp * cond", df, groups=df["subj"]).fit(reml=False, method=method)
    out[method] = f.params[terms].to_dict() | {
        "Var(RE)": float(f.cov_re.iloc[0, 0] / f.scale), "converged": f.converged}
res = pd.DataFrame(out).T
print(res.to_string(formatters={"Var(RE)": "{:.3e}".format}))

Intercept  10.264184
grp[T.B]   -2.064992
cond[T.c2] 27.193142
grp[T.B]:cond[T.c2] 27.249206

## Possible cause

`statsmodels/regression/mixed_linear_model.py`, `get_fe_params` (line 1269 in 0.14.4):


w, v = np.linalg.eigh(cov_re)
if w.min() < tol:                              # tol = 1e-10
    sing = True
    ii = np.flatnonzero(w >= tol)
    if len(ii) == 0:
        cov_re_inv = np.zeros_like(cov_re)     # <-- here
    ...
solver = _smw_solver(1., ex_r, ex2_r, cov_re_inv, vc_vari)
...
if sing:
    fe_params = np.dot(np.linalg.pinv(xtxy[:, 0:-1]), xtxy[:, -1])


Writing `D` for the random-effect covariance, the Woodbury identity used by `_smw_solver`
is


(I + Z D Z')^-1  =  I - Z (D^-1 + Z'Z)^-1 Z'


* the limit for `D -> 0` is `D^-1 -> inf`, the bracket vanishes, the solver returns `I`,
  and the GLS reduces to OLS;
* setting `D^-1 := 0` is instead the limit `D -> inf`. The solver then returns
  `I - Z (Z'Z)^-1 Z'`, the **within-group projection**.

Any column constant within a group is annihilated by that projection, `xtxy` loses rank in
exactly those columns, and `pinv` returns zero for them. Consistent with this, the `lbfgs`
coefficients above are numerically identical to a fixed-effects (within) regression on
subject-demeaned data, with the person-constant terms simply absent.

* statsmodels **0.14.4**, numpy 2.2.6, scipy 1.15.2, pandas 2.2.3, Python 3.11.12
* statsmodels **0.14.5**, numpy 1.26.4, Python 3.10.18

macOS 15 (arm64), both environments conda-forge.

INSTALLED VERSIONS
------------------
Python: 3.11.12.final.0
OS: Darwin 25.6.0 Darwin Kernel Version 25.6.0: Fri Jul 31 19:17:26 PDT 2026; root:xnu-12377.161.14~5/RELEASE_ARM64_T6041 arm64
byteorder: little
LC_ALL: None
LANG: C.UTF-8

statsmodels
===========

Installed: 0.14.4 (/Users/<user>/miniforge3/envs/pymc_env/lib/python3.11/site-packages/statsmodels)

Required Dependencies
=====================

cython: Not installed
numpy: 2.2.6 (/Users/<user>/miniforge3/envs/pymc_env/lib/python3.11/site-packages/numpy)
scipy: 1.15.2 (/Users/<user>/miniforge3/envs/pymc_env/lib/python3.11/site-packages/scipy)
pandas: 2.2.3 (/Users/<user>/miniforge3/envs/pymc_env/lib/python3.11/site-packages/pandas)
    dateutil: 2.9.0.post0 (/Users/<user>/miniforge3/envs/pymc_env/lib/python3.11/site-packages/dateutil)
...
pytest: Not installed
virtualenv: Not installed