# Gradient descent to optimize regularization parameter $\lambda$ instead of doing grid search?

Let's take a simple regression $y = ax$ with a $L^2$ regularization parameter $\lambda$. In Andrew Ng's coursera ML course, he uses a grid search to find $\lambda$.

So for example, he will find for lambdas $\lambda = 0.01,0.02,0.04,\dots,1$ the optimal $a$ $$\arg\min_a \sum_{(x,y) \in train} (ax-y)^2 + \lambda a^2$$ which gives us a function $opt_a(\lambda)$. Then he will return the $\lambda$ for which the error on the validation set is the lowest: $$\arg\min_\lambda \sum_{(x,y) \in valid} (opt_a(\lambda)x-y)^2$$ and then return $opt_a(\lambda)$ of this optimal $\lambda$.

But couldn't we use gradient descent (the same way we do to find $a$) to find $\lambda$ instead of doing a (seemingly primitive) grid search? Is there a reason I'm not seeing? Perhaps computationally too expensive? Or not enough data in the validation set so that there's a way to overtrain the hyperparameter $\lambda$, or something else?

• For gradient descent to "work" we need a convex function. Think about what you know about the function characteristics of "a".
– knk
Sep 2, 2017 at 0:43
• what's the gradient of $opt_a(\lambda)$? Sep 2, 2017 at 3:01
• Samlaf, can you explain more clearly the process you have have in mind. You have to calculate the lambda gradient numerically i believe, but you could certainly partly use gradient descent Sep 2, 2017 at 9:02
• I don't have a process in mind, but rather an intuition that seems wrong at this point. I just thought it was weird that we optimize the parameter $a$, but only heuristically choose the hyperparameter $\lambda$ (by testing a few points and taking the best one). Sep 2, 2017 at 13:59
• @Ben I just assumed that it was going to be the same as the gradient of any least squares with regularization, but I forgot that I am taking $\arg \min$. I guess this is where it breaks down, and there is no continuous analogue of softmax right? Sep 2, 2017 at 14:00

Short answer: Yes, gradient descent could conceivably be used in this situation because of the special form of ridge regression. Below, I derive an expression for the derivative of the validation set error w.r.t. $\lambda$. But, other optimization methods may be more efficient and more broadly applicable.

Longer answer: The validation set error as a function of $\lambda$ isn't convex but, empirically, it does tend to be well-behaved--e.g. smooth with a single minimum (perhaps pseudoconvex? If anyone could elaborate on its properties I'd be interested to hear). This means there's structure that can be exploited to optimize $\lambda$ more efficiently than grid search.

Gradient descent could conceivably be used to exploit some of this structure, although I've never seen anybody use it in this context. It would require differentiating the validation set error with respect to $\lambda$, which in turn requires differentiating the optimal weights (on the training set) with respect to $\lambda$. Because of the convenient form of ridge regression, it turns out that a closed form expression for the derivative exists (see below). A further tweak is needed: $\lambda$ must be constrained to be nonnegative. This could be done trivially using projected gradient descent (in this case, simply thresholding $\lambda$ after each step).

But, we can probably do better than gradient descent. For example, methods incorporating quadratic interpolation can exploit curvature information, which gradient descent can't. We could use a root finding algorithm like Brent's method to solve for $\lambda$ where the derivative is zero. Because this is a 1d optimization problem, we could also use methods like golden section search (or again, Brent's method) to directly optimize $\lambda$. These methods can exploit local structure in the objective function, but don't require an expression for the gradient. This means we can also apply them in cases where one doesn't exist. For example, there's almost certainly no closed-form expression for the derivative of the validation set error w.r.t. the lasso regularization parameter.

All of these methods only find local minima so, unlike grid search, they may not be good choices when multiple local minima exist. But, there are global search techniques that can still beat grid search in this situation by exploiting structure in the objective function.

Derivative of validation set error w.r.t. $\lambda$:

Say the training set contains inputs stored in matrix $X$ (with points on the rows, features on the columns) and outputs stored in column vector $y$. The validation set contains inputs $\tilde{X}$ and outputs $\tilde{y}$. Assume that the columns of $X$ and $\tilde{X}$ have been standardized (using the mean and standard deviation of $X$), and that $y$ and $\tilde{y}$ have been centered (using the mean of $y$). This is done to avoid including a constant term in the model.

Given $\lambda$, the weights are learned on the training set:

$$w_\lambda = \underset{w}{\arg \min} \| y - X w \|^2 + \lambda \|w\|^2$$

Ridge regression has a well known closed form solution:

$$w_\lambda = (X^T X + \lambda I)^{-1} X^T y$$

We want to find $\lambda$ that minimizes the error on the validation set: $$\min_\lambda \| \tilde{y} - \tilde{X} w_\lambda \|^2$$

Let $L$ denote the validation set error. Its derivative is:

$$\frac{d}{d\lambda} L = \frac{d}{d\lambda} \| \tilde{y} - \tilde{X} w_\lambda \|^2$$

After some matrix wrangling, we have:

$$\frac{d}{d\lambda} L = -2(\tilde{y} - \tilde{X} w_\lambda)^T \tilde{X} \frac{d}{d\lambda} w_\lambda$$

Now we must differentiate $w_\lambda$:

$$\frac{d}{d\lambda} w_\lambda = \frac{d}{d\lambda} (X^T X + \lambda I)^{-1} X^T y$$

If $A$ is an invertible square matrix that's a function of $t$ then $\frac{d}{dt} A^{-1} = -A^{-1} (\frac{d}{dt} A) A^{-1}$ (proof here). Using this identity and more matrix wrangling, we have:

$$\frac{d}{d\lambda} w_\lambda = -(X^T X + \lambda I)^{-2} X^T y$$

Plug this back into the derivative of the validation set error, along with the expression for $w_\lambda$. Let $M = X^T X + \lambda I$. Then:

$$\frac{d}{d\lambda} L = 2 (\tilde{y} - \tilde{X} M^{-1} X^T y)^T \tilde{X} M^{-2} X^T y$$

Naively computing this expression for different values of $\lambda$ would be computationally inefficient because it would require many matrix inversions. Matrix inversion may also be numerically unstable for small values of $\lambda$. To avoid this, re-express using the eigenvalue decomposition of $X^T X$. Let matrix $V$ contain the eigenvectors on the columns and $\gamma$ be a vector containing the corresponding eigenvalues. So $X^T X = V \text{diag}(\gamma) V^T$. The eigenvectors of $M = X^T X + \lambda I$ are also given by $V$, and its eigenvalues are $\gamma + \lambda$. So $M = V \text{diag}(\gamma + \lambda) V^T$. The inverse of a matrix can be found by inverting the eigenvalues so, letting $D = \text{diag}((\gamma + \lambda)^{-1})$, $M^{-1} = V D V^T$. Because $V$ is orthonormal, $M^{-2}$ simplifies to $V D^2 V^T$. Plug these expressions in:

$$\frac{d}{d\lambda} L = 2 (\tilde{y} - \tilde{X} V D V^T X^T y)^T \tilde{X} V D^2 V^T X^T y$$

Computational efficiency can be increased by pre-computing parts of the expression that don't depend on $\lambda$. Let:

$$A = 2 V^T X^T y$$ $$B = V^T \tilde{X}^T \tilde{X} V$$ $$q=\tilde{y}^T \tilde{X} V$$ $$r=y^T X V$$

Then:

$$\frac{d}{d\lambda} L = (q - r D B) D^2 A$$

Further computational gains can be had by replacing matrix operations involving $D$ with element-wise operations involving its diagonal.

• there's a discussion developing at stats.stackexchange.com/questions/298509 about the convexity (wrt $\lambda$) of the test error Sep 3, 2017 at 4:31
• Do you mean the training set loss as a function of $\lambda$? In that case, the loss should decrease monotonically as $\lambda$ increases, because the hypothesis space for larger $\lambda$ always contains the hypothesis space for smaller $\lambda$ as a subset (so the loss on the training set can't be worse). The minimal loss would occur when $\lambda = \|\beta_{OLS}\|^2$, where $\beta_{OLS}$ is the ordinary least squares (unconstrained) solution. And, after this point, all solutions for larger $\lambda$ would be equivalent. Jun 21, 2018 at 16:08