Support orthogonal polynomial features (via QR decomposition) in `PolynomialFeatures`

### Describe the workflow you want to enable

I want to introduce support for orthogonal polynomial features via QR decomposition in `PolynomialFeatures`, closely mirroring the behavior of R's `poly()` function.

In regression modeling, using orthogonal polynomials can often lead to improved numerical stability and reduced multi-collinearity among polynomial terms

As an example of what the difference looks like in R,
<pre>
#fits raw polynomial data without an orthogonal basis
model_raw <- lm(y ~ I(x) + I(x^2) + I(x^3), data = data)
#model_raw <- lm(y ~poly(x,3,raw=TRUE), data = data)

#fits the same degree-3 polynomial using an orthogonal basis
model_poly <- lm(y ~ poly(x, 3), data = data)
</pre>

This behavior cannot currently be replicated with `scikit-learn`'s `PolynomialFeatures`, which only produces the raw monomial terms. As a result transitioning from R to Python often leads to discrepancies in model behavior and performance.


### Describe your proposed solution

I propose extending `PolynomialFeatures` with a new parameter:
<pre>
PolynomialFeatures(..., method="raw")
</pre>
Accepted values:
- `"raw"` (default): retains existing behavior, returning standard raw terms
- `"qr"`: applies QR decomposition to each feature to generate orthogonal polynomial features.

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.

This feature would interact with other parameters as follows:

- `include_bias`: When `method="qr"`, The orthogonal polynomial basis inherently includes a transformed first column. However, this column is not a plain column of ones. Therefore, the concept of `include_bias=True` (which appends a column of ones) becomes redundant or misleading in this context. One option is to always set  `include_bias=False` if `method=qr` and always return orthogonal columns only, or raise a warning.

- `interaction_only`: This would be incompatible with `method="qr"` since the QR-based transformation does not naturally support selective inclusion of interaction terms.

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

Currently, users must implement QR decomposition manually when orthogonal polynomials are needed. This is a common pattern in statistical workflows but lacks "off the shelf" support in any major python library. This feature would eliminate the need to do this decomposition manually and would improve workflows for researchers who are used to R's statistical tools.

### Additional context

This idea stemmed from a broader effort to convert statistical modeling pipelines from R to python, where discrepencies in regression results were traced to the lack of orthogonal polynomial support in `PolynomialFeatures`.

I have drafted and tested a 1D implementation of this feature but wanted feedback on whether this idea aligns with `scikit-learn`'s scope before moving on. In particular, I'd appreciate input on

- Acceptability of feature-wise orthogonalization for multi-feature input.
- Preferred parameter naming (e.g., `method="qr"` vs. `orthogonal=True`).
- Compatibility decisions around parameters like `include_bias` and `interaction_only`.

COMMENT:
I also noticed that [`formulaic`](https://matthewwardrop.github.io/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.

Labels: New Feature, Needs Decision - Include Feature