Logit

Binary logistic regression

1 Where it fits

Group: Regression

Logit models a binary outcome through

\[ \Pr(Y_i=1\mid X_i=x_i)=\Lambda(\alpha+x_i'\beta), \]

with optional L2 regularization controlled by alpha. predict(x) returns probabilities for label 1; predict_label(x) returns class labels.

2 Likelihood and solver

With \(\eta_i=\alpha_0+x_i'\beta\) and \(p_i=\{1+\exp(-\eta_i)\}^{-1}\), the native fit minimizes the unnormalized penalized negative log-likelihood

\[ Q(\alpha_0,\beta) = \sum_{i=1}^n \left[ \log\{1+\exp(\eta_i)\}-y_i\eta_i \right] +\frac{\lambda}{2}\|\beta\|_2^2, \]

where the constructor argument named alpha is \(\lambda\) and the intercept is unpenalized. The implementation evaluates \(\log(1+\exp\eta)\) with a branch-stable softplus and evaluates the logistic function without overflowing for large negative indices. It supplies the analytic gradient to ten-vector-memory L-BFGS with a More-Thuente line search, starts all parameters at zero, and stops only when Argmin reports solver convergence or a target cost.

Reaching max_iterations is not convergence. fit() raises ValueError and clears any previous fitted state if the budget is exhausted or optimization otherwise fails. A successful summary() records converged, iterations, termination_reason, and the final penalized objective. The linear index and probability always refer to label 1.

3 Implementation walkthrough

This is a direct likelihood implementation, not a call to Linfa and not iteratively reweighted least squares (IRLS). The fit proceeds as follows.

  1. fit() first clears the stored model and training arrays, converts the NumPy inputs into owned Rust ndarray arrays, and checks finiteness, dimensions, the penalty, the iteration controls, and that both integer labels 0 and 1 occur. A failed refit therefore cannot leave a stale fitted model behind.
  2. The optimizer vector stores the \(p\) slopes first and the intercept last. The objective callback walks over observations, computes \(\eta_i=x_i'\beta+\alpha_0\), and accumulates softplus(eta) - y * eta. The softplus is evaluated as \(\eta+\log(1+e^{-\eta})\) for positive \(\eta\) and \(\log(1+e^\eta)\) otherwise. The sigmoid likewise branches on the sign of \(\eta\). These branches are what keep the likelihood and gradient finite when the linear index is large in magnitude.
  3. The analytic score callback accumulates \(x_i(p_i-y_i)\) for the slopes and \(p_i-y_i\) for the intercept, then adds \(\lambda\beta\) only to the slopes. It does not construct an \(n\times p\) working design or solve a weighted least-squares problem.
  4. Argmin’s L-BFGS stores ten correction pairs and obtains its step length from a More-Thuente line search. The initial vector is exactly zero, so the initial fitted probability is \(1/2\) for every row. gradient_tolerance is passed directly to the L-BFGS gradient stopping test.
  5. The wrapper accepts Argmin’s best parameter only when its termination status is a genuine solver convergence or target-cost condition. Hitting the iteration budget is converted into an exception. Only then are the coefficients and original training arrays installed on the Python object.
  6. predict_lin() performs one dense matrix-vector product and adds the fitted intercept. predict() applies the same stable sigmoid elementwise. predict_label() compares those probabilities with the requested cutoff, including equality on the label-1 side.

The deliberate tradeoff is simplicity and a low-memory first-order iteration. IRLS would exploit the exact Hessian through a weighted least-squares solve at every iteration; this implementation instead spends \(O(np)\) per objective/gradient call and lets L-BFGS approximate curvature. That avoids repeatedly forming and factorizing a \(p\times p\) matrix, but it can take more iterations and has no special separation detector.

4 Inference

Analytic inference is available only when \(\lambda=0\). Let \(\tilde X=[\mathbf1,X]\) and \(W=\operatorname{diag}\{\hat p_i(1-\hat p_i)\}\). The returned model-based covariance is

\[ \widehat V = (\tilde X'W\tilde X)^{-1}. \]

The inverse is computed through QR of \(W^{1/2}\tilde X\), without an inference ridge. Rank-deficient information leaves the point estimate available but returns inference_available=False, null covariance/SEs, and an inference_reason; Wald inference raises. There is no robust logit sandwich. The pairs bootstrap refits the same penalty and aborts on an unestimable resample.

5 Performance and numerical behavior

Each L-BFGS objective and gradient evaluation is \(O(np)\), with solver state linear in \(p\). Fisher covariance adds \(O(np^2+p^3)\) QR/triangular-solve work. Complete or near separation can still make inference unreliable even if numerical rank checks pass. Bootstrap multiplies the optimization cost by the draw count; indices are streamed one replicate at a time.

6 Python API

Constructor: cm.Logit

Fit with integer labels using fit(x, y_int32). Labels must be exactly 0 and 1, with both represented. summary() returns intercept and coefficient estimates plus fit diagnostics. Fisher-information standard errors are available only for the unpenalized fit (alpha=0); penalized fits are explicitly prediction-only. bootstrap() resamples observations and refits the classifier.

Code 
print(inspect.signature(cm.Logit))
(alpha=0.0, max_iterations=100, gradient_tolerance=0.0001)
Code 
cls = cm.Logit
display(HTML(html_table(["Public method"], public_methods(cls))))
Public method
Logit(alpha=0.0, max_iterations=100, gradient_tolerance=0.0001)
bootstrap(self, /, n_bootstrap, seed=None)
fit(self, /, x, y)
predict(self, /, x)
predict_label(self, /, x, cutoff=0.5)
predict_lin(self, /, x)
summary(self, /)
wald_test(self, /, r, q=None)

7 Minimal example

rng = np.random.default_rng(5)
x = rng.normal(size=(220, 3))
eta = -0.2 + x @ np.array([0.7, -0.4, 0.9])
p = 1 / (1 + np.exp(-eta))
y = rng.binomial(1, p, size=220).astype(np.int32)
model = cm.Logit(max_iterations=200)
model.fit(x, y)
fit = model.summary()
print({key: fit[key] for key in ['converged', 'iterations', 'termination_reason', 'objective']})
print(fit['coef'])
print(model.predict(x[:5]))
{'converged': True, 'iterations': 7, 'termination_reason': 'Solver converged', 'objective': 129.69250448711378}
[ 0.56568906 -0.58550069  0.85818987]
[0.5239276  0.4143503  0.68494357 0.42404209 0.2112137 ]

8 summary() contract

The table below is generated by fitting the live class in this repository and then inspecting summary(). Shapes are shown because most values are plain NumPy arrays or scalars.

Code 
rng = np.random.default_rng(105)
x = rng.normal(size=(90, 3))
p = 1 / (1 + np.exp(-(x @ np.array([0.7, -0.4, 0.9]))))
y = rng.binomial(1, p, size=90).astype(np.int32)
model = cm.Logit(max_iterations=200)
model.fit(x, y)
summary = model.summary()
display(HTML(html_table(["summary() key", "shape"], summary_shape_rows(summary))))
summary() key shape
intercept ()
coef (3,)
penalty ()
inference_available ()
converged ()
iterations ()
termination_reason ()
objective ()
intercept_se ()
coef_se (3,)
vcov (4, 4)