Efficiency of diag() - MATLAB

Question

Motivation:

In writing out a matrix operation that was to be performed over tens of thousands of vectors I kept coming across the warning:

Requested 200000x200000 (298.0GB) array exceeds maximum array size preference. Creation of arrays greater than this limit may take a long time and cause MATLAB to become unresponsive. See array size limit or preference panel for more information.

The reason for this was my use of diag() to get the values down the diagonal of an matrix inner product. Because MATLAB is generally optimized for vector/matrix operations, when I first write code, I usually go for the vectorized form. In this case, however, MATLAB has to build the entire matrix in order to get the diagonal which causes the memory and speed issues.

Experiment:

I decided to test the use of diag() vs a for loop to see if at any point it was more efficient to use diag():

num = 200000; % Matrix dimension
x = ones(num, 1);
y = 2 * ones(num, 1);

% z = diag(x*y'); % Expression to solve

% Loop approach
tic
z = zeros(num,1);
for i = 1 : num
   z(i) =  x(i)*y(i);
end
toc

% Dividing the too-large matrix into process-able chunks
fraction = [10, 20, 50, 100, 500, 1000, 5000, 10000, 20000];
time = zeros(size(fraction)); 
for k = 1 : length(fraction)
    f = fraction(k);

    % Operation to time
    tic
    z = zeros(num,1);
    for i = 1 : k
        first = (i-1) * (num / f);
        last = first + (num / f);
        z(first + 1 : last) = diag(x(first + 1: last) * y(first + 1 : last)');
    end
    time(k) = toc;
end

% Plot results
figure;
hold on
plot(log10(fraction), log10(chunkTime));
plot(log10(fraction), repmat(log10(loopTime), 1, length(fraction)));
plot(log10(fraction), log10(chunkTime), 'g*'); % Plot points along time
legend('Partioned Running Time', 'Loop Running Time');
xlabel('Log_{10}(Fractional Size)'), ylabel('Log_{10}(Running Time)'), title('Running Time Comparison');

This is the result of the test: (NOTE: The red line represents the loop time as a threshold--it's not to say that the total loop time is constant regardless of the number of loops)

From the graph it is clear that it takes breaking the operations down into roughly 200x200 square matrices to be faster to use diag than to perform the same operation using loops.

Question:

Can someone explain why I'm seeing these results? Also, I would think that with MATLAB's ever-more optimized design, there would be built-in handling of these massive matrices within a diag() function call. For example, it could just perform the i = j indexed operations. Is there a particular reason why this might be prohibitive?

I also haven't really thought of memory implications for diag using the partition method, although it's clear that as the partition size decreases, memory requirements drop.

Answer 1

Test of speed of diag vs. a loop.

Initialization:

n = 10000;
M = randn(n, n);  %create a random matrix.

Test speed of diag:

tic;
d = diag(M);
toc;

Test speed of loop:

tic;
d = zeros(n, 1);
for i=1:n
   d(i) = M(i,i);
end;
toc;

This would test diag. Your code is not a clean test of diag...

Comment on where there might be confusion

Diag only extracts the diagonal of a matrix. If x and y are vectors, and you do d = diag(x * y'), MATLAB first constructs the n by n matrix x*y' and calls diag on that. This is why, you get the error, "cannot construct 290GB matrix..." Matlab interpreter does not optimize in a crazy way, realize you only want the diagonal and construct just a vector (rather than full matrix with x*y', that does not happen.

Not sure if you're asking this, but the fastest way to calculate d = diag(x*y') where x and y are n by 1 vectors would simply be: d = x.*y