Matrix Guru Needed : How to compute a weighted covariance matrix fast

Hi,

I'm trying to speed up the function below and have hit a wall. It's a bog standard covariance weighting calculation but is incredibly slow because the calculation is O(n^4).

If this helps for optimisation, I actually call this at every time step of a model such that both matricies below are really of dimension [num_time, num_series, num_series] but I do a for loop over the function. I've often used the great tool below but can't see how to make this one faster

https://au.mathworks.com/matlabcentral/fileexchange/8773-multiple-matrix-multiplications-with-array-expansion-enabled

Many math thanks,

Lyle

function Vw = covWeight(V,weights)
% Vw = covWeight(V,weights)
%
% Recompute a covariance matrix based on a set of weightings for each
% row,column
%
% For example, say you want to blend a set of time series X with known 
% covariance V into a set of new time series Y such that
%
%           Y1 = aX1 + bX2 + cX3
%           Y2 = dX1 + eX2 + fX3
%           Y3 = gX1 + hX2 + iX3
%
% then the weights would be
% 
%                [a b c]
%     weights =  [d e f]
%                [g h i]
%
% and the covariance of the new time series would be given by Vw.
%
% DEBUG
%           V - covWeight(V, eye(size(V))) == 0
%
% because the weighting is simply the original time series X
%
% ***** This is VERY computationally expensive *********

Vw     = nan(size(V));
inplay = abs(weights) >eps('single');

n_a = size(V,1);
for j=1:n_a
    for k=j:n_a    
        a   = weights(j,inplay(j,:));
        b   = weights(k,inplay(k,:));
        tmp = a'*b;
        
        Vs = V(inplay(j,:), inplay(k,:));
                
        Vw(j,k) = tmp(:)'*Vs(:);
    end
end

% Reflect upper to lower 
Vw = triu(Vw,0) + triu(Vw,1)';

end