search

SciML / SciMLOperators.jl

SciMLOperators.jl: Matrix-Free Operators for the SciML Scientific Machine Learning Common Interface in Julia
https://docs.sciml.ai/SciMLOperators/stable
MIT License
44 stars 10 forks source link

How do I use this package to pass analytical jvp? #223

Open rveltz opened 1 year ago

rveltz commented 1 year ago

Hi,

Say I have the following vector field:

f(du,u,p,t) = @. du = -u^2

I know that the jacobian is given by

myjac(out,in,u,p,t) = @. out = -2 * u * in

How do I pass this information to ODEFunction or ODEProblem?

SciMLOperators.FunctionOperator is just asking for the linear op in the form op(out, in, p, t), but where is the information about u?

ChrisRackauckas commented 11 months ago

@avik-pal you just took a look through this could you comment?

avik-pal commented 11 months ago

myjac(out,in,u,p,t) = @. out = -2 * u * in

What is in? This looks like the JVP / VJP and not the Jacobian itself.

rveltz commented 11 months ago

yes that's the jvp

avik-pal commented 11 months ago

Right, see how the JacVec operator is defined. Instead of using ad for the internal op, you will have to specify it analytically for every u. The other option is the use the kwargs interface for FunctionOperator and pass in in as a kwarg, for this see the VecJac implementation. Both are in https://github.com/JuliaDiff/SparseDiffTools.jl/blob/master/src/differentiation/jaches_products.jl & https://github.com/JuliaDiff/SparseDiffTools.jl/blob/master/src/differentiation/vecjac_products.jl

rveltz commented 11 months ago

I thought DE would provide a built in solution like DiffEqOperator used to.

ChrisRackauckas commented 11 months ago

This seems like a mistake. You could use parameters here, but it feels off. The function operator is L(u,p,t) so that it defines mul!(y,L,x) to be y = L(u,p,t)*x. The u seems to just be lopped off here. You can do L((u,p),t)*x but that seems like a mistake. input isn't u: you're not necessarily always doing L(u,p,t)*u, so the state of the operator needs to update independent of the vector that is being multiplied.

Wouldn't this make JacVec easier to implement too?

ChrisRackauckas commented 11 months ago

We talked, FunctionOperator is just wrong. We need to fix that.

rveltz commented 11 months ago

👍

vpuri3 commented 11 months ago

We have come across this issue before when trying to form nonlinear operators eg, f(u) * u, df(u)/du * u, etc with ComposedOperator (https://github.com/SciML/SciMLOperators.jl/issues/159). This seemed like a design constraint for SciMLOperators, and so when writing SparseDiffTools.VecJac, we created the kwarg interface for FunctionOperator (https://github.com/SciML/SciMLOperators.jl/issues/125). That would allow calls to VecJac like L(v, u, p, t; VJP_input = w) which is equivalent to v .= df(w)/dw * u (https://github.com/JuliaDiff/SparseDiffTools.jl/pull/245). This interface wasn't extended to SparseDiffTools.JacVec but maybe it should have. I don't know what the current problem you're facing is, but maybe a similar solution for JacVec could help alleviate the problem?

ChrisRackauckas commented 11 months ago

That would allow calls to VecJac like L(v, u, p, t; VJP_input = w) which is equivalent to v .= df(w)/dw * u (https://github.com/JuliaDiff/SparseDiffTools.jl/pull/245).

Let's standardize our notation.

w = L(u,p,t)*v
ChrisRackauckas commented 11 months ago

The L someone would define here is just L = FunctionOperator((w,v,u,p,t)-> ..., u, p, t), where update_coefficients!(L,u,p,t) updates the internal states.

ChrisRackauckas commented 11 months ago

L should have two calls:

rveltz commented 11 months ago

Let's standardize our notation.

If I understand you well, L(u,p,t) is a linear operator representing the differential of say F(u,p,t) wrt to the first variable at position u.

rveltz commented 11 months ago

How is the progress on this issue?

vpuri3 commented 11 months ago

I'm tied down with school work at the moment. I have a long flight this weekend when I'm gonna start working on this. although it's a small change it's gonna require releasing v0.4 and subsequent modifications in downstream packages.

rveltz commented 10 months ago

Sorry to be pushy... did you make any progress?

ChrisRackauckas commented 10 months ago

No. But, CI in a few places has been blocked by https://github.com/bifurcationkit/BifurcationKit.jl/pull/127 so while I have your attention if that can get fixed that would help a lot.

rveltz commented 9 months ago

bump ;)

avik-pal commented 9 months ago

The next release of NonlinearSolve will have a JacobianOperator which could be a reference for the new FunctionOperator re-design https://github.com/SciML/NonlinearSolve.jl/blob/ap/rework_internals/src/internal/operators.jl

rveltz commented 5 months ago

Have we made progress on this side? it is still unusable for me.

rveltz commented 4 months ago

I see that you added the operators. Can you show me an example how to solve the current issue please?

ChrisRackauckas commented 4 months ago

No progress has been made. This is probably the package in the worst state right now. I need to find time for it.

HumHongeKamyaab commented 1 month ago

I am posting it here because this issue may be related to my following query I want to use matrix-free Newton Krylov to solve the ODE system using implicit methods (for e.g. KenCarp47(linsolve = KrylovJL_GMRES())) https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#Using-Jacobian-Free-Newton-Krylov I want to directly specify user-defined Jacobian vector product (JVP) for my system jvp(Jv,v,u,p,t) The ode function has this user option here https://docs.sciml.ai/DiffEqDocs/stable/types/ode_types/#SciMLBase.ODEFunction

But I don't know how to tell Krylov.jl that its mul! function is defined using above jvp. I am not sure if these options are connected to each other.

HumHongeKamyaab commented 1 month ago

I want to directly specify user-defined Jacobian vector product (JVP) for my system jvp(Jv,v,u,p,t)

Hi @ChrisRackauckas, Is there any update on this issue, or Is there any alternative way without using FunctionalOperator to specify JVP of the system so that Krylov.jl can use it in mul! via DifferencialEquation.jl.

avik-pal commented 1 month ago

We now have https://docs.sciml.ai/NonlinearSolve/dev/api/SciMLJacobianOperators/ but currently only has bindings if you specify it is a nonlinearproblem

HumHongeKamyaab commented 1 month ago

@avik-pal Thanks. I am not sure if can solve for problem defined using NonlinearSolve.jl during an impilcit ODE step of DifferencialEquation.jl

There is something about it mentioned in docs https://docs.sciml.ai/DiffEqDocs/stable/features/linear_nonlinear/#Nonlinear-Solvers:-nlsolve-Specification but its not using NonlinearSolve.jl

It looks like DifferencialEquation.jl uses their own NonLinearSolve type https://github.com/SciML/OrdinaryDiffEq.jl/blob/90ba736442ce9da6f38f2db907d203cb65f61dc5/lib/OrdinaryDiffEqNonlinearSolve/src/type.jl#L2