Optimization of the regularized least squares with gradient descent

Question

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 )
% GRADIENT DESCENT

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?

Answer 1

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.

Answer 2

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;