Блочная матрица внутренних произведений в Matlab

Question

Я использовал следующую пользовательскую функцию, чтобы выполнить умножение вектора на матрицу, в которой каждый элемент вектора умножает блок 3x3 в матрице (3xN) x (3):

function [B] = BlockScalar(v,A)

 N=size(v,2);
 B=zeros(3*N,3);

 for i=1:N
     B(3*i-2:3*i,:) = v(i).*A(3*i-2:3*i,:);
 end

end

Точно так же, когда я хочу умножить коллекцию матриц 3x3 на коллекцию векторов 3x3, я использую следующее

function [B] = BlockMatrix(A,u)

 N=size(u,2);
 B=zeros(N,3);

 for i=1:N
     B(i,:) = A(3*i-2:3*i,:)*u(:,i);
 end

end

Поскольку я вызываю их очень часто, это, к сожалению, замедляет работумой код значительно.Мне было интересно, существует ли более эффективная (возможно, векторизованная) версия вышеуказанных операций.

Answer 1

Вот векторизованные решения:

function [B] = BlockScalar(v,A)
    N = size(v,2);
    B = reshape(reshape(A,3,N,3) .* v, 3*N, 3);
end


function [B] = BlockMatrix(A,u)
    N = size(u,2);
    A_r = reshape(A,3,N,3);
    B = (A_r(:,:,1) .* u(1,:) + A_r(:,:,2) .* u(2,:) + A_r(:,:,3) .* u(3,:)).';
end

function [B] = BlockMatrix(A,u)
    N = size(u,2);
    B = sum(reshape(A,3,N,3) .* permute(u, [3 2 1]) ,3).';
end

Answer 2

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

Для первой функции вы можете выполнитьэто выглядит следующим образом:

function [B] = BlockScalar(v,A)
% We create a vector N = [1,1,1,2,2,2,3,3,3,...,N,N,N]
N=ceil((1:size(A,1))/3); 

% Use N to index v, and let matlab do the expansion
B = v(N).*A;

end

Для второй функции мы можем сделать блок-диагональную матрицу.

function [B] = BlockMatrix(A,u)

 N=size(u,2)*3;
 % We use a little meshgrid+sparse magic to convert A to a block matrix
 [X,Y] = meshgrid(1:N,1:3);

% Use sparse matrices to speed up multiplication and save space
 B = reshape(sparse(Y+(ceil((1:N)/3)-1)*3,X,A) * (u(:)),3,size(u,2))';

end

Обратите внимание, что если вы сможете получить доступ к отдельным матрицам 3x3, выпотенциально может сделать это быстрее / проще, используя собственный blkdiag:

function [B] = BlockMatrix(a,b,c,d,...,u)
% Where A = [a;b;c;d;...];

% We make one of the input matrices sparse to make the whole block matrix sparse
% This saves memory and potentially speeds up multiplication by a lot
% For small enough values of N, however, using sparse may slow things down.
reshape(blkdiag(sparse(a),b,c,d,...) * (u(:)),3,size(u,2))';

end