Improve numerical stability of polynomial fitting by replacing normal equations with QR decomposition

[gbowne1]

Summary

The current implementation of PolynomialSurface::fit() uses normal equations to solve the least squares problem:

Eigen::VectorXd coeffs = (A.transpose() * A).ldlt().solve(A.transpose() * b);

While this works for well-conditioned data, it can be numerically unstable for poorly distributed or nearly collinear input points.
Suggested Fix

Replace the normal equations with QR decomposition for better numerical stability:

Eigen::VectorXd coeffs = A.colPivHouseholderQr().solve(b);

This change improves robustness without significantly impacting performance.

Impact

• More accurate results for edge cases

• Better behavior with noisy or ill-conditioned data

• This should not require major API changes