Intuition for nonmonotonicity of coefficient paths in ridge regression Intuitively, why may some of the slope coefficients in ridge regression increase in magnitude when the penalty parameter $\lambda$ is increased? Or in other words, why are the coefficient paths nonmonotonic? 

Above is an example of ridge coefficient paths taken from James et al. "An Introduction to Statistical Learning" (Figure 6.4 on p. 216). As you can see, Income and Rating are nonmonotonic w.r.t. $\lambda$.
Algebraic and geometric arguments will be appreciated.
 A: My reasoning is somewhat similar to Cagdas', but I'd like to look at how the things develop as $\lambda$ goes into extremes.
If I did my algebra right, the derivative of the coefficients is given by:
$$\frac{\partial \hat{\beta}}{\partial\lambda} = -(X'X + \lambda I)^{-1}\hat{\beta}$$
Now, for $\lambda \rightarrow 0$, ridge regression approaches ordinary linear regression. The penalty term becomes negligible and you can approximate the derivative by:
$$\frac{\partial \hat{\beta}}{\partial\lambda} \approx -(X'X)^{-1}\hat{\beta}$$
That is, the gradient of $\hat{\beta}$ is mostly determined by the distribution of your data. The direction of change of $\hat{\beta}$ with increasing $\lambda$ can be positive or negative for any of the coefficients, depending both on the coefficients and on your data.
On the other hand, when $\lambda \rightarrow \infty$, $\lambda I \gg X'X$ and you can approximate the derivative by:
$$\frac{\partial \hat{\beta}}{\partial\lambda} \approx -(\lambda I)^{-1}\hat{\beta}$$
Here, the gradient is determined almost entirely by the value of $\hat{\beta}$, and ridge regression tries to force it to a null vector. The path of $\hat{\beta}$ for large and increasing $\lambda$'s is almost a straight line towards the origin.
Edit:
To illustrate this behaviour, here a simple data set with correlated variables:
set.seed(0)
tb1 = tibble(
  x1 = seq(-1, 1, by=.01),
  x2 = -2*x1 + rnorm(length(x1)),
  y  = -3*x1+x2 + rnorm(length(x1))
)

The path of the coefficients $\beta$ with an increasing $\lambda$ looks like this:

A: Geometrical point of view
The ridge path is not a straight line. See for instance an image from this question The limit of "unit-variance" ridge regression estimator when $\lambda\to\infty$

Note that the path crosses the points where the circles $\Vert \beta \Vert = constant$ and ellipses $\Vert y - X\beta \Vert = constant$ are touching. It is this elliptical shape that allows the path to pass 0 and be decreasing/increasing in various parts.
An question with an image that shows this more dramatically is Why under joint least squares direction is it possible for some coefficients to decrease in LARS regression?
Plot of lasso path as coordinates of $\beta$

This plot is for Lasso instead of ridge regression, but as shown in the first image the principle is the same. Typically parameters decrease when we increase the penalty, but due to correlation it might be better that parameters decrease while simultaneously increasing another. This happens in the image with parameter $b_1$. Increasing $b_1$ makes that decreasing $b_2$ and $b_3$ coincides with less increase of the  squared error part of the loss function (the green surface).
Shrinking the parameters in a straight line from the OLS solution to 0, means that the sum of squared error becomes high. Taking a detour allows to shrink the parameters with less reduction of the sum of least squared error.
A: I think the key intuition to think of here is a situation where you are interested in modelling a single outcome of interest, $Y$, with three possible independent variables, $W$, $X$, $Z$. Imagine in a standard regression, $W$ and $X$ are found significant whereas the coefficient on $Z$ is found insignificant and close to zero, although $Z$ and $Y$ are correlated one-on-one.
Let's assume the reason that $Z$ is not significant in the full regression is because all the effect of $Z$ is subsumed in $W$ and $X$ (i.e. they jointly form a better predictor). Now in a ridge regression context, for some parameter values of $\lambda$ it might be the case that using $Z$ is preferred over having both $W$ and $X$ due to the inclusion of a parameter error term.
This example is also illustrative why ridge regression is used less within an inferential context, because even when controlling for potential confounders, they might still be 'the last coefficient standing'.
See below an example of an outcome $Y$ which is solely related to two variables, $Z$ and $X$. The variable $W$ is affected by both $Z$ and $X$ but does not affect $Y$. Therefore, the SSE (which is minimised in the regression contex) is much lower in a regression of the form $y_i = \beta_0 + \beta_x X + \beta_z Z + \beta_w W + u_i$ relative to one with only $W$: $y_i = \beta_0 + \beta_w W + u_i$. 
This is illustrated in the R-squared of the former (at 0.97) relative to the latter (at 0.35). However, the L1-norm is strictly higher in the former case than in the latter (namely $L_1 = \sqrt{\beta_0^2 + \beta_1^2 + \beta_2^2}$). As you drive up $\lambda$, the model will prioritise minimising the $L_1$-norm over the SSE. At some point, the model will prefer having just a single non-zero $\beta$ driving the $L_1$-norm, even at the cost of a much poorer model fit. Therefore, you can see that the coefficients on $\beta_x$ and $\beta_z$ start to decrease up to the point that only $\beta_w$ is non-zero. Also note that further increasing $\lambda$ makes the $L_1$-norm important enough to completely lose interest in minimising any SSE (and therefore generating any model fit) and drives all coefficients to zero.

library(glmnet)

x <- rnorm(100, 10, 2)
z <- rnorm(100, 5, 2)
w <- 0.8*x + 0.8*z + rnorm(100, 0, 1)
u <- rnorm(100, 0, 1)

y <- 2*x + 2*z + u

# R-squared of 0.97
summary(lm(y ~ 1 + x + z + w))

# R-squared of 0.35
summary(lm(y ~ 1 + z))

fit <- glmnet(as.matrix(cbind(x, z, w)), as.matrix(y))
plot(fit)

A: Ridge solution:
$$\hat{\beta_\lambda} = (X'X + \lambda I)^{-1}X'y$$
If my matrix algebra is right, derivative of Ridge with respect to $\lambda$:
$$\frac{\partial \hat{\beta_\lambda}}{\partial\lambda} = -(X'X + \lambda I)^{-2}X'y$$
which is:
$$\frac{\partial \hat{\beta_\lambda}}{\partial\lambda} = -(X'X + \lambda I)^{-1}\hat{\beta_\lambda} = -A_\lambda \hat{\beta_\lambda}$$
The derivative of each component of $\hat{\beta_\lambda}$ at each $\lambda$ depends on the values of other components of $\hat{\beta_\lambda}$. At this point it would be unreasonable to think that the derivative of a particular component of $\hat{\beta_\lambda}$ will not have a zero crossing when other components are changing.
Suppose that $X$s are not co-linear and are whitened before regression. In that case $\frac{X'X}{n}$ is $I$ where $n$ is the number of data points. Hence:
$$\frac{\partial \hat{\beta_\lambda}}{\partial\lambda} = -(n + \lambda)^{-1}\hat{\beta_\lambda} = -\frac{\hat{\beta_\lambda}}{n+\lambda} = 0$$
which will be $0$ when a component of $\hat{\beta_\lambda}$ becomes $0$. In this case we should have a monotonic scenario. Basically if we start with no correlation between regressors they will go to $0$ with increasing $\lambda$ without changing direction starting from their cross-covariance with $y$ as their initial values.
