# Optimization of the regularized least squares with gradient descent

I have this regularized least square formula: $$\sum\limits_{i=1}^N (\omega^T x_i - y_i)^2 + \lambda \left\|\omega\right\|^2$$

And the gradient: $$2 \sum\limits_{i=1}^N ((\sum\limits_{j=1}^d x_{ij}\omega_j)x_{ik} - x_{ik} y_i) + 2\lambda \omega_k$$

I want to use gradient descent to find the vector w. I am using matlab. I though I would be able to make two loops and calculate the ws but my solution is very unstable and I need to use very small learning term a (a=0.000000001) in order to get not NAN solution. But I thought the values of w should head towards 0 when the lambda is large but it does not happen... My data set is a matrix X (400x64) and y (400x1). This is two class problem where contains the class labels (+1 for class 1 and -1 for class 2).

Here is my matlab code:

function [ w ] = gradDecent( X, Y, a, lambda, iter )

w = zeros(size(X(1,:)))';

for it=1:iter % For each iteration
for k = 1:size(w,1)
s = 0;
for i = 1:size(X,1)
s = s + (X(i,:)*w - Y(i))*X(i,k);
end
w(k) = w(k) - a*(2*s+2*lambda*w(k));
end
end


Am I making some stupid mistakes?

• Because this is a quadratic function of $\omega$, with linear gradient, why don't you just use the usual least-squares solution?
– whuber
Feb 19, 2014 at 18:57
• You mean the analytic solution for ridge (Tikhonov regularized) regression? That seems to be what the question is looking for... not "usual least-squares" which is generally used to mean ordinary least squares.
– ely
Feb 21, 2014 at 18:06

1. Instead of using a constant learning rate, set your learning rate to a reasonably large value (try $\alpha=1$) and then use "backoff" to shrink your learning rate as appropriate using the following algorithm:

while $f(w^{(k)} - \alpha d) \gt f(w^{(k)})$: $\alpha := \alpha/2$

$w^{(k+1)} := w^{(k)} - \alpha d$

where $f$ is your objective function, $d$ is your descent direction, and $\alpha$ is the learning rate.

2. I don't think you've correctly implemented the gradient in your code. In particular, I believe the way you have it coded the s term evaluates to $s=\sum_{i=1}^N ( x_i \omega - y_i x_i )$ when you really want $s=\sum\limits_{i=1}^N ((\sum\limits_{j=1}^d x_{ij}\omega_j)x_{ik} - x_{ik} y_i)$.

3. Gradient Descent probably isn't the best solution here. Gradient descent is slow: you shouldn't be surprised that it's taking a long time to converge, because gradient descent usually does. It gets the job done, but it's generally a slow option. Try playing with other optimization alogorithms and see what happens. Try using newton-raphson, for instance. I'm sure if you poke around the literature you'll be able to determine the "state-of-the-art" algorithm choice for optimizing the ridge regression loss function. I'd put money on it not being simple gradient descent.

As whuber commented, there is no need to bother with optimization algorithms. Gradient in matrix notation is $$2X'(Xw-y)+2\lambda w$$ Equating it to zero in Matlab is a one-liner:

w_star = (X'*X+lambda*eye(d)) \ X'*y;