Hi,
Gradient descent is really bad for optimization 
In particular, you have a smooth convex function for which gd will have a bad convergence rate (one over t?). The original Adam does not always converge even in the convex case.
For the gradient descent methods, precise tuning of learning rate (and use of momentum) might be necessary to converge fast enough.
I am less familiar with LBFGS so I’m not sure why it converges to a different point.
I don’t know how LinearRegression from sklearn works. But does it really just minimizes the average l2 difference? Doesn’t it contain some regularization? Why do you optimize the l2 and compare the final l1 norms?