Cannot solve with `AdjointFactorization` of Cholesky factorizations and `Adjoint` right-hand side

[jonas-schulze]

Consider the following MWE:

using LinearAlgebra, SparseArrays

n = 10
A = sprand(n, n, 1/n)
A = A + A' + 10I

B = rand(n, 2)
Bt = Matrix(B')

F = factorize(A')
F \ B   # works just fine
F \ Bt' # breaks

ERROR: MethodError: no method matching ldiv!(::LinearAlgebra.AdjointFactorization{Float64, SparseArrays.CHOLMOD.Factor{Float64, Int64}}, ::Matrix{Float64})

Closest candidates are:
  ldiv!(::BunchKaufman{T, <:StridedMatrix{T} where T}, ::StridedVecOrMat{T}) where T<:Union{Float32, Float64}
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.5+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/bunchkaufman.jl:329
  ldiv!(::LinearAlgebra.QRCompactWY{T, M, C} where {M<:AbstractMatrix{T}, C<:AbstractMatrix{T}}, ::AbstractMatrix{T}) where T
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.5+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/qr.jl:523
  ldiv!(::UpperHessenberg, ::AbstractVecOrMat; shift)
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.5+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/hessenberg.jl:179
  ...

Stacktrace:
 [1] \(F::LinearAlgebra.AdjointFactorization{Float64, SparseArrays.CHOLMOD.Factor{Float64, Int64}}, B::Adjoint{Float64, Matrix{Float64}})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.5+0.x64.linux.gnu/share/julia/stdlib/v1.10/LinearAlgebra/src/factorization.jl:143
 [2] top-level scope
   @ REPL[27]:1
 [3] top-level scope
   @ none:1

A similar situation happens in one of my codes that implements the Sherman-Morrison-Woodbury formula. What is the proper way to fix this?

1. Define adjoint(F::SparseArrays.CHOLMOD.Factor) = F?
2. Make Adjoint{T, Matrix{T}} a StridedMatrix?
3. Define the missing method of ldiv!?
4. ...?

Option 1 would work, but seems out of place given this method of (\) (permalink created from release-1.10). I didn't test options 2 or 3. What would be the best option to have this work with the current LTS?

[dkarrasch]

You may need to catch the case where the rhs is a AdjOrTrans, and materialize the outer wrapper:

\(F::AdjointFactorization{...}, Bx::AdjOrTrans) = F \ copy(Bx)

julia> F \ copy(Bt') # breaks
10×2 Matrix{Float64}:
  0.0229374   0.0222694
 -0.00846876  0.0206852
 -0.00239363  0.00159422
  0.0441779   0.044418
  0.0681237   0.0617245
  0.0571472   0.032829
  0.0510035   0.0234235
  0.0745302   0.0976153
  0.0636938   0.0970002
  0.0142583   0.0903773

Edit: I guess this is your option 2.

[jonas-schulze]

Thanks for the very swift response! No, by option 2 I meant to somehow adjust the definition of StridedMatrix, such that F \ Bt' hits this method (which I did not try). There has been some effort in that direction:

• https://discourse.julialang.org/t/adjoint-of-matrix-is-not-a-stridedarray/43038
• https://discourse.julialang.org/t/should-slices-of-adjoint-transpose-be-stridedvectors/98297

Your solution is more related to option 3, but sounds reasonable to me as well.

Edit: removed mention of SparseArrays (fallacy / "Denkfehler") and added two links to Discourse.