Julia: ldiv!, rdiv!, and missing methods

[antoine-levitt]

(this is three issues in one, I can split if preferred; I just want to get feedback on whether the current behavior is by design or an oversight)

• ldiv! has a ldiv!(Y,A,B) form which rdiv! doesn't.

• There appears to be a number of missing methods, such as with Diagonal.

• I find it strange that ldiv!(A,A) when A is a matrix doesn't work while A\A does. I understand the rationale explained in the docstring, but presumably if an user is advanced enough to go to ldiv! they also know that factorizations are expensive and require memory. Not having this means that, for instance, a library that wants to support preconditioning can not just call ldiv! and expect to work on a wide range of types, and be performant when necessary.

[andreasnoack]

ldiv! has a ldiv!(Y,A,B) form which rdiv! doesn't.

The apparent symmetry between ldiv! and rdiv! is a bit misleading. In most other places ldiv! would have the name solve! and you probably wouldn't have expected versions of solve! that did the same as rdiv!. By far, the most common situation is ldiv! and this is reflected in the library support. It would generally be difficult to support (efficient) mutating versions of rdiv! simply because the libraries don't support them. If we wrapped the C API of LAPACK, we could probably completely cover rdiv! by utilizing the row major flag but I don't see how we could easily support rdiv! more generally. My opinion (which some don't agree with) is also that the mutating versions should be relatively thin and have as predictable memory behavior as possible. In a generic situation, where you can't reason about the memory usage, I think the right approach is to use the non-mutating versions.

There appears to be a number of missing methods, such as with Diagonal.

There are probably some three argument versions missing. Originally, dense ldiv! was an in-place two-argument business but SuiteSparse doesn't support in-place solves so the lowest common denominator was three-argument out-of-place where we copy! if the underlying solver is in-place. There are general versions for Factorizations but we'd have to add specific versions for the matrix types that support direct solves, which leads to

I find it strange that ldiv!(A,A) when A is a matrix doesn't work while AA does. I understand the rationale explained in the docstring, but presumably if an user is advanced enough to go to ldiv! they also know that factorizations are expensive and require memory.

As mentioned, I strongly prefer that the mutating versions have somewhat predictable memory behavior and you'd sell that if ldiv!(::AbstractArray,::AbstractArray,[::AbstractArray]) is allowed to factorize the matrix. Why not just callfactorize?

[antoine-levitt]

My problem is that this is a pretty exposed API. My motivation comes from using ldiv! as interface to preconditioning, e.g. in https://github.com/JuliaMath/IterativeSolvers.jl/blob/master/src/cg.jl#L75 and https://github.com/JuliaNLSolvers/Optim.jl/blob/master/src/multivariate/solvers/first_order/cg.jl#L132.
That's a pretty natural API to use in this case. This leads to either preconditioners not working easily, or plain old type piracy in https://github.com/JuliaNLSolvers/Optim.jl/blob/master/src/multivariate/precon.jl#L45. (cc @pkofod)

So it'd be good to either flesh out the API (including support for inefficient operations), or hide it more/discourage explicitly its use in generic code. It seems like in julia the ship of having generic fallbacks at the cost of potentially surprising performance has sailed a long time ago (sparse, symmetric wrappers, etc.), so might as well embrace it (dare I suggest ldiv!(Y,A,B) = copy!(Y, A\B)?). In cases where ldiv! not allocating really makes a difference, the user is probably already pretty good at finding allocations anyway.

[andreasnoack]

My problem is that this is a pretty exposed API. My motivation comes from using ldiv! as interface to preconditioning

I don't see why you couldn't just require that the matrix is passed in a form that allows for direct solves, i.e. a Factorization or Diagonal, [X]Triangular, or Q[X].

This leads to either preconditioners not working easily, or plain old type piracy

As mentioned, this is just an oversight and the definition can just be moved to LinearAlgebra.

[dkarrasch]

Regarding the second issue, there is indeed a two-arg version of ldiv!(::Diagonal,::Matrix). Which else are missing @antoine-levitt ? I could easily add those, based on the existing methods.

[antoine-levitt]

How about defining ldiv!(Y,A,B) = (copy!(Y,B); ldiv!(A,Y)), and the same for rdiv!? Then e.g. SuiteSparse can overload that, there's a consistent API and we don't need to keep chasing those three-arg forms.

[dkarrasch]

Wouldn't that mutate B, opposed to the expectation that A and B are left untouched and the result is returned in Y?

[antoine-levitt]

No, why? It copies B into Y (leaving B untouched) and then operates in place on Y.

edit: ah, sorry, I had a mistake in my previous post, I edited it

[andreasnoack]

How about defining ldiv!(Y,A,B) = (copy!(Y,B); ldiv!(A,Y))

We already have this definition when A is a Factorization. I'd be fine adding a similar definition for AbstractMatrix. (There is a risk that such a fallback definition causes some ambiguities).

[baggepinnen]

I found that both Optim and DiffEqBase pirate the same method of ldiv!, causing conflicts if both packages are loaded at the same time, could this method signature be added to base to avoid packages having to pirate it?

MethodError: ldiv!(::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}) is ambiguous. Candidates:
  ldiv!(Y::AbstractArray, A::AbstractArray, B::AbstractArray) in DiffEqBase at /home/fredrikb/.julia/packages/DiffEqBase/4V8I6/src/init.jl:5
  ldiv!(x, P::AbstractArray{T,2} where T, b) in Optim at /home/fredrikb/.julia/packages/Optim/EhyUl/src/multivariate/precon.jl:68
Possible fix, define
  ldiv!(::AbstractArray, ::AbstractArray{T,2} where T, ::AbstractArray)

[andreasnoack]

@baggepinnen So far the view has been that ldiv! shouldn't compute factorizations for dense inputs. The idea was that ldiv! shouldn't allocate to avoid a false sense of optimization when switching to a ! method, i.e. if you haven't factorized your matrix upfront anyway then you might as well just us \.

[dlfivefifty]

Not sure if this is the right place for this but I found out that ldiv! is missing adjoint versions:

julia> ldiv!(randn(5,5), factorize(randn(5,5)), randn(5,5))
5×5 Array{Float64,2}:
  1.28414   -1.45517      9.28033   12.1509     -3.98101  
 -0.940981   0.00740987  -2.06598   -2.5928      1.1498   
  0.349083   0.296984    -0.844023  -0.0261503   0.0485404
 -0.78573    0.102968    -3.71458   -3.50948     3.50485  
  0.707247  -1.91725      4.7809     8.19342    -0.918122 

julia> ldiv!(randn(5,5)', factorize(randn(5,5))', randn(5,5)')
ERROR: MethodError: no method matching ldiv!(::Adjoint{Float64,LU{Float64,Array{Float64,2}}}, ::Adjoint{Float64,Array{Float64,2}})
Closest candidates are:
  ldiv!(::Number, ::AbstractArray) at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/generic.jl:152
  ldiv!(::Adjoint{T<:Real,#s617} where #s617<:(LU{T<:Real,#s616} where #s616<:(Union{DenseArray{T,2}, Base.ReinterpretArray{T,2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T,2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}} where T)), ::Union{DenseArray{T<:Real,1}, DenseArray{T<:Real,2}, Base.ReinterpretArray{T<:Real,2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReinterpretArray{T<:Real,1,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T<:Real,1,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T<:Real,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T<:Real,1,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}, SubArray{T<:Real,2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}}) where T<:Real at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/lu.jl:359
  ldiv!(::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T, ::Factorization, ::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/factorization.jl:100
  ...
Stacktrace:
 [1] ldiv!(::Adjoint{Float64,Array{Float64,2}}, ::Adjoint{Float64,LU{Float64,Array{Float64,2}}}, ::Adjoint{Float64,Array{Float64,2}}) at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/factorization.jl:113
 [2] top-level scope at REPL[74]:1

julia> ldiv!(randn(5,5)', factorize(randn(5,5)'), randn(5,5)')
ERROR: MethodError: no method matching ldiv!(::Adjoint{Float64,LU{Float64,Array{Float64,2}}}, ::Adjoint{Float64,Array{Float64,2}})
Closest candidates are:
  ldiv!(::Number, ::AbstractArray) at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/generic.jl:152
  ldiv!(::Adjoint{T<:Real,#s617} where #s617<:(LU{T<:Real,#s616} where #s616<:(Union{DenseArray{T,2}, Base.ReinterpretArray{T,2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T,2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}} where T)), ::Union{DenseArray{T<:Real,1}, DenseArray{T<:Real,2}, Base.ReinterpretArray{T<:Real,2,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReinterpretArray{T<:Real,1,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T<:Real,1,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T<:Real,2,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T<:Real,1,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}, SubArray{T<:Real,2,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}}) where T<:Real at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/lu.jl:359
  ldiv!(::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T, ::Factorization, ::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/factorization.jl:100
  ...
Stacktrace:
 [1] ldiv!(::Adjoint{Float64,Array{Float64,2}}, ::Adjoint{Float64,LU{Float64,Array{Float64,2}}}, ::Adjoint{Float64,Array{Float64,2}}) at /Users/sabae/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/factorization.jl:113
 [2] top-level scope at REPL[75]:1

[antoine-levitt]

That's a different issue: in this case it's correctly dispatched to the 2arg ldiv!. I think the transpose types are mostly "read-only" and do not support fast copy to them, so it makes sense that this method is not supported.

Ie this works:
ldiv!(randn(5,5), factorize(randn(5,5))', randn(5,5)')

[dlfivefifty]

Usually BLAS/LAPACK support transposes (ie row major layout). I’ll have double check if gesv! is an exception...