How to transform large matrix calculations into smaller scale in matlab

I have to calculate fourth-order cumulants with large matrices in my current research project. The following is part of my code.

N=zeros((Kse*Ksf*Ksr)^2); 

for k=1:Ke   % Ke=3
for i=1:Kr  % Kr=3
  for j=1:Kf  % Kf=1;
      Mat=Xb(k:k+Kse-1,i:i+Ksr-1,j:j+Ksf-1); % size of Xb: (Ke+Kse-1)*(Kr+Ksr-1)*(Kf+Ksf-1)

        Mat_vec=Mat(:);

        snapshot_number=70;

        Mat_kr=kron(Mat_vec,conj(Mat_vec));

        N=N+Mat_kr*Mat_kr'/snapshot_number;
        N=N-Mat_kr*Mat_kr'/snapshot_number.^2;

        temp=Mat_vec*Mat_vec'/snapshot_number;
        N=N-kron(temp,conj(temp));
      end
    end
  end

Under the conditions of Kse=Ksr=5 and Ksf=9, the running memory is about 200GB. The size of N is 50625*50625. My available server has the maximum memory of 250GB. The key limit is memory. So I can not increase the parameters of Kse, Ksr and Ksf any more. My target is to set Kse=Ksr=Ksf=11. So I have thought of matlab distributed computing. Firstly, I shrinked the parameter to Kse=3, Ksr=5 and Ksf=8 just for test. I modified the code to the following editon:

myPool=parpool();
Xb_dis=distributed(Xb);
N_dis=distributed(zeros((Kse*Ksf*Ksr)^2)); 

for k=1:Ke   % Ke=3
  for i=1:Kr  % Kr=3
    for j=1:Kf  % Kf=1;
      Mat_dis=Xb_dis(k:k+Kse-1,i:i+Ksr-1,j:j+Ksf-1);
      Mat_vec_dis=Mat_dis(:);

        snapshot_number=70;

        Mat_kr_dis=kron(Mat_vec_dis,conj(Mat_vec_dis));

        N_dis=N_dis+Mat_kr_dis*Mat_kr_dis'/snapshot_number;
        N_dis=N_dis-Mat_kr_dis*Mat_kr_dis'/snapshot_number.^2;

        temp_dis=Mat_vec_dis*Mat_vec_dis'/snapshot_number;
        N_dis=N_dis-kron(temp_dis,conj(temp_dis)); 
      end
    end
  end 
  N=gather(N_dis);

The experimental distributed computing cluster consists of two computers. I found that there is little data communication between the two computers during the triple loops. So I thought that the calculations can be done in smaller scale by sequences instead of loading the whole matrices into the memory. In other words, the calculation of large scale matrix may equal to combinations of smaller ones. But after several days in searching the Internet, I still have no idea.

Anybody have some advices or suggestions for this question?

[David Goodmanson]

Hello ZZ,

Is there anything you are leaving out in the code above for simplicity? The matrix N depends only on Mat_vec and snapshot_number (sn). It appears that after one run through the for loop, N is incremented by the amount

Mat_kr*Mat_kr'*(1/sn - 2/sn^2)