[LinearALgebra] Should Diagonal specialize on AdjOrTransAbsVec?

[briochemc]

I was surprised that Diagonal turns 1×n adjoints into 1×1 diagonals:

julia> Diagonal(rand(2)')
1×1 Diagonal{Float64,Array{Float64,1}}:
 0.05075540276803281

That's because Diagonal considers the adjoint vector as a 2D array via this method:

https://github.com/JuliaLang/julia/blob/28ff641fa9130d760173122bb231ced2fedf839f/stdlib/LinearAlgebra/src/diagonal.jl#L36

Maybe we could instead add another specialized Diagonal method for adjoints/transposes of vectors (using AdjOrTransAbsVec)?

Specifically, I was thinking of essentially duplicating

https://github.com/JuliaLang/julia/blob/28ff641fa9130d760173122bb231ced2fedf839f/stdlib/LinearAlgebra/src/diagonal.jl#L13-L14

but using AdjOrTransAbsVec instead of AbstractVector.

What do you think?

[dkarrasch]

It seems this is a feature, and we have tests for that. See JuliaLang/julia#23706. From the discussion in #466, this was done on purpose, to make AdjOrTransAbsVec behave like a 1xn matrix.

[mcabbott]

I think it wouldn't be crazy to call this a bug, and fix it.

Adjoint vectors pretend to be 1 x N, but within LinearAlgebra the this first dimension is sometimes dropped. It's a slightly awkward compromise that lets many things work well. The issue giving them matrix-like indexing was I think from when this compromise was being adopted.

The tests that Diagonal gives a 1x1 matrix seem to have been added in the process of fixing a more serious bug. Not crashing is very desirable but is the 1x1 matrix actually wanted? Maybe there was some infinite loop between the then-new matrix-like indexing and earlier behaviour... what did Diagonal(::RowVector) do on Julia 0.6 or so?

[dkarrasch]

               _
   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: https://docs.julialang.org
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.6.4 (2018-07-09 19:09 UTC)
 _/ |\__'_|_|_|\__'_|  |  Official http://julialang.org/ release
|__/                   |  x86_64-apple-darwin14.5.0

julia> x = rand(2)'
1×2 RowVector{Float64,Array{Float64,1}}:
 0.658066  0.302965

julia> Diagonal(x)
ERROR: StackOverflowError:
Stacktrace:
 [1] steprange_last(::Int64, ::Int64, ::Int64) at ./range.jl:0
 [2] Type at ./range.jl:93 [inlined]
 [3] _range at ./range.jl:61 [inlined]
 [4] range at ./range.jl:60 [inlined]
 [5] diagind(::Int64, ::Int64, ::Int64) at ./linalg/dense.jl:206
 [6] diag(::RowVector{Float64,Array{Float64,1}}, ::Int64) at ./linalg/dense.jl:250
 [7] Diagonal(::RowVector{Float64,Array{Float64,1}}) at ./linalg/diagonal.jl:29 (repeats 79997 times)

So, it hasn't been fixed in v0.6.x. I don't have a real preference. I think I can see both points: AdjOrTransAbsVec are backed by a vector and should be regarded as such, vs. they are functionals on vectors, hence linear operators/row matrices, and should be regarded as a matrix. I suspect there is no right or wrong here, just a choice to be made and then documented.