Reputation: 1681
In Julia, what would be an efficient way of turning the diagonal of a matrix to zero?
Upvotes: 4
Views: 1230
Reputation: 2332
Performance wise, simple loop is faster (and more explicit, but it is taste dependent)
julia> @btime foreach(i -> $m[i, i] = 0, 1:20)
11.881 ns (0 allocations: 0 bytes)
julia> @btime setindex!.(Ref($m), 0.0, 1:20, 1:20);
50.923 ns (1 allocation: 240 bytes)
And it is faster then diagind
version, but not by much
julia> m = rand(1000, 1000);
julia> @btime foreach(i -> $m[i, i] = 0.0, 1:1000)
1.456 μs (0 allocations: 0 bytes)
julia> @btime foreach(i -> @inbounds($m[i, i] = 0.0), 1:1000)
1.338 μs (0 allocations: 0 bytes)
julia> @btime $m[diagind($m)] .= 0.0;
1.495 μs (2 allocations: 80 bytes)
Upvotes: 2
Reputation: 1450
Here is a more general way how to performantly use setindex!
on an Array by accessing custom indices:
Using an array for in array indexing
Linear indexing is best performing, that's why diagind
runs better than the Cartesian Indices.
Upvotes: 0
Reputation: 130
Przemyslaw Szufel's solutions bench-marked for 1_000 x 1_000 matrix size show that diagind
performs best:
julia> @btime setindex!.(Ref($m), 0.0, 1:1_000, 1:1_000);
2.222 μs (1 allocation: 7.94 KiB)
julia> @btime $m[diagind($m)] .= 0.0;
1.280 μs (2 allocations: 80 bytes)
Upvotes: 0
Reputation: 42194
Assuming m
is your matrix of size N x N
this could be done as:
setindex!.(Ref(m), 0.0, 1:N, 1:N)
Another option:
using LinearAlgebra
m[diagind(m)] .= 0.0
And some performance tests:
julia> using LinearAlgebra, BenchmarkTools
julia> m=rand(20,20);
julia> @btime setindex!.(Ref($m), 0.0, 1:20, 1:20);
55.533 ns (1 allocation: 240 bytes)
julia> @btime $m[diagind($m)] .= 0.0;
75.386 ns (2 allocations: 80 bytes)
Upvotes: 5