I'm looking for a good package to train a linear quantile regression model, i.e. $\hat y = \sum_{i=1}^n w_i \cdot X_i$. With $x_i$ are the input features, and $w_i$ are the bounded trainable weights. This should all be trained with quantile loss. Here's the catch: I want to optimize under the constraint that all weights $w_i$ are positive (or strictly positive, doesn't matter much). I just cannot seem to find a good package that does it. Preferably, I'd like to fit 50'ish quantiles (0.01, 0.03, ..., 0.99) simultaneously (both computationally and memory efficient) (I have access to multiple cpu's). Scikit-learn's quantile regressor doesn't allow bounded parameters. Statsmodels doesn't allow bounded parameters either as far as I understand? Although I do see something interesting about ``Estimate a quantile regression model using iterative reweighted least squares.'' What is this about? Maybe I can use bounded least squares in combination with their iterative algorithm? I can't seem to find that approach in their sources, does anyone know about this? I tried my own implementation using linear programming, which used way too much memory. def bounded_quantile_regression(X, y, quantile): n_samples, n_features = X.shape # Variables: [beta_1, ..., beta_p, u_1, ..., u_n, v_1, ..., v_n] # u_i, v_i are positive parts for quantile regression loss c = np.hstack([ np.zeros(n_features), # Coefficients (no cost in objective) quantile * np.ones(n_samples), # u_i (1-quantile) * np.ones(n_samples) # v_i ]) # For each sample: y_i - X_i beta = u_i - v_i => X_i*beta + v_i - u_i = y_i A_eq = np.hstack([X, -np.eye(n_samples), np.eye(n_samples)]) b_eq = y coef_bounds = [(0, None)] * n_features u_bounds = [(0, np.inf)] * n_samples v_bounds = [(0, np.inf)] * n_samples bounds = coef_bounds + u_bounds + v_bounds result = linprog( c, A_eq=A_eq, b_eq=b_eq, bounds=bounds, method='highs' ) if not result.success: raise ValueError("Optimization did not converge.") beta = result.x[:n_features] return beta Scipy.optimize was really slow, even when passing a gradient function. In the end I settled on a torch implementation which uses softplus activated variables with a custom QR loss function (autograd wasn't fast enough). But this feels a bit like unnecessary complexity, but in implementation as in package-usage (I don't need torch anywhere else in my codebase). Does anyone know any better alternatives? Because I've kinda ran out of ideas, but it feels like it shouldn't be such a hard problem to tackle. Continue reading...
