Why are there large coefficents for higher-order polynomial In Bishop's book on machine learning, it discusses the problem of curve-fitting a polynomial function to a set of data points.
Let M be the order of the polynomial fitted. It states as that 

We see that, as M increases, the magnitude of the coefficients
  typically gets larger. In particular for the M = 9 polynomial, the
  coefficients have become finely tuned to the data by developing large
  positive and negative values so that the corresponding polynomial
  function matches each of the data points exactly, but between data
  points (particularly near the ends of the range) the function exhibits
  the large oscillations.

I don't understand why large values implies more closely fitting the data points. I would think the values would become more precise after the decimal instead for better fitting.
 A: Abhishek,
you are right that improving precision of coefficients will improve accuracy.

We see that, as M increases, the magnitude of the coefficients typically gets larger. In particular for the M = 9 polynomial, the coefficients have become finely tuned to the data by developing large positive and negative values so that the corresponding polynomial function matches each of the data points exactly, but between data points (particularly near the ends of the range) the function exhibits large oscillations.

I think the magnitude issue is rather irrelevant to Bishop's overall point - that using a complicated model on limited data leads to 'overfitting'.
In his example 10 datapoints are used to estimate a 9 dimensional polynomial (ie 10 variables and 10 unknowns).
If we fit a sine wave (no noise), then the fit works perfectly, since sine waves [over a fixed interval] can be approximated with arbitrary accuracy using polynomials.  However, in Bishop's example we have a certain amount of 'noise' that we should not fit.  The way we do this is by keeping the number of datapoints to number of model variables (polynomial coefficents) large or by using regularisation.

Regularisation imposes 'soft' constraints on the model (eg in ridge regression) the cost function you try to minimise is a combination of 'fitting error' and model complexity : eg in ridge regression the complexity is measured by the sum of squared coefficients- in effect this imposes a cost on reducing error - increasing the coefficients will only be allowed if it has a large enough reduction in the fitting error [how large is large enough is specified by a multiplier on the model complexity term].  Therefore the hope is that by choosing the appropriate multiplier we will not fit to additional small noise term, since the improvement in fit does not justify the increase in coefficients.
You asked why large coefficients improve the quality of the fit.  Essentially, the reason is that the function estimated (sin + noise) is not a polynomial, and the large changes in curvature required to approximate the noise effect with polynomials require large coefficients.
Note that using orthogonal polynomials has no effect ( I have added an offset of 0.1 just so that the orthogonal and raw polynomials are not on top of each other)
require (penalized)
poly_order<-9
x_long <- seq(0, 1, length.out = 100)
nx <- 10
x <- seq(0, 1, length.out = nx)
noise <- rnorm(nx, 0, 1)
noise_scale <- 0.2
y <- sin(2*pi*x) + noise_scale*noise

training_data <- data.frame(x=x, y=y)
y_long <- sin(2*pi*x_long)

plot(x, y, col ='blue', ylim=c(-1.5,1.5))
lines(x_long, y_long, col='green')

polyfit_raw <- lm(y ~ poly(x, poly_order, raw=TRUE), 
                  data=training_data)
summary(polyfit_raw)

polyfit_raw_ridge1 <- penalized(y, ~poly(x, poly_order, raw=TRUE), 
           model='linear', data=training_data, lambda2=0.0001, 
           maxiter=10000, standardize=TRUE)

polyfit_orthog <- lm(y ~ poly(x, poly_order), data=training_data)
summary(polyfit_orthog)

pred_raw <- predict(polyfit_raw, data.frame(x=x_long))
pred_ortho <- predict(polyfit_orthog, data.frame(x=x_long))
pred_raw_ridge <- predict(polyfit_raw_ridge1, 
                   data=data.frame(x=x_long))[,'mu']
lines(x_long, pred_raw, col='red')
# add 0.1 offset to make visible
lines(x_long, pred_ortho+0.1, col='black')
lines(x_long, pred_raw_ridge, col='purple')
legend("bottomleft", legend=c('data sin(2 pi x) + 
             noise','sin(2 pi x)', 'raw poly', 'orthog poly + 0.1 
             offset', 'raw poly + ridge regression'),
       fill=c('blue', 'green', 'red', 'black', 'purple'))

A: This is a well known issue with high-order polynomials, known as Runge's phenomenon. Numerically it is associated with ill-conditioning of the Vandermonde matrix, which makes the coefficients very sensitive to small variations in the data and/or roundoff in the computations (i.e. the model is not stably identifiable). See also this answer on the SciComp SE.
There are many solutions to this problem, for example Chebyshev approximation, smoothing splines, and Tikhonov regularization.
Tikhonov regularization is a generalization of ridge regression, penalizing a norm $||\Lambda \theta]||$ of the coefficient vector $\theta$, where for smoothing the weight matrix $\Lambda$ is some derivative operator. To penalize oscillations, you might use $\Lambda \theta=p^{\prime\prime}[x]$, where $p[x]$ is the polynomial evaluated at the data.
EDIT: The answer by user hxd1011 notes that some of the numerical ill-conditioning problems can be addressed using orthogonal polynomials, which is a good point. I would note however that the identifiability issues with high-order polynomials still remain. That is, numerical ill-conditioning is associated with sensitivity to "infinitesimal" perturbations (e.g. roundoff), while "statistical" ill-conditioning concerns sensitivity to "finite" perturbations (e.g. outliers; the inverse problem is ill-posed).
The methods mentioned in my second paragraph are concerned with this outlier sensitivity. You can think of this sensitivity as violation of the standard linear regression model, which by using an $L_2$ misfit implicitly assumes the data is Gaussian. Splines and Tikhonov regularization deal with this outlier sensitivity by imposing a smoothness prior on the fit. Chebyshev approximation deals with this by using an $L_{\infty}$ misfit applied over the continuous domain, i.e. not just at the data points. Though Chebyshev polynomials are orthogonal (w.r.t. a certain weighted inner product), I believe that if used with an $L_2$ misfit over the data they would still have outlier sensitivity.
