Improving calculation of multivariate multiple regression¶

A prototype for multivariate multiple regression computation assuming fixed residual covariance matrix.

data = mvsusieR::mvsusie_sim2(1000,7000,50,0.8)

Here we do not use an existing R routine because lm is very slow (even using .lm.fit compared to providing pre-computed d as done below).

method_1 = function(d,X,y,residual_covar,p) {
    XtY=t(X)%*%y
    bhat = diag(1/d) %*% XtY
    s2 = diag(residual_covar)
    sbhat = sqrt(do.call(rbind, lapply(1:p, function(j) s2 / d[j])))
    return(list(b=bhat,s=sbhat))
}

method = method_1

res = with(data, method(d,X,y,diag(r),p))

The speed is decent but the code seems ugly. Must be a more elegant matrix representation (and hopefully faster)?