iterator for diagonal of a matrix ?

[jlapeyre]

This is so obvious that I guess it exists elsewhere, but I don't know about it. But, if not, maybe it belongs (after some modifications) in LinearAlgebra ? ... Or some other package ?

A prototype iterator over the diagonal of a matrix is in this package.

It's called diagonal. It's like diag but does very little allocation.

julia> N = 1000;

julia> m = rand(N, N);

julia> @btime sum(diag($m));         # copy
  1.718 μs (1 allocation: 7.94 KiB)

julia> @btime sum(diagonal($m));    # no copy
  936.793 ns (1 allocation: 32 bytes)

Off diagonals are not yet supported.

[yuyichao]

If you just want a iterator diagonal(m) = (m[i, i] for i in 1:checksquare(m))?

[jlapeyre]

diagonal(m) = (m[i, i] for i in 1:checksquare(m))

Yes. I'm thinking of something like that. But,

1. For users coming from, say, MATLAB, you can just say that diagonal(m) gets you the diagonal (maybe with some caveats in the fine print)

2. The prototype I pointed to is constructed slightly faster. It supports getindex, setindex!, ==, etc. (maybe you get some of that by subtyping on AbstractVector, I did not do that.) If it's in a package, it can be modified, improved, debugged, etc.

[jlapeyre]

I see (JuliaLang/LinearAlgebra.jl#542) that formerly diag sometimes returned a copy and sometimes a view. It was changed so that it always returns a copy. So, diagonal would be complementary; it always returns a view.

[ararslan]

To always return a view, you could just use view(A, diagind(A, k)) (for the kth diagonal).

julia> A = rand(5, 5);

julia> diagonal(A::AbstractMatrix, k::Integer=0) = view(A, diagind(A, k))
diagonal (generic function with 2 methods)

julia> diagonal(A)
5-element view(::Array{Float64,1}, 1:6:25) with eltype Float64:
 0.46312443951994453
 0.9667117622181021 
 0.08385355771707848
 0.2588063486420449 
 0.7049707946935773 

Since it's just a SubArray, it supports all of those things like indexing, equality, mutating elements, etc.

[jlapeyre]

That may be the best choice. I put a note about it in the README. I'll probably replace the existing code with your one liner.

UPDATE: Thanks @ararslan . Your submission is the winner

[mbauman]

So I guess the actionable question is if we want a view version of diag in base Julia or the LinearAlgebra stdlib?

[ararslan]

I almost kind of think that diag should always return a view, then people can call copy if they need a copy. That is potentially disruptive and breaking though.

[KristofferC]

I think the lazy adjoint has shown us that we should not underestimate how many new methods that need to be implemented for lazy wrappers to not hit the slow AbstractArray fallback. We still have multiple regressions from that, so I would personally think that holding off on more lazy wrappers for a while is perhaps for the best?

[jlapeyre]

I almost kind of think that diag should always return a view, then people can call copy if they need a copy. That is potentially disruptive and breaking though.

It also brings up again the question of how non-coders react to this. Julia the language vs. Julia the better MATLAB-like tool. I think (have no statistics) that many are upset with the disappearance of eye. It's too late now, but pre v1.0, I would have said that if eye goes, so does diag.

Having just copy(diagonal(m)) or copy(eye(n)) (or Eye, or whatever) is obviously more orthogonal. Maybe the question is whether 1) with enough documentation, examples, culture change, etc., everyone finds copy(x) natural. (umm, or another verb, but that's another question) or 2) One could have a package that gives people the matlab (or python or whatever) -like things that they miss. A high-quality package that people will trust.

The language is complicated enough that there is an argument for the core language and stdlib to have a set of composable pieces that is minimal, yet usable. For instance, Matrix{T}(I, n, n) is actually very confusing to some casual (and not so casual) coders. Maybe Eye should go in stdlib.

I don't have the breadth of perspective that others in this thread do. This summer, I heard Jeff say "whatever the question, putting it in base is not the solution".

Given all this, my answer to @mbauman's question is put it in stdlib.

And in answer to to @KristofferC, yeah. Leave it in a package for a while, and let things shake out. Then, maybe stdlib.

[jlapeyre]

There are 19 methods (just in base and stdlib) for diag. The default method, for AbstractMatrix, calls diagind. There are only generic methods for diagind. The current implementation of diagonal, the iterator, calls diagind. So copy(diagonal(m)) is currently not a generic replacement.

[mbauman]

One thing I've considered in the past is that this could (or even should) be expressed via indexing. I often just use A[1:size(A, 1)+1:end] for this computation without giving it a first-class name. But since I is an "array-like" object of booleans that expands to the correct size, we could even implement this as the logical indexing expression A[I] — that is, A[LinearAlgebra.I], not some generic placeholder I. That makes the view/copy dichotomy crystal clear. Easy to implement, too:

julia> Base.to_index(A::AbstractArray, ::UniformScaling{Bool}) = diagind(A, 0)

julia> A = reshape(1:25, 5, 5);

julia> A[I]
5-element Array{Int64,1}:
  1
  7
 13
 19
 25

julia> @view A[I]
5-element view(::UnitRange{Int64}, 1:6:25) with eltype Int64:
  1
  7
 13
 19
 25

If this is a common thing folks reach for, we should probably specialize diagind(::SparseMatrix, …) to return a range of appropriate cartesian indices — it'd take a new type, but that's simple enough.

[ararslan]

I like that, but then you don't have a way to specify other diagonals. (Though I guess you could just use good ol' diagind(A, k) there.)

[cbartondock]

I think the fastest way in 2021 to get the diagonal of m is m[CartesianIndex.(axes(m,1),axes(m,2))]

[jishnub]

Ideally, the indices may be a StepRangeLen of CartesianIndexes if #50457 is approved. Currently, this allocates a Vector:

julia> m = Diagonal(1:4)
4×4 Diagonal{Int64, UnitRange{Int64}}:
 1  ⋅  ⋅  ⋅
 ⋅  2  ⋅  ⋅
 ⋅  ⋅  3  ⋅
 ⋅  ⋅  ⋅  4

julia> CartesianIndex.(axes(m)...)
4-element Vector{CartesianIndex{2}}:
 CartesianIndex(1, 1)
 CartesianIndex(2, 2)
 CartesianIndex(3, 3)
 CartesianIndex(4, 4)