MATLAB中矩阵和“矩阵列表”的有效乘法

Question

我想计算以下矩阵

M(i,j)=sum_(k,l) c(k,l) kron(N(:,:,k),N(:,:,l))

其中N 的前两个索引是行数和列数（在这种情况下相等），最后一个索引表示我有一个矩阵族N，而矩阵c 只是一个系数矩阵。指数 k 和 l 从 1 到 50（在我的具体情况下）。

为此，我实现了以下算法：

c=randn(50,50);
N=randn(26,26,50);
M=zeros(size(N,1)^2);
for k=1:size(N,3)
for l=1:size(N,3)
M=M+c(k,l).*kron(N(:,:,k),N(:,:,l));
end
end

这段代码有一些有趣的地方，我不知道为什么会发生：

假设我截断了k 中的总和，直到一个有限的数字n。 k 从 1 变为 n 与从 1 变为 n+1 时的计算时间差远高于（大约和数量级）仅计算由索引的矩阵所需的时间k=n+1 并将其与计算到 k=n 的总矩阵相加。

我的问题是：是否有更有效的方法来计算我正在寻找的矩阵？这个算法出了什么问题？

非常感谢。

Answer 1

TL;DR：这比你的代码慢，我只是想分享我尝试过的......

---

下面是我在 cmets (funcAmro) 中提到的想法的实现，与您的原始代码 (funcAlex) 相比：

function [t,v] = testKron()
    s1 = 26;
    s3 = 50;
    c = randn(s3,s3);
    N = randn(s1,s1,s3);

    funcs = {
        @() funcAlex(N,c)
        @() funcAmro(N,c)
    };
    t = cellfun(@timeit, funcs);
    v = cellfun(@feval, funcs, 'Uniform',false);
    norm(v{1}-v{2})
end

function M = funcAlex(N,c)
    [s1,~,s3] = size(N);
    M = zeros(s1*s1);
    for i=1:s3
        for j=1:s3
            M = M + c(i,j) .* kron(N(:,:,i), N(:,:,j));
        end
    end
end

function M = funcAmro(N,c)
    [s1,~,s3] = size(N);
    % tall matrix, NN = cat(1, N(:,:,1), ..., N(:,:,s3))
    NN = reshape(permute(N, [1 3 2]), s1*s3, s1);
    % accumulate kron results
    M = zeros(s1*s1*s3, s1*s1);
    for i=1:s3
        M = M + bsxfun(@times, repelem(c(:,i), s1*s1), kron(NN, N(:,:,i)));
    end
    % split M back into s3 slices each of size s1^2-by-s1^2
    M = permute(reshape(M.', s1*s1, s1*s1, s3), [2 1 3]);
    % sum along slices
    M = sum(M,3);
end

使用您提到的大小（N 是一个 26×26×50 数组），我在运行 R2016b 的机器上得到以下时序：

>> t = testKron
t =
    3.2770    % funcAlex
    5.5119    % funcAmro

所以我的方法比较慢（更不用说可读性差了！）。我想所有的置换/重塑都有不可忽略的开销..

当然，得到的矩阵 M 是等价的（差异是 4.5801e-12 接近机器 epsilon）。