Create sparse matrix with zig-zag diagonal

Is there a way to shift every column c by (c-1)*n rows down, where n is size(A,1) in the following example, without a loop and without excessive memory usage? So to get

   A = [
      a11 a12 a13 a14 ...
      a21 a22 a23 a24 ...
      ... ... ... ... ...
   ]

to read

    B = [
       a11   0   0   0 ...
       a21   0   0   0 ...
       ... ... ... ... ...
         0 a12   0   0 ...
         0 a22   0   0 ...
       ... ... ... ... ...
         0   0 a13   0 ...
         0   0 a23   0 ...
       ... ... ... ... ...
         0   0   0 a14 ...
         0   0   0 a24 ...
       ... ... ... ... ...
   ]

I tried diag and sparse with no luck so far. I need Matrix A in an operation X=J\A afterwards. It is crucial that the zero elements do not consume memory as A can have up to 20e6 columns and probably 10 rows.

Thanks in advance.

Answer 1

[Cedric]

You could use sparse(), building appropriate vectors of indices and values, e.g.:

 nCol = size(A, 2) ;
 nRow = nCol * size(A, 1) ;
 I = (1:nRow).' ;
 J = reshape(repmat(1:nCol, size(A, 1), 1), [], 1) ;
 B = sparse(I, J, A(:), nRow, nCol) ;

(that you can write in two lines actually, but I wanted to keep it clear)

[Florian]

Outstanding! It works!

I still need to figure out the magic behind this, though... sparse and especially reshape are still a closed book to me :(. I didn't manage to get behind the syntax, not to mention the usage of []...

[Cedric]

Well, one way to use sparse() is to provide 5 args: a vector or row indices, a vector of column indices, a vector of values, the number of rows and the number of columns. Here we build these "things", using one "trick" that simplifies a bit the problem: we create the vector of values using linear indexing on A (in this case, it means lining up columns of A and reading that as a vector). Example:

 >> A = [1, 2, 3; 4, 5, 6]
 A =
     1     2     3
     4     5     6
 >> A(:)
 ans =
     1
     4
     2
     5
     3
     6

Knowing that we will pass A(:) to sparse() as a vector of values, we have to build appropriate vectors of row and column indices..

Looking at what you depicted, row indices will just increase linearly, so we build I as 1:nRow, transposed so it is a column vector (same dim. as A(:)).

For column indices, it is a little more tricky as it will be something like

1 1 1 ... 2 2 2 ... etc, with the number of each integer being the number of rows of A. We can obtain this by creating a matrix

 1 2 3 4 ...
 1 2 3 4 ...
 . . . . 
 . . . .

that we ultimately index linearly. This matrix is actually the vertical repetition of (1:nCol) a number of times that matches the number of rows of A. You obtain that with repmat(). Example:

 >> repmat(1:4, 3, 1)
 ans =
     1     2     3     4
     1     2     3     4
     1     2     3     4

Now I could have written

 ...
 J = repmat(1:nCol, size(A, 1), 1)) ;
 B = sparse(I, J(:), A(:), nRow, nCol) ;

and avoid the reshape, but you would have had I given as a vector, and J as a matrix indexed linearly, which I thought would have brought confusion. The reshape transforms the output of repmat() into a column vector (last arg = 1). As reshape() is flexible, it allows to specify only one size (the other is left to an empty array), and computes the other. The reshape does therefore essentially what J(:) would have done.

[Florian]

Pretty cool, thank's a lot!

Answer 2

[bym]

s = spdiags(A',[0,-1],c*(n-1),c);

[Florian]

Thank you for your answer.

Unfortunately it didn't do the trick. Actually no specific value is assign to "c". In my question it was meant to be a variable to describe that the first column (c=1) should be moved (c-1)*n=0 rows down, the second column (c=2) should be moved by (c-2)*n=n rows down, the third column (c=3) should be moved (c-2)*n=2n rows down, etc.

From what I see, diags or spdiags places its elements only on true diagonals, shifted or not. In my case, however, the elements are not placed on a straight diagonal line but in a rather zig-zag way.