Why doesn't penalized cubic regression reduce the number of knots in a GAM? As far as I understand, cubic regression penalization prevents overfitting by reducing the number of knots by penalizing wiggliness. The supplied parameter k serves only as a starting points for choosing the knots.
From ?mgcv::choose.k:

In practice k-1 (or k) sets the upper limit on the degrees of freedom associated with an s smooth [...] However the actual effective degrees of freedom are controlled by the degree of penalization selected during fitting, by GCV, AIC, REML or whatever is specified. 

However, with my model specifications, penalization does not seem to have an effect. The demo below starts and ends with 114 knots. Obviously I'm mistaken somewhere. What am I missing?
> require(mgcv)
> data(beavers)
> mod <- gam(temp ~ s(time, bs="cr", k = length(beaver1$time)), data = beaver1)
> length(mod$smooth[[1]]$xp)
[1] 114
> plot(mod, scheme = 1, residuals = T)


 A: You do have a slight misunderstanding of how penalization works in the spline regression situation.  Regression with (cubic, in this case) splines does not implement penalization by reducing the number of knots, but instead by reducing the magnitude of the coefficients of the spline terms.  "Degrees of freedom" doesn't refer to the number of knots, but rather to "actual effective degrees of freedom", which take into account the fact that the regression wasn't allowed to fully optimize on the coefficients (unlike, say, least squares), therefore the effective loss in degrees of freedom is somewhat less than the number of coefficients estimated.  
One way of writing out the regression problem is:
$$\min_{\beta} \sum(y_i- Z_i^T\gamma - X_i^T\beta)^2 + \lambda||\beta||_2$$
The $X_i$ are the cubic spline terms, the $Z_i$ are terms whose coefficients aren't being regularized.  The last term is a penalty on the magnitude of the $\beta$ coefficients; the larger $\lambda$ is, the more the coefficients are shrunk towards zero.  It's this last term that penalizes "wiggliness"; to see this, imagine what happens as $\lambda \rightarrow \infty$ - the fitted curve comes closer and closer to ignoring the $X_i$ altogether, which is about as non-wiggly (with respect to the $X_i$) as you can get - it's saying that the estimated function is constant with respect to changes in the $X_i$.
You certainly could regularize by reducing the number of knots, but, from an optimization perspective, this is very difficult.  If you are optimizing with respect to knots, not only do you have to choose the number of knots but also their location.  As far as I know, there is no good, fast solution to the latter problem.  With regularization of the coefficients, there are high-quality approximations and shortcuts that are very effective at speeding up the computations.  
