Хорошие пакеты для ограниченной линейной квантильной регрессии? ⇐ Python

[Anonymous]

Я ищу хороший пакет для обучения модели линейной квантильной регрессии, т. е. $\hat y = \sum_{i=1}^n w_i \cdot X_i$. $x_i$ — это входные признаки, а $w_i$ — ограниченные обучаемые веса. Все это должно быть обучено с потерей квантилей.
Вот в чем загвоздка: я хочу оптимизировать при условии, что все веса $w_i$ положительны (или строго положительны, не имеет большого значения). Кажется, я просто не могу найти хороший пакет, который бы это сделал.
Предпочтительно, я бы хотел одновременно разместить 50 квантилей (0,01, 0,03, ..., 0,99) (как с точки зрения вычислений, так и с точки зрения памяти) (у меня есть доступ к нескольким процессорам).

• Квантильный регрессор Scikit-learn не позволяет ограничить параметры.
• Насколько я понимаю, Statsmodels также не допускает ограниченных параметров? Хотя я вижу кое-что интересное в статье «Оценка модели квантильной регрессии с использованием итеративного перевзвешенного метода наименьших квадратов». О чем идет речь? Может быть, я могу использовать ограниченный метод наименьших квадратов в сочетании с их итеративным алгоритмом? Кажется, я не могу найти этот подход в их источниках, кто-нибудь знает об этом?
• Я попробовал свою собственную реализацию с использованием линейного программирования, которое использовало слишком много памяти.

Код: Выделить всё

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 работал очень медленно, даже при передаче функции градиента.
• В конце концов я остановился на реализации факела, которая использует активируемые переменные softplus с настраиваемой функцией потери QR (autograd был недостаточно быстрым). Но это немного похоже на ненужную сложность, как в реализации, так и в использовании пакета (мне не нужно больше нигде в моей кодовой базе).

Кто-нибудь знает лучшие альтернативы? Потому что у меня закончились идеи, но мне кажется, что решить эту проблему не составит большого труда.

Подробнее здесь: https://stackoverflow.com/questions/798 ... regression