Populating off-diagonal blocks of diagonal matrix and making it sparse

Hi,

I've been thinking about what is the best way to populate a block diagonal matrix, make it sparse and then update those matrices. I've scoured the answers section and realize that everyone's matrices somehow appear to be special and performance is not acconted for in many of the answers.

Here is the matrix that I am trying to create:

B_j are (n x n) matrices. The whole matrix is T*n for an arbitrary T. "p" is given. EDIT: Both n and p are somewhat small, say n=15, p = 3. T is maybe 200.

The main diagonal is simple: kron(eye(n),T)

I know that the diag function can populate off-diagonals with diag(A,-1) for example. the blkdiag does not have such option.

What is important here is that the B matrices will be updated many, many times (say 20 000). So I want the matrix to be sparse in the end and I think it is good to avoid loops. So my question is what is the best way to populate and use the matrix such that performance is not impacted? I want to avoid reassigning memory for the matrix on every iteration. Should I try to create the indices and then only change the numbers and populate it this way?

e.g. for n=2 and T = 10 the I matrices would be assigned as

[1 2; 
11 12;
3 4;
15 16; ... ]

or are there better ways to do this?

Here is an example code that generates p=5 B (3x3) matrices (for n=3) both in a column format and in a 3D format, such that Bmat is (3x3x5) and B is 15x3 matrix stacking the transposed Bmat matrices in a column.

n = 3;     % # number of variables
p = 5;     % # number of lags

b0 = ones(n,1)*0.01;
B1  = -0.2 + (0.2+0.2).*rand(n,n); % off-diagonal elements are U(-0.2 0.2)
B1(1:(n+1):end) = 0 + (0.5-0.0).*rand(1,n); % diagonal elements are U(0, 0.5)

Bmat = zeros(n,n,p);        % This is the form y_t_(nx1) = B1 y_{t-1} + B2 y_{t-2} + ... + B_p y_{t-p} + e
Bmat(:,:,1) = B1;
B = zeros(n*p,n);           % This is the form y_t_(1xn) = [y_tm1' y_tm2',... y_tmp']*[B1';B2';...;Bp'] + e
B(1:n,:) = B1';
for ll = 2:p
    B_ll = randn(n,n)*((0.05/ll)^2);
    Bmat(:,:,ll) = B_ll;
    B(1+(ll-1)*n:n+(ll-1)*n,:) = B_ll';
end

Thank you for your time!

[ANTHONY REID]

Use cell arrays for block materices. Then use "cell2mat" to build actual matrix.

[Boris Blagov]

Thank you, but the answer from Matt J is superior in terms of speed and I need to do this tens of thousands of times.

Answer 1

[Matt J]

Perhaps as below. Note that the first 3 lines are just fixed precomputations.

n=2; p=3; T=5;
Matrix=kron(eye(n),eye(T)); 
[L,R]=makeIndices(n,p,T);


Bcell=arrayfun(@(z)ones(n)*z,2:p+1,'uni',0);  %create B matrices
B=cat(3,Bcell{:})

B = 
B(:,:,1) =

     2     2
     2     2


B(:,:,2) =

     3     3
     3     3


B(:,:,3) =

     4     4
     4     4

Matrix(L)=-B(R)  %Update matrix

Matrix = 10×10

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

function [L,R]=makeIndices(n,p,T)

  
Icell=mat2cell(reshape(1:n^2*p,n,[]), n,n*ones(1,p) );

Icell=[zeros(n), Icell, zeros(n)];

idx=toeplitz(min(p+2,1:T), [1,repelem(p+2,T-1)]  );

 Q = cell2mat(Icell(idx));

 L=logical(Q);
 R=nonzeros(Q);

end

[Boris Blagov]

Hey, thank you for your answer!

Either I have been misunderstood or, more probably, I have worded my question poorly.

I have not been reading at all about sparse matrices lately. I have read all the posts I could find about populating the off-diagonal blocks of block matrices.

My matrix is not 3d. My matrix actually resembles a banded matrix with very few non-zero elements but my band is below the main diagonal. Forgive me if my terminology is off, I do not know if it classifies as a banded matrix if it is not symmetric.

I have included a picture in my post exactly how it is supposed to look, maybe it doesn't show. Here is a wording of the matrix I will use:

• have a full matrix H which is Tn*Tn
• I have T identity matrices on the block diagonal, which are each is nxn. "n" will typically be small number, but still something like, say 15 or 20.
• On the block diagonals below (the bands), I would have "p" number of B matrices, that are all nxn, p will be small, 3,4,5,6, max 7
• Everything else will be fixed 0s

Another way to try to summarize it: the main block diagonal is Identity matrices and I have "p" number of bands below with, respectively B1 band, B2 band until Bp band. p is small. Then everyhting else is 0. The matrix is square!

Alternatively If the picture in the original post doesn't show, here is a visualisation with the code editor (not code!):

H = [I    O O O O ...  O    O  O
    -B1   I O O O ...  O    O  O
    -B2  B1 O O O ...  O    O  O
      .   . . . . ...  .    .  .
    -Bp -Bp-1 . . ...  O    .  .
     O  -Bp   . . ...  I    O  O
     .    .   . . ... -B1   I  O
     0    0   . . ... -B2 -B1 I]_(Tn x Tn)
I is (nxn), B1 ... Bp are all (nxn)
p is say 3, n = 15, T = 180 (all example numbers), O are nxn matrices with zeros, 

Take the following example: for a dataset of T = 180, n=15 and p = 3,

• I have a full matrix H which is Tn*Tn so 900x900.
• Consider each "n column block", so each column that is 900x15. It has an nxn identity matrix on top, then I have three (p=3) 15x15 dense matrices (my B matrices in the picture above)
• That is in total 15+15*15*3 = 690 elements in each block of 900x15 = 13500 elements. In the last three 900x15 blocks the dense elements will be less than that (specificall in the last column I have only 15 out of 13500 elements that are non-zero), but let's use the "average column block" as a rough guide: 690/13500 ~ 5% will be non-zeros. Everything else will be 0.

My code will do the following: it will create the H matrix and then use it in a calculation. The B matrices will then change and I need to (1) construct it again or (2) use the fact that I know exactly where those B matrices go. That is why I thought that a sparse matrix should be a good idea for that. All other elements remain 0s. I don't have experience with sparse matrices, so that may not be wise, hence my question. However, the scientific paper that already uses the above matrix implements sparse matrices and claim they help a lot with the computation (but haven't shared their code yet, so I am coding this myself).

The only "3d" thing that I have mentioned above was that i have currently stored the B_j = {B1, B2, ... Bp} in a 3d matrix and as a stacked matrix [B1'; B2';...Bp]'. That doesn't matter at all for the H matrix, I can store the B matrices in any form it is required.

My question in general is how to construct such a matrix in the fastest way possible and then do so again and agian.

[Boris Blagov]

Hmm, that is interesting, I saw the toepliz trick here https://nl.mathworks.com/matlabcentral/answers/1654045-block-sub-diagonals-matrix but wasn*t sure how to integrate it.

You've refined your answer and I think it should be exactly what I am looking for, thank you!

If anyone has other ideas, I would love to compare them to this one for speed. I will keep you posted, thanks!

[Matt J]

You've refined your answer

Yes, Matrix(L)=-B(R) should be much faster than what I had initially.

[Boris Blagov]

I tried that in my code and have to come back and say this, the makeIndices function is genious! I would have never came up with it on my own. And I only need to call it once! Thanks @Matt J!

Answer 2

[John D'Errico]

There are several mistakes you seem to be making here, surprising if you have been reading a lot about sparse matrices lately.

You mention a 3-d matrix. Sorry. There is no 3-d sparse matrix capability in MATLAB.

Frequent updates to a sparse matrix? A seriously TERRIBLE idea, certainly if the number or location of the elements in the matrix will be changing. The worst thing to do is to delete an element (making it zero) or inserting a new non-zero element. That forces MATLAB to reorder the entire set of elements in memory. If you are doing this often, then it would be a good idea to reconsider using a sparse matrix.

Next, how sparse is your matrix? Typically, it is reasonable to use a sparse form for linear algebra if you have only a few non-zero elements per row. Perhaps 2 to 5 non-zeros is common. But if there are too many non-zeros, then anything you do with that matrix in the form of factorizations will cause serious computational problems. You may still be able to gain in some cases, perhaps if you are using iterative methods to solve equations, which will not require a matrix factorization.

A = sprand(2000,2000,0.002); % 0.2% sparse, so roughly only 4 non-zeros per row 
timeit(@() lu(A))

ans = 0.2420

B = full(A);
timeit(@() lu(B))

ans = 0.0576

So computing an LU factorization was actually 4 times SLOWER for the sparse matrix. Admittedly, this was not a banded matrix, where the sparse matrix would SOMETIMES shine more. But the matrix you show seems to have many bands.

How large is your matrix? In the example you show, the matrix is trivially small. Computations with full matrices are extremely fast, and sparsity often will not gain you.

Why are you hoping to use sparse matrices? If only half of the elements in the matrix are non-zero, then you barely save any memory. And the computations done with the full versions of those matrices will be considerably faster. For example:

A = sprand(1000,1000,0.5);
B = full(A);
whos A B

  Name         Size                Bytes  Class     Attributes

  A         1000x1000            6296360  double    sparse    
  B         1000x1000            8000000  double              

timeit(@() A*A)

ans = 0.0800

timeit(@() B*B)

ans = 0.0131

As you can see, the 50% full sparse matrix was barely a memory gain, and a multiply was 6 times slower with the sparse matrix.

So, are you sure you understand sparse matrices, and the reason for using them?

Answer 3

[Chris Turnes]

Here's another option to consider. Depending on your problem quantities, it may or may not be faster than Matt J's suggestion above -- the main problem is really the size of the sparse kernel you have to create.

Since you've mentioned that you're updating B frequently, I'll make the following assumptions:

• You have freedom to store B how you want, and so when we're looking at timing we can mostly consider the cost of assembling the matrix you want given some fixed B
• You want things to be "sparse" in the sense that this matrix will have a lot of zeros, but not the actual sparse storage MATLAB provides. I'm inferring this from the examples.

Since you have a Toeplitz structure, this begs to use convolution, especially since there are a lot of zeros here which can hit some optimized code paths. You can do this by setting up an N-D kernel and playing games with reshaping. Using the example from above:

n=2; p=3; T=5;
Bcell=arrayfun(@(z)ones(n)*z,2:p+1,'uni',0);  %create B matrices
B=cat(3,Bcell{:});

% Just time the creation of the matrix, assuming B will be known a priori
% so we don't want to time it here.
tic;
Matrix=kron(eye(n),eye(T)); 
[L,R]=makeIndices(n,p,T);
Matrix(L)=-B(R);  %Update matrix;
toc

Elapsed time is 0.016943 seconds.

% Generate B as a series of matrices.
Bset = ones(n,p+1,n,1).*(-(1:(p+1)));
% This is set up so:
%  Bset(:,1,:) = reshape(I, n, 1, n);
%  Bset(:,1+k,:) = -reshape(Bk, n, 1, n);
Bset(:,1,:) = reshape(eye(n), n, 1, n);
% Generate a kernel to convolve against.
kern = reshape(eye(T), 1, T, 1, T);
tic;
Matrix2 = convn(kern, Bset);
Matrix2 = reshape(Matrix2(1:n, 1:T, :, :), T*n, []);
toc;

Elapsed time is 0.005491 seconds.

assert(isequal(Matrix2, Matrix))

function [L,R]=makeIndices(n,p,T)
  
Icell=mat2cell(reshape(1:n^2*p,n,[]), n,n*ones(1,p) );
Icell=[zeros(n), Icell, zeros(n)];
idx=toeplitz(min(p+2,1:T), [1,repelem(p+2,T-1)]  );
 Q = cell2mat(Icell(idx));
 L=logical(Q);
 R=nonzeros(Q);
end

You can play around with scaling the the quantities up and down. The relative values of the quantities may make this a more or less appealing solution.

I'd also note that depending on what you actually want to do with this matrix, you may not actually need to explicitly construct it, as you may be able to do whatever manipulations you're trying to do via calls to convn.

Edit: I reworked the arrangment of B to help me keep the kernel small, which helps a lot with the scaling.