Implementation of Order Recursive Matching Pursuit

### Describe the workflow you want to enable

The sparse linear regression problem is the NP-hard problem of performing ordinary L² linear regression with the catch that at most a fixed number of the resulting coefficients can be non-zero (or the variation in which one wants as few non-zero columns as possible while reaching a given score). Scikit-learn includes a greedy approximation algorithm for solving this problem in the form of [`sklearn.linear_model.OrthogonalMatchingPursuit`](https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.OrthogonalMatchingPursuit.html#sklearn.linear_model.OrthogonalMatchingPursuit) (OMP).

Many variations on this algorithm exist; one of them (which dates back [at least to the 70s](https://www.jstor.org/stable/1267353)) is "Order Recursive Matching Pursuit" (ORMP). Like OMP, ORMP is an iterative algorithm which picks one column from the input at a time, while maintaining the current residual. In OMP, one picks the column whose angle with the current residual is minimal, whereas in ORMP, one chooses the column resulting in the minimal residual. That is, while the approach is still greedy, it's somewhat more to the point, choosing the column leading to the minimal objective in each step. A priori, the cost of this approach is increased computational time, but with proper care, it is possible to make a [simple implementation](https://github.com/Kvantify/ormp/blob/master/src/ormp/impl_np.py) that rivals the runtime performance of the existing OMP implementation, and is faster on the instances that motivated the work.

### Describe your proposed solution

In our company, we needed fit accuracy greater than what could be achieved with OMP and thus implemented ORMP in a handful of different ways. The implementations are available on https://github.com/Kvantify/ormp. As I believe it fits the inclusion criteria, these implementations could be lifted into scikit-learn more or less directly, some parts of it more relevant than others: I imagine we don't want a JAX dependency in scikit-learn, for example, but note that even the pure NumPy implementation appears to be generally faster than scikit-learn's existing OMP, while delivering the promised score improvements on typical (but far from all) problem instances. See the README in that repository for [a few crude performance benchmarks](https://github.com/Kvantify/ormp#choosing-an-implementation) and [notes on a result quality comparison with OMP](https://github.com/Kvantify/ormp#comparison-with-orthogonal-matching-pursuit).

We can submit a pull request if you are interested.

COMMENT:
> To allow using both backends in a single implementation, somehow? If so, good question. If you diff the [NumPy implementation](https://github.com/Kvantify/ormp/blob/d2fc645ab3db636bc2971a22ac5f50e7550a3514/src/ormp/impl_np.py) to the [two JAX implementations](https://github.com/Kvantify/ormp/blob/d2fc645ab3db636bc2971a22ac5f50e7550a3514/src/ormp/impl_jax.py) you'll find that we did have to introduce some JAX-specific tweaks. To keep array sizes static, and to avoid unrolling the potentially long loop when short compile times are important. But maybe it's possible to get the best of both worlds in a single implementation?

maybe @thomasjpfan or @betatim would know how to tackle the issue?

Labels: Enhancement