diag() 的效率 - MATLAB

Efficiency of diag() - MATLAB

动机:

在编写要对数万个向量执行的矩阵运算时,我不断遇到警告:

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.

原因是我使用 diag() 来获取矩阵内积对角线上的值。因为 MATLAB 通常针对 vector/matrix 操作进行了优化,所以当我第一次编写代码时,我通常会选择矢量化形式。然而,在这种情况下,MATLAB 必须构建整个矩阵才能获得导致内存和速度问题的对角线。

实验:

我决定测试 diag()for 循环的使用,看看在任何时候使用 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');

这是测试结果: (注意:红线表示循环时间作为一个阈值——并不是说无论循环次数多少,总循环时间是恒定的)

从图中可以清楚地看出,将操作分解为大约 200x200 的方矩阵,使用 diag 比使用循环执行相同的操作要快。

问题:

有人可以解释为什么我会看到这些结果吗?此外,我认为随着 MATLAB 不断优化的设计,将在 diag() 函数调用中内置处理这些庞大的矩阵。例如,它可以只执行 i = j 索引操作。是否有特殊原因导致这可能令人望而却步?

我也没有真正考虑过使用分区方法 diag 的内存影响,尽管很明显随着分区大小的减小,内存需求也会下降。

诊断与循环的速度测试。

初始化:

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

诊断测试速度:

tic;
d = diag(M);
toc;

循环测试速度:

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

这将测试诊断。您的代码不是诊断的干净测试...

评论可能存在混淆的地方

Diag 仅提取矩阵的对角线。如果 xy 是向量,而您执行 d = diag(x * y'),MATLAB 首先构造 n x n 矩阵 x*y' 并对其调用 diag。这就是为什么您会收到错误消息“无法构造 290GB 矩阵...”Matlab 解释器 不会 以疯狂的方式进行优化,意识到您只需要对角线并只构造一个向量(而不是 x*y' 的完整矩阵, 不会 发生。

不确定您是否在问这个问题,但是计算 d = diag(x*y') 的最快方法是 xy 是 n x 1 个向量:d = x.*y