[np-ml] Linear Regression

算法描述

A weighted linear least-squares regression model.

Notes
In weighted linear least-squares regression, a real-valued target vector, \(\mathbf{y}\), is modeled as a linear combination of covariates \(\mathbf{X}\), and model coefficients, \(\beta\):

\[y_i = \beta^{\top} \mathbf{x}_{i} + \epsilon_i \]

where \(\epsilon_i\sim\mathcal{N}(0, \sigma^2_i )\) is the error term associated with \(i\)-th example, and \(\sigma^{2}_i\) is the variance of the corresponding example.

Under the model, the maximum-likelihood estimate for the regression coefficients, \(\beta\), is:

\[\hat{\beta} = \Sigma^{-1} \mathbf{X}^{\top} \mathbf{Wy} \]

where \(\Sigma^{-1} = (\mathbf{X}^{\top} \mathbf{WX})^{-1}\) and \(\mathbf{W}\) is a diagonal matrix of weights, with each entry inversely proportional to the variance of the corresponding measurement.
When \(\mathbf{W}\) is the identity matrix the examples are weighted equally and the model reduces to standard linear least squares.

算法实现

直接按照解析解进行更新：

import numpy as np


def fit(X, y, fit_intercept=True, weights=None):
    N = X.shape[0]

    weights = np.ones(N) if weights is None else np.atleast_1d(weights)
    weights = np.squeeze(weights) if weights.size > 1 else weights
    err_str = f"`weights` must have shape ({N},), but got {weights.shape}"
    assert weights.shape == (N,), err_str

    # scale X and y by the weight associated with each examples
    W = np.diag(np.sqrt(weights))
    X, y = W @ X, W @ y

    # convert X to a design matrix if we are fitting an intercept
    if fit_intercept:
        X = np.c_[np.sqrt(weights), X]

    sigma_inv = np.linalg.pinv(X.T @ X)
    beta = np.atleast_2d(sigma_inv @ X.T @ y)

    return sigma_inv, beta

---

SGD方式进行模型参数更新：

import numpy as np
from numpy import ndarray


def main(X: ndarray, y: ndarray, num_iteration=1000, learning_rate=1e-2, weights=None):
    _, input_dim = X.shape

    # initialize weights
    weights = weights or np.random.randn(input_dim, 1)

    for iteration in range(num_iteration):
        # (batch_size, 1)
        y_pred = X @ weights
        # (batch_size, 1)
        diff = y_pred - y
        # (input_dim, 1)
        grad = X.T @ diff

        weights -= learning_rate * grad

    return weights


def test_main():
    A = np.array([[1, 1, 0]]).T
    X = np.random.rand(5, 3)
    y = X @ A
    weights = main(X, y, num_iteration=10000, learning_rate=1e-2)
    assert np.allclose(A, weights, rtol=1e-4), f"{weights=}\n but expect {A=}"


if __name__ == "__main__":
    test_main()

参考资料

code
doc