Thanks @ptrblck! I fixed the regularization but the error did not get much better. I was in fact using RidgeCV from sklearn, but the idea is same. I looked at the source code (thanks, that is always a good research strategy advice) and I think that the only difference now is that the code uses SVD for finding the coefficients: we find the derivative with respect to weights and find the solution in terms of SVD decomposition matrices.
As far as I understand, backprop does that for me in a network anyways and I don’t care about how to update weights other than my loss function. I can use pytorch to implement this but if I use SVD then it won’t be a network anymore.
For interested future reader, as can be seen in Bishop’s PRML book:

The particular case of a quadratic regularizer is called ridge
regression (Hoerl and Kennard, 1970). In the context of neural
networks, this approach is known as weight decay.

So we can either set weight_decay parameter in optimizer to a non-zero value or we can regularize the loss function with L2 as follows:
criterion = torch.nn.MSELoss()
lmbd = 1e-8  # for custom L2 regularization
loss = criterion(y_pred, y_train)

reg_loss = None
for param in model.parameters():
    if reg_loss is None:
        reg_loss = torch.sum(param ** 2)
    else:
        reg_loss = reg_loss + param.norm(2) ** 2

loss += lmbd * reg_loss