linreg.py

#Python 3.6.5 (default, Jun 17 2018, 12:13:06) 
#Type 'copyright', 'credits' or 'license' for more information
#IPython 6.4.0 -- An enhanced Interactive Python. Type '?' for help.
#Using matplotlib backend: MacOSX

In [6]: import scipy

In [7]: from scipy import linalg

In [1]: n, p = 10**5, 10**2

In [2]: X, y = np.random.randn(n,p), np.random.randn(n)

In [3]: XX, Xy = X.T @ X, X.T @ y

In [4]: %timeit np.linalg.inv(XX) @ Xy
728 µs ± 46.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [5]: %timeit np.linalg.solve(XX, Xy)
413 µs ± 42 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [8]: %timeit scipy.linalg.solve(XX, Xy)
673 µs ± 55.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [9]: %timeit scipy.linalg.solve(XX, Xy, sym_pos=True)
344 µs ± 25.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [10]: %timeit scipy.linalg.cho_solve(scipy.linalg.cho_factor(XX), Xy, check_finite=False)
222 µs ± 14 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Note: these times can vary significantly between samples and different processes running in the background, but np.linalg.solve() and scipy.linalg.cho_solve(scipy.linalg.cho_factor( ... )) tend to be consistently the best performing.