Replicate vectors shifting them to the right

Question

In Matlab, I have two single row (1x249) vectors in a 2x249 matrix and I have to create a matrix A by replicating them many times, each time shifting the vectors of 2 positions to the right. I would like to fill the entries on the left with zeros. Is there a smart way to do this? Currently, I am using a for loop and circshift, and I add at each iteration I add the new row to A, but probably this is highly inefficient.

Code (myMat is the matrix I want to shift):

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

@cholo14 thanks for listening to the comments and editing your question, this question is now well formed and really answerable... unfortunately I've not got my linear algebra hat on to work out the matrix operation which I think must exist and would give you this result simply. However, you could at least improve the code by pre-allocating the size of `A` using `zeros` (or similar), then assign `aux` to specific rows of `A` instead of appending each loop. — Wolfie, May 02 '18 at 13:10
Prelocating as Wolfie suggests will make a big difference. Also, the second line inside the loop replaces zeros with zeros? Don’t you get the same result without? — Cris Luengo, May 02 '18 at 13:29
@CrisLuengo preallocating is indeed one option to speed up the code. The second line inside the loop replaces zeros with zeros, but in general `myMat` could have elements different from 0 in the last columns, so I need to force them to 0 — cholo14, May 02 '18 at 13:39
@cholo14 so you are not really doing a circular shift, but a right shift with zero padding. It probably is more efficient to just copy over the relevant bit of the matrix. — Cris Luengo, May 02 '18 at 13:44

score 2 · Accepted Answer · answered May 02 '18 at 17:14

As you are probably aware, loops in Matlab are not so efficient. I know that the Mathworks keep saying this is no longer so with JIT compilation, but I haven't experienced the fast loops yet.

I put your method for constructiong the matrix A in a function:

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

Your matrix-operation looks a lot like a kronecker product, but the block-matrixces have overlapping column ranges so a direct kronecker product will not work. Instead, I constructed the following function:

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));

You can follow the algorithm by removing semicolons and looking at intermediate results.

The following main program runs your example, and can easily be modified to run similar examples:

% 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)))

"I haven't experienced the fast loops yet" -> check out my answer. :) -- sure, the difference is peanuts, but 10 years ago that loop would have been 100 times slower than your code. — Cris Luengo, May 02 '18 at 19:25

Cris Luengo · Answer 2 · 2018-05-02T19:31:15.483

I've taken Nathan's code (see his answer to this question), and added another possible implementation (replvector3).

My idea here stems from you not really needing a circular shift. You need to right-shift and add zeros to the left. If you start with a pre-allocated array (this is really where the big wins in time are for you, the rest is peanuts), then you already have the zeros. Now you just need to copy over myMat to the right locations.

These are the times I see (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

This is the full copy, copy-paste into an M-file and run:

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

Replicate vectors shifting them to the right

2 Answers2