diff --git a/.github/workflows/pre-commit.yml b/.github/workflows/pre-commit.yml index b13f7c95..4db2ca03 100644 --- a/.github/workflows/pre-commit.yml +++ b/.github/workflows/pre-commit.yml @@ -10,5 +10,5 @@ jobs: - uses: actions/checkout@v2 - uses: actions/setup-python@v2 with: - python-version: '3.9' + python-version: '3.10' - uses: pre-commit/action@v2.0.3 \ No newline at end of file diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 78f361ba..5ff24e13 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -7,7 +7,7 @@ jobs: timeout-minutes: 5 strategy: matrix: - python-version: [3.8, 3.7] + python-version: ['3.10'] steps: - uses: actions/checkout@v2 - name: Setup Python diff --git a/CHANGELOG.md b/CHANGELOG.md index 1580391f..c2774e4a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,4 +1,6 @@ # Changelog +### 2.0.0 - 2023-09-22 +- feat: adding ordinary least squares regression solver. Updating quantile regression to solve dual [#11](https://github.com/washingtonpost/elex-solver/pull/11) ### 1.1.0 - 2023-04-21 - fix: Not regularizing intercept coefficient + better warning handling [#8](https://github.com/washingtonpost/elex-solver/pull/8) diff --git a/README.md b/README.md index a59d0739..fcd98ee0 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,7 @@ # elex-solver This packages includes solvers for: +* Ordinary least squares regression * Quantile regression * Transition matrices @@ -9,8 +10,11 @@ This packages includes solvers for: * We recommend that you set up a virtualenv and activate it (IE ``mkvirtualenv elex-solver`` via http://virtualenvwrapper.readthedocs.io/en/latest/). * Run ``pip install elex-solver`` +## Ordinary least squares +We have our own implementation of ordinary least squares in Python because this let us optimize it towards the bootstrap by storing and re-using the normal equations. This allows for significant speed up. + ## Quantile Regression -Since we did not find any implementations of quantile regression in Python that fit our needs, we decided to write one ourselves. This uses [`cvxpy`](https://www.cvxpy.org/#) and sets up quantile regression as a normal optimization problem. We use quantile regression for our election night model. +Since we did not find any implementations of quantile regression in Python that fit our needs, we decided to write one ourselves. At the moment this uses two libraries, the version that solves the non-regularized problem uses `numpy`and solves the dual based on [this](https://arxiv.org/pdf/2305.12616.pdf) paper. The version that solves the regularized problem uses [`cvxpy`](https://www.cvxpy.org/#) and sets up the problem as a normal optimization problem. Eventually, we are planning on replacing the regularized version with the dual also. ## Transition matrices We also have a solver for transition matrices. While this works arbitrarily, we have used this in the past for our primary election night model. We can still use this to create the sankey diagram coefficients. diff --git a/setup.py b/setup.py index a1e5824e..681c4df9 100644 --- a/setup.py +++ b/setup.py @@ -3,7 +3,7 @@ from setuptools import find_packages, setup -INSTALL_REQUIRES = "cvxpy<=1.2.0" +INSTALL_REQUIRES = ["cvxpy<=1.2.0", "scipy<1.11.0"] THIS_FILE_DIR = os.path.dirname(__file__) @@ -13,7 +13,7 @@ LONG_DESCRIPTION = f.read() # The full version, including alpha/beta/rc tags -RELEASE = "1.1.0" +RELEASE = "2.0.0" # The short X.Y version VERSION = ".".join(RELEASE.split(".")[:2]) @@ -29,7 +29,7 @@ "Intended Audience :: Developers", "License :: OSI Approved :: MIT License", "Programming Language :: Python", - "Programming Language :: Python :: 3.7", + "Programming Language :: Python :: 3.10", ], description="A package for optimization solvers", long_description=LONG_DESCRIPTION, diff --git a/src/elexsolver/LinearSolver.py b/src/elexsolver/LinearSolver.py new file mode 100644 index 00000000..949cb5bd --- /dev/null +++ b/src/elexsolver/LinearSolver.py @@ -0,0 +1,79 @@ +import logging +import warnings +from abc import ABC + +import numpy as np + +from elexsolver.logging import initialize_logging + +initialize_logging() + +LOG = logging.getLogger(__name__) + + +class LinearSolverException(Exception): + pass + + +class IllConditionedMatrixException(LinearSolverException): + pass + + +class LinearSolver(ABC): + """ + An abstract base class for a linear solver + """ + + CONDITION_WARNING_MIN = 50 # arbitrary + CONDITION_ERROR_MIN = 1e8 # based on scipy + + def __init__(self): + self.coefficients = None + + @classmethod + def fit(self, x: np.ndarray, y: np.ndarray, weights: np.ndarray | None = None, lambda_: float = 0.0, *kwargs): + """ + Fits model + """ + raise NotImplementedError + + def predict(self, x: np.ndarray) -> np.ndarray: + """ + Use coefficients to predict + """ + self._check_any_element_nan_or_inf(x) + + return x @ self.coefficients + + def get_coefficients(self) -> np.ndarray: + """ + Returns model coefficients + """ + return self.coefficients + + def _check_matrix_condition(self, x): + """ + Check condition number of the design matrix as a check for multicolinearity. + This is equivalent to the ratio between the largest and the smallest singular value of the design matrix. + """ + condition_number = np.linalg.cond(x) + if condition_number >= self.CONDITION_ERROR_MIN: + raise IllConditionedMatrixException( + f"Ill-conditioned matrix detected. Matrix condition number >= {self.CONDITION_ERROR_MIN}" + ) + elif condition_number >= self.CONDITION_WARNING_MIN: + warnings.warn("Warning: Ill-conditioned matrix detected. result is not guaranteed to be accurate") + + def _check_any_element_nan_or_inf(self, x): + """ + Check whether any element in a matrix or vector is NaN or infinity + """ + if np.any(np.isnan(x)) or np.any(np.isinf(x)): + raise ValueError("Array contains NaN or Infinity") + + def _check_intercept(self, x): + """ + Check whether the first column is all 1s (normal intercept) otherwise raises a warning. + """ + if ~np.all(x[:, 0] == 1): + warnings.warn("Warning: fit_intercept=True and not all elements of the first columns are 1s") diff --git a/src/elexsolver/OLSRegressionSolver.py b/src/elexsolver/OLSRegressionSolver.py new file mode 100644 index 00000000..ca055354 --- /dev/null +++ b/src/elexsolver/OLSRegressionSolver.py @@ -0,0 +1,143 @@ +import logging + +import numpy as np + +from elexsolver.LinearSolver import LinearSolver +from elexsolver.logging import initialize_logging + +initialize_logging() + +LOG = logging.getLogger(__name__) + + +class OLSRegressionSolver(LinearSolver): + """ + A class for Ordinary Least Squares Regression optimized for the bootstrap + """ + + # OLS setup: + # X \beta = y + # since X might not be square, we multiply the above equation on both sides by X^T to generate X^T X, which is guaranteed + # to be square + # X^T X \beta = X^T y + # Since X^T X is square we can invert it + # \beta = (X^T X)^{-1} X^T y + # Since our version of the model bootstraps y, but keeps X constant we can + # pre-compute (X^T X)^{-1} X^T and then re-use it to compute \beta_b for every bootstrap sample + + def __init__(self): + super().__init__() + self.normal_eqs = None + self.hat_vals = None + + def _get_regularizer( + self, lambda_: float, dim: int, fit_intercept: bool, regularize_intercept: bool, n_feat_ignore_reg: int + ) -> np.ndarray: + """ + Returns the regularization matrix + """ + # lambda_I is the matrix for regularization, which need to be the same shape as R and + # have the regularization constant lambda_ along the diagonal + lambda_I = lambda_ * np.eye(dim) + + # we don't want to regularize the coefficient for intercept + # but we also might not want to fit the intercept + # for some number of features + # so set regularization constant to zero for intercept + # and the first n_feat_ignore_reg features + for i in range(fit_intercept + n_feat_ignore_reg): + # if we are fitting an intercept and want to regularize intercept then we don't want + # to set the regularization matrix at lambda_I[0, 0] to zero + if fit_intercept and i == 0 and regularize_intercept: + continue + lambda_I[i, i] = 0 + + return lambda_I + + def _compute_normal_equations( + self, + x: np.ndarray, + L: np.ndarray, + lambda_: float, + fit_intercept: bool, + regularize_intercept: bool, + n_feat_ignore_reg: int, + ) -> np.ndarray: + """ + Computes the normal equations for OLS: (X^T X)^{-1} X^T + """ + # Inverting X^T X directly is computationally expensive and mathematically unstable, so we use QR factorization + # which factors x into the sum of an orthogonal matrix Q and a upper tringular matrix R + # L is a diagonal matrix of weights + Q, R = np.linalg.qr(L @ x) + + # get regularization matrix + lambda_I = self._get_regularizer(lambda_, R.shape[0], fit_intercept, regularize_intercept, n_feat_ignore_reg) + + # substitute X = QR into the normal equations to get + # R^T Q^T Q R \beta = R^T Q^T y + # R^T R \beta = R^T Q^T y + # \beta = (R^T R)^{-1} R^T Q^T y + # since R is upper triangular it is eqsier to invert + # lambda_I is the regularization matrix + return np.linalg.inv(R.T @ R + lambda_I) @ R.T @ Q.T + + def fit( + self, + x: np.ndarray, + y: np.ndarray, + weights: np.ndarray | None = None, + lambda_: float = 0.0, + normal_eqs: np.ndarray | None = None, + fit_intercept: bool = True, + regularize_intercept: bool = False, + n_feat_ignore_reg: int = 0, + ): + self._check_any_element_nan_or_inf(x) + self._check_any_element_nan_or_inf(y) + + if fit_intercept: + self._check_intercept(x) + + # if weights are none, all rows should be weighted equally + if weights is None: + weights = np.ones((y.shape[0],)) + + # normalize weights and turn into diagional matrix + # square root because will be squared when R^T R happens later + L = np.diag(np.sqrt(weights.flatten() / weights.sum())) + + # if normal equations are provided then use those, otherwise compute them + # in the bootstrap setting we can now pass in the normal equations and can + # save time re-computing them + if normal_eqs is None: + self.normal_eqs = self._compute_normal_equations( + x, L, lambda_, fit_intercept, regularize_intercept, n_feat_ignore_reg + ) + else: + self.normal_eqs = normal_eqs + + # compute hat matrix: X (X^T X)^{-1} X^T + self.hat_vals = np.diag(x @ self.normal_eqs @ L) + + # compute coefficients: (X^T X)^{-1} X^T y + self.coefficients = self.normal_eqs @ L @ y + + def residuals(self, y: np.ndarray, y_hat: np.ndarray, loo: bool = True, center: bool = True) -> np.ndarray: + """ + Computes residuals for the model + """ + # compute standard residuals + residuals = y - y_hat + + # if leave one out is True, inflate by (1 - P) + # in OLS setting inflating by (1 - P) is the same as computing the leave one out residuals + # the un-inflated training residuals are too small, since training covariates were observed during fitting + if loo: + residuals /= (1 - self.hat_vals).reshape(-1, 1) + + # centering removes the column mean + if center: + residuals -= np.mean(residuals, axis=0) + + return residuals diff --git a/src/elexsolver/QuantileRegressionSolver.py b/src/elexsolver/QuantileRegressionSolver.py index 3fe12560..cfe4c525 100644 --- a/src/elexsolver/QuantileRegressionSolver.py +++ b/src/elexsolver/QuantileRegressionSolver.py @@ -1,9 +1,10 @@ import logging -import warnings import cvxpy as cp import numpy as np +from scipy.optimize import linprog +from elexsolver.LinearSolver import LinearSolver from elexsolver.logging import initialize_logging initialize_logging() @@ -11,133 +12,115 @@ LOG = logging.getLogger(__name__) -class QuantileRegressionSolverException(Exception): - pass +class QuantileRegressionSolver(LinearSolver): + def __init__(self): + super().__init__() + self.coefficients = [] - -class IllConditionedMatrixException(QuantileRegressionSolverException): - pass - - -class QuantileRegressionSolver: - - VALID_SOLVERS = {"SCS", "ECOS", "MOSEK", "OSQP", "CVXOPT", "GLPK"} - KWARGS = {"ECOS": {"max_iters": 10000}} - - CONDITION_WARNING_MIN = 50 # arbitrary - CONDITION_ERROR_MIN = 1e8 # based on scipy - - def __init__(self, solver="ECOS"): - if solver not in self.VALID_SOLVERS: - raise ValueError(f"solver must be in {self.VALID_SOLVERS}") - self.tau = cp.Parameter() - self.coefficients = None - self.problem = None - self.solver = solver - - def _check_matrix_condition(self, x): - """ - Check condition number of the design matrix as a check for multicolinearity. - This is equivalent to the ratio between the largest and the smallest singular value of the design matrix. - """ - condition_number = np.linalg.cond(x) - if condition_number >= self.CONDITION_ERROR_MIN: - raise IllConditionedMatrixException( - f"Ill-conditioned matrix detected. Matrix condition number >= {self.CONDITION_ERROR_MIN}" - ) - elif condition_number >= self.CONDITION_WARNING_MIN: - warnings.warn("Warning: Ill-conditioned matrix detected. result is not guaranteed to be accurate") - - def _check_any_element_nan_or_inf(self, x): - """ - Check whether any element in a matrix or vector is NaN or infinity - """ - if np.any(np.isnan(x)) or np.any(np.isinf(x)): - raise ValueError("Array contains NaN or Infinity") - - def _check_intercept(self, x): + def _fit(self, x: np.ndarray, y: np.ndarray, weights: np.ndarray, tau: float) -> np.ndarray: """ - Check whether the first column is all 1s (normal intercept) otherwise raises a warning. + Fits the dual problem of a quantile regression, for more information see appendix 6 here: https://arxiv.org/pdf/2305.12616.pdf """ - if ~np.all(x[:, 0] == 1): - warnings.warn("Warning: fit_intercept=True and not all elements of the first columns are 1s") + S = y + Phi = x + zeros = np.zeros((Phi.shape[1],)) + N = y.shape[0] + bounds = weights.reshape(-1, 1) * np.asarray([(tau - 1, tau)] * N) - def get_loss_function(self, x, y, coefficients, weights): - """ - Get the quantile regression loss function - """ - y_hat = x @ coefficients - residual = y - y_hat - return cp.sum(cp.multiply(weights, 0.5 * cp.abs(residual) + (self.tau.value - 0.5) * residual)) + # A_eq are the equality constraint matrix + # b_eq is the equality constraint vector (ie. A_eq @ x = b_eq) + # bounds are the (min, max) possible values of every element of x + res = linprog(-1 * S, A_eq=Phi.T, b_eq=zeros, bounds=bounds, method="highs", options={"presolve": False}) + # marginal are the dual values, since we are solving the dual this is equivalent to the primal + return -1 * res.eqlin.marginals - def get_regularizer(self, coefficients, fit_intercept): + def _get_regularizer( + self, coefficients: cp.expressions.variable.Variable, regularize_intercept: bool, n_feat_ignore_reg: int + ) -> cp.atoms.norm: """ Get regularization component of the loss function. Note that this is L2 (ridge) regularization. """ - # if we are fitting an intercept in the model, then that coefficient should not be regularized. - # NOTE: assumes that if fit_intercept=True, that the intercept is in the first column - coefficients_to_regularize = coefficients - if fit_intercept: - coefficients_to_regularize = coefficients[1:] + # this assumes that if regularize_intercept=True that the intercept is the first column + # also note that even if regularize_intercept is True BUT n_feat_ignore_req > 0 and fit_intercept + # is true that we are NOT regularizing the intercept + coefficients_to_regularize = coefficients[n_feat_ignore_reg:] + if not regularize_intercept: + coefficients_to_regularize = coefficients[1 + n_feat_ignore_reg :] # noqa: E203 return cp.pnorm(coefficients_to_regularize, p=2) ** 2 - def __solve(self, x, y, weights, lambda_, fit_intercept, verbose): + def _fit_with_regularization( + self, + x: np.ndarray, + y: np.ndarray, + weights: np.ndarray, + tau: float, + lambda_: float, + regularize_intercept: bool, + n_feat_ignore_reg: int, + ): """ - Sets up the optimization problem and solves it + Fits quantile regression with regularization + TODO: convert this problem to use the dual like in the non regularization case """ - self._check_matrix_condition(x) + arguments = {"ECOS": {"max_iters": 10000}} coefficients = cp.Variable((x.shape[1],)) - loss_function = self.get_loss_function(x, y, coefficients, weights) - loss_function += lambda_ * self.get_regularizer(coefficients, fit_intercept) + y_hat = x @ coefficients + residual = y - y_hat + loss_function = cp.sum(cp.multiply(weights, 0.5 * cp.abs(residual) + (tau - 0.5) * residual)) + loss_function += lambda_ * self._get_regularizer(coefficients, regularize_intercept, n_feat_ignore_reg) objective = cp.Minimize(loss_function) problem = cp.Problem(objective) - problem.solve(solver=self.solver, verbose=verbose, **self.KWARGS.get(self.solver, {})) - return coefficients, problem + problem.solve(solver="ECOS", **arguments.get("ECOS", {})) + return coefficients.value def fit( self, - x, - y, - tau_value=0.5, - weights=None, - lambda_=0, - fit_intercept=True, - verbose=False, - save_problem=False, - normalize_weights=True, + x: np.ndarray, + y: np.ndarray, + taus: list | float = 0.5, + weights: np.ndarray | None = None, + lambda_: float = 0.0, + fit_intercept: bool = True, + regularize_intercept: bool = False, + n_feat_ignore_reg: int = 0, + normalize_weights: bool = True, ): """ - Fit the (weighted) quantile regression problem. - Weights should not sum to one. - If fit_intercept=True then intercept is assumed to be the first column in `x` + Fits quantile regression """ - self._check_any_element_nan_or_inf(x) self._check_any_element_nan_or_inf(y) if fit_intercept: self._check_intercept(x) - if weights is None: # if weights are none, give unit weights - weights = [1] * x.shape[0] + # if weights are none, all rows should be weighted equally + if weights is None: + weights = np.ones((y.shape[0],)) + + # normalize weights. default to true if normalize_weights: weights_sum = np.sum(weights) if weights_sum == 0: - # This should not happen raise ZeroDivisionError weights = weights / weights_sum - self.tau.value = tau_value - coefficients, problem = self.__solve(x, y, weights, lambda_, fit_intercept, verbose) - self.coefficients = coefficients.value - if save_problem: - self.problem = problem - else: - self.problem = None + # _fit assumes that taus is list, so if we want to do one value of tau then turn into a list + if isinstance(taus, float): + taus = [taus] + + for tau in taus: + if lambda_ > 0: + coefficients = self._fit_with_regularization( + x, y, weights, tau, lambda_, regularize_intercept, n_feat_ignore_reg + ) + else: + coefficients = self._fit(x, y, weights, tau) + self.coefficients.append(coefficients) - def predict(self, x): + def predict(self, x: np.ndarray) -> np.ndarray: """ - Returns predictions + Use coefficients to predict """ self._check_any_element_nan_or_inf(x) diff --git a/tests/test_linear_solver.py b/tests/test_linear_solver.py new file mode 100644 index 00000000..98e1f163 --- /dev/null +++ b/tests/test_linear_solver.py @@ -0,0 +1,10 @@ +import numpy as np +import pytest + +from elexsolver.LinearSolver import LinearSolver + + +def test_fit(): + solver = LinearSolver() + with pytest.raises(NotImplementedError): + solver.fit(np.ndarray((5, 3)), np.ndarray((1, 3))) diff --git a/tests/test_ols.py b/tests/test_ols.py new file mode 100644 index 00000000..935b2e63 --- /dev/null +++ b/tests/test_ols.py @@ -0,0 +1,275 @@ +import numpy as np +import pytest + +from elexsolver.OLSRegressionSolver import OLSRegressionSolver + +# relatively high tolerance, since different implementation. +TOL = 1e-3 + +# the outputs are compared against lm + +############### +# Basic tests # +############### + + +def test_basic1(): + lm = OLSRegressionSolver() + x = np.asarray([[1], [1], [1], [2]]) + y = np.asarray([3, 8, 9, 15]) + lm.fit(x, y, fit_intercept=False) + preds = lm.predict(x) + assert all(np.abs(preds - [7.142857, 7.142857, 7.142857, 14.285714]) <= TOL) + + +def test_basic2(): + lm = OLSRegressionSolver() + x = np.asarray([[1, 1], [1, 1], [1, 1], [1, 2]]) + y = np.asarray([3, 8, 9, 15]) + lm.fit(x, y, fit_intercept=True) + preds = lm.predict(x) + assert all(np.abs(preds - [6.666667, 6.666667, 6.666667, 15]) <= TOL) + + +###################### +# Intermediate tests # +###################### + + +def test_random_no_intercept(random_data_no_weights): + lm = OLSRegressionSolver() + x = random_data_no_weights[["x0", "x1", "x2", "x3", "x4"]].values + y = random_data_no_weights["y"].values.reshape(-1, 1) + lm.fit(x, y, fit_intercept=False) + predictions = lm.predict(x) + # check coefficients + assert all(np.abs(lm.coefficients - [[1.037], [7.022], [4.794], [4.776], [4.266]]) <= TOL) + + # check hat values + assert len(lm.hat_vals) == 100 + assert lm.hat_vals[0] == pytest.approx(0.037102687) + assert lm.hat_vals[-1] == pytest.approx(0.038703403) + + # check predictions + assert predictions[0] == pytest.approx(10.743149) + assert predictions[-1] == pytest.approx(9.878154) + + +def test_random_intercept(random_data_no_weights): + lm = OLSRegressionSolver() + random_data_no_weights["intercept"] = 1 + x = random_data_no_weights[["intercept", "x0", "x1", "x2", "x3", "x4"]].values + y = random_data_no_weights["y"].values.reshape(-1, 1) + lm.fit(x, y, fit_intercept=True) + predictions = lm.predict(x) + # check coefficients + assert all(np.abs(lm.coefficients - [[0.08111], [1.01166], [6.98917], [4.77003], [4.73864], [4.23325]]) <= TOL) + + # check hat values + assert len(lm.hat_vals) == 100 + assert lm.hat_vals[0] == pytest.approx(0.03751295) + assert lm.hat_vals[-1] == pytest.approx(0.03880323) + + # check predictions + assert predictions[0] == pytest.approx(10.739329) + assert predictions[-1] == pytest.approx(9.880039) + + +###################### +# Tests with weights # +###################### + + +def test_random_weights_no_intercept(random_data_weights): + lm = OLSRegressionSolver() + x = random_data_weights[["x0", "x1", "x2", "x3", "x4"]].values + y = random_data_weights["y"].values.reshape(-1, 1) + weights = random_data_weights["weights"].values + lm.fit(x, y, weights=weights, fit_intercept=False) + predictions = lm.predict(x) + + # coefficients + assert all(np.abs(lm.coefficients - [[1.455], [2.018], [4.699], [3.342], [9.669]]) <= TOL) + + # check hat values + assert len(lm.hat_vals) == 100 + assert lm.hat_vals[0] == pytest.approx(0.013961893) + assert lm.hat_vals[-1] == pytest.approx(0.044913885) + + # check predictions + assert predictions[0] == pytest.approx(16.090619) + assert predictions[-1] == pytest.approx(12.538442) + + +def test_random_weights_intercept(random_data_weights): + lm = OLSRegressionSolver() + random_data_weights["intercept"] = 1 + x = random_data_weights[["intercept", "x0", "x1", "x2", "x3", "x4"]].values + y = random_data_weights["y"].values.reshape(-1, 1) + weights = random_data_weights["weights"].values + lm.fit(x, y, weights=weights, fit_intercept=True) + predictions = lm.predict(x) + + # coefficients + assert all(np.abs(lm.coefficients - [[0.1151], [1.4141], [1.9754], [4.6761], [3.2803], [9.6208]]) <= TOL) + + # check hat values + assert len(lm.hat_vals) == 100 + assert lm.hat_vals[0] == pytest.approx(0.014940744) + assert lm.hat_vals[-1] == pytest.approx(0.051229900) + + # check predictions + assert predictions[0] == pytest.approx(16.069887) + assert predictions[-1] == pytest.approx(12.565045) + + +######################## +# Test regularization # +######################## + + +def test_regularizer(): + lm = OLSRegressionSolver() + + lambda_I = lm._get_regularizer(7, 10, fit_intercept=True, regularize_intercept=True, n_feat_ignore_reg=2) + + assert lambda_I.shape == (10, 10) + assert lambda_I[0, 0] == pytest.approx(7) + assert lambda_I[1, 1] == pytest.approx(0) + assert lambda_I[2, 2] == pytest.approx(0) + assert lambda_I[3, 3] == pytest.approx(7) + assert lambda_I[4, 4] == pytest.approx(7) + assert lambda_I[5, 5] == pytest.approx(7) + assert lambda_I[6, 6] == pytest.approx(7) + assert lambda_I[7, 7] == pytest.approx(7) + + lambda_I = lm._get_regularizer(7, 10, fit_intercept=True, regularize_intercept=False, n_feat_ignore_reg=2) + + assert lambda_I.shape == (10, 10) + assert lambda_I[0, 0] == pytest.approx(0) + assert lambda_I[1, 1] == pytest.approx(0) + assert lambda_I[2, 2] == pytest.approx(0) + assert lambda_I[3, 3] == pytest.approx(7) + assert lambda_I[4, 4] == pytest.approx(7) + assert lambda_I[5, 5] == pytest.approx(7) + assert lambda_I[6, 6] == pytest.approx(7) + assert lambda_I[7, 7] == pytest.approx(7) + + +def test_regularization_with_intercept(random_data_no_weights): + x = random_data_no_weights[["x0", "x1", "x2", "x3", "x4"]].values + x[:, 0] = 1 + y = random_data_no_weights["y"].values + + lm = OLSRegressionSolver() + lambda_ = 1e6 + lm.fit(x, y, lambda_=lambda_, fit_intercept=True, regularize_intercept=False) + coefficients_w_reg = lm.coefficients + assert all(np.abs(coefficients_w_reg[1:] - [0, 0, 0, 0]) <= TOL) + assert np.abs(coefficients_w_reg[0]) > TOL + + +def test_regularization_with_intercept_and_unreg_feature(random_data_no_weights): + x = random_data_no_weights[["x0", "x1", "x2", "x3", "x4"]].values + x[:, 0] = 1 + y = random_data_no_weights["y"].values + + lm = OLSRegressionSolver() + lambda_ = 1e6 + lm.fit(x, y, lambda_=lambda_, fit_intercept=True, regularize_intercept=False, n_feat_ignore_reg=2) + coefficients_w_reg = lm.coefficients + assert all(np.abs(coefficients_w_reg[3:] - [0, 0]) <= TOL) + assert np.abs(coefficients_w_reg[0]) > TOL + assert np.abs(coefficients_w_reg[1]) > TOL + assert np.abs(coefficients_w_reg[2]) > TOL + + +def test_regularization_with_intercept_but_no_intercept_reg_and_unreg_feature(random_data_no_weights): + x = random_data_no_weights[["x0", "x1", "x2", "x3", "x4"]].values + x[:, 0] = 1 + y = random_data_no_weights["y"].values + + lm = OLSRegressionSolver() + lambda_ = 1e6 + lm.fit(x, y, lambda_=lambda_, fit_intercept=True, regularize_intercept=True, n_feat_ignore_reg=2) + coefficients_w_reg = lm.coefficients + assert all(np.abs(coefficients_w_reg[3:] - [0, 0]) <= TOL) + assert np.abs(coefficients_w_reg[0]) < TOL + assert np.abs(coefficients_w_reg[1]) > TOL + assert np.abs(coefficients_w_reg[2]) > TOL + + +################## +# Test residuals # +################## + + +def test_residuals_no_weights(random_data_no_weights): + lm = OLSRegressionSolver() + x = random_data_no_weights[["x0", "x1", "x2", "x3", "x4"]].values + y = random_data_no_weights["y"].values.reshape(-1, 1) + lm.fit(x, y, fit_intercept=False) + predictions = lm.predict(x) + + residuals = lm.residuals(y, predictions, loo=False, center=False) + + assert residuals[0] == pytest.approx(0.885973530) + assert residuals[-1] == pytest.approx(0.841996302) + + residuals = lm.residuals(y, predictions, loo=True, center=False) + + assert residuals[0] == pytest.approx(0.920112164) + assert residuals[-1] == pytest.approx(0.875896477) + + residuals = lm.residuals(y, predictions, loo=True, center=True) + assert np.sum(residuals) == pytest.approx(0) + + +def test_residuals_weights(random_data_weights): + lm = OLSRegressionSolver() + x = random_data_weights[["x0", "x1", "x2", "x3", "x4"]].values + y = random_data_weights["y"].values.reshape(-1, 1) + weights = random_data_weights["weights"].values + + lm.fit(x, y, weights=weights, fit_intercept=False) + predictions = lm.predict(x) + + residuals = lm.residuals(y, predictions, loo=False, center=False) + + assert residuals[0] == pytest.approx(-1.971798590) + assert residuals[-1] == pytest.approx(-1.373951578) + + residuals = lm.residuals(y, predictions, loo=True, center=False) + + assert residuals[0] == pytest.approx(-1.999718445) + assert residuals[-1] == pytest.approx(-1.438563033) + + residuals = lm.residuals(y, predictions, loo=True, center=True) + assert np.sum(residuals) == pytest.approx(0) + + +################################ +# Test saving normal equations # +################################ + + +def test_saving_normal_equations(random_data_no_weights): + lm = OLSRegressionSolver() + x = random_data_no_weights[["x0", "x1", "x2", "x3", "x4"]].values + y = random_data_no_weights["y"].values.reshape(-1, 1) + lm.fit(x, y, fit_intercept=False) + + normal_equations = lm.normal_eqs + + # passing normal equations so should stay the same + x_new = np.ones_like(x) + y_new = np.zeros_like(y) + lm.fit(x_new, y_new, normal_eqs=normal_equations) + np.testing.assert_array_equal(lm.normal_eqs, normal_equations) + + # not passing normal equations so should now change + x_new = random_data_no_weights[["x0", "x1"]].values + y_new = np.zeros_like(y) + lm.fit(x_new, y_new, fit_intercept=False) + with np.testing.assert_raises(AssertionError): + np.testing.assert_array_equal(lm.normal_eqs, normal_equations) diff --git a/tests/test_quantile.py b/tests/test_quantile.py index 72c44e28..792c6bd5 100644 --- a/tests/test_quantile.py +++ b/tests/test_quantile.py @@ -1,7 +1,8 @@ import numpy as np import pytest -from elexsolver.QuantileRegressionSolver import IllConditionedMatrixException, QuantileRegressionSolver +from elexsolver.LinearSolver import IllConditionedMatrixException +from elexsolver.QuantileRegressionSolver import QuantileRegressionSolver # relatively high tolerance, since different implementation. TOL = 1e-3 @@ -22,7 +23,7 @@ def test_basic_median_1(): preds = quantreg.predict(x) # you'd think it would be 8 instead of 7.5, but run quantreg in R to confirm # has to do with missing intercept - assert all(np.abs(preds - [7.5, 7.5, 7.5, 15]) <= TOL) + np.testing.assert_array_equal(preds, [[7.5, 7.5, 7.5, 15]]) def test_basic_median_2(): @@ -32,7 +33,7 @@ def test_basic_median_2(): y = np.asarray([3, 8, 9, 15]) quantreg.fit(x, y, tau) preds = quantreg.predict(x) - assert all(np.abs(preds - [8, 8, 8, 15]) <= TOL) + np.testing.assert_array_equal(preds, [[8, 8, 8, 15]]) def test_basic_lower(): @@ -42,7 +43,7 @@ def test_basic_lower(): y = np.asarray([3, 8, 9, 15]) quantreg.fit(x, y, tau) preds = quantreg.predict(x) - assert all(np.abs(preds - [3, 3, 3, 15]) <= TOL) + np.testing.assert_array_equal(preds, [[3, 3, 3, 15]]) def test_basic_upper(): @@ -52,7 +53,7 @@ def test_basic_upper(): y = np.asarray([3, 8, 9, 15]) quantreg.fit(x, y, tau) preds = quantreg.predict(x) - assert all(np.abs(preds - [9, 9, 9, 15]) <= TOL) + np.testing.assert_array_equal(preds, [[9, 9, 9, 15]]) ###################### @@ -67,7 +68,7 @@ def test_random_median(random_data_no_weights): y = random_data_no_weights["y"].values quantreg.fit(x, y, tau, fit_intercept=False) quantreg.predict(x) - assert all(np.abs(quantreg.coefficients - [1.57699, 6.74906, 4.40175, 4.85346, 4.51814]) <= TOL) + np.testing.assert_allclose(quantreg.coefficients, [[1.57699, 6.74906, 4.40175, 4.85346, 4.51814]], rtol=TOL) def test_random_lower(random_data_no_weights): @@ -77,7 +78,7 @@ def test_random_lower(random_data_no_weights): y = random_data_no_weights["y"].values quantreg.fit(x, y, tau, fit_intercept=False) quantreg.predict(x) - assert all(np.abs(quantreg.coefficients - [0.17759, 6.99588, 4.18896, 4.83906, 3.22546]) <= TOL) + np.testing.assert_allclose(quantreg.coefficients, [[0.17759, 6.99588, 4.18896, 4.83906, 3.22546]], rtol=TOL) def test_random_upper(random_data_no_weights): @@ -87,7 +88,25 @@ def test_random_upper(random_data_no_weights): y = random_data_no_weights["y"].values quantreg.fit(x, y, tau, fit_intercept=False) quantreg.predict(x) - assert all(np.abs(quantreg.coefficients - [1.85617, 6.81286, 6.05586, 5.51965, 4.19864]) <= TOL) + np.testing.assert_allclose(quantreg.coefficients, [[1.85617, 6.81286, 6.05586, 5.51965, 4.19864]], rtol=TOL) + + +################# +# Test multiple # +################# + + +def test_multiple(random_data_no_weights): + quantreg = QuantileRegressionSolver() + taus = [0.1, 0.5, 0.9] + x = random_data_no_weights[["x0", "x1", "x2", "x3", "x4"]].values + y = random_data_no_weights["y"].values + quantreg.fit(x, y, taus, fit_intercept=False) + quantreg.predict(x) + assert len(quantreg.coefficients) == 3 + np.testing.assert_allclose(quantreg.coefficients[0], [0.17759, 6.99588, 4.18896, 4.83906, 3.22546], rtol=TOL) + np.testing.assert_allclose(quantreg.coefficients[1], [1.57699, 6.74906, 4.40175, 4.85346, 4.51814], rtol=TOL) + np.testing.assert_allclose(quantreg.coefficients[2], [1.85617, 6.81286, 6.05586, 5.51965, 4.19864], rtol=TOL) ###################### @@ -103,7 +122,7 @@ def test_basic_median_weights(): weights = np.asarray([1, 1, 100, 3]) quantreg.fit(x, y, tau, weights) preds = quantreg.predict(x) - assert all(np.abs(preds - [9, 9, 9, 15]) <= TOL) + np.testing.assert_array_equal(preds, [[9, 9, 9, 15]]) def test_random_median_weights(random_data_weights): @@ -114,7 +133,7 @@ def test_random_median_weights(random_data_weights): weights = random_data_weights["weights"].values quantreg.fit(x, y, tau, weights=weights, fit_intercept=False) quantreg.predict(x) - assert all(np.abs(quantreg.coefficients - [1.59521, 2.17864, 4.68050, 3.10920, 9.63739]) <= TOL) + np.testing.assert_allclose(quantreg.coefficients, [[1.59521, 2.17864, 4.68050, 3.10920, 9.63739]], rtol=TOL) def test_random_lower_weights(random_data_weights): @@ -125,7 +144,7 @@ def test_random_lower_weights(random_data_weights): weights = random_data_weights["weights"].values quantreg.fit(x, y, tau, weights=weights, fit_intercept=False) quantreg.predict(x) - assert all(np.abs(quantreg.coefficients - [0.63670, 1.27028, 4.81500, 3.08055, 8.69929]) <= TOL) + np.testing.assert_allclose(quantreg.coefficients, [[0.63670, 1.27028, 4.81500, 3.08055, 8.69929]], rtol=TOL) def test_random_upper_weights(random_data_weights): @@ -136,62 +155,7 @@ def test_random_upper_weights(random_data_weights): weights = random_data_weights["weights"].values quantreg.fit(x, y, tau, weights=weights, fit_intercept=False) quantreg.predict(x) - assert all(np.abs(quantreg.coefficients - [3.47742, 2.07360, 4.51754, 4.15237, 9.58856]) <= TOL) - - -######################## -# Test changing solver # -######################## - - -def test_changing_solver(random_data_no_weights): - tau = 0.5 - x = random_data_no_weights[["x0", "x1", "x2", "x3", "x4"]].values - y = random_data_no_weights["y"].values - - quantreg_scs = QuantileRegressionSolver(solver="SCS") - quantreg_ecos = QuantileRegressionSolver(solver="ECOS") - quantreg_scs.fit(x, y, tau, save_problem=True, fit_intercept=False) - quantreg_ecos.fit(x, y, tau, save_problem=True, fit_intercept=False) - - assert quantreg_scs.problem.value == pytest.approx(quantreg_ecos.problem.value, TOL) - - -def test_changing_solver_weights(random_data_weights): - tau = 0.5 - x = random_data_weights[["x0", "x1", "x2", "x3", "x4"]].values - y = random_data_weights["y"].values - weights = random_data_weights["weights"].values - - quantreg_scs = QuantileRegressionSolver(solver="SCS") - quantreg_ecos = QuantileRegressionSolver(solver="ECOS") - quantreg_scs.fit(x, y, tau, weights=weights, save_problem=True, fit_intercept=False) - quantreg_ecos.fit(x, y, tau, weights=weights, save_problem=True, fit_intercept=False) - - assert quantreg_scs.problem.value == pytest.approx(quantreg_ecos.problem.value, TOL) - - -####################### -# Test saving problem # -####################### - - -def test_saving_problem(random_data_no_weights): - tau = 0.5 - x = random_data_no_weights[["x0", "x1", "x2", "x3", "x4"]].values - y = random_data_no_weights["y"].values - - quantreg = QuantileRegressionSolver(solver="ECOS") - - quantreg.fit(x, y, tau, save_problem=False, fit_intercept=False) - assert quantreg.problem is None - - quantreg.fit(x, y, tau, save_problem=True, fit_intercept=False) - assert quantreg.problem is not None - - # testing whether overwrite works - quantreg.fit(x, y, tau, save_problem=False, fit_intercept=False) - assert quantreg.problem is None + np.testing.assert_allclose(quantreg.coefficients, [[3.47742, 2.07360, 4.51754, 4.15237, 9.58856]], rtol=TOL) ############################# @@ -203,9 +167,9 @@ def test_weight_normalization_divide_by_zero(random_data_no_weights): tau = 0.5 x = random_data_no_weights[["x0", "x1", "x2", "x3", "x4"]].values y = random_data_no_weights["y"].values - weights = [0] * x.shape[0] # all zero weights + weights = np.asarray([0] * x.shape[0]) # all zero weights - quantreg = QuantileRegressionSolver(solver="ECOS") + quantreg = QuantileRegressionSolver() # Will succeed without weight normalization quantreg.fit(x, y, tau, normalize_weights=False, weights=weights, fit_intercept=False) @@ -216,20 +180,21 @@ def test_weight_normalization_divide_by_zero(random_data_no_weights): def test_weight_normalization_same_fit(random_data_weights): - quantreg = QuantileRegressionSolver() tau = 0.5 x = np.asarray([[1, 1], [1, 1], [1, 1], [1, 2]]) y = np.asarray([3, 8, 9, 15]) weights = np.asarray([1, 1, 100, 3]) # Test predictions are right when normalize_weights is True and False + quantreg = QuantileRegressionSolver() quantreg.fit(x, y, tau, weights, normalize_weights=True) preds = quantreg.predict(x) - assert all(np.abs(preds - [9, 9, 9, 15]) <= TOL) + np.testing.assert_allclose(preds, [[9, 9, 9, 15]], rtol=TOL) + quantreg = QuantileRegressionSolver() quantreg.fit(x, y, tau, weights, normalize_weights=False) preds = quantreg.predict(x) - assert all(np.abs(preds - [9, 9, 9, 15]) <= TOL) + np.testing.assert_allclose(preds, [[9, 9, 9, 15]], rtol=TOL) ######################## @@ -237,46 +202,47 @@ def test_weight_normalization_same_fit(random_data_weights): ######################## -def test_regularization_with_intercept(random_data_no_weights): +def test_regularization_without_intercept(random_data_no_weights): tau = 0.5 x = random_data_no_weights[["x0", "x1", "x2", "x3", "x4"]].values - x[:, 0] = 1 y = random_data_no_weights["y"].values quantreg = QuantileRegressionSolver() lambda_ = 1e6 - quantreg.fit(x, y, tau, lambda_=lambda_, fit_intercept=True, save_problem=True) - coefficients_w_reg = quantreg.coefficients - assert all(np.abs(coefficients_w_reg[1:] - [0, 0, 0, 0]) <= TOL) - assert np.abs(coefficients_w_reg[0]) > TOL - - objective_w_reg = quantreg.problem.value - quantreg.fit(x, y, tau, save_problem=True) - assert quantreg.problem.value < objective_w_reg + quantreg.fit(x, y, tau, lambda_=lambda_, fit_intercept=False, regularize_intercept=True) + np.testing.assert_allclose( + quantreg.coefficients, [[0, 0, 0, 0, 0]], atol=TOL + ) # using absolute tolerance since comparing to zero -def test_regularization_with_intercept_warning(random_data_no_weights, caplog): - caplog.clear() +def test_regularization_with_intercept(random_data_no_weights): tau = 0.5 x = random_data_no_weights[["x0", "x1", "x2", "x3", "x4"]].values + x[:, 0] = 1 y = random_data_no_weights["y"].values quantreg = QuantileRegressionSolver() lambda_ = 1e6 - with pytest.warns(UserWarning): - quantreg.fit(x, y, tau, lambda_=lambda_, fit_intercept=True, save_problem=True) + quantreg.fit(x, y, tau, lambda_=lambda_, fit_intercept=True, regularize_intercept=False) + coefficients_w_reg = quantreg.coefficients + np.testing.assert_allclose(quantreg.coefficients[0][1:], [0, 0, 0, 0], atol=TOL) + assert np.abs(coefficients_w_reg[0][0]) > TOL -def test_regularization_without_intercept(random_data_no_weights): +def test_regularization_with_intercept_and_features(random_data_no_weights): tau = 0.5 x = random_data_no_weights[["x0", "x1", "x2", "x3", "x4"]].values + x[:, 0] = 1 y = random_data_no_weights["y"].values quantreg = QuantileRegressionSolver() lambda_ = 1e6 - quantreg.fit(x, y, tau, lambda_=lambda_, fit_intercept=False, save_problem=True) + quantreg.fit(x, y, tau, lambda_=lambda_, fit_intercept=True, regularize_intercept=False, n_feat_ignore_reg=2) coefficients_w_reg = quantreg.coefficients - assert all(np.abs(coefficients_w_reg - [0, 0, 0, 0, 0]) <= TOL) + np.testing.assert_allclose(quantreg.coefficients[0][3:], [0, 0], atol=TOL) + assert np.abs(coefficients_w_reg[0][0]) > TOL + assert np.abs(coefficients_w_reg[0][1]) > TOL + assert np.abs(coefficients_w_reg[0][2]) > TOL ######################## diff --git a/tox.ini b/tox.ini index ed82e997..c19a5e08 100644 --- a/tox.ini +++ b/tox.ini @@ -1,5 +1,5 @@ [tox] -envlist=py{3} +envlist=py3.10 skipdist=True [base] @@ -17,6 +17,7 @@ commands= [testenv] deps= {[base]deps} +basepython = python3.10 commands= {[base]commands} pytest --cov-report term-missing --cov=elexsolver \ No newline at end of file