Yes, the function of R is to rotate F to left.
julia> R
8×8 SparseMatrixCSC{Float64, Int64} with 7 stored entries:
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ 1.0 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ 1.0 ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ 1.0 ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ 1.0 ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 1.0 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 1.0 ⋅
Ah I see now, can you specify how Ik changes?
If Ik .+= i is true than you don’t need a for loop, even if it is in the time-domain. hmm actually not entirely sure
As for reducing the number of times you need to create sparse matrices you can operate on it in place like how it was suggested earlier in this discussion:
for iter = 1:100
@. F.nzval = iter * F.nzval + rand() # update nonzero entries in-place
end
Changing Ik relates to many parameters in the code (long equation with several parameters).
I see, then yeah the for loop stays haha
Yes, but after introducing the term F*R, I am not able to determine how to make in-place update.
So now you have two matrices that you want to update in place, one that is the final F matrix you want to return (which gets rotated every timestep), and one that only gets its values updated.
Yes. I am not sure if dividing it into two lines can help. However, I am still not able to make the correct action.
So you would introduce another matrix, say O, that keeps track of the nonrotated portions to compute at each timestep like:
F .= F * R .+ O
And you would update O’s non zero values using the .nzsval field like so: O.nzval .+= i
And the original value of O? O = sparse(1:6, Nh, Ik) which is defined only once outside the for loop.
An optimized version (with my suggestions up above) would look something like this:
R = sparse(2:8, 1:7, 1.0, 8, 8)
Nh = [6, 6, 6, 8, 7, 7] # fixed
Ik = [0, 1, -4, -0, 4, 1.1]
O = sparse(1:6, Nh, Ik)
F = sparse(
[
3 3 3 4 4 4 0.0 0.0
1 2 1 7 9 3 0.0 0.0
7 8 1 4 8 8 0.0 0.0
1 1 2 8 5 1 0.0 0.0
5 6 1 4 2 5 4.6 0.0
1 6 8 1 1 4 6.8 0.0
],
)
for i = 1:100
O.nzval .+= i
F .= F * R .+ O
end
Profiling the code:
julia> include("sparse.jl") # Original code
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min … max): 127.978 μs … 8.009 ms ┊ GC (min … max): 0.00% … 91.08%
Time (median): 282.130 μs ┊ GC (median): 0.00%
Time (mean ± σ): 330.082 μs ± 556.243 μs ┊ GC (mean ± σ): 17.04% ± 9.76%
█ ▆
███▇▄▃▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▂▂▂▁▂▂ ▂
128 μs Histogram: frequency by time 4.84 ms <
Memory estimate: 489.38 KiB, allocs estimate: 1910.
julia> include("sparse.jl") # Improved code
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min … max): 94.017 μs … 1.798 ms ┊ GC (min … max): 0.00% … 82.30%
Time (median): 131.914 μs ┊ GC (median): 0.00%
Time (mean ± σ): 148.972 μs ± 102.671 μs ┊ GC (mean ± σ): 5.10% ± 7.64%
▇▇▅▅▅██▆▅▄▃▂▂▁ ▁▁▁▂▂▂▂▂▂▂▁▁ ▃
███████████████████████████████▇▇▇▆▆▆▆▆▁▅▅▅▅▅▆▃▅▅▅▅▅▅▄▅▅▅▁▁▄▄ █
94 μs Histogram: log(frequency) by time 508 μs <
Memory estimate: 377.45 KiB, allocs estimate: 729.
Also is there any particular reason why you are using sparse matrices? For the dimension size you have given us, dense matrices might give a faster result. At the very least for the F matrix, which seems to be a dense matrix instead of a sparse matrix.
That also gives me a considerable speed up in the computation.
Thank you very much for your efforts. It worked in this example, but it did not in my real code. I will check where is the problem.
In may real code, R is 33x33 and F is 102x33. O is 102x33
You can update your example with the required dimensions, by making fake matrices with random values. For example, the below works and has the same dimensions as your original problem.
row_dim = 102
col_dim = 33
R = sparse(2:col_dim, 1:col_dim-1, 1, col_dim, col_dim)
size = row_dim * col_dim
O_density = 6 / size
F_density = 38 / size
O = sprand(row_dim, col_dim, O_density)
F = sprand(row_dim, col_dim, F_density)
for i = 1:100
O.nzval .+= i
F .= F * R .+ O
end
Note, I’m using density values based on the nonzeros values in your above examples.
What error are you getting?
These are probably too small for sparse matrices to be worthwhile.
If you have banded structure, it totally will be worth it, but CSC has a minimum size to be worthwhile of roughly 128x128
I believe that the density is the key, even if you have a low enough matrix density you can still use sparse matrices at a more performant level than dense matrices. Even below the 128x128 range.
In this problem, if the number of nonzeros is the same as in the example code provided and in larger dimensions case, then using sparse matrices is still faster at newer dimensions (102 x 33).
However, if you use the same density (notice that in the example code F is quite dense) than using dense matrices is faster.
@Oscar_Smith @stevengj @Jeff_Emanuel @acxz Thank you very much for your efforts. Could you help me here please?
I found the causing problem. In my real code, I calculated O then update F2. However, as can be seen O is not as S. In O.nzval .= Ik; the order of filling is based on column-based while it is row-based in S . Consequently, F1 is not as F2.
Is there a way to add indices to O.nzval to correct the filling in O.nzval .= Ik;?
R = sparse([zeros(1,8); Matrix(1.0I, 8-1, 8-1) zeros(8-1,1)]);
Nh = [6, 6, 6, 8, 7, 7]; # fixed
Ik = [7, 1, -4, 8, 4, 1];
F = sparse([3 3 3 4 4 4 0.0 0.0;
1 2 1 7 9 3 0.0 0.0;
7 8 1 4 8 8 0.0 0.0;
1 1 2 8 5 1 0.0 0.0;
5 6 1 4 2 5 4.6 0.0;
1 6 8 1 1 4 6.8 0.0]);
F1 = copy(F);
F2 = copy(F);
O = sparse(1:6, Nh, Ik)
Ik = [10, 20, 30, 40, 50, 60];
O.nzval .= Ik;
S = sparse(1:6, Nh, Ik);
F1 .= F1 * R .+ O;
F2 .= F2 * R .+ S;
julia> O
6×8 SparseMatrixCSC{Int64, Int64} with 6 stored entries:
⋅ ⋅ ⋅ ⋅ ⋅ 10 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 20 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 30 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 60
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 40 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 50 ⋅
julia> S
6×8 SparseMatrixCSC{Int64, Int64} with 6 stored entries:
⋅ ⋅ ⋅ ⋅ ⋅ 10 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 20 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ 30 ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 40
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 50 ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ 60 ⋅
julia> F1
6×8 SparseMatrixCSC{Float64, Int64} with 38 stored entries:
3.0 3.0 4.0 4.0 4.0 10.0 ⋅ ⋅
2.0 1.0 7.0 9.0 3.0 20.0 ⋅ ⋅
8.0 1.0 4.0 8.0 8.0 30.0 ⋅ ⋅
1.0 2.0 8.0 5.0 1.0 ⋅ ⋅ 60.0
6.0 1.0 4.0 2.0 5.0 4.6 40.0 ⋅
6.0 8.0 1.0 1.0 4.0 6.8 50.0 ⋅
jjulia> F2
6×8 SparseMatrixCSC{Float64, Int64} with 38 stored entries:
3.0 3.0 4.0 4.0 4.0 10.0 ⋅ ⋅
2.0 1.0 7.0 9.0 3.0 20.0 ⋅ ⋅
8.0 1.0 4.0 8.0 8.0 30.0 ⋅ ⋅
1.0 2.0 8.0 5.0 1.0 ⋅ ⋅ 40.0
6.0 1.0 4.0 2.0 5.0 4.6 50.0 ⋅
6.0 8.0 1.0 1.0 4.0 6.8 60.0 ⋅
julia> F1==F2
false
The issue is the ordering:
Change this to Nh = [6, 6, 6, 7, 7, 8] When using .nzval it returns the nonzeros values of the array based on the order of the row and col indices, not based on the original order you passed in.
Instead of using the .nzval field directly, it may be safer to use the nzrange function and iterate. See: Sparse Arrays · The Julia Language
Thanks for your reply!
I made Nh = sort([6, 6, 6, 8, 7, 7]); so it worked in this example but failed in others. I am not sure if there is another way?
Thanks for this. However, using this iteration in my time-domain for-loop is expensive. Is there another way?
Using sort is incorrect, sort will not generalize.
Is there another way?
If you want to use the .nzval field (or the nonzeros function) directly you need to make sure that Ik is in the order of how it is represented internally, which is in the CSC ordering, instead of the input CSR ordering you used to create the sparse matrix. (I think this is true, you can make sure by checking the source code: https://github.com/JuliaLang/julia/blob/742b9abb4dd4621b667ec5bb3434b8b3602f96fd/stdlib/SparseArrays/src/sparsematrix.jl or just try assigning the nzval with CSC ordering and if it works it work!)
using this iteration in my time-domain
for-loopis expensive.
Can you be more quantitative? What are the time differences? I would be surprised if the performance is drastically hurt bottled.
Can you also comment on the density of the matrices that you are working with, specifically the F and O matrices? This information along with the dimensions will help us understand which is the appropriate data structure to use.