Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add diagview to obtain a view along a diagonal #56175

Open
wants to merge 1 commit into
base: master
Choose a base branch
from

Conversation

jishnub
Copy link
Contributor

@jishnub jishnub commented Oct 15, 2024

A function to obtain a view of a diagonal of a matrix is useful, and this is clearly being used widely within LinearAlgebra.

The implementation here iterates according to the IndexStyle of the array:

julia> using LinearAlgebra

julia> A = reshape(1:9, 3, 3)
3×3 reshape(::UnitRange{Int64}, 3, 3) with eltype Int64:
 1  4  7
 2  5  8
 3  6  9

julia> diagview(A,1)
2-element view(::UnitRange{Int64}, 4:4:8) with eltype Int64:
 4
 8

julia> T = Tridiagonal(1:3, 3:6, 4:6)
4×4 Tridiagonal{Int64, UnitRange{Int64}}:
 3  4    
 1  4  5  
   2  5  6
     3  6

julia> diagview(T,1)
3-element view(::Tridiagonal{Int64, UnitRange{Int64}}, StepRangeLen(CartesianIndex(1, 2), CartesianIndex(1, 1), 3)) with eltype Int64:
 4
 5
 6

Closes #30250

@jishnub jishnub added linear algebra Linear algebra arrays [a, r, r, a, y, s] labels Oct 15, 2024
@mcabbott
Copy link
Contributor

mcabbott commented Oct 15, 2024

Seems OK, it's used a lot internally, and this is usually the only reason to want diagind.

However, it seems like the one-arg form should understand structured matrices, so that e.g. diagview(T) isa UnitRange. And likewise diagview(UpperTriangular(rand(4,4))) should probably act on the parent (which has linear indexing).

What diagview(UpperTriangular(rand(4,4)), 2) should do is less obvious... it could make an efficient view(::Matrix) in the useful case, but not for the -2 case, hence at the cost of type instability.

@nsajko nsajko added the feature Indicates new feature / enhancement requests label Oct 15, 2024
@jishnub
Copy link
Contributor Author

jishnub commented Oct 16, 2024

Type-stability is indeed the reason I didn't specialize this for banded matrices. This parallels how diag doesn't copy the bands, but writes to a similar vector. While the single-argument method may be specialized, I wonder if this is particularly useful? Accessing any band other than the principal diagonal would require the 2-argument method. We would also want to ensure that diagview(A) and diagview(A,0) return the same type, so whether we specialize the one-argument method goes along with whether we want type-stability in the two-argument method.

@mcabbott
Copy link
Contributor

We would also want to ensure that diagview(A) and diagview(A,0) return the same type

Why do we want this?

@jishnub
Copy link
Contributor Author

jishnub commented Oct 16, 2024

Mainly for ease of use. Having code that works for diagview(A) but not for diagview(A,0) doesn't sound ideal.

@mcabbott
Copy link
Contributor

I think you're suggesting someone may pass diagview(A) to a function which accepts only ::Vector, and may be surprised if they later change to diagview(A, 0). That seems unlikely to me, and sort-of backwards? If you ask for a view of some unknown AbstractMatrix, you generically expect some AbstractVector whose type is complicated... and if you accept that, then when someone passes in A::Diagonal, getting a simple Vector seems like it should be fine.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
arrays [a, r, r, a, y, s] feature Indicates new feature / enhancement requests linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

iterator for diagonal of a matrix ?
3 participants