The best thing to do would be to exploit conjugacy directly and avoid evaluating the normal altogether.
You should put this in transformed data becaue it’s a constant. I’d also recommend using compound declare-define, e.g…, matrix[N, K] w_mat = make_w_matrix(N, K, w);.
It’d be even more efficient if you could work at the level of a Cholesky factor. The Jacobian looks daunting for that, but it’d be the way to go for underlying speed.