Update a sparse matrix efficiently

Question

mafoEV 2020-3-5

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/509170-update-a-sparse-matrix-efficiently

编辑： Bruno Luong 2023-1-26

Dear MATLAB community,

I have a large sparse matrix (~ 100,000 by 100,000, with about 100,000 non-zero elements), which I need to update regularly (at each solver step).

So far, I've done it by creating its elements first (rows, columns and values) as column vectors, then assembling them using the sparse function:

A = sparse(iA, jA, vA, m, n, nnzA);

However this line is now taking 13% of my total solve time (called 355 times, taking about 0.29s each time).

The matrix structure does not change from one iteration to the next, only about half of the values need updating. Is there a more efficient way to do it? I haven't been able to find any solution so far looking at the forums.

I'm using Matlab 2019b Update 3 (9.7.0.1261785).

Many thanks.

2 个评论
显示无隐藏无

Bruno Luong 2020-8-25

Are you updated only non-zero elements? only zero elements? Not known?

In three cases, the required time are differents.

Joel Lynch 2023-1-14

Are you absolutely certain you need to regenerate A every iteration? A basic strategy in nonlinear systems is to transform the Ax=b problem to solve for perturbations in x, allowing you to recycle old A matrices (Jacobians) until they no longer converge. More importantly, you can also recycle decompositions of A as well, which is where the real time-savings are, especially if the changes in A are relatively small.

请先登录，再进行评论。

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

Answer 1

sjhstone 2020-8-25

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/509170-update-a-sparse-matrix-efficiently#answer_484745

You can refer to the following blog post. It seems that you have reached the optimal efficiency.

https://blogs.mathworks.com/loren/2007/03/01/creating-sparse-finite-element-matrices-in-matlab/

If the sparsity pattern does not change, and there are nonzero entries on the last column and the last row, I think it might be redundant to pass arguments nnzA. Because MATLAB seems to use MEX routine, we do need to do numerical experiments to see whether not passing them into sparse can accelerate your code.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Christine Tobler 2023-1-16

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/509170-update-a-sparse-matrix-efficiently#answer_1149995

在 MATLAB Online 中打开

The fastest way to construct a sparse matrix will be when the inputs are sorted, first by columns and then by rows. You can verify that your inputs iA, jA match this by calling find on the resulting matrix A: The first two outputs of find(A) should match the order of indices in iA and jA.

I tried running sparse for just a random set of indices of the size you mentioned. As you can see, using sorted indices gives a bit of speed-up:

iA = randi(1e5, 1e5, 1);
jA = randi(1e5, 1e5, 1);
vA = randi(1e5, 1e5, 1);
tic; A = sparse(iA, jA, vA, 1e5, 1e5); toc
Elapsed time is 0.006209 seconds.
[siA, sjA, svA] = find(A);
tic; sA = sparse(siA, sjA, svA, 1e5, 1e5); toc
Elapsed time is 0.002676 seconds.
isequal(A, sA)
ans = logical
   1

Both my calls here are faster than what you mentioned (0.006 seconds instead of 0.29 seconds), which might come down to the machine used.

2 个评论
显示无隐藏无

Joel Lynch 2023-1-25

A very minor comment, but sparse does sum co-located points in iA/jA, so the second snippet using find(A) may represent a smaller reduced dataset which skews the timing. This is not noticable for the sizes used here, but if you were to try this with a different set of I/J/V data, you may get inaccurate timing results.

Bruno Luong 2023-1-26

编辑：Bruno Luong 2023-1-26

在 MATLAB Online 中打开

@Joel Lynch Christine code returns sIA, sjA without duplication and svA accumulated of vA.

The purpose to is show that calling sparse with sorted indexes is faster.

This code might be better (fairer) since it is also accumulate with eventually duplication for both cases:

iA = randi(1e5, 1e5, 1);
jA = randi(1e5, 1e5, 1);
vA = randi(1e5, 1e5, 1);
tic; A = sparse(iA, jA, vA, 1e5, 1e5); toc
Elapsed time is 0.005231 seconds.
[~, is] = sortrows([jA iA]);
siA = iA(is);
sjA = jA(is);
svA = vA(is);
tic; sA = sparse(siA, sjA, svA, 1e5, 1e5); toc
Elapsed time is 0.002319 seconds.