The Hidden Engine Behind LinearRegression

1. A Small Framework Hiding in One File

This module lives in scikit-learn/linear_model/_base.py, alongside Ridge, coordinate‑descent, and logistic implementations. It looks like “the place where LinearRegression lives,” but it actually defines a tiny framework that many estimators depend on.

scikit-learn/
  linear_model/
    _base.py  <-- this module
    _ridge.py
    _coordinate_descent.py
    _logistic.py

Key relationships in _base.py:

  +------------------+          +------------------------+
  |  LinearModel     |<---------|  LinearRegression      |
  |  (base class)    |          |  (concrete regressor)  |
  +------------------+          +------------------------+
             ^                             ^
             |                             |
   +--------------------+         +---------------------+
   | LinearClassifier   |         | SparseCoefMixin     |
   | Mixin              |         | (optional mixin)    |
   +--------------------+         +---------------------+

  +----------------------+    +------------------------------+
  |  _preprocess_data    |    |  _pre_fit / _check_gram      |
  |  (shared helper)     |    |  (L1/L0 pre-fit helpers)     |
  +----------------------+    +------------------------------+

  +------------------+
  |  make_dataset    |
  |  (SeqDataset)    |
  +------------------+

This file:

• Defines base behavior for linear estimators via LinearModel.
• Adds classification behavior with LinearClassifierMixin.
• Handles sparse coefficient storage with SparseCoefMixin.
• Implements the concrete LinearRegression estimator.
• Offers shared preprocessing and Gram‑matrix helpers reused across linear models.

Conceptually, this is less “one model implementation” and more “an engine room for linear models”. LinearRegression, Lasso, ElasticNet, and friends all plug into the same centering, weighting, and prediction machinery instead of rebuilding it.

2. One Preprocessing Pipeline for All Linear Models

The framework starts with a shared contract for turning user input into a clean optimization problem. The core helper, _preprocess_data, encodes the rules for centering, intercepts, sample weights, and sparse vs. dense handling.

2.1. Centralizing validation and intercept semantics

def _preprocess_data(
    X,
    y,
    *,
    fit_intercept,
    copy=True,
    copy_y=True,
    sample_weight=None,
    check_input=True,
    rescale_with_sw=True,
):
    """Common data preprocessing for fitting linear models."""
    xp, _, device_ = get_namespace_and_device(X, y, sample_weight)
    n_samples, n_features = X.shape
    X_is_sparse = sp.issparse(X)

    if check_input:
        X = check_array(
            X,
            copy=copy,
            accept_sparse=["csr", "csc"],
            dtype=supported_float_dtypes(xp),
        )
        y = check_array(y, dtype=X.dtype, copy=copy_y, ensure_2d=False)

    if fit_intercept:
        if X_is_sparse:
            X_offset, X_var = mean_variance_axis(
                X, axis=0, weights=sample_weight
            )
        else:
            X_offset = _average(X, axis=0, weights=sample_weight, xp=xp)
            X_offset = xp.astype(X_offset, X.dtype, copy=False)
            X -= X_offset

        y_offset = _average(y, axis=0, weights=sample_weight, xp=xp)
        y -= y_offset
    else:
        X_offset = xp.zeros(n_features, dtype=X.dtype, device=device_)
        y_offset = 0.0

    X_scale = xp.ones(n_features, dtype=X.dtype, device=device_)

    if sample_weight is not None and rescale_with_sw:
        X, y, sample_weight_sqrt = _rescale_data(
            X, y, sample_weight, inplace=copy
        )
    else:
        sample_weight_sqrt = None

    return X, y, X_offset, y_offset, X_scale, sample_weight_sqrt

Critical design choices here:

1. Validation is standardized. All callers go through check_array and the array‑API utilities, so types and shapes are consistent before optimization.
2. Intercept handling is explicit. If fit_intercept=True, the function always computes and returns X_offset and y_offset. For dense data it actually centers X and y in place; for sparse data it keeps X uncentered and pushes the centering burden to the solver logic.
3. Sample weights are normalized into the data. Instead of special‑casing weights in each solver, the pipeline optionally rescales samples with sqrt(sample_weight) via _rescale_data, turning weighted least squares into ordinary least squares.

After this step, solvers can assume they operate on a consistent, centered, optionally rescaled dataset with known offsets. All of the “API surface” complexity—intercepts, weights, sparse formats—has been concentrated into one helper instead of leaking into every algorithm.

2.2. Making sample weights disappear

Weighted problems often complicate both centering and loss computation. This module uses a standard reduction: convert the weighted quadratic form into an unweighted one by folding weights into the data.

def _rescale_data(X, y, sample_weight, inplace=False):
    """Rescale data sample-wise by square root of sample_weight."""
    xp, _ = get_namespace(X, y, sample_weight)
    n_samples = X.shape[0]
    sample_weight_sqrt = xp.sqrt(sample_weight)

    if sp.issparse(X) or sp.issparse(y):
        sw_matrix = sparse.dia_matrix(
            (sample_weight_sqrt, 0), shape=(n_samples, n_samples)
        )

    if sp.issparse(X):
        X = safe_sparse_dot(sw_matrix, X)
    else:
        if inplace:
            X *= sample_weight_sqrt[:, None]
        else:
            X = X * sample_weight_sqrt[:, None]

    # y is rescaled similarly (omitted for brevity)
    return X, y, sample_weight_sqrt

Once _rescale_data runs, the rest of the pipeline no longer needs to “know” about sample weights. The weighted loss has been turned into a standard ||y_rescaled - X_rescaled w||² objective. The entire weight behavior is isolated here, which simplifies reasoning and testing.

The underlying idea generalizes: express complex options as early data transformations so that downstream components can remain simple and generic.

3. Strategy Dispatch Inside LinearRegression.fit

With preprocessing in place, LinearRegression just needs to choose how to solve the least‑squares problem. The fit method is essentially a strategy dispatcher: it selects among three solving approaches—non‑negative constraints, sparse least squares, and dense least squares—behind one API.

3.1. The high‑level fit structure

@_fit_context(prefer_skip_nested_validation=True)
def fit(self, X, y, sample_weight=None):
    n_jobs_ = self.n_jobs
    accept_sparse = False if self.positive else ["csr", "csc", "coo"]

    X, y = validate_data(
        self,
        X,
        y,
        accept_sparse=accept_sparse,
        y_numeric=True,
        multi_output=True,
        force_writeable=True,
    )

    has_sw = sample_weight is not None
    if has_sw:
        sample_weight = _check_sample_weight(
            sample_weight, X, dtype=X.dtype, ensure_non_negative=True
        )

    copy_X_in_preprocess_data = self.copy_X and not sp.issparse(X)

    X, y, X_offset, y_offset, _, sample_weight_sqrt = _preprocess_data(
        X,
        y,
        fit_intercept=self.fit_intercept,
        copy=copy_X_in_preprocess_data,
        sample_weight=sample_weight,
    )

    if self.positive:
        ...  # non-negative branch
    elif sp.issparse(X):
        ...  # sparse lsqr branch
    else:
        ...  # dense lstsq branch

    if y.ndim == 1:
        self.coef_ = np.ravel(self.coef_)
    self._set_intercept(X_offset, y_offset)
    return self

The flow is deliberately simple:

1. Standardize inputs with validate_data and _preprocess_data.
2. Use configuration (positive) and structure (sparse vs. dense) to pick a solving strategy.
3. Normalize the learned parameters’ shapes and compute intercept_ via _set_intercept.

The branches themselves show how to plug specialized solvers into this common framework without changing the outward API.

3.2. Non‑negative least squares via nnls

When positive=True, coefficients must satisfy w_i ≥ 0. Instead of teaching the dense or sparse solvers about constraints, the code switches to scipy.optimize.nnls and keeps the rest of the pipeline unchanged.

if self.positive:
    if y.ndim < 2:
        self.coef_ = optimize.nnls(X, y)[0]
    else:
        outs = Parallel(n_jobs=n_jobs_)(
            delayed(optimize.nnls)(X, y[:, j])
            for j in range(y.shape[1])
        )
        self.coef_ = np.vstack([out[0] for out in outs])

Each target dimension is solved independently, which parallelizes naturally. Upstream and downstream logic (preprocessing, intercept setting, predict) doesn’t change at all; only the solving strategy swaps out.

3.3. Sparse least squares with a centered LinearOperator

For sparse X, materializing a dense centered matrix would destroy sparsity and memory efficiency. Instead, the code builds a LinearOperator that represents “apply X and subtract the mean contribution” without ever allocating the centered matrix.

elif sp.issparse(X):
    if has_sw:
        def matvec(b):
            return X.dot(b) - sample_weight_sqrt * b.dot(X_offset)

        def rmatvec(b):
            return X.T.dot(b) - X_offset * b.dot(sample_weight_sqrt)
    else:
        def matvec(b):
            return X.dot(b) - b.dot(X_offset)

        def rmatvec(b):
            return X.T.dot(b) - X_offset * b.sum()

    X_centered = sparse.linalg.LinearOperator(
        shape=X.shape, matvec=matvec, rmatvec=rmatvec
    )

    if y.ndim < 2:
        self.coef_ = lsqr(X_centered, y, atol=self.tol, btol=self.tol)[0]
    else:
        outs = Parallel(n_jobs=n_jobs_)(
            delayed(lsqr)(
                X_centered, y[:, j].ravel(), atol=self.tol, btol=self.tol
            )
            for j in range(y.shape[1])
        )
        self.coef_ = np.vstack([out[0] for out in outs])

Two key ideas:

• Preserve sparsity. The design avoids ever building a dense centered representation of X.
• Keep semantics aligned. The solver still operates on “centered” data in the conceptual sense, matching the dense branch’s assumptions.

3.4. Dense ordinary least squares

For dense data without positivity constraints, LinearRegression is a thin wrapper around SciPy’s lstsq with a condition‑number cutoff.

else:
    cond = max(X.shape) * np.finfo(X.dtype).eps
    self.coef_, _, self.rank_, self.singular_ = linalg.lstsq(X, y, cond=cond)
    self.coef_ = self.coef_.T

After any branch, the method normalizes coef_’s shape and calls _set_intercept(X_offset, y_offset), which applies the same intercept logic that all linear estimators share.

4. Mixins as Capability Sockets

The design becomes even more reusable when you look beyond LinearRegression. The module defines mixins that act as “capability sockets”: drop them into a class, implement a few attributes in fit, and you get a lot of behavior for free.

4.1. LinearModel: shared prediction and intercept math

LinearModel is an abstract base class that centralizes linear prediction and intercept computation.

class LinearModel(BaseEstimator, metaclass=ABCMeta):
    """Base class for Linear Models"""

    @abstractmethod
    def fit(self, X, y):
        """Fit model."""

    def _decision_function(self, X):
        check_is_fitted(self)
        X = validate_data(
            self, X, accept_sparse=["csr", "csc", "coo"], reset=False
        )
        coef_ = self.coef_
        if coef_.ndim == 1:
            return X @ coef_ + self.intercept_
        else:
            return X @ coef_.T + self.intercept_

    def predict(self, X):
        return self._decision_function(X)

    def _set_intercept(self, X_offset, y_offset, X_scale=None):
        xp, _ = get_namespace(X_offset, y_offset, X_scale)

        if self.fit_intercept:
            self.coef_ = xp.astype(self.coef_, X_offset.dtype, copy=False)
            if X_scale is not None:
                self.coef_ = xp.divide(self.coef_, X_scale)

            if self.coef_.ndim == 1:
                self.intercept_ = y_offset - X_offset @ self.coef_
            else:
                self.intercept_ = y_offset - X_offset @ self.coef_.T
        else:
            self.intercept_ = 0.0

Any subclass that, in its fit, sets coef_, intercept_ (implicitly via _set_intercept), and fit_intercept automatically gets a correct predict and consistent intercept computation. This removes a common source of subtle bugs and keeps intercept logic in one place.

4.2. LinearClassifierMixin: from scores to classes and probabilities

Classification adds interpretation on top of linear scores. LinearClassifierMixin packages that behavior: given coef_, intercept_, and classes_, it knows how to compute decision scores, labels, and logistic probabilities.

class LinearClassifierMixin(ClassifierMixin):
    """Mixin for linear classifiers."""

    def decision_function(self, X):
        check_is_fitted(self)
        xp, _ = get_namespace(X)

        X = validate_data(self, X, accept_sparse="csr", reset=False)
        coef_T = self.coef_.T if self.coef_.ndim == 2 else self.coef_
        scores = safe_sparse_dot(X, coef_T, dense_output=True) + self.intercept_
        return (
            xp.reshape(scores, (-1,))
            if (scores.ndim > 1 and scores.shape[1] == 1)
            else scores
        )

    def predict(self, X):
        xp, _ = get_namespace(X)
        scores = self.decision_function(X)
        if len(scores.shape) == 1:
            indices = xp.astype(scores > 0, indexing_dtype(xp))
        else:
            indices = xp.argmax(scores, axis=1)
        return xp.take(self.classes_, indices, axis=0)

    def _predict_proba_lr(self, X):
        prob = self.decision_function(X)
        expit(prob, out=prob)
        if prob.ndim == 1:
            return np.vstack([1 - prob, prob]).T
        else:
            prob /= prob.sum(axis=1).reshape((prob.shape[0], -1))
            return prob

Any linear classifier that mixes this in, learns coef_, intercept_, and classes_ in fit, and delegates prediction to the mixin instantly gains consistent behavior, including sparse support and probability normalization.

4.3. SparseCoefMixin: memory‑aware coefficients

L1‑regularized models often learn mostly‑zero coefficient vectors. SparseCoefMixin gives estimators an opt‑in way to store those coefficients sparsely and to switch representations when needed.

class SparseCoefMixin:
    """Mixin for converting coef_ to and from CSR format."""

    def densify(self):
        msg = "Estimator, %(name)s, must be fitted before densifying."
        check_is_fitted(self, msg=msg)
        if sp.issparse(self.coef_):
            self.coef_ = self.coef_.toarray()
        return self

    def sparsify(self):
        msg = "Estimator, %(name)s, must be fitted before sparsifying."
        check_is_fitted(self, msg=msg)
        self.coef_ = sp.csr_matrix(self.coef_)
        return self

The unifying pattern is simple but powerful: design small, orthogonal mixins that each provide one focused capability (“can predict linearly”, “can classify”, “can sparsify coefficients”) and compose them as needed. This avoids deep inheritance hierarchies while maximizing reuse.

5. Design Patterns to Steal

Stepping back, this module offers a compact playbook for designing reusable, maintainable model code. The same patterns apply far beyond linear regression.

5.1. Normalize options into data early

• Centralize cross‑cutting concerns like fit_intercept, centering, and sample_weight in a helper such as _preprocess_data.
• Have solvers work on a standard form of the problem (centered, rescaled) instead of raw user input.

This keeps solver implementations smaller and easier to test; the tricky semantics live in a single, well‑documented place.

5.2. Use strategy dispatch in fit

• In your public fit (or equivalent), treat configuration and input structure as a strategy selector: “dense vs sparse”, “constrained vs unconstrained”, “precomputed Gram vs raw features”.
• Delegate each strategy to a focused helper or external library, then normalize outputs back into a shared contract (coef_, intercept_).

This pattern scales much better than embedding all logic in a single monolithic method filled with intertwined conditionals.

5.3. Design mixins as reusable capability sockets

• Extract repeating behaviors—prediction formulas, class decoding, sparse conversions—into mixins like LinearModel, LinearClassifierMixin, and SparseCoefMixin.
• Make each mixin small, with clear expectations about which attributes a class must define to use it.

This approach lets new estimators snap into a shared ecosystem of behavior by implementing a minimal fit and a few attributes, rather than re‑implementing boilerplate around them.

5.4. Guard against subtle misuse of advanced paths

The Gram‑matrix helpers in this module illustrate defensive design around expensive or fragile operations:

• _check_precomputed_gram_matrix validates user‑supplied Gram matrices by recomputing a single entry and comparing, catching silent inconsistencies.
• _pre_fit reconciles centering options with Gram matrices and can recompute safely when they conflict.

Whenever you expose advanced, performance‑oriented options (like precomputed structures), add cheap consistency checks to avoid numerically wrong results that are hard to trace.

5.5. Keep the API small, the internals structured, and performance observable

From the outside, the contract is tiny:

• fit(X, y, sample_weight=None),
• well‑defined coef_ and intercept_,
• standard predict / decision_function semantics.

Internally, the module orchestrates a shared preprocessing pipeline, multiple solving strategies, sparse and dense paths, intercept logic, and mixin‑based capabilities. That combination—a small public surface over a well‑structured internal framework—is exactly what you want when building reusable libraries.

For real systems, you also need to see how this engine behaves under load. Simple metrics like overall fit time, prediction throughput, Gram‑related memory usage in Gram‑based models, and iteration counts for iterative solvers (such as lsqr) turn performance questions into concrete numbers you can act on.

The core lesson from _base.py is this: treat your model implementations as instances of a framework you’re designing. Build a shared preprocessing pipeline, express variations as strategy choices, and package common behavior into mixins. Once you do, adding a new estimator stops being an exercise in copy‑paste and becomes a matter of wiring the right solver into a well‑designed engine.