Fit ordinary least squares in closed form. With an intercept column of ones added to X, the coefficients solve
Implement fit_linear(X, y) returning (intercept, coefs) and predict_linear(X, intercept, coefs). Use np.linalg.solve, not an explicit inverse.