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.
