Write a code for linear regression in bare NumPy.
Implement ordinary least squares using only NumPy, automatically adding an intercept term and supporting multiple features. The implementation should remain defined for rank-deficient or underdetermined design matrices by returning the minimum-norm least-squares solution. Input: X, a numeric 2D array of shape (n_samples, n_features), and y, a numeric 1D array of shape (n_samples,). Output: a 1D array of shape (n_features + 1,), ordered as intercept followed by feature coefficients.
def linear_regression(X, y):