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

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

Всего 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 .


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

 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 я получаю:

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

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


Есть идеи?

10000