沿对角线子集对称矩阵

Question

我有一个8x8的矩阵，例如A=rand(8,8)。我需要做的是沿对角线子集所有2x2的矩阵。这意味着我需要保存矩阵A(1:2,1:2)、A(3:4,3:4)、A(5:6,5:6)和A(7:8,7:8)。为了更好地解释我的需求，我目前使用的版本如下：

 AA = rand(8,8);
 BB = zeros(8,2);
 for i = 1:4
     BB(2*i-1:2*i,:) = AA(2*i-1:2*i,2*i-1:2*i);
 end

这对于小的AA矩阵和小的AA子矩阵工作得很好，但是随着尺寸显著增长（甚至可以达到50,000x50,000），像上面那样使用for循环是不可行的。有没有一种方法可以在不使用循环的情况下实现上述内容？我考虑过其他方法，也许可以利用上三角和下三角矩阵，但即使这些矩阵在某些时候似乎也需要循环。任何帮助都将不胜感激！

- J.Avra

---

AA 的大小总是可以被2整除的吗？ - Luis Mendo

不一定。我们需要对子集进行分割的矩阵大小可能会发生变化，但它们始终是原始矩阵的完美分割器。例如，对于当前矩阵，我们也可以沿对角线子集所有4x4矩阵（但不能是3x3或5x5）。对于一个9x9矩阵，我们将沿对角线子集3x3矩阵等。 - J.Avra

Answer 1

这里有一种替代方法，不像Luis Mendo的答案那样生成一个完整的矩阵来选择块对角线元素，而是直接生成这些元素的索引。对于非常大的矩阵来说，创建索引矩阵会很昂贵，因此这种方法可能会更快。

AA = rand(8,8); % example matrix. Assumed square
n = 2; % submatrix size. Assumed to divide the size of A

m=size(AA,1);
bi = (1:n)+(0:m:n*m-1).'; % indices for elements of one block
bi = bi(:);               % turn into column vector
di = 1:n*(m+1):m*m;       % indices for first element of each block
BB = AA(di+bi-1);         % extract the relevant elements
BB = reshape(BB,n,[]).'   % put these elements in the desired order

---

基准测试

AA = rand(5000); % couldn't do 50000x50000 because that's too large!
n = 2;

BB1 = method1(AA,n);
BB2 = method2(AA,n);
BB3 = method3(AA,n);
assert(isequal(BB1,BB2))
assert(isequal(BB1,BB3))

timeit(@()method1(AA,n))
timeit(@()method2(AA,n))
timeit(@()method3(AA,n))

% OP's loop
function BB = method1(AA,n)
m = size(AA,1);
BB = zeros(m,n);
for i = 1:m/n
   BB(n*(i-1)+1:n*i,:) = AA(n*(i-1)+1:n*i,n*(i-1)+1:n*i);
end
end

% Luis' mask matrix
function BB = method2(AA,n)
mask = repelem(logical(eye(size(AA,1)/n)), n, n);
BB = reshape(permute(reshape(AA(mask), n, n, []), [1 3 2]), [], n);
end

% Cris' indices
function BB = method3(AA,n)
m = size(AA,1);
bi = (1:n)+(0:m:n*m-1).';
bi = bi(:);
di = 0:n*(m+1):m*m-1;
BB = reshape(AA(di+bi),n,[]).';
end

在我的电脑上，使用MATLAB R2017a，我得到以下结果：

• method1（OP的循环）：0.0034秒
• method2（Luis的掩码矩阵）：0.0599秒
• method3 （Cris的索引）：1.5617e-04秒

请注意，对于一个5000x5000的数组，这个答案中的方法比一个循环快约20倍，而循环又比Luis的解决方案快约20倍。

对于较小的矩阵，情况略有不同，对于50x50的矩阵，Luis的方法几乎比循环代码快两倍（尽管这个方法仍然比它快约3倍）。

- Cris Luengo

Answer 2

这是一种方法：

AA = rand(8,8); % example matrix. Assumed square
n = 2; % submatrix size. Assumed to divide the size of A
mask = repelem(logical(eye(size(AA,1)/n)), n, n);
BB = reshape(permute(reshape(AA(mask), n, n, []), [1 3 2]), [], n);

这将生成一个逻辑掩码，选择所需的元素，然后使用reshape和permute按需要重新排列它们。

- Luis Mendo

---

LM，真是太棒了。 - greengrass62

@greengrass62 谢谢 :-) - Luis Mendo

太棒了！我自己永远也想不出如此优雅的东西。非常感谢LM！！ - J.Avra

@J.Avra 很高兴你喜欢它。 - Luis Mendo

你会讨厌我的。 你可能已经讨厌我了。 在我的电脑上，OP的循环代码比这个解决方案快约20倍。 :) - Cris Luengo

@Cris 我就猜到了 :-) - Luis Mendo