6

I would like to repeat a vector A of length n on a diagonal m times to obtain a (n+m-1) x m matrix B. As an example let's say A = [a;b;c;d], m = 4. This should result in

B = 
[a 0 0 0;
 b a 0 0;
 c b a 0;
 d c b a;
 0 d c b;
 0 0 d c;
 0 0 0 d]

Any suggestions for a fast way to achieve this? blkdiag(repmat(A,1,m)) does not help me in this case as it creates a (n*m) x m matrix.

In the end I am actually just interested in the matrix product D of a third matrix C with B:

D=C*B

If you see another option to obtain D without having to generate B, I would appreciate it. But a solution for the problem above would make me very happy as well! n and m will be quite big by the way.

Thanks!

Robert Seifert
  • 25,078
  • 11
  • 68
  • 113
julian
  • 63
  • 5
  • 3
    That looks like a transposed toeplitz matrix. I believe that matlab has a function to generate toeplitz matrix. It might be of use to you. – mathematician1975 Jan 06 '14 at 21:05
  • @julian Is my mothod works for you? If m and n are high values, it may save you the time of creating matrix `B`, and multiply all the zeros over there... I just corious to know if it works... – Adiel Jan 07 '14 at 02:47

3 Answers3

7

Because @mathematician1975 is too lazy to write a proper answer.

Matlab has a function for this, called toeplitz

You would call it like this:

c=[1;2;3;4;0;0;0];
r=[0, 0, 0, 0];
toeplitz(c,r)

ans =

   1   0   0   0
   2   1   0   0
   3   2   1   0
   4   3   2   1
   0   4   3   2
   0   0   4   3
   0   0   0   4

You can play with the zeroes to shape the matrix the way you want it.

rubenvb
  • 74,642
  • 33
  • 187
  • 332
  • 1
    It wasn't laziness I just could not recall what the arguments were - I always find its better to write a comment than post an incorrect answer! You get +1 from me anyway – mathematician1975 Jan 07 '14 at 17:21
5

A clumsy, but generic one-liner

n = 3;      %number of elements in A;
m = 5;      %repetitions
A = (1:n);  

B = full( spdiags( repmat(A(:),1,m)' , 1-(1:n) , n+m-1, m) )

returns:

B =

     1     0     0     0     0
     2     1     0     0     0
     3     2     1     0     0
     0     3     2     1     0
     0     0     3     2     1
     0     0     0     3     2
     0     0     0     0     3

Alternatively an improved, generic version of rubenvb's solution

B = toeplitz( [A(:);zeros(m-1,1)] , zeros(1,m) )

in both cases A can be either a row or column vector.

The faster solution (factor 2x) is the first one with spdiags!


Edit: even clumsier, but up to 10x faster (it depends on n,m) than the toeplitz-approach:

B = reshape( [repmat([A(:);zeros(m,1)],m-1,1) ; A3(:)] ,[],m ) 
Community
  • 1
  • 1
Robert Seifert
  • 25,078
  • 11
  • 68
  • 113
2

Total solution, without matrix B, is to do convolution of each row of C with A. you can do it by for loop:

for k=1:size(C,1)
   D(k,:)=conv(C(k,:),A');
end
D=D(:,length(A)-1:end-length(A)+1);    % elliminate the convolution edges

I think it can be done also without the loop, by arrayfun:

k=1:size(C,1);
D=arrayfun(@(x) conv(C(x,:),A'), k);
D=D(:,length(A)-1:end-length(A)+1);    % elliminate the convolution edges
Adiel
  • 3,071
  • 15
  • 21
  • As far as I see it the convolution only gives the result as wanted above for symmetric `A` (`[a b b a]` for instance). My `A` is actually symmetric, so that's not a problem. Your for loop is a lot slower than the solution by rubenvb but actually a bit faster than your `arrayfun` version, which must be epanded to: `colNoise = cell2mat(arrayfun(@(x) conv(C(x,:), A', 'valid')', k, 'UniformOutput', 0))';` – julian Jan 08 '14 at 16:36
  • The `valid` saves your last line of eliminating the edges. The `UniformOutput` is necessary as the output is a vector and not a scalar. Even if I would change my data format and could therefore get rid of all "`'`" this is still slower. But still thank you very much for the suggestions – julian Jan 08 '14 at 16:41
  • And of course I meant `D` and not colNoise in the first comment... copy and paste error... – julian Jan 08 '14 at 16:43
  • This is interest, that it slower than building `B` as a different matrix and multiply it. Perhaps it becuase the `conv` function also pad zeros, so it doesn't save anything. I don't understand why it doesn't work for asymmetric `A`. And I know that usually `arrayfun`s are slower than for loops, but people like it sometimes, maybe it seems more "clean", or clear, or smart... :) – Adiel Jan 08 '14 at 18:36
  • `conv` mutltiplies elements from one side of a vector with an element from the other side of the other vector. What I suggested above was an algorithm for multiplying elements from the same side. For a symmetric vector this doesn't matter. I actually just came across 'conv2' which can solve my problem in a very short way: `D=conv2(C, A, 'valid')`. Interestingly even that takes longer than the method from thewaywewalk... – julian Jan 08 '14 at 18:44
  • You right, it should be written `D(k,:)=conv(C(k,:),flipud(A)');`. Or with `corr`, that similar to `conv`, but does not reverse the order of elements. Does `conv2` works as well? I think you need 1-D convolution on each row, it's not like 2-D convolution... – Adiel Jan 08 '14 at 19:02