Generalize ReHLine to solve quadratic problems with linear constraints

[aorazalin]

ReHLine is designed to solve problems of the kind
$$\min_{\mathbf{\beta} \in \mathbb{R}^d} \frac{1}{2}||\beta||^2 + \sum_{i=1}^n L_i(\mathbf{x}_i^T \mathbf{\beta}) (*)$$
s.t. $\mathbf{Aw} + \mathbf{b} \geq \mathbf{0}$ where $L_i$ are PLQ convex loss functions.

This PR aims to add a straightforward generalization (for details, see the attachment) to the ReHLine algorithm to solve problems the kind
$$\min_{\mathbf{\beta} \in \mathbb{R}^d} \frac{1}{2}\mathbf{\beta}^T \mathbf{P} \mathbf(\beta) - \mathbf{\mu}^T \mathbf{\beta} + \sum_{i=1}^n L_i(\mathbf{x}_i^T \mathbf{\beta})$$
s.t. $\mathbf{Aw} + \mathbf{b} \geq \mathbf{0}$ where $L_i$ are PLQ convex loss functions, and $\mathbf{P}$ is a PSD matrix.

It is of course possible to transform the above problem to the original problem (*) using linear transformations. However, there are cases where that might not be desirable. For example, instead of running cholesky decomposition one might have a very fast way to solve $\mathbf{Px} = \mathbf{y}$ (which is intrinsic to how the modified ReHLine algorithm runs), e.g. portfolio optmization with factor modelling.

PR_details.pdf

[statmlben]

As its performance appears to be slightly suboptimal compared to variable substitution approach, I will put on hold this push for further investigation and improvement.