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:
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?

Labels: New Feature, Needs Decision - Include Feature