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

37 次查看（过去 30 天）
Boris Blagov
Boris Blagov 2023-3-21
评论： Boris Blagov 2023-3-21
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!

请先登录，再回答此问题。

采纳的回答

Matt J
Matt J 2023-3-21
编辑：Matt J 2023-3-21
Ran in:
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
  3 个评论
显示 1更早的评论
Matt J
Matt J 2023-3-21
You've refined your answer
Yes, Matrix(L)=-B(R) should be much faster than what I had initially.
Boris Blagov
Boris Blagov 2023-3-21
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!

请先登录，再进行评论。

更多回答（1 个）

John D'Errico
John D'Errico 2023-3-21
编辑：John D'Errico 2023-3-21
Ran in:
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?
  1 个评论
Boris Blagov
Boris Blagov 2023-3-21
编辑：Boris Blagov 2023-3-21
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Operating on Diagonal Matrices 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by