Julia: Why ridge regression not working (optim)

Question

I am trying to implement ridge-regression from scratch in Julia but something is going wrong.

# Imports
using DataFrames
using LinearAlgebra: norm, I
using Optim: optimize, LBFGS, minimizer

# Read Data
out = CSV.read(download("https://raw.githubusercontent.com/jbrownlee/Datasets/master/housing.csv"), DataFrame, header=0)

# Separate features and response
y = Vector(out[:, end])
X = Matrix(out[:, 1:(end-1)])
λ = 0.1

# Functions
loss(beta) = norm(y - X * beta)^2 + λ*norm(beta)^2
function grad!(G, beta) 
    G = -2*transpose(X) * (y - X * beta) + 2*λ*beta
end
function hessian!(H, beta)
   H =  X'X + λ*I
end

# Optimization
start = randn(13)
out = optimize(loss, grad!, hessian!, start, LBFGS())

However, the result of this is terrible and we essentially get back start since it is not moving. Of course, I know I could simply use (X'X + λ*I) \ X'y or IterativeSolvers.lmsr(X, y) but I would like to implement this myself.

user_17153821 · Answer

The problem is with the implementation of the grad! and hessian! functions: you should use dot assignment to change the content of the G and H matrices:

G .= -2*transpose(X) * (y - X * beta) + 2*λ*beta
H .= X'X + λ*I

Without the dot you replace the matrix the function parameter refers to, but the matrix passed to the function (which will then be used by the optimizer) remains unchanged (presumably a zero matrix, that's why you got back the start vector).

Julia: Why ridge regression not working (optim)

1 Answers