From b7fbb25347b737d097cb25025a3a74c1c1d0ddcd Mon Sep 17 00:00:00 2001 From: lbvienna Date: Wed, 6 Sep 2023 11:29:16 -0400 Subject: [PATCH 01/20] added ols regression and started adding unit tests --- src/elexsolver/LinearSolver.py | 75 +++++++++++++++++ src/elexsolver/OLSRegressionSolver.py | 116 ++++++++++++++++++++++++++ tests/test_ols.py | 57 +++++++++++++ 3 files changed, 248 insertions(+) create mode 100644 src/elexsolver/LinearSolver.py create mode 100644 src/elexsolver/OLSRegressionSolver.py create mode 100644 tests/test_ols.py diff --git a/src/elexsolver/LinearSolver.py b/src/elexsolver/LinearSolver.py new file mode 100644 index 00000000..3798aa0f --- /dev/null +++ b/src/elexsolver/LinearSolver.py @@ -0,0 +1,75 @@ +from abc import ABC +import logging +import warnings + +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 LinearSovler(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 + """ + 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") \ No newline at end of file diff --git a/src/elexsolver/OLSRegressionSolver.py b/src/elexsolver/OLSRegressionSolver.py new file mode 100644 index 00000000..19bd1021 --- /dev/null +++ b/src/elexsolver/OLSRegressionSolver.py @@ -0,0 +1,116 @@ +import logging + +import numpy as np + +from elexsolver.LinearSolver import LinearSovler +from elexsolver.logging import initialize_logging + +initialize_logging() + +LOG = logging.getLogger(__name__) + +class OLSRegressionSolver(LinearSovler): + """ + 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_matrix = None + + def _get_regularizer(self, lambda_: float, dim: int, fit_intercept: bool, n_feat_ignore_req: 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_req features + for i in range(fit_intercept + n_feat_ignore_req): + lambda_I[i, i] = 0 + + return lambda_I + + def _compute_normal_equations(self, x: np.ndarray, L: np.ndarray, lambda_: float, fit_intercept: bool, n_feat_ignore_req: 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, n_feat_ignore_req) + + # 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, n_feat_ignore_req: 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, n_feat_ignore_req) + 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 \ No newline at end of file diff --git a/tests/test_ols.py b/tests/test_ols.py new file mode 100644 index 00000000..d3a4e19e --- /dev/null +++ b/tests/test_ols.py @@ -0,0 +1,57 @@ +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_basic_median_1(): + 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_basic_median_2(): + 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_median(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) + lm.predict(x) + assert all(np.abs(lm.coefficients - [[1.037], [7.022], [4.794], [4.776], [4.266]]) <= TOL) + + +###################### +# Tests with weights # +###################### + +def test_random_median_weights(random_data_weights): + lm = OLSRegressionSolver() + tau = 0.5 + 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) + lm.predict(x) + assert all(np.abs(lm.coefficients - [[1.455], [2.018], [4.699], [3.342], [9.669]]) <= TOL) \ No newline at end of file From 47a48c991736cd58bba4400467169aa59d0bfd40 Mon Sep 17 00:00:00 2001 From: lbvienna Date: Wed, 6 Sep 2023 15:45:42 -0400 Subject: [PATCH 02/20] made progress on unit tests --- src/elexsolver/OLSRegressionSolver.py | 2 +- tests/test_ols.py | 115 ++++++++++++++++++++++++-- 2 files changed, 108 insertions(+), 9 deletions(-) diff --git a/src/elexsolver/OLSRegressionSolver.py b/src/elexsolver/OLSRegressionSolver.py index 19bd1021..426e7df4 100644 --- a/src/elexsolver/OLSRegressionSolver.py +++ b/src/elexsolver/OLSRegressionSolver.py @@ -27,7 +27,7 @@ class OLSRegressionSolver(LinearSovler): def __init__(self): super().__init__() self.normal_eqs = None - self.hat_matrix = None + self.hat_vals = None def _get_regularizer(self, lambda_: float, dim: int, fit_intercept: bool, n_feat_ignore_req: int) -> np.ndarray: """ diff --git a/tests/test_ols.py b/tests/test_ols.py index d3a4e19e..38fb698e 100644 --- a/tests/test_ols.py +++ b/tests/test_ols.py @@ -12,7 +12,7 @@ # Basic tests # ############### -def test_basic_median_1(): +def test_basic1(): lm = OLSRegressionSolver() x = np.asarray([[1], [1], [1], [2]]) y = np.asarray([3, 8, 9, 15]) @@ -20,7 +20,7 @@ def test_basic_median_1(): preds = lm.predict(x) assert all(np.abs(preds - [7.142857, 7.142857, 7.142857, 14.285714]) <= TOL) -def test_basic_median_2(): +def test_basic2(): lm = OLSRegressionSolver() x = np.asarray([[1, 1], [1, 1], [1, 1], [1, 2]]) y = np.asarray([3, 8, 9, 15]) @@ -33,25 +33,124 @@ def test_basic_median_2(): ###################### -def test_random_median(random_data_no_weights): +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) - lm.predict(x) + 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_median_weights(random_data_weights): +def test_random_weights_no_intercept(random_data_weights): lm = OLSRegressionSolver() - tau = 0.5 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) - lm.predict(x) - assert all(np.abs(lm.coefficients - [[1.455], [2.018], [4.699], [3.342], [9.669]]) <= TOL) \ No newline at end of file + 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_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) + 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, n_feat_ignore_req=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 # +################## + +################################ +# Test saving normal equations # +################################ \ No newline at end of file From 248585a5d2a46cb052edb1cc8189bb2e8bc869c6 Mon Sep 17 00:00:00 2001 From: lbvienna Date: Wed, 6 Sep 2023 18:06:46 -0400 Subject: [PATCH 03/20] added unit tests --- tests/test_ols.py | 79 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 78 insertions(+), 1 deletion(-) diff --git a/tests/test_ols.py b/tests/test_ols.py index 38fb698e..4b4c4f05 100644 --- a/tests/test_ols.py +++ b/tests/test_ols.py @@ -120,6 +120,20 @@ def test_random_weights_intercept(random_data_weights): # Test regularization # ######################## +def test_regularizer(): + lm = OLSRegressionSolver() + + lambda_I = lm._get_regularizer(7, 10, fit_intercept=True, n_feat_ignore_req=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 @@ -151,6 +165,69 @@ def test_regularization_with_intercept_and_unreg_feature(random_data_no_weights) # 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 # -################################ \ No newline at end of file +################################ + +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) \ No newline at end of file From 761580741cd0ac579353a41389ae9615e545e27d Mon Sep 17 00:00:00 2001 From: lbvienna Date: Wed, 6 Sep 2023 19:03:39 -0400 Subject: [PATCH 04/20] finalized ols tests --- src/elexsolver/LinearSolver.py | 2 +- src/elexsolver/OLSRegressionSolver.py | 2 +- tests/test_linear_solver.py | 10 ++++++++++ 3 files changed, 12 insertions(+), 2 deletions(-) create mode 100644 tests/test_linear_solver.py diff --git a/src/elexsolver/LinearSolver.py b/src/elexsolver/LinearSolver.py index 3798aa0f..eb8675cb 100644 --- a/src/elexsolver/LinearSolver.py +++ b/src/elexsolver/LinearSolver.py @@ -17,7 +17,7 @@ class LinearSolverException(Exception): class IllConditionedMatrixException(LinearSolverException): pass -class LinearSovler(ABC): +class LinearSolver(ABC): """ An abstract base class for a linear solver """ diff --git a/src/elexsolver/OLSRegressionSolver.py b/src/elexsolver/OLSRegressionSolver.py index 426e7df4..7d0223c5 100644 --- a/src/elexsolver/OLSRegressionSolver.py +++ b/src/elexsolver/OLSRegressionSolver.py @@ -9,7 +9,7 @@ LOG = logging.getLogger(__name__) -class OLSRegressionSolver(LinearSovler): +class OLSRegressionSolver(LinearSolver): """ A class for Ordinary Least Squares Regression optimized for the bootstrap """ diff --git a/tests/test_linear_solver.py b/tests/test_linear_solver.py new file mode 100644 index 00000000..87d49552 --- /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))) + From ad22afcc4479060dd4e55dbfb517b163cc2f308f Mon Sep 17 00:00:00 2001 From: lbvienna Date: Wed, 6 Sep 2023 19:37:24 -0400 Subject: [PATCH 05/20] fixed typo --- src/elexsolver/OLSRegressionSolver.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/elexsolver/OLSRegressionSolver.py b/src/elexsolver/OLSRegressionSolver.py index 7d0223c5..3a70029b 100644 --- a/src/elexsolver/OLSRegressionSolver.py +++ b/src/elexsolver/OLSRegressionSolver.py @@ -2,7 +2,7 @@ import numpy as np -from elexsolver.LinearSolver import LinearSovler +from elexsolver.LinearSolver import LinearSolver from elexsolver.logging import initialize_logging initialize_logging() From a453fcaa8354b3d7272741fdde0e6806f275a4e3 Mon Sep 17 00:00:00 2001 From: lbvienna Date: Wed, 6 Sep 2023 22:47:52 -0400 Subject: [PATCH 06/20] possibility to regularize intercept or not now --- src/elexsolver/OLSRegressionSolver.py | 13 ++++++----- tests/test_ols.py | 32 ++++++++++++++++++++++++--- 2 files changed, 37 insertions(+), 8 deletions(-) diff --git a/src/elexsolver/OLSRegressionSolver.py b/src/elexsolver/OLSRegressionSolver.py index 3a70029b..23cc3755 100644 --- a/src/elexsolver/OLSRegressionSolver.py +++ b/src/elexsolver/OLSRegressionSolver.py @@ -29,7 +29,7 @@ def __init__(self): self.normal_eqs = None self.hat_vals = None - def _get_regularizer(self, lambda_: float, dim: int, fit_intercept: bool, n_feat_ignore_req: int) -> np.ndarray: + def _get_regularizer(self, lambda_: float, dim: int, fit_intercept: bool, regularize_intercept: bool, n_feat_ignore_req: int) -> np.ndarray: """ Returns the regularization matrix """ @@ -43,11 +43,14 @@ def _get_regularizer(self, lambda_: float, dim: int, fit_intercept: bool, n_feat # so set regularization constant to zero for intercept # and the first n_feat_ignore_req features for i in range(fit_intercept + n_feat_ignore_req): + # 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, n_feat_ignore_req: int) -> np.ndarray: + def _compute_normal_equations(self, x: np.ndarray, L: np.ndarray, lambda_: float, fit_intercept: bool, regularize_intercept: bool, n_feat_ignore_req: int) -> np.ndarray: """ Computes the normal equations for OLS: (X^T X)^{-1} X^T """ @@ -57,7 +60,7 @@ def _compute_normal_equations(self, x: np.ndarray, L: np.ndarray, lambda_: float Q, R = np.linalg.qr(L @ x) # get regularization matrix - lambda_I = self._get_regularizer(lambda_, R.shape[0], fit_intercept, n_feat_ignore_req) + lambda_I = self._get_regularizer(lambda_, R.shape[0], fit_intercept, regularize_intercept, n_feat_ignore_req) # substitute X = QR into the normal equations to get # R^T Q^T Q R \beta = R^T Q^T y @@ -67,7 +70,7 @@ def _compute_normal_equations(self, x: np.ndarray, L: np.ndarray, lambda_: float # 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, n_feat_ignore_req: int = 0): + 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_req: int = 0): self._check_any_element_nan_or_inf(x) self._check_any_element_nan_or_inf(y) @@ -86,7 +89,7 @@ def fit(self, x: np.ndarray, y: np.ndarray, weights: np.ndarray | None = None, l # 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, n_feat_ignore_req) + self.normal_eqs = self._compute_normal_equations(x, L, lambda_, fit_intercept, regularize_intercept, n_feat_ignore_req) else: self.normal_eqs = normal_eqs diff --git a/tests/test_ols.py b/tests/test_ols.py index 4b4c4f05..1e3663cb 100644 --- a/tests/test_ols.py +++ b/tests/test_ols.py @@ -123,7 +123,19 @@ def test_random_weights_intercept(random_data_weights): def test_regularizer(): lm = OLSRegressionSolver() - lambda_I = lm._get_regularizer(7, 10, fit_intercept=True, n_feat_ignore_req=2) + lambda_I = lm._get_regularizer(7, 10, fit_intercept=True, regularize_intercept=True, n_feat_ignore_req=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_req=2) assert lambda_I.shape == (10, 10) assert lambda_I[0, 0] == pytest.approx(0) @@ -142,7 +154,7 @@ def test_regularization_with_intercept(random_data_no_weights): lm = OLSRegressionSolver() lambda_ = 1e6 - lm.fit(x, y, lambda_=lambda_, fit_intercept=True) + 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 @@ -154,13 +166,27 @@ def test_regularization_with_intercept_and_unreg_feature(random_data_no_weights) lm = OLSRegressionSolver() lambda_ = 1e6 - lm.fit(x, y, lambda_=lambda_, fit_intercept=True, n_feat_ignore_req=2) + lm.fit(x, y, lambda_=lambda_, fit_intercept=True, regularize_intercept=False, n_feat_ignore_req=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_req=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 # ################## From c29027337087aa6a2210816f572a233e2b9424af Mon Sep 17 00:00:00 2001 From: lbvienna Date: Fri, 8 Sep 2023 09:30:36 -0400 Subject: [PATCH 07/20] pinned scipy dependency to avoid unit test issue --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index a1e5824e..07365d2c 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__) From 7c8f6cdaaf7397a6dfa91af1464b5d8ec8d29f8a Mon Sep 17 00:00:00 2001 From: lbvienna Date: Fri, 8 Sep 2023 17:30:23 -0400 Subject: [PATCH 08/20] added reg, working on unit tests --- src/elexsolver/LinearSolver.py | 2 + src/elexsolver/QuantileRegressionSolver.py | 221 +++++++++------------ tests/test_quantile.py | 125 ++---------- 3 files changed, 117 insertions(+), 231 deletions(-) diff --git a/src/elexsolver/LinearSolver.py b/src/elexsolver/LinearSolver.py index eb8675cb..5bbc424c 100644 --- a/src/elexsolver/LinearSolver.py +++ b/src/elexsolver/LinearSolver.py @@ -39,6 +39,8 @@ 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: diff --git a/src/elexsolver/QuantileRegressionSolver.py b/src/elexsolver/QuantileRegressionSolver.py index 3fe12560..6ddb6762 100644 --- a/src/elexsolver/QuantileRegressionSolver.py +++ b/src/elexsolver/QuantileRegressionSolver.py @@ -1,9 +1,10 @@ import logging -import warnings +from scipy.optimize import linprog import cvxpy as cp import numpy as np +from elexsolver.LinearSolver import LinearSolver from elexsolver.logging import initialize_logging initialize_logging() @@ -11,134 +12,110 @@ LOG = logging.getLogger(__name__) -class QuantileRegressionSolverException(Exception): - pass - - -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): - """ - 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") - - def get_loss_function(self, x, y, coefficients, weights): +class QuantileRegressionSolver(LinearSolver): + + def __init__(self): + super().__init__() + self.coefficients = [] + + def _fit( + self, x: np.ndarray, y: np.ndarray, weights: np.ndarray, tau: float + ) -> np.ndarray: + """ + Fits the dual problem of a quantile regression, for more information see appendix 6 here: https://arxiv.org/pdf/2305.12616.pdf + """ + 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) + + # 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 _fit_with_regularization(self, x: np.ndarray, y: np.ndarray, weights: np.ndarray, tau: float, lambda_: float): + S = cp.Constant(y.reshape(-1, 1)) + N = y.shape[0] + Phi = cp.Constant(x) + radius = 1 / lambda_ + gamma = cp.Variable() + C = radius / (N + 1) + eta = cp.Variable(name="weights", shape=N) + constraints = [ + C * (tau - 1) <= eta, + C * tau >= eta, + eta.T @ Phi == gamma + ] + prob = cp.Problem( + cp.Minimize(0.5 * cp.sum_squares(Phi.T @ eta) - cp.sum(cp.multiply(eta, cp.vec(S)))), + constraints + ) + prob.solve() + + return prob.constraints[-1].dual_value + """ + S -> scores + tau -> quantile + eta -> optimization_variables + + eta = cp.Variable(name="weights", shape=n_calib) + + scores = cp.Constant(scores_calib.reshape(-1,1)) + + Phi = cp.Constant(phi_calib) + + radius = 1 / infinite_params.get('lambda', FUNCTION_DEFAULTS['lambda']) + + C = radius / (n_calib + 1) + + constraints = [ + C * (quantile - 1) <= eta, + C * quantile >= eta, + eta.T @ Phi == 0] + prob = cp.Problem( + cp.Minimize(0.5 * cp.sum_squares(eta) - cp.sum(cp.multiply(eta, cp.vec(scores)))), + constraints + ) + + coefficients = prob.constraints[-1].dual_value + """ + + def fit(self, x: np.ndarray, y: np.ndarray, taus: list | float = 0.5, weights: np.ndarray | None = None, lambda_: float = 0.0, fit_intercept: bool = True) -> np.ndarray: + """ + Fits quantile regression """ - 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)) - - def get_regularizer(self, coefficients, fit_intercept): - """ - 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:] - return cp.pnorm(coefficients_to_regularize, p=2) ** 2 - - def __solve(self, x, y, weights, lambda_, fit_intercept, verbose): - """ - Sets up the optimization problem and solves it - """ - self._check_matrix_condition(x) - 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) - objective = cp.Minimize(loss_function) - problem = cp.Problem(objective) - problem.solve(solver=self.solver, verbose=verbose, **self.KWARGS.get(self.solver, {})) - return coefficients, problem - - def fit( - self, - x, - y, - tau_value=0.5, - weights=None, - lambda_=0, - fit_intercept=True, - verbose=False, - save_problem=False, - normalize_weights=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` - """ - 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 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 - - def predict(self, x): + # if weights are none, all rows should be weighted equally + if weights is None: + weights = np.ones((y.shape[0], )) + + # normalize weights + weights = weights / np.sum(weights) + + # _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_) + else: + coefficients = self._fit(x, y, weights, tau) + self.coefficients.append(coefficients) + + def predict(self, x: np.ndarray) -> np.ndarray: """ - Returns predictions + Use coefficients to predict """ self._check_any_element_nan_or_inf(x) - return self.coefficients @ x.T + return self.coefficients @ x.T \ No newline at end of file diff --git a/tests/test_quantile.py b/tests/test_quantile.py index 72c44e28..d06006e7 100644 --- a/tests/test_quantile.py +++ b/tests/test_quantile.py @@ -1,7 +1,7 @@ import numpy as np import pytest -from elexsolver.QuantileRegressionSolver import IllConditionedMatrixException, QuantileRegressionSolver +from elexsolver.QuantileRegressionSolver import QuantileRegressionSolver # relatively high tolerance, since different implementation. TOL = 1e-3 @@ -22,7 +22,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 +32,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 +42,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 +52,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,8 +67,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): quantreg = QuantileRegressionSolver() @@ -77,7 +76,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 +86,7 @@ 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) ###################### @@ -103,7 +102,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 +113,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 +124,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,100 +135,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 - - -############################# -# Test weight normalization # -############################# - - -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 - - quantreg = QuantileRegressionSolver(solver="ECOS") - - # Will succeed without weight normalization - quantreg.fit(x, y, tau, normalize_weights=False, weights=weights, fit_intercept=False) - - # Will fail with weight normalization - with pytest.raises(ZeroDivisionError): - quantreg.fit(x, y, tau, normalize_weights=True, weights=weights, fit_intercept=False) - - -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.fit(x, y, tau, weights, normalize_weights=True) - preds = quantreg.predict(x) - assert all(np.abs(preds - [9, 9, 9, 15]) <= TOL) - - 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(quantreg.coefficients, [[3.47742, 2.07360, 4.51754, 4.15237, 9.58856]], rtol=TOL) ######################## @@ -245,8 +151,9 @@ def test_regularization_with_intercept(random_data_no_weights): quantreg = QuantileRegressionSolver() lambda_ = 1e6 - quantreg.fit(x, y, tau, lambda_=lambda_, fit_intercept=True, save_problem=True) + quantreg.fit(x, y, tau, lambda_=lambda_, fit_intercept=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 @@ -264,7 +171,7 @@ def test_regularization_with_intercept_warning(random_data_no_weights, caplog): 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) def test_regularization_without_intercept(random_data_no_weights): @@ -274,7 +181,7 @@ def test_regularization_without_intercept(random_data_no_weights): 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=False) coefficients_w_reg = quantreg.coefficients assert all(np.abs(coefficients_w_reg - [0, 0, 0, 0, 0]) <= TOL) From 296a62401def82f4d4e43be5bfdbac0d2b8f3c68 Mon Sep 17 00:00:00 2001 From: lbvienna Date: Mon, 18 Sep 2023 22:11:49 -0400 Subject: [PATCH 09/20] updated all tests --- src/elexsolver/OLSRegressionSolver.py | 14 ++-- src/elexsolver/QuantileRegressionSolver.py | 74 ++++++++-------------- tests/test_quantile.py | 51 +++++++++------ 3 files changed, 67 insertions(+), 72 deletions(-) diff --git a/src/elexsolver/OLSRegressionSolver.py b/src/elexsolver/OLSRegressionSolver.py index 23cc3755..41d7e694 100644 --- a/src/elexsolver/OLSRegressionSolver.py +++ b/src/elexsolver/OLSRegressionSolver.py @@ -29,7 +29,7 @@ def __init__(self): 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_req: int) -> np.ndarray: + 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 """ @@ -41,8 +41,8 @@ def _get_regularizer(self, lambda_: float, dim: int, fit_intercept: bool, regula # 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_req features - for i in range(fit_intercept + n_feat_ignore_req): + # 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 @@ -50,7 +50,7 @@ def _get_regularizer(self, lambda_: float, dim: int, fit_intercept: bool, regula 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_req: int) -> np.ndarray: + 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 """ @@ -60,7 +60,7 @@ def _compute_normal_equations(self, x: np.ndarray, L: np.ndarray, lambda_: float 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_req) + 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 @@ -70,7 +70,7 @@ def _compute_normal_equations(self, x: np.ndarray, L: np.ndarray, lambda_: float # 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_req: int = 0): + 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) @@ -89,7 +89,7 @@ def fit(self, x: np.ndarray, y: np.ndarray, weights: np.ndarray | None = None, l # 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_req) + self.normal_eqs = self._compute_normal_equations(x, L, lambda_, fit_intercept, regularize_intercept, n_feat_ignore_reg) else: self.normal_eqs = normal_eqs diff --git a/src/elexsolver/QuantileRegressionSolver.py b/src/elexsolver/QuantileRegressionSolver.py index 6ddb6762..30343b91 100644 --- a/src/elexsolver/QuantileRegressionSolver.py +++ b/src/elexsolver/QuantileRegressionSolver.py @@ -37,54 +37,34 @@ def _fit( # marginal are the dual values, since we are solving the dual this is equivalent to the primal return -1 * res.eqlin.marginals - def _fit_with_regularization(self, x: np.ndarray, y: np.ndarray, weights: np.ndarray, tau: float, lambda_: float): - S = cp.Constant(y.reshape(-1, 1)) - N = y.shape[0] - Phi = cp.Constant(x) - radius = 1 / lambda_ - gamma = cp.Variable() - C = radius / (N + 1) - eta = cp.Variable(name="weights", shape=N) - constraints = [ - C * (tau - 1) <= eta, - C * tau >= eta, - eta.T @ Phi == gamma - ] - prob = cp.Problem( - cp.Minimize(0.5 * cp.sum_squares(Phi.T @ eta) - cp.sum(cp.multiply(eta, cp.vec(S)))), - constraints - ) - prob.solve() - - return prob.constraints[-1].dual_value + def _get_regularizer(self, coefficients: cp.expressions.variable.Variable, regularize_intercept: bool, n_feat_ignore_reg: int) -> cp.atoms.norm: """ - S -> scores - tau -> quantile - eta -> optimization_variables - - eta = cp.Variable(name="weights", shape=n_calib) - - scores = cp.Constant(scores_calib.reshape(-1,1)) - - Phi = cp.Constant(phi_calib) - - radius = 1 / infinite_params.get('lambda', FUNCTION_DEFAULTS['lambda']) - - C = radius / (n_calib + 1) - - constraints = [ - C * (quantile - 1) <= eta, - C * quantile >= eta, - eta.T @ Phi == 0] - prob = cp.Problem( - cp.Minimize(0.5 * cp.sum_squares(eta) - cp.sum(cp.multiply(eta, cp.vec(scores)))), - constraints - ) - - coefficients = prob.constraints[-1].dual_value - """ + Get regularization component of the loss function. Note that this is L2 (ridge) regularization. + """ + # this assumes that if regularize_intercept=True that the intercept is the first column + coefficients_to_regularize = coefficients[n_feat_ignore_reg: ] + if not regularize_intercept: + coefficients_to_regularize = coefficients[1 + n_feat_ignore_reg: ] + return cp.pnorm(coefficients_to_regularize, p=2) ** 2 - def fit(self, x: np.ndarray, y: np.ndarray, taus: list | float = 0.5, weights: np.ndarray | None = None, lambda_: float = 0.0, fit_intercept: bool = True) -> np.ndarray: + 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): + """ + Fits quantile regression with regularization + TODO: convert this problem to use the dual like in the non regularization case + """ + arguments = {"ECOS": {"max_iters": 10000}} + coefficients = cp.Variable((x.shape[1], )) + 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='ECOS', **arguments.get('ECOS', {})) + return coefficients.value + + + def fit(self, 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) -> np.ndarray: """ Fits quantile regression """ @@ -107,7 +87,7 @@ def fit(self, x: np.ndarray, y: np.ndarray, taus: list | float = 0.5, weights: n for tau in taus: if lambda_ > 0: - coefficients = self._fit_with_regularization(x, y, weights, tau, lambda_) + 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) diff --git a/tests/test_quantile.py b/tests/test_quantile.py index d06006e7..383a2095 100644 --- a/tests/test_quantile.py +++ b/tests/test_quantile.py @@ -2,6 +2,7 @@ import pytest from elexsolver.QuantileRegressionSolver import QuantileRegressionSolver +from elexsolver.LinearSolver import IllConditionedMatrixException # relatively high tolerance, since different implementation. TOL = 1e-3 @@ -89,6 +90,22 @@ def test_random_upper(random_data_no_weights): 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) + ###################### # Tests with weights # ###################### @@ -143,47 +160,45 @@ def test_random_upper_weights(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) - 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) + 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) + 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 ######################## From 936b2dd3d7d3371b62ecb544b87fade4cf844d6c Mon Sep 17 00:00:00 2001 From: lbvienna Date: Mon, 18 Sep 2023 22:12:41 -0400 Subject: [PATCH 10/20] ran pre-commit hook --- src/elexsolver/LinearSolver.py | 14 +++--- src/elexsolver/OLSRegressionSolver.py | 58 +++++++++++++++------- src/elexsolver/QuantileRegressionSolver.py | 53 +++++++++++++------- tests/test_linear_solver.py | 2 +- tests/test_ols.py | 40 ++++++++++----- tests/test_quantile.py | 9 +++- 6 files changed, 121 insertions(+), 55 deletions(-) diff --git a/src/elexsolver/LinearSolver.py b/src/elexsolver/LinearSolver.py index 5bbc424c..949cb5bd 100644 --- a/src/elexsolver/LinearSolver.py +++ b/src/elexsolver/LinearSolver.py @@ -1,6 +1,6 @@ -from abc import ABC import logging import warnings +from abc import ABC import numpy as np @@ -10,6 +10,7 @@ LOG = logging.getLogger(__name__) + class LinearSolverException(Exception): pass @@ -17,6 +18,7 @@ class LinearSolverException(Exception): class IllConditionedMatrixException(LinearSolverException): pass + class LinearSolver(ABC): """ An abstract base class for a linear solver @@ -34,7 +36,7 @@ def fit(self, x: np.ndarray, y: np.ndarray, weights: np.ndarray | None = None, l Fits model """ raise NotImplementedError - + def predict(self, x: np.ndarray) -> np.ndarray: """ Use coefficients to predict @@ -42,13 +44,13 @@ def predict(self, x: np.ndarray) -> np.ndarray: self._check_any_element_nan_or_inf(x) return x @ self.coefficients - - def get_coefficients(self) -> np.ndarray: + + 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. @@ -74,4 +76,4 @@ 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") \ No newline at end of file + 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 index 41d7e694..ca055354 100644 --- a/src/elexsolver/OLSRegressionSolver.py +++ b/src/elexsolver/OLSRegressionSolver.py @@ -9,19 +9,20 @@ LOG = logging.getLogger(__name__) + class OLSRegressionSolver(LinearSolver): """ A class for Ordinary Least Squares Regression optimized for the bootstrap """ - # OLS setup: + # 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 + # 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): @@ -29,15 +30,17 @@ def __init__(self): 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: + 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 + + # 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 @@ -45,12 +48,21 @@ def _get_regularizer(self, lambda_: float, dim: int, fit_intercept: bool, regula 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 + 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: + + 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 """ @@ -70,16 +82,26 @@ def _compute_normal_equations(self, x: np.ndarray, L: np.ndarray, lambda_: float # 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): + 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], )) + 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 @@ -89,22 +111,24 @@ def fit(self, x: np.ndarray, y: np.ndarray, weights: np.ndarray | None = None, l # 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) + 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) + 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 @@ -116,4 +140,4 @@ def residuals(self, y: np.ndarray, y_hat: np.ndarray, loo: bool = True, center: if center: residuals -= np.mean(residuals, axis=0) - return residuals \ No newline at end of file + return residuals diff --git a/src/elexsolver/QuantileRegressionSolver.py b/src/elexsolver/QuantileRegressionSolver.py index 30343b91..1ab9a64b 100644 --- a/src/elexsolver/QuantileRegressionSolver.py +++ b/src/elexsolver/QuantileRegressionSolver.py @@ -1,8 +1,8 @@ import logging -from scipy.optimize import linprog import cvxpy as cp import numpy as np +from scipy.optimize import linprog from elexsolver.LinearSolver import LinearSolver from elexsolver.logging import initialize_logging @@ -13,14 +13,11 @@ class QuantileRegressionSolver(LinearSolver): - def __init__(self): super().__init__() self.coefficients = [] - def _fit( - self, x: np.ndarray, y: np.ndarray, weights: np.ndarray, tau: float - ) -> np.ndarray: + def _fit(self, x: np.ndarray, y: np.ndarray, weights: np.ndarray, tau: float) -> np.ndarray: """ Fits the dual problem of a quantile regression, for more information see appendix 6 here: https://arxiv.org/pdf/2305.12616.pdf """ @@ -37,46 +34,66 @@ def _fit( # 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: cp.expressions.variable.Variable, regularize_intercept: bool, n_feat_ignore_reg: int) -> cp.atoms.norm: + 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. """ # this assumes that if regularize_intercept=True that the intercept is the first column - coefficients_to_regularize = coefficients[n_feat_ignore_reg: ] + coefficients_to_regularize = coefficients[n_feat_ignore_reg:] if not regularize_intercept: - coefficients_to_regularize = coefficients[1 + n_feat_ignore_reg: ] + coefficients_to_regularize = coefficients[1 + n_feat_ignore_reg :] # noqa: E203 return cp.pnorm(coefficients_to_regularize, p=2) ** 2 - 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): + 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, + ): """ Fits quantile regression with regularization TODO: convert this problem to use the dual like in the non regularization case """ arguments = {"ECOS": {"max_iters": 10000}} - coefficients = cp.Variable((x.shape[1], )) + coefficients = cp.Variable((x.shape[1],)) 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='ECOS', **arguments.get('ECOS', {})) + problem.solve(solver="ECOS", **arguments.get("ECOS", {})) return coefficients.value - - def fit(self, 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) -> np.ndarray: + def fit( + self, + 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, + ) -> np.ndarray: """ 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 are none, all rows should be weighted equally if weights is None: - weights = np.ones((y.shape[0], )) + weights = np.ones((y.shape[0],)) # normalize weights weights = weights / np.sum(weights) @@ -87,7 +104,9 @@ def fit(self, x: np.ndarray, y: np.ndarray, taus: list | float = 0.5, weights: n for tau in taus: if lambda_ > 0: - coefficients = self._fit_with_regularization(x, y, weights, tau, lambda_, regularize_intercept, n_feat_ignore_reg) + 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) @@ -98,4 +117,4 @@ def predict(self, x: np.ndarray) -> np.ndarray: """ self._check_any_element_nan_or_inf(x) - return self.coefficients @ x.T \ No newline at end of file + return self.coefficients @ x.T diff --git a/tests/test_linear_solver.py b/tests/test_linear_solver.py index 87d49552..98e1f163 100644 --- a/tests/test_linear_solver.py +++ b/tests/test_linear_solver.py @@ -3,8 +3,8 @@ 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 index 1e3663cb..02d4da60 100644 --- a/tests/test_ols.py +++ b/tests/test_ols.py @@ -12,6 +12,7 @@ # Basic tests # ############### + def test_basic1(): lm = OLSRegressionSolver() x = np.asarray([[1], [1], [1], [2]]) @@ -20,6 +21,7 @@ def test_basic1(): 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]]) @@ -28,6 +30,7 @@ def test_basic2(): preds = lm.predict(x) assert all(np.abs(preds - [6.666667, 6.666667, 6.666667, 15]) <= TOL) + ###################### # Intermediate tests # ###################### @@ -41,7 +44,7 @@ def test_random_no_intercept(random_data_no_weights): 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) @@ -51,16 +54,17 @@ def test_random_no_intercept(random_data_no_weights): 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 + 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) @@ -75,6 +79,7 @@ def test_random_intercept(random_data_no_weights): # Tests with weights # ###################### + def test_random_weights_no_intercept(random_data_weights): lm = OLSRegressionSolver() x = random_data_weights[["x0", "x1", "x2", "x3", "x4"]].values @@ -82,8 +87,8 @@ def test_random_weights_no_intercept(random_data_weights): weights = random_data_weights["weights"].values lm.fit(x, y, weights=weights, fit_intercept=False) predictions = lm.predict(x) - - #coefficients + + # coefficients assert all(np.abs(lm.coefficients - [[1.455], [2.018], [4.699], [3.342], [9.669]]) <= TOL) # check hat values @@ -95,6 +100,7 @@ def test_random_weights_no_intercept(random_data_weights): 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 @@ -103,23 +109,25 @@ def test_random_weights_intercept(random_data_weights): weights = random_data_weights["weights"].values lm.fit(x, y, weights=weights, fit_intercept=True) predictions = lm.predict(x) - - #coefficients + + # 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 ) + 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() @@ -134,7 +142,7 @@ def test_regularizer(): 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_req=2) assert lambda_I.shape == (10, 10) @@ -147,6 +155,7 @@ def test_regularizer(): 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 @@ -159,6 +168,7 @@ def test_regularization_with_intercept(random_data_no_weights): 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 @@ -173,6 +183,7 @@ def test_regularization_with_intercept_and_unreg_feature(random_data_no_weights) 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 @@ -187,10 +198,12 @@ def test_regularization_with_intercept_but_no_intercept_reg_and_unreg_feature(ra 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 @@ -202,7 +215,7 @@ def test_residuals_no_weights(random_data_no_weights): 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) @@ -211,6 +224,7 @@ def test_residuals_no_weights(random_data_no_weights): 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 @@ -222,7 +236,7 @@ def test_residuals_weights(random_data_weights): residuals = lm.residuals(y, predictions, loo=False, center=False) - assert residuals[0] == pytest.approx(-1.971798590 ) + assert residuals[0] == pytest.approx(-1.971798590) assert residuals[-1] == pytest.approx(-1.373951578) residuals = lm.residuals(y, predictions, loo=True, center=False) @@ -233,10 +247,12 @@ def test_residuals_weights(random_data_weights): 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 @@ -256,4 +272,4 @@ def test_saving_normal_equations(random_data_no_weights): 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) \ No newline at end of file + np.testing.assert_array_equal(lm.normal_eqs, normal_equations) diff --git a/tests/test_quantile.py b/tests/test_quantile.py index 383a2095..a172d38c 100644 --- a/tests/test_quantile.py +++ b/tests/test_quantile.py @@ -1,8 +1,8 @@ import numpy as np import pytest -from elexsolver.QuantileRegressionSolver import QuantileRegressionSolver from elexsolver.LinearSolver import IllConditionedMatrixException +from elexsolver.QuantileRegressionSolver import QuantileRegressionSolver # relatively high tolerance, since different implementation. TOL = 1e-3 @@ -70,6 +70,7 @@ def test_random_median(random_data_no_weights): quantreg.predict(x) 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): quantreg = QuantileRegressionSolver() tau = 0.1 @@ -94,6 +95,7 @@ def test_random_upper(random_data_no_weights): # Test multiple # ################# + def test_multiple(random_data_no_weights): quantreg = QuantileRegressionSolver() taus = [0.1, 0.5, 0.9] @@ -106,6 +108,7 @@ def test_multiple(random_data_no_weights): 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) + ###################### # Tests with weights # ###################### @@ -168,7 +171,9 @@ def test_regularization_without_intercept(random_data_no_weights): quantreg = QuantileRegressionSolver() lambda_ = 1e6 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 + 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(random_data_no_weights): From 21cdfbca815894c3d3090cdd66f9005268b9690e Mon Sep 17 00:00:00 2001 From: lbvienna Date: Mon, 18 Sep 2023 22:16:17 -0400 Subject: [PATCH 11/20] updated readme --- README.md | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) 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. From fde1457f6334d9cc82073bd21816b89213eec8b3 Mon Sep 17 00:00:00 2001 From: lbvienna Date: Mon, 18 Sep 2023 22:26:35 -0400 Subject: [PATCH 12/20] updated python versions the unit tests run for --- tests/test_ols.py | 8 ++++---- tox.ini | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/test_ols.py b/tests/test_ols.py index 02d4da60..935b2e63 100644 --- a/tests/test_ols.py +++ b/tests/test_ols.py @@ -131,7 +131,7 @@ def test_random_weights_intercept(random_data_weights): def test_regularizer(): lm = OLSRegressionSolver() - lambda_I = lm._get_regularizer(7, 10, fit_intercept=True, regularize_intercept=True, n_feat_ignore_req=2) + 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) @@ -143,7 +143,7 @@ def test_regularizer(): 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_req=2) + 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) @@ -176,7 +176,7 @@ def test_regularization_with_intercept_and_unreg_feature(random_data_no_weights) lm = OLSRegressionSolver() lambda_ = 1e6 - lm.fit(x, y, lambda_=lambda_, fit_intercept=True, regularize_intercept=False, n_feat_ignore_req=2) + 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 @@ -191,7 +191,7 @@ def test_regularization_with_intercept_but_no_intercept_reg_and_unreg_feature(ra lm = OLSRegressionSolver() lambda_ = 1e6 - lm.fit(x, y, lambda_=lambda_, fit_intercept=True, regularize_intercept=True, n_feat_ignore_req=2) + 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 diff --git a/tox.ini b/tox.ini index ed82e997..69c14fac 100644 --- a/tox.ini +++ b/tox.ini @@ -1,5 +1,5 @@ [tox] -envlist=py{3} +envlist=py3.10,py3.11 skipdist=True [base] From db8cbee9de88b5103f2f764537c43a43c4587ffb Mon Sep 17 00:00:00 2001 From: lbvienna Date: Mon, 18 Sep 2023 22:43:04 -0400 Subject: [PATCH 13/20] updated python version in GH also --- .github/workflows/test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 78f361ba..57d5ed18 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, 3.11] steps: - uses: actions/checkout@v2 - name: Setup Python From bb52d6b3484eea3a9f8fc45f58ba3c9173066c3d Mon Sep 17 00:00:00 2001 From: lbvienna Date: Mon, 18 Sep 2023 22:45:23 -0400 Subject: [PATCH 14/20] added quotes to 3.1 --- .github/workflows/test.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 57d5ed18..64d19e3b 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -7,7 +7,7 @@ jobs: timeout-minutes: 5 strategy: matrix: - python-version: [3.10, 3.11] + python-version: ['3.10', 3.11] steps: - uses: actions/checkout@v2 - name: Setup Python From c149b04279328ed1f805bdaeea906a28f3b346be Mon Sep 17 00:00:00 2001 From: lbvienna Date: Tue, 19 Sep 2023 09:53:03 -0400 Subject: [PATCH 15/20] updated comments --- src/elexsolver/QuantileRegressionSolver.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/elexsolver/QuantileRegressionSolver.py b/src/elexsolver/QuantileRegressionSolver.py index 1ab9a64b..aab610ff 100644 --- a/src/elexsolver/QuantileRegressionSolver.py +++ b/src/elexsolver/QuantileRegressionSolver.py @@ -41,6 +41,8 @@ def _get_regularizer( Get regularization component of the loss function. Note that this is L2 (ridge) regularization. """ # 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 From 5ca87f767b396d6c121bf12e6d38643eb23b6c2e Mon Sep 17 00:00:00 2001 From: lbvienna Date: Tue, 19 Sep 2023 11:06:48 -0400 Subject: [PATCH 16/20] only using python 3.10 for now --- .github/workflows/pre-commit.yml | 2 +- .github/workflows/test.yml | 2 +- tox.ini | 3 ++- 3 files changed, 4 insertions(+), 3 deletions(-) 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 64d19e3b..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.10', 3.11] + python-version: ['3.10'] steps: - uses: actions/checkout@v2 - name: Setup Python diff --git a/tox.ini b/tox.ini index 69c14fac..c19a5e08 100644 --- a/tox.ini +++ b/tox.ini @@ -1,5 +1,5 @@ [tox] -envlist=py3.10,py3.11 +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 From 02bff49e9bd29f6bc6df422babe28b0996715ca5 Mon Sep 17 00:00:00 2001 From: lbvienna Date: Tue, 19 Sep 2023 14:15:16 -0400 Subject: [PATCH 17/20] removed accidental return object --- src/elexsolver/QuantileRegressionSolver.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/elexsolver/QuantileRegressionSolver.py b/src/elexsolver/QuantileRegressionSolver.py index aab610ff..69d77cdb 100644 --- a/src/elexsolver/QuantileRegressionSolver.py +++ b/src/elexsolver/QuantileRegressionSolver.py @@ -83,7 +83,7 @@ def fit( fit_intercept: bool = True, regularize_intercept: bool = False, n_feat_ignore_reg: int = 0, - ) -> np.ndarray: + ): """ Fits quantile regression """ From a49644b7748793fbb4076a778a1aa5d164dc80a1 Mon Sep 17 00:00:00 2001 From: lbvienna Date: Tue, 19 Sep 2023 14:29:05 -0400 Subject: [PATCH 18/20] readded normalizing weights --- src/elexsolver/QuantileRegressionSolver.py | 9 +++-- tests/test_quantile.py | 39 ++++++++++++++++++++++ 2 files changed, 46 insertions(+), 2 deletions(-) diff --git a/src/elexsolver/QuantileRegressionSolver.py b/src/elexsolver/QuantileRegressionSolver.py index 69d77cdb..cfe4c525 100644 --- a/src/elexsolver/QuantileRegressionSolver.py +++ b/src/elexsolver/QuantileRegressionSolver.py @@ -83,6 +83,7 @@ def fit( fit_intercept: bool = True, regularize_intercept: bool = False, n_feat_ignore_reg: int = 0, + normalize_weights: bool = True, ): """ Fits quantile regression @@ -97,8 +98,12 @@ def fit( if weights is None: weights = np.ones((y.shape[0],)) - # normalize weights - weights = weights / np.sum(weights) + # normalize weights. default to true + if normalize_weights: + weights_sum = np.sum(weights) + if weights_sum == 0: + raise ZeroDivisionError + weights = weights / weights_sum # _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): diff --git a/tests/test_quantile.py b/tests/test_quantile.py index a172d38c..792c6bd5 100644 --- a/tests/test_quantile.py +++ b/tests/test_quantile.py @@ -158,6 +158,45 @@ def test_random_upper_weights(random_data_weights): np.testing.assert_allclose(quantreg.coefficients, [[3.47742, 2.07360, 4.51754, 4.15237, 9.58856]], rtol=TOL) +############################# +# Test weight normalization # +############################# + + +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 = np.asarray([0] * x.shape[0]) # all zero weights + + quantreg = QuantileRegressionSolver() + + # Will succeed without weight normalization + quantreg.fit(x, y, tau, normalize_weights=False, weights=weights, fit_intercept=False) + + # Will fail with weight normalization + with pytest.raises(ZeroDivisionError): + quantreg.fit(x, y, tau, normalize_weights=True, weights=weights, fit_intercept=False) + + +def test_weight_normalization_same_fit(random_data_weights): + 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) + 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) + np.testing.assert_allclose(preds, [[9, 9, 9, 15]], rtol=TOL) + + ######################## # Test regularization # ######################## From 086688ff88b53a0a88242f89705323b50f7033b6 Mon Sep 17 00:00:00 2001 From: lbvienna Date: Wed, 20 Sep 2023 19:23:17 -0400 Subject: [PATCH 19/20] updated python version in classifier --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 07365d2c..8a2236d7 100644 --- a/setup.py +++ b/setup.py @@ -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, From 37d995be0d1f4fd26c99fdc56f00b828b856d6e6 Mon Sep 17 00:00:00 2001 From: lbvienna Date: Fri, 22 Sep 2023 18:09:04 -0400 Subject: [PATCH 20/20] release 2.0.0 --- CHANGELOG.md | 2 ++ setup.py | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) 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/setup.py b/setup.py index 8a2236d7..681c4df9 100644 --- a/setup.py +++ b/setup.py @@ -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])