(This could be implemnted via glm(... family=gaussian) but the implementation with lsfit is much faster
fit_linear(y, X)
y
numeric vector
X
numeric matrix (covariates)