Cholesky of a view of a BandedMatrix forgets the structure

[DanielVandH]

julia> using BandedMatrices, LinearAlgebra

julia> L = brand(10, 10, 2, 0); B = BandedMatrix(Symmetric(L * L'));

julia> cholesky(Symmetric(B)) # good
Cholesky{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}
U factor:
10×10 UpperTriangular{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}:
 0.959058  0.698054  0.15382    ⋅         ⋅          ⋅         ⋅         ⋅         ⋅         ⋅
  ⋅        0.184258  0.432745  0.686227   ⋅          ⋅         ⋅         ⋅         ⋅         ⋅
  ⋅         ⋅        0.311502  0.977993  0.642646    ⋅         ⋅         ⋅         ⋅         ⋅
  ⋅         ⋅         ⋅        0.148787  0.285363   0.155638   ⋅         ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅         ⋅        0.0493197  0.259591  0.55198    ⋅         ⋅         ⋅
  ⋅         ⋅         ⋅         ⋅         ⋅         0.293538  0.124905  0.113931   ⋅         ⋅
  ⋅         ⋅         ⋅         ⋅         ⋅          ⋅        0.769331  0.620623  0.106919   ⋅
  ⋅         ⋅         ⋅         ⋅         ⋅          ⋅         ⋅        0.462947  0.723893  0.75173
  ⋅         ⋅         ⋅         ⋅         ⋅          ⋅         ⋅         ⋅        0.714449  0.32878
  ⋅         ⋅         ⋅         ⋅         ⋅          ⋅         ⋅         ⋅         ⋅        0.215385

julia> cholesky(Symmetric(B[1:end-1, 1:end-1])) # good
Cholesky{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}
U factor:
9×9 UpperTriangular{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}:
 0.959058  0.698054  0.15382    ⋅         ⋅          ⋅         ⋅         ⋅         ⋅
  ⋅        0.184258  0.432745  0.686227   ⋅          ⋅         ⋅         ⋅         ⋅
  ⋅         ⋅        0.311502  0.977993  0.642646    ⋅         ⋅         ⋅         ⋅ 
  ⋅         ⋅         ⋅        0.148787  0.285363   0.155638   ⋅         ⋅         ⋅
  ⋅         ⋅         ⋅         ⋅        0.0493197  0.259591  0.55198    ⋅         ⋅
  ⋅         ⋅         ⋅         ⋅         ⋅         0.293538  0.124905  0.113931   ⋅
  ⋅         ⋅         ⋅         ⋅         ⋅          ⋅        0.769331  0.620623  0.106919
  ⋅         ⋅         ⋅         ⋅         ⋅          ⋅         ⋅        0.462947  0.723893
  ⋅         ⋅         ⋅         ⋅         ⋅          ⋅         ⋅         ⋅        0.714449

julia> cholesky(Symmetric(@view B[1:end-1, 1:end-1])) # bad
Cholesky{Float64, Matrix{Float64}}
U factor:
9×9 UpperTriangular{Float64, Matrix{Float64}}:
 0.959058  0.698054  0.15382   0.0       0.0        0.0       0.0       0.0       0.0
  ⋅        0.184258  0.432745  0.686227  0.0        0.0       0.0       0.0       0.0
  ⋅         ⋅        0.311502  0.977993  0.642646   0.0       0.0       0.0       0.0
  ⋅         ⋅         ⋅        0.148787  0.285363   0.155638  0.0       0.0       0.0
  ⋅         ⋅         ⋅         ⋅        0.0493197  0.259591  0.55198   0.0       0.0
  ⋅         ⋅         ⋅         ⋅         ⋅         0.293538  0.124905  0.113931  0.0
  ⋅         ⋅         ⋅         ⋅         ⋅          ⋅        0.769331  0.620623  0.106919
  ⋅         ⋅         ⋅         ⋅         ⋅          ⋅         ⋅        0.462947  0.723893
  ⋅         ⋅         ⋅         ⋅         ⋅          ⋅         ⋅         ⋅        0.714449

Presumably (?) it would be possible for the latter to also return something like a

Cholesky{Float64, SubArray{Float64, 2, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}}

which knows the structure. Not 100% sure how views propagate in all the code.

[dlfivefifty]

This is probably an ArrayLayouts.jl issue: we need to override cholesky for SubArray of LayoutMatrix to call cholesky_layout or something

[DanielVandH]

It does go through that code, I think the issue is when it calls cholcopy since

julia> LinearAlgebra.cholcopy(B) # bad
10×10 Matrix{Float64}:
 0.617164  0.308407  0.170248  0.0       0.0       0.0       0.0       0.0        0.0       0.0
 0.308407  1.00335   0.332173  0.605791  0.0       0.0       0.0       0.0        0.0       0.0
 0.170248  0.332173  0.214511  0.236248  0.159116  0.0       0.0       0.0        0.0       0.0
 0.0       0.605791  0.236248  0.550771  0.27541   0.018336  0.0       0.0        0.0       0.0
 0.0       0.0       0.159116  0.27541   1.04702   0.591206  0.56485   0.0        0.0       0.0
 0.0       0.0       0.0       0.018336  0.591206  0.976661  1.05175   0.184048   0.0       0.0
 0.0       0.0       0.0       0.0       0.56485   1.05175   1.24      0.510496   0.100412  0.0
 0.0       0.0       0.0       0.0       0.0       0.184048  0.510496  0.983012   0.314483  0.0262631
 0.0       0.0       0.0       0.0       0.0       0.0       0.100412  0.314483   0.245317  0.266145
 0.0       0.0       0.0       0.0       0.0       0.0       0.0       0.0262631  0.266145  0.532169

julia> LinearAlgebra.cholcopy(Symmetric(B)) # good
10×10 Symmetric{Float64, BandedMatrix{Float64, Matrix{Float64}, Base.OneTo{Int64}}}:
 0.617164  0.308407  0.170248   ⋅         ⋅         ⋅         ⋅         ⋅          ⋅         ⋅
 0.308407  1.00335   0.332173  0.605791   ⋅         ⋅         ⋅         ⋅          ⋅         ⋅
 0.170248  0.332173  0.214511  0.236248  0.159116   ⋅         ⋅         ⋅          ⋅         ⋅
  ⋅        0.605791  0.236248  0.550771  0.27541   0.018336   ⋅         ⋅          ⋅         ⋅
  ⋅         ⋅        0.159116  0.27541   1.04702   0.591206  0.56485    ⋅          ⋅         ⋅
  ⋅         ⋅         ⋅        0.018336  0.591206  0.976661  1.05175   0.184048    ⋅         ⋅
  ⋅         ⋅         ⋅         ⋅        0.56485   1.05175   1.24      0.510496   0.100412   ⋅
  ⋅         ⋅         ⋅         ⋅         ⋅        0.184048  0.510496  0.983012   0.314483  0.0262631
  ⋅         ⋅         ⋅         ⋅         ⋅         ⋅        0.100412  0.314483   0.245317  0.266145
  ⋅         ⋅         ⋅         ⋅         ⋅         ⋅         ⋅        0.0262631  0.266145  0.532169

julia> LinearAlgebra.cholcopy(Symmetric(view(B, :, :))) # bad
10×10 Symmetric{Float64, Matrix{Float64}}:
 0.617164  0.308407  0.170248  0.0       0.0       0.0       0.0       0.0        0.0       0.0
 0.308407  1.00335   0.332173  0.605791  0.0       0.0       0.0       0.0        0.0       0.0
 0.170248  0.332173  0.214511  0.236248  0.159116  0.0       0.0       0.0        0.0       0.0
 0.0       0.605791  0.236248  0.550771  0.27541   0.018336  0.0       0.0        0.0       0.0
 0.0       0.0       0.159116  0.27541   1.04702   0.591206  0.56485   0.0        0.0       0.0
 0.0       0.0       0.0       0.018336  0.591206  0.976661  1.05175   0.184048   0.0       0.0
 0.0       0.0       0.0       0.0       0.56485   1.05175   1.24      0.510496   0.100412  0.0
 0.0       0.0       0.0       0.0       0.0       0.184048  0.510496  0.983012   0.314483  0.0262631
 0.0       0.0       0.0       0.0       0.0       0.0       0.100412  0.314483   0.245317  0.266145
 0.0       0.0       0.0       0.0       0.0       0.0       0.0       0.0262631  0.266145  0.532169

[DanielVandH]

https://github.com/JuliaLinearAlgebra/BandedMatrices.jl/blob/c31e007e642962ec737649328568c43e57bc0c9f/src/symbanded/BandedCholesky.jl#L99

Probably just needs a similar method for SubArrays here. I'll make a pr