Thanks for your proposal.

Because R's poly() only operates on 1D input vectors, my thought was to apply QR decomposition feature by feature when the input is multi-dimensional. Each column is processed independently, mirroring R's approach.

That sounds quite different from the current behavior of PolynomialFeatures which always considers interactions between input features. I wonder if trying to implement the two approaches into the same class makes sense. Maybe we could implement that as a new transformer but before doing so, I would rather make sure that people actually need this feature before implementing it in scikit-learn and adding the maintenance of a new estimator and cognitive overhead choosing from too many estimators and options.

Could you please publish prototype code to a gist.github.com or other small repo or notebook? It would be great if you could highlight a simple regression tasks for which this kind of feature engineering is a game changer (either in terms of predictive performance, model size, fitting speed, numerical stability of the fit) compared to what is already available in scikit-learn.

It would be great to compare to:

• SplineTransformer()
• PolynomialFeatures()
• make_pipeline(SplineTransformer(), PolynomialFeatures(interaction_only=True)

and all of the above followed by a PCA step.

BTW: maybe @lorentzenchr would like to share his views on such a new feature in scikit-learn.

I also noticed that formulaic provides an implementation of the poly function that seems to follow the R-style orthogonalization convention by default:

• https://matthewwardrop.github.io/formulaic/latest/guides/splines/#poly

Note that it is possible to wrap formulaic as a scikit-learn transformer but it requires copy-pasting some boilerplate code as explained here:

• https://matthewwardrop.github.io/formulaic/latest/guides/integration/#scikit-learn

Rather than replicating such features in scikit-learn, it might be more fruitful to see if formulaic authors be open to the idea of adding first-class integration with scikit-learn API/pipelines in directly into formulaic or via a new extension package. If so, we could extend on of scikit-learn examples on polynomial / spline features to show how to use formulaic with scikit-learn a go beyond the options implemented by default in scikit-learn.

@cottnich What is your motivation?

As a result transitioning from R to Python often leads to discrepancies in model behavior and performance.

Without penalty, the predicted values are the same for all design matrices related to orthogonal transformations.
B-Splines are almost always preferred to pure polynomials of features.
If you want to be as close to R's lm and glm, I recommend to use https://www.statsmodels.org .

There are existing orthogonal polynomial bases, such as the Legendre basis, which are already orthogonal and don't need any transformation at all. We can also support interactions using tensor-product bases.
Do we really need the heavy computational overhead of "orthogonalizing" the coefficients?

Again, what is your motivation? Why do you want orthogonal polynomials?

I don't. The OP does. But I can guess why.

You can fit high degree polynomials without the well-known "bad effects" of overfitting. Moreover and you can easily prune models by simply removing the tail of the coefficients, and staying with low-degree polynomials. This is because they are a sort of a spectrum, where higher degree polynomials are like frequency components - this pruning is just denoising.

What I don't understand is why it belongs to scikit-learn, and not their own research repository, where the OP does their research on polynomials. This library is easily extensible - just inherit the right transformer base-class, and you can put your new super-duper polynomial basis into your pipeline.

Thank you for all the responses and my apologies for not checking on this thread. I actually wasn't aware that the formulaic library handles this. I assumed it was like statsmodels.formula which doesn't do any additional computations.