Populating off-diagonal blocks of diagonal matrix and making it sparse
27 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1331670/image.png)
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!
0 commentaires
Réponse acceptée
Matt J
le 21 Mar 2023
Modifié(e) : Matt J
le 21 Mar 2023
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{:})
Matrix(L)=-B(R) %Update 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
3 commentaires
Matt J
le 21 Mar 2023
You've refined your answer
Yes, Matrix(L)=-B(R) should be much faster than what I had initially.
Plus de réponses (1)
John D'Errico
le 21 Mar 2023
Modifié(e) : John D'Errico
le 21 Mar 2023
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))
B = full(A);
timeit(@() lu(B))
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
timeit(@() A*A)
timeit(@() B*B)
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 commentaire
Voir également
Catégories
En savoir plus sur Sparse Matrices dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!