`cholsigmainv` is not a matrix, but a `nobs` array, and so element-by-element with broadcasting is needed.

When `nobs` is very large and `sigma` is a vector or scalar, creating a `nobs` by `nobs` array is problematic.  For example, when `nobs=100000` then this array requires about 64GB of memory. 

This tries to avoid this penalty.  It is also algorithmically superior since multiplying a diagonal `nobs` array requires _O_(`nobs**2`) operations while this is _O_(`nobs`) 

Or maybe I am misunderstanding?