Умножьте 5D матрицу на 2D матрицу - PullRequest
2 голосов
/ 19 декабря 2011

У меня есть 5D матрица Cij (3,3, Nx, Ny, Nz), где Nx, Ny и Nz приведены в качестве входных данных.

Мне нужно выполнить что-то вроде этого:

for ikx=1:Nx,
    for iky=1:Ny,
        for ikz=1:Nz,

            %Random simulation of fourier components
            n=zeros((3),'double');
            for j=1:9,
                ncomponent=randn(2);
                n(j)=complex(ncomponent(1),ncomponent(2));
                %Calculation of H
                H(:,ikx,iky,ikz)=dot(Cij(:,:,ikx,iky,ikz),n);
            end;
        end;
    end;
end;

Проблема в том, что для увеличения Nx, Ny, Nz в цикле требуется очень много времени для вычисления матрицы H.

Кто-нибудь знает более быстрый способ получить H-матрицу?

Ответы [ 2 ]

3 голосов
/ 19 декабря 2011

Прежде всего следует отметить, что в самом внутреннем цикле вы выполняете скалярное произведение 9 раз, каждый раз перезаписывая H(:,ikx,iky,ikz). Там нет никакого смысла в этом. Вы должны просто заполнить случайные значения для n в цикле и вычислить H(:,ikx,iky,ikz) один раз после этого цикла.

Однако все циклы не нужны, поскольку вы можете воспользоваться тем фактом, что функция DOT векторизована и может обрабатывать 5-мерные массивы (т. Е. Она будет автоматически выполнять операцию с точками на первом одно измерение). Все, что вам нужно сделать, это сделать n матрицей комплексных значений 3 на 3 на Nx на Ny на Nz. Эти две строки должны дать вам тот же результат, что и код выше:

n = complex(rand([3 3 Nx Ny Nz]), rand([3 3 Nx Ny Nz]));
H = squeeze(dot(Cij, n));

Функция SQUEEZE используется для удаления одноэлементных измерений из H, что сделает ее матрицей 3 на Nx на Ny на Nz.

1 голос
/ 19 декабря 2011

Вы можете сделать это, используя reshapepermute)

C=rand(3,3,Nx,Ny,Nz);
n=rand(3,3);

Если вы хотите умножение mtrix между n и каждым элементом C:

H=reshape(n*reshape(C,3,3*Nx*Ny*Nz),[3,3,Nx,Ny,Nz])

Если вы хотите получить скалярное произведение между n и каждым элементом C:

H=reshape(reshape(n,1,[])*reshape(C,3*3,Nx*Ny*Nz),[Nx,Ny,Nz])
...