Raw or orthogonal polynomial regression? I want to regress a variable $y$ onto $x,x^2,\ldots,x^5$. Should I do this using raw or orthogonal polynomials? I looked at the question on the site that deals with these, but I don't really understand what's the difference between using them.
Why can't I just do a "normal" regression to get the coefficients $\beta_i$ of $y=\sum_{i=0}^5 \beta_i x^i$ (along with p-values and all the other nice stuff) and instead have to worry whether using raw or orthogonal polynomials? This choice seems to me to be outside the scope of what I want to do.
In the stat book I'm currently reading (ISLR by Tibshirani et al) these things weren't mentioned. Actually, they were downplayed in a way.
The reason is, AFAIK, that in the lm() function in R, using y ~ poly(x, 2) amounts to using orthogonal polynomials and using y ~ x + I(x^2) amounts to using raw ones. But on pp. 116 the authors say that we use the first option because the latter is "cumbersome" which leaves no indication that these commands actually do two completely different things (and have different outputs as a consequence).
(third question) Why would the authors of ISLR confuse their readers like that?
 A: 
Why can't I just do a "normal" regression to get the coefficients?

Because it is not numerically stable. Remember computer is using fixed number of bits to represent a float number. Check IEEE754 for details, you may surprised that even simple number $0.4$, computer need to store it as $0.4000000059604644775390625$. You can try other numbers here
Using raw polynomial will cause problem because we will have huge number. Here is a small proof: we are comparing matrix condition number with raw and orthogonal polynomial.
> kappa(model.matrix(mpg~poly(wt,10),mtcars))
[1] 5.575962
> kappa(model.matrix(mpg~poly(wt,10, raw = T),mtcars))
[1] 2.119183e+13

You can also check my answer here for an example.
Why are there large coefficents for higher-order polynomial
A: 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)
cor(raw.poly)
cor(orthogonal.poly)

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. This is typically referred to as the problem of multicollinearity. 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 inference routines. This will manifest as larger standard errors (and thus smaller t-stats) 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))
summary(raw.mod)
summary(orthogonal.mod)

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. 
A: I would have just commented to mention this, but I do not have enough rep, so I'll try to expand into an answer. You might be interested to see that in Lab Section 7.8.1 in "Introduction to Statistical Learning" (James et. al., 2017, corrected 8th printing), they do discuss some differences between using orthogonal polynomials or not, which is using the raw=TRUE or raw=FALSE in the poly() function. For example, the coefficient estimates will change, but the fitted values do not:
# using the Wage dataset in the ISLR library
fit1 <- lm(wage ~ poly(age, 4, raw=FALSE), data=Wage)
fit2 <- lm(wage ~ poly(age, 4, raw=TRUE), data=Wage)
print(coef(fit1)) # coefficient estimates differ
print(coef(fit2))
all.equal(predict(fit1), predict(fit2)) #returns TRUE    

The book also discusses how when orthogonal polynomials are used, the p-values obtained using the anova() nested F-test (to explore what degree polynomial might be warranted) are the same as those obtained when using the standard t-test, output by summary(fit). This illustrates that the F-statistic is equal to the square of the t-statistic in certain situations. 
