How to extract an array of diagonal entries from a diagonal matrix in Julia?

Viewed 178

Given some diagonal matrix in Julia like A = Diagonal(rand(3,3)), is there any way I can create a one-dimensional array whose elements are the diagonal entries of this Diagonal matrix A?

2 Answers

There is diag(A, k::Integer=0) that returns the kth diagonal of a matrix A, as a vector.

julia> A = Diagonal(rand(3,3))
3×3 Diagonal{Float64, Vector{Float64}}:
 0.213159   ⋅         ⋅
  ⋅        0.034186   ⋅
  ⋅         ⋅        0.539693

julia> diag(A)
3-element Vector{Float64}:
 0.21315894297089488
 0.03418604147090787
 0.5396925608269262

In case speed of the essence, it would make sense not to allocate/copy any memory and access the diagonal of a matrix directly out of the stored matrix. Ideally with an abstract vector interface to this diagonal hiding the necessary memory location calculations.

This can be achieved using the Strided package. To install the package, the usual:

using Pkg
Pkg.add("Strided")

should work.

The example in the question would be solved with:

julia> using Strided

julia> A = rand(3,3)
3×3 Matrix{Float64}:
 0.627885   0.996657  0.5304
 0.290007   0.19639   0.0311003
 0.0931983  0.912228  0.227603

julia> D = Strided.UnsafeStridedView(pointer(A), (3,), (4,), 0)
3-element Strided.UnsafeStridedView{Float64, 1, Float64, typeof(identity)}:
 0.6278853842183714
 0.1963898010179035
 0.22760324615233707

This example has 3x3 matrix size hard-coded. The 3 and 4 in StridedView definition should be changed to matrix size (and size+1). It is important that original matrix is a dense normal layout matrix and not sparse or any other of the many matrix implementations in Julia.

The benefit of using this approach is about 10x in speed. Benchmarking yielded the following:

julia> @btime Strided.UnsafeStridedView(pointer($A), (3,), (4,), 0)
  1.885 ns (0 allocations: 0 bytes)
Related