CODEWITHSUNDEEP
arraysmatlabmatrixsparse-matrix
Michael A
Michael A

Reputation: 4613

"Index exceeds matrix dimensions" error with spdiags when creating band diagonal matrix?

I have this sample code that creates a band diagonal matrix

T = 6;
d1 = ones(T-2, 1);
d2 = 2*ones(T-1, 1);
d3 = 3*ones(T, 1);

f = sparse(diag(d1, -2) + diag(d2, -1) + diag(d3) + diag(d2, 1) + diag(d1, 2));

but I want to avoid creating the full TxT matrix, so I want to use spdiags, like this:

f2 = spdiags(d1, -2, T, T) + spdiags(d2, -1, T, T) + spdiags(d3, 0, T, T) + ...
    spdiags(d2, 1, T, T) + spdiags(d1, 2, T, T);

matlab tells me that "index exceeds matrix dimensions" and the problem comes from these commands

spdiags(d2, 1, T, T)
spdiags(d3, 2, T, T)

but these commands work normally:

spdiags(d1, -2, T, T)
spdiags(d2, -1, T, T)

What's going on here? The final matrix should look like the sparse form of this:

f =
     3     2     1     0     0     0
     2     3     2     1     0     0
     1     2     3     2     1     0
     0     1     2     3     2     1
     0     0     1     2     3     2
     0     0     0     1     2     3

Also these are sample matrices that Im using as examples only.

This code works too:

T = 6;
d1 = ones(T-2, 1);
d2 = 2*ones(T-1, 1);
d3 = 3*ones(T, 1);

f = sparse(diag(d1, -2) + diag(d2, -1) + diag(d3) + diag(d2, 1) + diag(d1, 2));

B = [[d1;0;0], [d2;0], d3, [0;d2], [0;0;d1]];
f2 = spdiags(B, -2:2, T, T);

Upvotes: 1

Views: 974

Answers (1)

Luis Mendo
Luis Mendo

Reputation: 112699

The documentation is not very clear. It looks like you need your d vector to be of length T, even if some values will be ignored (namely, the first values are ignored for positive diagonals, and the last are ignored for negative diagonals). But somehow Matlab only actually complains for positive diagonals; for negative diagonals it does accept shorter vectors.

So: use all d vectors of length T:

T = 6;
d1 = ones(T, 1);
d2 = 2*ones(T, 1);
d3 = 3*ones(T, 1);
f2 = spdiags(d1, -2, T, T) + spdiags(d2, -1, T, T) + spdiags(d3, 0, T, T) + ...
    spdiags(d2, 1, T, T) + spdiags(d1, 2, T, T);

By the way, you can build a matrix containing all d vectors as columns (now that they all have the same length) and call spdiags just once:

f2 = spdiags([d1 d2 d3 d2 d1], -2:2, T, T);

Upvotes: 1

Related Questions