I'm not 100% certain here, but I believe the answer is less about numeric stability (though that plays a role) and more about reducing correlation.
In essence -- the issue boils down to the fact that when we regress against a bunch of high order polynomials, the covariates we are regressing against become highly correlated. Example code below:
x = rnorm(1000)
raw.poly = poly(x,6,raw=T)
orthogonal.poly = poly(x,6)
This is tremendously important. As the covariates become more correlated, our ability to determine which are important (and what the size of their effects are) erodes rapidly. At the limit, if we had two variables that were fully correlated, when we regress them against something, its impossible to distinguish between the two -- you can think of this as an extreme version of the problem, but this problem affects our estimates for lesser degrees of correlation as well. Thus in a real sense -- even if numerical instability wasn't a problem -- the correlation from higher order polynomials does tremendous damage to our estimation routines. This will manifest as larger standard errors that you would otherwise see (see example regression below). For this reason, we might choose to orthogonalize our polynomials before regressing them.
y = x*2 + 5*x**3 - 3*x**2 + rnorm(1000)
raw.mod = lm(y~poly(x,6,raw=T))
orthogonal.mod = lm(y~poly(x,6))
If you run this code, interpretation is a touch hard because the coefficients all change and so things are hard to compare. Looking at the T-stats though, we can see that the ability to determine the coefficients was MUCH larger with orthogonal polynomials. For the 3 relevant coefficients, I got t-stats of (560,21,449) for the orthogonal model, and only (28,-38,121) for the raw polynomial model. This is a huge difference for a simple model with only a few relatively low order polynomial terms that mattered.
That is not to say that this comes without costs. There are two primary costs to bear in mind. 1) we lose some interpretability with orthogonal polynomials. We might understand what the coefficient on
x**3 means, but interpreting the coefficient on
x**3-3x (the third hermite poly -- not necessarily what you will use) can be much harder.
Second -- when we say that these are polynomials are orthogonal -- we mean that they are orthogonal with respect to some measure of distance. Picking a measure of distance that is relevant to your situation can be difficult. However, having said that, I believe that the
poly function is designed to choose such that it is orthogonal with respect to covariance -- which is useful for linear regressions.