Q:

Lots of nested loops - How can I make it faster?

This part of my code, shown below, is taking long time to run due to many nested loops. Is there a way these nested loops can be obviated to make it run faster?
for k = 1:numel(UpscaledZLen.pq)
    for j = 1:numel(UpscaledRowLen.pq)
        for i = 1:numel(UpscaledColLen.pq)

            iZ = 1;
            while iZ <= UpscaledZLen.pq(k)
                for ZLag = 1:iZ

                    ihRow = 1;
                    while ihRow <= UpscaledRowLen.pq(j)
                        for hRowLag = 1:ihRow

                            ihCol = 1;
                            while ihCol <= UpscaledColLen.pq(i)
                                for hColLag = 1:ihCol
                                    temp1(hColLag) = trapz(AnalGamma(hRowLag,...
                                        1:hColLag,...
                                        ZLag));
                                end
                                temp2(ihCol) = trapz(temp1);
                                ihCol = ihCol + 1;
                            end

                            temp3(hRowLag) = trapz(temp2);
                        end
                        temp4(ihRow) = trapz(temp3);
                        ihRow = ihRow + 1;
                    end

                    temp5(ZLag) = trapz(temp4);
                end
                temp6(iZ) = trapz(temp5);
                iZ = iZ + 1;
            end

            NormVariance_AnalCorrAvg(i, j, k) = (2/((UpscaledRowLen.pq(j)*...
                UpscaledColLen.pq(i)*UpscaledZLen.pq(k))^2))*trapz(temp6);
        end
    end
end

This piece of code is trying to implement the following integral expression:

EDIT (SSCCE): For a brief example, I took following variable sizes and then used the above piece of code to see how much time it takes to run for these specific variable sizes:
nRow = 4;
nCol = 4;
nZ = 4;

RowLenScale.pq = 1:nRow;
ColLenScale.pq = 1:nCol;
UpscaledRowLen.pq = RowLenScale.pq(rem(RowLenScale.pq(end), RowLenScale.pq) == 0);
UpscaledColLen.pq = RowLenScale.pq(rem(ColLenScale.pq(end), ColLenScale.pq) == 0);

ZLenScale.pq = 1:nZ;
UpscaledZLen.pq = ZLenScale.pq(rem(ZLenScale.pq(end), ZLenScale.pq) == 0);

AnalGamma = rand(nRow, nCol, nZ);

For this example it takes only 0.321976 seconds, whereas, for original case where nRow = 100;, nCol = 100; and nZ = 20; it took more than 24 hours and was still running with outermost loop index at k = 4. 

A:

First of all, your code is faulty, because temporary variables are not initialized and values left from the last iteration are still present in the next. After correcting this and replacing while by for loops, the code looks like this:
NormVariance_AnalCorrAvg = nan(numel(UpscaledColLen.pq), numel(UpscaledRowLen.pq), numel(UpscaledZLen.pq));
for k = 1 : numel(UpscaledZLen.pq)
    for j = 1 : numel(UpscaledRowLen.pq)
        for i = 1 : numel(UpscaledColLen.pq)

            temp6 = nan(1, UpscaledZLen.pq(k));
            for iZ = 1 : UpscaledZLen.pq(k)
                temp5 = nan(1, iZ);
                for ZLag = 1 : iZ
                    temp4 = nan(1, UpscaledRowLen.pq(j));
                    for ihRow = 1 : UpscaledRowLen.pq(j)
                        temp3 = nan(1, ihRow);
                        for hRowLag = 1 : ihRow
                            temp2 = nan(1, UpscaledColLen.pq(i));
                            for ihCol = 1 : UpscaledColLen.pq(i)
                                temp1 = nan(1, ihCol);
                                for hColLag = 1:ihCol
                                    temp1(hColLag) = trapz(AnalGamma(hRowLag,...
                                        1:hColLag,...
                                        ZLag));
                                end
                                temp2(ihCol) = trapz(temp1);
                            end
                            temp3(hRowLag) = trapz(temp2);
                        end
                        temp4(ihRow) = trapz(temp3);
                    end
                    temp5(ZLag) = trapz(temp4);
                end
                temp6(iZ) = trapz(temp5);
            end
            NormVariance_AnalCorrAvg(i, j, k) = (2/((UpscaledRowLen.pq(j)*...
                UpscaledColLen.pq(i)*UpscaledZLen.pq(k))^2))*trapz(temp6);

        end
    end
end

The initialization here also leads to preallocation of arrays, which is one of the basic recommended things to speed up loops. In this case, however, it doesn't have an effect.
In order to benchmark the code I use your snippet using rand, but with
nRow = 20;
nCol = 20;
nZ = 2;

On my computer, the cleaned-up code needs 25.5 seconds to complete.

How to make it faster? Good old vectorization:
The innermost loop over hColLag can be replaced by a cumtrapz, descreasing the runtime to 8.2 seconds.
temp1 = cumtrapz(AnalGamma(hRowLag, 1 : ihCol, ZLag), 2);

The same holds for the loop over ihCol, decreasing the runtime to 1.85 seconds.
temp2 = cumtrapz(cumtrapz(AnalGamma(hRowLag, 1 : UpscaledColLen.pq(i), ZLag), 2), 2);

The next outer loop over hRowLag just serves to compute several analogous trapzs. That's not necessary, since trapz is fully vectorized; the loop can be replaced by one call, which decreases the runtime to 0.34 seconds.
temp3 = trapz(cumtrapz(cumtrapz(AnalGamma(1 : ihRow, 1 : UpscaledColLen.pq(i), ZLag), 2), 2), 2);

The loop over ihRow computes a cumulative integral; using cumtrapz brings the runtime down to about 0.07 seconds.
temp4 = cumtrapz(trapz(cumtrapz(cumtrapz(AnalGamma(1 : UpscaledRowLen.pq(j), 1 : UpscaledColLen.pq(i), ZLag), 2), 2), 2), 1);

Applying the same logic to the loops over iZ and ZLag, we arrive at the following code
NormVariance_AnalCorrAvg = nan(numel(UpscaledColLen.pq), numel(UpscaledRowLen.pq), numel(UpscaledZLen.pq));
for k = 1 : numel(UpscaledZLen.pq)
    for j = 1 : numel(UpscaledRowLen.pq)
        for i = 1 : numel(UpscaledColLen.pq)
            temp6 = cumtrapz(trapz(cumtrapz(trapz(cumtrapz(cumtrapz(...
                AnalGamma(1 : UpscaledRowLen.pq(j), 1 : UpscaledColLen.pq(i), 1 : UpscaledZLen.pq(k)), ...
                2), 2), 2), 1), 1), 3);
            NormVariance_AnalCorrAvg(i, j, k) = (2/((UpscaledRowLen.pq(j)*...
                UpscaledColLen.pq(i)*UpscaledZLen.pq(k))^2))*trapz(temp6);

        end
    end
end

with a runtime of 0.043 seconds, an improvement by a factor of about 600.
Writing all these nested trapzs and cumtrapzs is of course very error-prone. I used a checksum applied to the result saved to a file in order to make sure it is not changed by my changes to the code. Still, you should double check in order to be completely sure.

It is probably possible to go even further. The first three for loops again do something cumulative-y, so maybe more cumtrapz are in order. trapz is not very effectively implemented because it has a lot of error-checking code. Since what it does is basically a sum over averaged values (see trapezoidal rule), precomputing multidimensional averages on the full grid would allow to translate these calls into calls to sum. Simpler than the trapezoidal rule is the rectangle method, for which you would compute values of your function not at the edges but in the center of the grid cells. And then of course you could try to partially compute the integral analytically; and you can decrease the resolution of your grid.
Hope this helps!

