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:
Thanks for checking it out!

@glemaitre wrote:
> I assume that we would not be against the inclusion of a pure NumPy implementation. We would only need to have a proper benchmark against the current OMP solver.

Yep, we can expand on what's in the other repository and collect that information in a suitable place.

> I would probably use a `solver` parameter in the current `OrthogonalMatchingPursuit` instead of creating a new class.

Makes sense as an approach to keep the API minimal. A couple of things do speak for adding a new class IMO:

* `OrthogonalMatchingPursuit` is a well-established term, and even if ORMP is quite similar to OMP, it's not OMP; rather, both are variants on the same theme, with none of them being instances of the other, so if a user comes looking for ORMP specifically, chances are OMP wouldn't be the first place they would look.
* As implemented right now, ORMP comes with one additional parameter, a [`greediness`](https://github.com/Kvantify/ormp/blob/d2fc645/src/ormp/ormp.py#L42C1-L44C62) which can be used to weigh execution time against expected accuracy; exposing such a parameter when it only works for one method always becomes a bit awkward.
* But the two have similarities; much of the pre-processing can be shared, either by subclassing or if you would rather avoid deep(er) class hierarchies, by just collecting what they have in common in appropriate functions.

> We would also need a set of tests.

Right, we have a handful in the other repo, and can adapt to the existing OMP test machinery as well.

@adrinjalali wrote:
> We could also probably use array namespace here to enable using JAX? But I agree adding it as a different solver is probably an easier path forward.

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?

Labels: Enhancement