
function w = reg_least_sq(X,y,lambda,w_0)

n = size(X,1); % number of (x,y) samples
M = size(X,2); % number of features

% precompute a_k's (doesn't vary w/ w)
a = sum(X.*X)'/n;

w = w_0;
niter = 0;
while (1)

  maxerr = 0.0;
  for k = 1:M

    w_diff_k = w;
    w_diff_k(k) = 0.0;
    e_diff_k = y - X*w_diff_k; % prediction error w/o feature k
    c_k = X(:,k)'*e_diff_k/n;  % inner product of feature k w/ prediction error

    if (c_k < -lambda)
      w_hat_k = (c_k+lambda)/a(k);
    elseif (c_k > lambda)
      w_hat_k = (c_k-lambda)/a(k);
    else
      w_hat_k = 0.0; % l1-regularization forces weights of uninformative features to zero
    end
    maxerr = max(maxerr,abs(w(k) - w_hat_k));
    w(k) = w_hat_k;

  end

  niter = niter + 1;
  if (maxerr < 1e-8) return; end

end
