just a thought: 
using linalg.solve looks good, In similar quadratic forms, I still use linalg.inv
I think scipy.linalg.cho_factor and cho_solve could be used in a similar way, but we could reuse cho_factor, for example in line 424 vinv_resid