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)