Копируем векторы, сдвигая их вправо - PullRequest
Новая
       35

Копируем векторы, сдвигая их вправо

0 голосов
/ 02 мая 2018

В Matlab у меня есть два однорядных (1x249) вектора в матрице 2x249, и мне нужно создать матрицу A, реплицируя их много раз, каждый раз сдвигая векторы на 2 позиции вправо. Я хотел бы заполнить записи слева нулями. Есть ли умный способ сделать это? В настоящее время я использую цикл for и circshift, и на каждой итерации я добавляю новую строку в A, но, вероятно, это крайне неэффективно.

Код (myMat - это матрица, которую я хочу сдвинуть): 

A = [];
myMat = [1 0 -1 zeros(1,246); 0 2 0 -2 zeros(1,245)];
N = 20;
for i=1:N-1
    aux = circshift(myMat,[0,2*(i-1)]);
    aux(:,1:2*(i-1)) = 0;
    A =[A; aux];
end

Ответы [ 2 ]

0 голосов
/ 02 мая 2018

Я взял код Натана ( см. Его ответ на этот вопрос ) и добавил еще одну возможную реализацию (replvector3).

Моя идея здесь проистекает из того, что вам не нужен круговой сдвиг. Вам нужно сдвинуть вправо и добавить нули влево. Если вы начинаете с предварительно выделенного массива (это действительно то, где большие выигрыши по времени для вас, а остальное - арахис), то у вас уже есть нули. Теперь вам просто нужно скопировать myMat в нужное место.

Это время, которое я вижу (MATLAB R2017a): 

OP's, with pre-allocation: 1.1730e-04
Nathan's:                   5.1992e-05
Mine:                       3.5426e-05
                            ^ shift by one on purpose, to make comparison of times easier

Это полная копия, скопируйте и вставьте в M-файл и запустите: 

function so

shift_right = 2;
width       = 249;
myMat       =  [ 1 0 -1  0 ;
                 0 2  0 -2 ];
N           = 20;

A = replvector1(myMat,shift_right,width,N);
B = replvector2(myMat,shift_right,width,N);
norm(B(:)-A(:))
C = replvector3(myMat,shift_right,width,N);
norm(C(:)-A(:))

timeit(@()replvector1(myMat,shift_right,width,N))
timeit(@()replvector2(myMat,shift_right,width,N))
timeit(@()replvector3(myMat,shift_right,width,N))

% Original version, modified to pre-allocate
function A = replvector1(myMat,shift_right,width,N)
    % Assuming width > shift_right * (N-1) + size(myMat,2)
    myMat(1,width) = 0;
    M = size(myMat,1);
    A = zeros(M*(N-1),width);
    for i = 1:N-1
        aux = circshift(myMat,[0,shift_right*(i-1)]);
        aux(:,1:shift_right*(i-1)) = 0;
        A(M*(i-1)+(1:M),:) = aux;
    end

% Nathan's version
function A = replvector2(myMat,shift_right,width,N)
    [i,j,a] = find(myMat);
    i = kron(ones(N-1,1),i) + kron((0:N-2)',ones(size(i))) * size(myMat,1);
    j = kron(ones(N-1,1),j) + kron((0:N-2)',ones(size(j))) * shift_right;
    a = kron(ones(N-1,1),a);
    ok = j<=width;
    A = full(sparse(i(ok),j(ok),a(ok),(N-1)*size(myMat,1),width));

% My trivial version with loops
function A = replvector3(myMat,shift_right,width,N)
    % Assuming width > shift_right * (N-1) + size(myMat,2)
    [M,K] = size(myMat);
    A = zeros(M*(N-1),width);
    for i = 1:N-1
        A(M*(i-1)+(1:M),shift_right*(i-1)+(1:K)) = myMat;
    end
0 голосов
/ 02 мая 2018

Как вы, наверное, знаете, циклы в Matlab не так эффективны. Я знаю, что Mathworks постоянно говорит, что это не так с JIT компиляция, но я еще не испытывал быстрых циклов.

Я поместил ваш метод построения матрицы A в функцию: 

function A = replvector1(myMat,shift_right,width,N)

    pre_alloc = true; % make implementation faster using pre-allocation yes/no

    % Pad myMat with zeros to make it wide enough
    myMat(1,width)=0;

    % initialize A
    if pre_alloc
       A = zeros(size(myMat,1)*(N-1),width);
    else
       A = [];
    end

    % Fill A
    for i=1:N-1
        aux = circshift(myMat,[0,shift_right*(i-1)]);
        aux(:,1:min(width,shift_right*(i-1))) = 0;

        A(size(myMat,1)*(i-1)+1:size(myMat,1)*i,:) =aux;
    end

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

function A = replvector2(myMat,shift_right,width,N)
    [i,j,a] = find(myMat);
    i = kron(ones(N-1,1),i) + kron([0:N-2]',ones(size(i))) * size(myMat,1);
    j = kron(ones(N-1,1),j) + kron([0:N-2]',ones(size(j))) * shift_right;
    a = kron(ones(N-1,1),a);
    ok = j<=width;
    A = full(sparse(i(ok),j(ok),a(ok),(N-1)*size(myMat,1),width));

Вы можете следовать алгоритму, удалив точку с запятой и посмотрев на промежуточный Результаты.

Следующая основная программа запускает ваш пример и может быть легко изменена на запустите похожие примеры: 

% inputs (you may vary them to see that it always works)
shift_right = 2;
width       = 249;
myMat1      =  [ 1 0 -1  0 ;
                 0 2  0 -2 ];
N           = 20;

% Run your implementation
tic;
A = replvector1(myMat,shift_right,width,N);
disp(sprintf('\n   original implementation took %e sec',toc))

% Run the new implementation
tic;
B = replvector2(myMat,shift_right,width,N);
disp(sprintf('   new implementation took %e sec',toc))

disp(sprintf('\n   norm(B-A)=%e\n',norm(B-A)))
...