Issues with Layout of inv(Triangular{ApplyLayout})

[DanielVandH]

From JuliaLinearAlgebra/ArrayLayouts.jl#239

---

Consider

julia> using LazyArrays, ArrayLayouts, LinearAlgebra

julia> A = UpperTriangular(ApplyArray(inv, rand(5, 5)))
5×5 UpperTriangular{Float64, ApplyArray{Float64, 2, typeof(inv), Tuple{Matrix{Float64}}}}:
 2.35135   1.68449   -1.84848   -0.955078  -1.85021
  ⋅       -0.271777  -1.25844    0.94071   -1.27413
  ⋅         ⋅        -0.124992  -1.4524     2.35783
  ⋅         ⋅          ⋅         0.669905  -0.609163
  ⋅         ⋅          ⋅          ⋅         0.840739

julia> inv(A)
inv(5×5 UpperTriangular{Float64, ApplyArray{Float64, 2, typeof(inv), Tuple{Matrix{Float64}}}}):
  0.425287   2.63595  -32.8287   -74.2701    43.1852
 -0.0       -3.67949   37.0457    85.4846   -47.5317
 -0.0       -0.0       -8.00051  -17.3457     9.86934
  0.0        0.0        0.0        1.49275    1.08158
  0.0        0.0        0.0        0.0        1.18943

julia> MemoryLayout(inv(A))
LazyArrays.InvLayout{TriangularLayout{'U', 'N', LazyArrays.LazyLayout}}()

The LazyArrays.InvLayout{TriangularLayout{'U', 'N', LazyArrays.LazyLayout}}() could maybe be TriangularLayout{'L', 'N', LazyArrays.LazyLayout}()?

One issue that the current situation causes is e.g.

julia> using BandedMatrices

julia> colsupport(inv(A), 1)
Base.OneTo(5)

which is of course more problematic for infinite matrices.

My suggestion to just simplify the InvLayout using

applylayout(::Type{typeof(inv)}, ::TriangularLayout{'U', UNIT, ML}) where {UNIT, ML} = TriangularLayout{'L', UNIT, ML}()
applylayout(::Type{typeof(inv)}, ::TriangularLayout{'L', UNIT, ML}) where {UNIT, ML} = TriangularLayout{'U', UNIT, ML}()

doesn't seem to be good either since

julia> inv(A)
inv(5×5 UpperTriangular{Float64, ApplyArray{Float64, 2, typeof(inv), Tuple{Matrix{Float64}}}}):
Error showing value of type ApplyArray{Float64, 2, typeof(inv), Tuple{UpperTriangular{Float64, ApplyArray{Float64, 2, typeof(inv), Tuple{Matrix{Float64}}}}}}:
ERROR: StackOverflowError:
Stacktrace:
     [1] getindex(A::ApplyArray{Float64, 2, typeof(inv), Tuple{UpperTriangular{Float64, ApplyArray{…}}}}, kr::Colon, jr::Int64)
       @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\70nvu\src\ArrayLayouts.jl:156
     [2] getindex(A::ApplyArray{Float64, 2, typeof(inv), Tuple{UpperTriangular{Float64, ApplyArray{…}}}}, k::Int64, j::Int64)
       @ LazyArrays c:\Users\djv23\.julia\dev\LazyArrays.jl\src\linalg\inv.jl:178
     [3] getindex(V::SubArray{Float64, 1, ApplyArray{Float64, 2, typeof(inv), Tuple{…}}, Tuple{Base.Slice{…}, Int64}, false}, I::Int64)
       @ Base .\subarray.jl:290
     [4] iterate(A::SubArray{Float64, 1, ApplyArray{…}, Tuple{…}, false}, state::Tuple{Base.OneTo{…}})
       @ Base .\abstractarray.jl:1217
     [5] iterate
       @ .\abstractarray.jl:1215 [inlined]
     [6] copyto_unaliased!(deststyle::IndexLinear, dest::Vector{…}, srcstyle::IndexCartesian, src::SubArray{…})
       @ Base .\abstractarray.jl:1093
     [7] copyto!(dest::Vector{Float64}, src::SubArray{Float64, 1, ApplyArray{Float64, 2, typeof(inv), Tuple{…}}, Tuple{Base.Slice{…}, Int64}, false})
       @ Base .\abstractarray.jl:1068
     [8] copyto!_layout
       @ C:\Users\djv23\.julia\packages\ArrayLayouts\70nvu\src\ArrayLayouts.jl:259 [inlined]
     [9] copyto!_layout
       @ C:\Users\djv23\.julia\packages\ArrayLayouts\70nvu\src\ArrayLayouts.jl:265 [inlined]
    [10] copyto!
       @ C:\Users\djv23\.julia\packages\ArrayLayouts\70nvu\src\ArrayLayouts.jl:274 [inlined]
    [11] copyto_axcheck!
       @ .\abstractarray.jl:1177 [inlined]
    [12] Array
       @ .\array.jl:673 [inlined]
    [13] Array
       @ .\boot.jl:500 [inlined]
    [14] sub_materialize_axes
       @ C:\Users\djv23\.julia\packages\ArrayLayouts\70nvu\src\ArrayLayouts.jl:129 [inlined]
    [15] sub_materialize
       @ C:\Users\djv23\.julia\packages\ArrayLayouts\70nvu\src\ArrayLayouts.jl:130 [inlined]
    [16] sub_materialize
       @ C:\Users\djv23\.julia\packages\ArrayLayouts\70nvu\src\ArrayLayouts.jl:131 [inlined]
    [17] sub_materialize
       @ C:\Users\djv23\.julia\packages\ArrayLayouts\70nvu\src\ArrayLayouts.jl:132 [inlined]
    [18] layout_getindex(::ApplyArray{Float64, 2, typeof(inv), Tuple{UpperTriangular{Float64, ApplyArray{…}}}}, ::Function, ::Int64)
       @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\70nvu\src\ArrayLayouts.jl:138
--- the last 18 lines are repeated 3103 more times ---
 [55873] getindex(A::ApplyArray{Float64, 2, typeof(inv), Tuple{UpperTriangular{Float64, ApplyArray{…}}}}, kr::Colon, jr::Int64)
       @ ArrayLayouts C:\Users\djv23\.julia\packages\ArrayLayouts\70nvu\src\ArrayLayouts.jl:156
 [55874] getindex
       @ c:\Users\djv23\.julia\dev\LazyArrays.jl\src\linalg\inv.jl:178 [inlined]
 [55875] isassigned(::ApplyArray{Float64, 2, typeof(inv), Tuple{UpperTriangular{Float64, ApplyArray{…}}}}, ::Int64, ::Int64)
       @ Base .\multidimensional.jl:1578
 [55876] alignment(io::IOContext{…}, X::AbstractVecOrMat, rows::Vector{…}, cols::Vector{…}, cols_if_complete::Int64, cols_otherwise::Int64, sep::Int64, ncols::Int64)
       @ Base .\arrayshow.jl:68
 [55877] _print_matrix(io::IOContext{…}, X::AbstractVecOrMat, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64, rowsA::UnitRange{…}, colsA::UnitRange{…})
       @ Base .\arrayshow.jl:207
 [55878] print_matrix(io::IOContext{…}, X::ApplyArray{…}, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64)
       @ Base .\arrayshow.jl:171
 [55879] print_matrix
       @ .\arrayshow.jl:171 [inlined]
 [55880] print_array
       @ .\arrayshow.jl:358 [inlined]
 [55881] show(io::IOContext{Base.TTY}, ::MIME{Symbol("text/plain")}, X::ApplyArray{Float64, 2, typeof(inv), Tuple{UpperTriangular{…}}})
       @ Base .\arrayshow.jl:399
 [55882] (::REPL.var"#55#56"{REPL.REPLDisplay{REPL.LineEditREPL}, MIME{Symbol("text/plain")}, Base.RefValue{Any}})(io::Any)
       @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:273
 [55883] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
       @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:569
 [55884] display(d::REPL.REPLDisplay, mime::MIME{Symbol("text/plain")}, x::Any)
       @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:259
 [55885] display(d::REPL.REPLDisplay, x::Any)
       @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:278
 [55886] display(x::Any)
       @ Base.Multimedia .\multimedia.jl:340
 [55887] #invokelatest#2
       @ .\essentials.jl:892 [inlined]
 [55888] invokelatest
       @ .\essentials.jl:889 [inlined]
 [55889] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{Nothing, AbstractDisplay})
       @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:315
 [55890] (::REPL.var"#57#58"{REPL.LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
       @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:284
 [55891] with_repl_linfo(f::Any, repl::REPL.LineEditREPL)
       @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:569
 [55892] print_response(repl::REPL.AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
       @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:282
 [55893] (::REPL.var"#do_respond#80"{…})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
       @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:911
 [55894] (::VSCodeServer.var"#103#106"{REPL.var"#do_respond#80"{…}})(mi::REPL.LineEdit.MIState, buf::IOBuffer, ok::Bool)
       @ VSCodeServer c:\Users\djv23\.vscode\extensions\julialang.language-julia-1.83.2\scripts\packages\VSCodeServer\src\repl.jl:122
 [55895] #invokelatest#2
       @ .\essentials.jl:892 [inlined]
 [55896] invokelatest
       @ .\essentials.jl:889 [inlined]
 [55897] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
       @ REPL.LineEdit C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\LineEdit.jl:2656
 [55898] run_frontend(repl::REPL.LineEditREPL, backend::REPL.REPLBackendRef)
       @ REPL C:\Users\djv23\.julia\juliaup\julia-1.10.4+0.x64.w64.mingw32\share\julia\stdlib\v1.10\REPL\src\REPL.jl:1312
Some type information was truncated. Use `show(err)` to see complete types

[DanielVandH]

One issue that the current situation causes is e.g.

Of course for this example I could make a PR extending rowsupport and colsupport, but maybe there is a better approach for this somehow to help more broadly.

[dlfivefifty]

I think overloading col and row support for inverses of triangular matrices would indeed be a good idea

[DanielVandH]

I'll make a PR