Frisch-Newton Interior Point Solver for Quantile Regression

### Describe the workflow you want to enable

Hi @ scikit-learn devs! 

Over at [pyfixest](https://github.com/py-econometrics/pyfixest), we have implemented a Frisch-Newton Interior Point solver to fit quantile regressions. The algorithm goes back to work from Koenker. In practice, we have followed Koenker and Ng ["A Frisch-Newton Algorithm for Sparse Quantile Regression". ](https://link.springer.com/article/10.1007/s10255-005-0231-1)

The code is licensed under MIT and available [here](https://github.com/py-econometrics/pyfixest/blob/master/pyfixest/estimation/quantreg/frisch_newton_ip.py#L70). 

We (@apoorvalal) have collected some benchmarks [here](https://gist.github.com/apoorvalal/3e18eea79c6e9e8e8ee380e0fc0bab1f) - the FN solver seems to outperform the scikit default solver by an order of a magnitude.  

<img width="1362" height="534" alt="Image" src="https://github.com/user-attachments/assets/f25eb315-60d8-464f-80f8-bf1c6aedce3b" />

Would you be interested in a PR that adds the FN solver as a new estimation method to the quantile regression class? 

We've also implemented algorithms from [Chernozhukov et al ](https://arxiv.org/abs/1909.05782)that can drastically speed up estimation of the entire **quantile regression process**. 

All the best, Alex

### Describe your proposed solution

I open a PR and add a new solver "fn" to `QuantileRegressor`.

### Describe alternatives you've considered, if relevant

_No response_

### Additional context

Related to https://github.com/scikit-learn/scikit-learn/issues/20132

COMMENT:
Hi all, 

I have created a reference implementation here, where I integrate the FN solver into the `QuantileRegressor` class: https://github.com/s3alfisc/scikit-learn/blob/main/sklearn/linear_model/_quantile.py. Happy to open a PR if this will help you to inspect the code quickly (there is still a lot of cleanup to do, this is more a proof of concept). 

Below I provide point estimates using the highs solver from scipy and the frisch newtwon solver from pyfixest. 

Point estimates and pinball loss are identical in my tests, and the FN solver appears to be significantly faster. I have attached the code to replicate the numbers based on my fork below. 

<img width="2672" height="800" alt="Image" src="https://github.com/user-attachments/assets/489ad23f-ee60-4a59-a623-f052e488f6b0" />

<img width="550" height="307" alt="Image" src="https://github.com/user-attachments/assets/58c5901e-7dfb-4c8e-91a6-bc1e18cbc4c4" />

Unfortunately I cannot compute the sub-modality plots, as the highs-solver does not allow to track values of the objective function; and if the algorithm stops before convergence, no parameters are returned - so the trick that we discussed @ogrisel is unfortunately not possible. The only option I see to provide you with this would be to call the HighS solver directly (and not via scipy). 

Next steps for me would be to clean up my implementation, open a PR, and make the FN algo work with sparse matrices. 

Let me know what you think! 

Best, Alex 

## Code 

Replace the `_quantile.py` file with https://github.com/s3alfisc/scikit-learn/blob/main/sklearn/linear_model/_quantile.py. 

<details>

```python
import numpy as np
import time
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
from sklearn.linear_model import QuantileRegressor
from sklearn.metrics import mean_pinball_loss


class QuantileRegressorBenchmark(QuantileRegressor):

    def __init__(self, N, k, rng, solver, solver_options= None, n_iter = 5):
        super().__init__(solver=solver, solver_options=solver_options, alpha=0)

        self.N = N
        self.k = k
        self.n_iter = n_iter
        self.solver = solver
        # make sure options are passed correctly
        if solver_options is not None:
            solver_options["maxiter"] = solver_options.get("max_iter", None)
        self.rng = rng
        self.y, self.X = self.get_data(N = self.N, k = self.k)

        self.start_time = time.time()
        for i in range(self.n_iter):
            self.fit(X=self.X, y=self.y)
        self.end_time = time.time()
        self.get_pinball_loss()
        self.collect()

    def get_data(self, N, k):

        X = self.rng.normal(size=(N, k))
        y = 1 + 2 * X[:, 0] + 3 * X[:, 1] + self.rng.normal(0, 1, N)
        return y, X

    def get_pinball_loss(self):

        y = self.y
        yhat = self.predict(X = self.X)
        return mean_pinball_loss(y, yhat, alpha=self.quantile)

    def collect(self):

        return {
            "solver":self.solver,
            "N": self.N,
            "k": self.k,
            "coef": self.coef_[0],
            "nit": self.n_iter_,
            "pinball_loss": self.get_pinball_loss(),
            "time": (self.end_time - self.start_time) / self.n_iter,
            "estimation_n_iter": self.n_iter
        }


all_list = []

for N in [100, 1000, 10_000]:
    for k in [10, 50, 100]:

        qr_sp = QuantileRegressorBenchmark(
            solver = "highs",
            N = N,
            k = k,
            rng = np.random.default_rng(1),
            n_iter = 5
        ).collect()

        qr_fn = QuantileRegressorBenchmark(
            solver = "frisch-newton",
            N = N,
            k = k,
            rng = np.random.default_rng(1),
            n_iter = 5
        ).collect()

        all_list.append(qr_sp)
        all_list.append(qr_fn)


res = pd.DataFrame(all_list)


## create plot

df = res.sort_values(["k", "N", "solver"]).reset_index(drop=True)

# Unique k panels
ks = sorted(df["k"].unique())
fig, axes = plt.subplots(1, len(ks), figsize=(4.5 * len(ks), 3.8), sharey=True)

# If only one k, axes may not be an array
if len(ks) == 1:
    axes = [axes]

for ax, k in zip(axes, ks):
    d = df[df["k"] == k]
    # plot a line per solver
    for solver in d["solver"].unique():
        dd = d[d["solver"] == solver].sort_values("N")
        ax.plot(dd["N"], dd["time"], marker="o", label=solver)
    ax.set_title(f"k = {k}")
    ax.set_xlabel("N")
    ax.grid(True, alpha=0.3)

axes[0].set_ylabel("runtime (seconds)")
axes[-1].legend(title="solver", loc="best")
fig.suptitle("Quantile Regression Runtime vs N by Number of Covariates k", y=1.04)
fig.tight_layout()
# Optional: if the scales differ wildly, try log-scale on y:
# for ax in axes: ax.set_yscale("log")

if "agg" in matplotlib.get_backend().lower():
    fig.savefig("benchmark.png", dpi=200, bbox_inches="tight")
    print("Saved: benchmark")
else:
    plt.show()

res
```

</details>

Labels: New Feature, Needs Decision - Include Feature