Подстановка симметричных матриц по диагонали

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 000 x 50 000), используя цикл for, как показано выше, в нежизнеспособных.Есть ли способ достичь вышеупомянутого без цикла?Я думал о других подходах, которые, возможно, могли бы использовать верхнюю и нижнюю треугольные матрицы, однако даже они, кажется, нуждаются в петле в некоторой точке.Любая помощь приветствуется !!

Answer 1

Вот альтернатива, которая не генерирует полную матрицу для выбора диагоналей блока, как в Luis Mendo 'answer , но вместо этого напрямую генерирует индексы для этих элементов.Вполне вероятно, что это будет быстрее для очень больших матриц, поскольку в этом случае создание матрицы индексации будет дорогостоящим.

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 (матрица маски Луиса): 0,0599 с
• method3 (индексы Cris): 1,5617e-04 с

Обратите внимание, что для массива 5000x5000 метод в этом ответе в ~ 20 раз быстрее, чем цикл, тогда как цикл ~ 20быстрее, чем решение Луиса.

Для матриц меньшего размера метод немного отличается, метод Луиса почти в два раза быстрее, чем код цикла для матрицы 50x50 (хотя этот метод все еще превосходит его в ~ 3 раза).

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.