Create sparse matrix with zig-zag diagonal

5 views (last 30 days)
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.

Accepted Answer

Cedric
Cedric on 3 Feb 2013
Edited: Cedric on 3 Feb 2013
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)
  3 Comments
Cedric
Cedric on 3 Feb 2013
Edited: Cedric on 4 Feb 2013
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.

Sign in to comment.

More Answers (1)

bym
bym on 3 Feb 2013
s = spdiags(A',[0,-1],c*(n-1),c);
  1 Comment
Florian
Florian on 3 Feb 2013
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.

Sign in to comment.

Categories

Find more on Creating and Concatenating Matrices in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!