length(1:25e-6:0.6) = 0
is it 0 instead of 1?
Also it is important to note that after a certain point of shifting (col_dim < iteration) youβll end up with just a matrix of 0βs, unless I presume you are doing other operations to populate A in the loop.
a slight variant of your script that appears to be even faster (at least for arrays of this size)
begin
A[:, 1] .= 0
dropzeros!(A)
A[1:size(A, 1), vcat(2:size(A, 2), 1)]
end
Sorry, yes 0:25e-6:0.6.
Yes, I am doing operation later on A.
julia> @bprofile include("Code_withMatrixMul.jl")
8.476085 seconds (46.76 M allocations: 4.289 GiB, 9.84% gc time, 32.81% compilation time)
5.236861 seconds (38.98 M allocations: 3.900 GiB, 11.59% gc time)
4.828120 seconds (38.98 M allocations: 3.900 GiB, 11.54% gc time)
4.647583 seconds (38.98 M allocations: 3.900 GiB, 7.84% gc time)
BenchmarkTools.Trial: 1 sample with 1 evaluation.
Single result which took 5.780 s (6.58% GC) to evaluate,
with a memory estimate of 4.11 GiB, over 41170470 allocations.
julia> @bprofile include("Code_withMatrixPermute.jl")
8.411843 seconds (48.63 M allocations: 4.943 GiB, 10.01% gc time, 28.52% compilation time)
5.438983 seconds (40.90 M allocations: 4.558 GiB, 12.14% gc time)
5.436992 seconds (40.90 M allocations: 4.558 GiB, 11.77% gc time)
5.102955 seconds (40.90 M allocations: 4.558 GiB, 8.03% gc time)
BenchmarkTools.Trial: 1 sample with 1 evaluation.
Single result which took 6.177 s (7.08% GC) to evaluate,
with a memory estimate of 4.76 GiB, over 43090389 allocations.
BTW, I should take the average as the computation time or the one after the sentence βSingle result which tookβ
Run @bprofile on your function not the include method.
Nice!, I guess permute! doesnβt work as well as simple slicing. Looks like it is fast enough over repeated operations to beat the multiplication as well!
I didnβt quite understand what you have to do with this matrix, but it occurs to me that you could change the calculation algorithm a bit, using only a slice of the original matrix A [2: end], without making any shifts, and then take into account that the last column is all of zeros.
The exchange could be convenient. it depends on the operations you have to do with it.
Updated MWE:
using BenchmarkTools
using LinearAlgebra
using SparseArrays
using Random
function trie()
row_dim = 44
col_dim = 33
A_density = 0.2679063360881543
rng = MersenneTwister(1)
A = sprand(rng, row_dim, col_dim, A_density)
B = sprand(rng, row_dim, col_dim, A_density)
# Matrix mult
shift = spdiagm(-1 => ones(size(A, 2) - 1))
for _ in 1:24001
# Matrix mult
A .= A * shift .+ B
# Slicing
A[:, 1] .= 0
dropzeros!(A)
A[1:size(A, 1), vcat(2:size(A, 2), 1)]
A .+= B
end
end
Results in order of method:
julia> include("shift_sparse.jl")
BenchmarkTools.Trial: 16 samples with 1 evaluation.
Range (min β¦ max): 301.004 ms β¦ 324.754 ms β GC (min β¦ max): 4.78% β¦ 4.60%
Time (median): 317.756 ms β GC (median): 4.59%
Time (mean Β± Ο): 317.456 ms Β± 5.221 ms β GC (mean Β± Ο): 4.60% Β± 0.13%
β β β ββ βββ ββββ ββ β
βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ β
301 ms Histogram: frequency by time 325 ms <
Memory estimate: 1.13 GiB, allocs estimate: 144044.
julia> include("shift_sparse.jl")
BenchmarkTools.Trial: 28 samples with 1 evaluation.
Range (min β¦ max): 175.798 ms β¦ 184.465 ms β GC (min β¦ max): 5.95% β¦ 5.29%
Time (median): 179.539 ms β GC (median): 5.99%
Time (mean Β± Ο): 180.077 ms Β± 2.321 ms β GC (mean Β± Ο): 5.94% Β± 0.34%
β β β
βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ β
176 ms Histogram: frequency by time 184 ms <
Memory estimate: 509.46 MiB, allocs estimate: 168029.
Excuse me, this line should be:
A .= A[1:size(A, 1), vcat(2:size(A, 2), 1)]
right?
and if so, the speed is slower than multiplication.
I also did as below but gave slower speed:
A[:, 1] .= 0
dropzeros!(A)
A .= A[1:size(A, 1), vcat(2:size(A, 2), 1)] .+ B
Please post your results of bprofile so we can quantify.
A .= A[1:size(A, 1), vcat(2:size(A, 2), 1)]
Why would u say that?
If you donβt mind accessing internals, you can use this version
function shiftleft!(A::SparseMatrixCSC, n)
sz = size(A)
n < sz[2] || throw(DimensionMismatch("$(n) bigger than the number of columns $(size[2])"))
colnz = A.colptr[n+1] - 1
@view(A.colptr[begin:end-n]) .= @view(A.colptr[begin+n:end])
A.colptr[end-n+1:end] .= A.colptr[end]
A.colptr .-= colnz
deleteat!(A.rowval, 1:colnz)
deleteat!(A.nzval, 1:colnz)
A
end
shifleft(A::SparseMatrixCSC, n) = shiftleft!(copy(A), n)
Is in place (shifleft!) and maybe the fastest you can get.
@suavesito while a clever idea using your function turned out to be a slower than the slicing method.
Results in order of mult, slice, shiftleft
julia> include("shift_sparse.jl")
BenchmarkTools.Trial: 16 samples with 1 evaluation.
Range (min β¦ max): 281.928 ms β¦ 360.109 ms β GC (min β¦ max): 4.87% β¦ 4.85%
Time (median): 312.091 ms β GC (median): 4.79%
Time (mean Β± Ο): 316.609 ms Β± 25.090 ms β GC (mean Β± Ο): 4.92% Β± 0.52%
β β β β β β β β β β β β ββ β
βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ β
282 ms Histogram: frequency by time 360 ms <
Memory estimate: 1.13 GiB, allocs estimate: 144045.
julia> include("shift_sparse.jl")
BenchmarkTools.Trial: 23 samples with 1 evaluation.
Range (min β¦ max): 144.470 ms β¦ 301.354 ms β GC (min β¦ max): 7.42% β¦ 6.76%
Time (median): 234.378 ms β GC (median): 6.32%
Time (mean Β± Ο): 223.732 ms Β± 46.080 ms β GC (mean Β± Ο): 6.66% Β± 0.57%
ββ ββββ β β β β βββ β βββββ β ββ
βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ β
144 ms Histogram: frequency by time 301 ms <
Memory estimate: 509.46 MiB, allocs estimate: 168029.
julia> include("shift_sparse.jl")
BenchmarkTools.Trial: 18 samples with 1 evaluation.
Range (min β¦ max): 212.972 ms β¦ 322.766 ms β GC (min β¦ max): 4.63% β¦ 4.78%
Time (median): 281.343 ms β GC (median): 4.62%
Time (mean Β± Ο): 283.081 ms Β± 25.233 ms β GC (mean Β± Ο): 4.74% Β± 0.49%
β β β β βββ β ββ β ββ β
βββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ β
213 ms Histogram: frequency by time 323 ms <
Memory estimate: 487.23 MiB, allocs estimate: 96030.
According to my benchmarks, it looks like it is much faster for the in-place shiftleft!, but still x3 faster for the out-of-place shiftleft. Also, my approach looks like is almost not affected by the percentage of non-zero values, vs. the other that do are affected.
julia> using SparseArrays, LinearAlgebra, BenchmarkTools
julia> function shiftleft!(A::SparseMatrixCSC, n)
sz = size(A)
n <= sz[2] || throw(DimensionMismatch("$(n) bigger than the number of columns $(size[2])"))
colnz = A.colptr[n+1] - 1
@view(A.colptr[begin:end-n]) .= @view(A.colptr[begin+n:end])
A.colptr[end-n+1:end] .= A.colptr[end]
A.colptr .-= colnz
deleteat!(A.rowval, 1:colnz)
deleteat!(A.nzval, 1:colnz)
A
end
shiftleft! (generic function with 1 method)
julia> shiftleft(A, n) = shiftleft!(copy(A), n)
shiftleft (generic function with 1 method)
julia> shiftleft_mul(A, n) = A * spdiagm(-n => ones(Int, size(A, 2)-n))
shiftleft_mul (generic function with 1 method)
julia> function shiftleft_slide(A::SparseMatrixCSC, n)
A[:, begin:begin+n-1] .= zero(eltype(A))
dropzeros!(A)
A[:, vcat(begin+n:end, begin:begin+n-1)]
end
shiftleft_slide (generic function with 1 method)
Results in
julia> @benchmark shiftleft!(A, 10) setup=(A=sprand(100, 100, 0.1))
BenchmarkTools.Trial: 10000 samples with 803 evaluations.
Range (min β¦ max): 150.450 ns β¦ 1.189 ΞΌs β GC (min β¦ max): 0.00% β¦ 66.28%
Time (median): 164.403 ns β GC (median): 0.00%
Time (mean Β± Ο): 193.630 ns Β± 89.107 ns β GC (mean Β± Ο): 3.87% Β± 8.50%
βββββββββββββββββββ ββ β
ββββββββββββββββββββββββββββββββββββββββββββββββββ
βββ
β
ββββββ β
150 ns Histogram: log(frequency) by time 697 ns <
Memory estimate: 816 bytes, allocs estimate: 1.
julia> @benchmark shiftleft!(A, 10) setup=(A=sprand(100, 100, 0.9))
BenchmarkTools.Trial: 10000 samples with 800 evaluations.
Range (min β¦ max): 145.001 ns β¦ 2.450 ΞΌs β GC (min β¦ max): 0.00% β¦ 84.63%
Time (median): 167.564 ns β GC (median): 0.00%
Time (mean Β± Ο): 214.467 ns Β± 130.594 ns β GC (mean Β± Ο): 4.92% Β± 8.98%
βββββββββββββββββββ β
ββββββββββββββββββββββββββ
ββββββ
β
β
β
βββββββββ
β
β
ββ
β
β
βββββββ
ββββ β
145 ns Histogram: log(frequency) by time 860 ns <
Memory estimate: 816 bytes, allocs estimate: 1.
julia> @benchmark shiftleft(A, 10) setup=(A=sprand(100, 100, 0.1))
BenchmarkTools.Trial: 10000 samples with 9 evaluations.
Range (min β¦ max): 1.769 ΞΌs β¦ 297.860 ΞΌs β GC (min β¦ max): 0.00% β¦ 96.98%
Time (median): 2.857 ΞΌs β GC (median): 0.00%
Time (mean Β± Ο): 4.596 ΞΌs Β± 10.683 ΞΌs β GC (mean Β± Ο): 11.58% Β± 5.31%
ββββββ
ββββββββ β βββ β
βββββββββββββββββββββ
β
β
β
β
ββββββββββββββββββββββ
βββββ
β
β
ββ
ββ
β
β
1.77 ΞΌs Histogram: log(frequency) by time 20.2 ΞΌs <
Memory estimate: 15.80 KiB, allocs estimate: 4.
julia> @benchmark shiftleft(A, 10) setup=(A=sprand(100, 100, 0.9))
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min β¦ max): 11.614 ΞΌs β¦ 808.363 ΞΌs β GC (min β¦ max): 0.00% β¦ 94.41%
Time (median): 14.659 ΞΌs β GC (median): 0.00%
Time (mean Β± Ο): 18.201 ΞΌs Β± 35.230 ΞΌs β GC (mean Β± Ο): 10.51% Β± 5.41%
ββ
βββββββββββββ ββββββ β
ββββββββββββββββββββββββββββββββ
ββββββββ
βββββββββββ
βββββββββ β
11.6 ΞΌs Histogram: log(frequency) by time 48.1 ΞΌs <
Memory estimate: 140.52 KiB, allocs estimate: 6.
julia> @benchmark shiftleft_mul(A, 10) setup=(A=sprand(100, 100, 0.1))
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min β¦ max): 17.919 ΞΌs β¦ 1.694 ms β GC (min β¦ max): 0.00% β¦ 95.68%
Time (median): 22.169 ΞΌs β GC (median): 0.00%
Time (mean Β± Ο): 25.795 ΞΌs Β± 36.031 ΞΌs β GC (mean Β± Ο): 3.66% Β± 2.65%
β
ββββββ
ββββ
ββββββββββ
βββββββββββββββββββββββββββββββββββββββββββββ β
17.9 ΞΌs Histogram: frequency by time 41.7 ΞΌs <
Memory estimate: 25.09 KiB, allocs estimate: 19.
julia> @benchmark shiftleft_mul(A, 10) setup=(A=sprand(100, 100, 0.9))
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min β¦ max): 51.419 ΞΌs β¦ 938.537 ΞΌs β GC (min β¦ max): 0.00% β¦ 86.56%
Time (median): 55.294 ΞΌs β GC (median): 0.00%
Time (mean Β± Ο): 62.974 ΞΌs Β± 43.053 ΞΌs β GC (mean Β± Ο): 4.06% Β± 5.75%
βββββ
β
ββββββββββββββββββββββββ β
βββββββββββββββββββββββββββββββββββββββββ
ββ
β
β
ββββ
βββββββ
ββββ β
51.4 ΞΌs Histogram: log(frequency) by time 127 ΞΌs <
Memory estimate: 246.69 KiB, allocs estimate: 23.
julia> @benchmark shiftleft_slide(A, 10) setup=(A=sprand(100, 100, 0.1))
BenchmarkTools.Trial: 10000 samples with 7 evaluations.
Range (min β¦ max): 4.884 ΞΌs β¦ 230.295 ΞΌs β GC (min β¦ max): 0.00% β¦ 88.43%
Time (median): 5.887 ΞΌs β GC (median): 0.00%
Time (mean Β± Ο): 6.586 ΞΌs Β± 8.344 ΞΌs β GC (mean Β± Ο): 5.40% Β± 4.16%
βββββββββ
β
βββββββββββββββββββ
β
βββββββββββββββββββββββββββββββββββββββ β
4.88 ΞΌs Histogram: frequency by time 9.59 ΞΌs <
Memory estimate: 14.62 KiB, allocs estimate: 4.
julia> @benchmark shiftleft_slide(A, 10) setup=(A=sprand(100, 100, 0.9))
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min β¦ max): 27.227 ΞΌs β¦ 927.842 ΞΌs β GC (min β¦ max): 0.00% β¦ 90.21%
Time (median): 30.218 ΞΌs β GC (median): 0.00%
Time (mean Β± Ο): 37.314 ΞΌs Β± 37.987 ΞΌs β GC (mean Β± Ο): 4.89% Β± 4.92%
βββββ
β
β
ββββββββββββββββββββββββ ββ β
ββββββββββββββββββββββββββββββββββββββββββββββββββββ
β
βββββ
β
β
β
27.2 ΞΌs Histogram: log(frequency) by time 79 ΞΌs <
Memory estimate: 126.72 KiB, allocs estimate: 6.
I had also tried to go down this path, but I had problems with the modification of the structure which is immutable.
Now (seeing @suavesito 's example) I understand better how to do it.
Here is a version (of @suavesito 's solution) that saves some allocations.
using SparseArrays, ShiftedArrays
function spshiftl(A,n)
cnz=A.colptr[n+1]-1
SparseMatrixCSC{Float64,Int64}(size(A)..., ShiftedArray(A.colptr,-n, default=A.colptr[end]).-cnz,
A.rowval[cnz+1:end], A.nzval[cnz+1:end])
end
for example, if we were to multiply spshiftl (A, 1) by a matrix B, we could use the following relations
# A(mXn) = A'| A'' = A[:,2:n]|spzeros(n)
# B(nXk) = B' | B'' = B[1:n-1,:]|B[n,:]
# C[i, j] = dot(A[i, :], B[:, j]) = dot(Aβ[i, :], Bβ[:, j]) + 0*B[n,j] => A*B == A'*B'
I ran my benchmarks with your tests and I got similar results that agree with them.
So I think this is the fastest method we have come up so far for shifting a matrix, but I donβt think it will be faster for @Amro 's workload.
Notice that @Amro is shifting the same matrix length(0:25e-6:0.6) times and that he is doing an operation on A. I tried to recreate this scenario with an updated MWE here: What is the fastest method to shift matrix? - #27 by acxz
And under this test, shiftleft! falls short. Iβm not sure why, but I guess doing an operation ends up being slower on the shiftleft!'d matrix than the slideβd matrix? I really donβt understand why tho.
Do you mind running your benchmarks with the MWE linked?
I think this is a really useful observation!
Because:
# Slicing
A[1:size(A, 1), vcat(2:size(A, 2), 1)]
is not updating A. This is what I got.
Good one to save allocations.
# A(mXn) = A'| A'' = A[:,2:n]|spzeros(n)
# B(nXk) = B' | B'' = B[1:n-1,:]|B[n,:]
# C[i, j] = dot(A[i, :], B[:, j]) = dot(Aβ[i, :], Bβ[:, j]) + 0*B[n,j] => A*B == A'*B'
I am not sure if I got it well here.
Do you have any doubts about the correctness of the thesis or is it not clear what the thesis affirms?
It is not clear for me what the thesis affirms.