Speed up repeated kronecker product

I need to run several times (around 70-80 times) few lines of code, that include a kronecker product with a huge matrix.

In the code below this huge matrix is the sparse nanYt of dimensions 15k x 15k or greater, depending on t it can reach 50kx50k.

   
    denom = zeros(n*r,n*r);
    nom = zeros(n,r);
    y_b_nan = ~isnan(y_b(:,t))
    y_b(isnan(y_b(:,t))) = 0
    
    for t=1:T
        
        nanYt = sparse(diag(y_b_nan));
        denom = denom + kron(Z(:,t+1)*Z(:,t+1)',nanYt);
        nom = nom + y_b(:,t)*Zsmooth(r_cb_ind,t+1)';
       
    end
   
    vec_C = denom\nom(:);
    C_par = reshape(vec_C,n,r);
 

The three lines within the cycle need to be called 100 times. And, in particular, the one udating denom takes quite a while (4.5 second) and for repeated calculiations computational time increases too much.

Parallelizing is not an option, because this is a third-level cycle and I need to parallelize at the first (the more external) level.

Any idea or suggestion to speed up here?