It looks pretty good, but I’d love to see you make better use of the routines already in numpy. In particular, I see you are solving the OLS problem using direct inversion when you already have QR and SVD available to just call. There are other simple things that can do a lot for your results, centering, scaling, etc, too. I guess it works well enough for small well behaved problems as is though.