Overload operator application for ComposedOperator

[vpuri3]

For ComposedOperator the present update_coeffs(L, u, p, t) interface will update each op in L.ops with u, which is not it's input.

SciMLOperators.jl/src/interface.jl

Lines 22 to 27 in 31275d4

	function update_coefficients!(L::AbstractSciMLOperator, u, p, t)
	for op in getops(L)
	update_coefficients!(op, u, p, t)
	end
	nothing
	end

For example, with

N = 4
A = DiagonalOperator(rand(N); update_func = (d, u, p, t) -> copy!(d, u)) # A * u == u .^2
B = DiagonalOperator(rand(N); update_func = (d, u, p, t) -> copy!(d, u)) # A * u == u .^2
C = DiagonalOperator(rand(N); update_func = (d, u, p, t) -> copy!(d, u)) # A * u == u .^2

L = A ∘ B ∘ C
u = 0.1 * ones(N)

v = L(u, nothing, 0.0)

julia> v = L(u, nothing, 0.0)
4-element Vector{Float64}:
 0.00010000000000000003
 0.00010000000000000003
 0.00010000000000000003
 0.00010000000000000003

In this case, what we want to get by composing three squaring operations is u .^ 8, but we will get u .^ 4 because A,B,C will be updated by u, u, u, not u, u^2, u^4 respectively. What we want is to update C with u, B with its input C * u = u .^ 2, and A with B * C * u = u .^ 4. The operator update procedure must necessarily be interlaced with multiplication.

function (L::ComposedOperator)(u, p, t)
    v = u
    for op in reverse(L.ops)
        update_coefficients!(op, v, p, t)
        v = op * v
    end
    v
end

So we need to overload update-coeffs and operator application for ComposedOperator. It would be wise to take a second look at other composite operators to ensure such a thing isn't happening anywhere else. Maybe we need to define update_coeffs overload for AdjointOperator, TransposedOperator to pass in u'. Should also take a second look at FunctionOperator.

[vpuri3]

@gaurav-arya , would you be able to help with this?

[gaurav-arya]

Hey, I'd be happy to help, but: I don't think I agree this is a mistake? Shouldn't u be the state of current solution in OrdinaryDiffEq? And in general, something fixed, not changing for different operators? I don't really understand this use case. (A*u = u.^2 is a non-linear op, so I don't think it makes sense to try to enforce this.)

In a similar vein (but different context), I was also confused by lines such as https://github.com/SciML/SciMLOperators.jl/blob/31275d4e2b2f3671fb4c2811656dd555c4e1f24a/test/matrix.jl#LL174-L176C24 in tests. The b should be used as a buffer there, but it doesn't make sense to me for the state update to actually depend on what was in b.

Maybe we have different understandings of what update_coefficients means?

[vpuri3]

It is a mistake for the simple reason that

julia> A(B(C(u, p, t), p, t), p, t) ≈ (A ∘ B ∘ C)(u, p, t)
false

julia> A(B(C(u, p, t), p, t), p, t)
4-element Vector{Float64}:
 1.0000000000000008e-8
 1.0000000000000008e-8
 1.0000000000000008e-8
 1.0000000000000008e-8

julia> (A ∘ B ∘ C)(u, p, t)
4-element Vector{Float64}:
 0.00010000000000000003
 0.00010000000000000003
 0.00010000000000000003
 0.00010000000000000003

This problem occurs specifically when you want to compose nonlinear operators, i.e L depends on u. One application of SciMLOperators is to be the ODEFunction passed to OrdinaryDiffEq. If you want A to be the operator that represents pointwise multiplication with the ODE solution, you would write it as follows:

N = 4

uODE = rand(N)

A = DiagonalOperator(rand(N); update_func = (d, u, p, t) -> copy!(d, uODE)) # A * u = uODE .* u
B = DiagonalOperator(rand(N); update_func = (d, u, p, t) -> copy!(d, uODE))
C = DiagonalOperator(rand(N); update_func = (d, u, p, t) -> copy!(d, uODE))

L = A ∘ B ∘ C
u = rand(N)

v = L(u, nothing, 0.0) # ≈ uODE .^3 .* u

In my applications till now, I haven't needed to compose multiple nonlinear operations. I have only needed to compose nonlinear operations with linear operations, thus I never came across this error. For example, to solve 1D Burgers, du/dt = -u * du/dx - nu * d2u/d2x. So my SciMLOp is

A = nu * D2x
C =  DiagonalOperator(u0, update_func = (v, u, p, t) -> copy!(v, u)) * Dx
dudt = A - C

So when I pass dudt to OrdinaryDiffEq, things work correctly as update_func correctly forms the nonlinear operation. Of course, I'd have to go back and fix that code later.

However, we should do things correctly. L = A ∘ B ∘ C means you apply C, then B, and then A, even if they are nonlinear. If we want to make an operator that multiplies (pointwise) with the current ODE solution, which is what you have in mind, we should write it using your kwargs approach.

uODE 
A = DiagonalOperator(zeros(N), update_func = (d, u, p, t; uODE) -> copy!(d, uODE)))

[gaurav-arya]

I think I kind of see where you're coming from, but it still does seem weird to me to represent arbitrary nonlinear operators through this weird back and forth between update_coefficients! and a linear *, where the action of the former depends on the input to the latter. Seems like for the arbitrary nonlinear case, update_coefficients and * should both never be used as the distinction doesn't really make sense, only op(u,p,t) etc. directly like a usual DiffEq function.

[vpuri3]

Correct

[vpuri3]

We should always use L(du, u, p, t), L(u, p, t). The *, mul! are for convenience.