Binding Internals: OLS

From NumPy inputs to a native least-squares fit

OLS illustrates the package’s basic binding pattern: PyO3 exposes a Rust class, NumPy arrays cross the boundary, Rust performs the numerical work, and results return as arrays in a Python dictionary. See the class reference for covariance options and signatures.

Ownership and fitting

src/lib.rs registers OLS; its implementation is in src/estimators/linear.rs. Shared linear algebra and array conversion live in src/utils.rs.

  1. fit(x, y) clears the previous fitted state and validates finite arrays, compatible row counts, and a usable design.
  2. NumPy inputs are copied into owned Rust arrays. The public OLS() constructor includes an intercept; it has no fit_intercept argument.
  3. A QR least-squares solve obtains the parameters. The wrapper separates the intercept and slopes and retains data needed for inference and resampling.
  4. predict(x_new) checks the fitted feature count and finite values before forming the linear prediction.

Do not interpret a failed refit as leaving an earlier fit available. Python argument-conversion errors occur before Rust entry and are a separate boundary.

Covariance and weights

The point estimate minimizes squared residuals. summary() selects classical, HC0–HC3, Newey-West, or cluster covariance without refitting coefficients. The inverse Gram factor used in covariance is obtained from QR rather than explicitly forming and inverting the normal equations.

fit_weighted() uses square-root transformed rows and analytic weights. Zero-weight rows are excluded from inference counts; zero-only clusters do not increase the cluster count. HAC lags refer to the active-row sequence. These weights are not frequency counts or a survey-design specification.

Code 
import numpy as np
import crabbymetrics as cm

rng = np.random.default_rng(17)
x = rng.normal(size=(400, 3))
y = 0.5 + x @ np.array([1.0, -0.7, 0.25]) + rng.normal(size=400)
model = cm.OLS()
model.fit(x, y)
result = model.summary(vcov="hc1")
print("Coefficients:", result["coef"])
print("HC1 SE:", result["coef_se"])
print("Predictions:", model.predict(x[:3]))
print("Bootstrap shape:", model.bootstrap(20, seed=17).shape)

weights = np.ones(len(y))
weights[:20] = 0
weighted = cm.OLS()
weighted.fit_weighted(x, y, weights)
active = cm.OLS()
active.fit(x[20:], y[20:])
np.testing.assert_allclose(
    weighted.summary(vcov="hc1")["coef_se"],
    active.summary(vcov="hc1")["coef_se"],
)
Coefficients: [ 1.03016958 -0.83282979  0.33450749]
HC1 SE: [0.05136313 0.05123947 0.04863709]
Predictions: [1.20890514 0.82276708 0.35398366]
Bootstrap shape: (20, 4)

Resampling and approximation

Bootstrap draws resample stored rows and solve the same least-squares problem. Indices are generated one draw at a time, so they do not require a full n_bootstrap by n index matrix.

fit_sketch() is a separate approximate fitting path. Its summary reports the fit method, sketch size, and seed. A reported covariance is not a proof that sketching preserves interval coverage; see the sketching ablation.

For the next level of numerical complexity, read the native Poisson solver. For custom Python objectives, read the MEstimator callback bridge.