How do I use this package to pass analytical jvp?

[rveltz]

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]

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

[avik-pal]

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]

yes that's the jvp

[avik-pal]

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]

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

[ChrisRackauckas]

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]

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

[rveltz]

👍

[vpuri3]

We have come across this issue before when trying to form nonlinear operators eg, f(u) * u, df(u)/du * u, etc with ComposedOperator (#159). This seemed like a design constraint for SciMLOperators, and so when writing SparseDiffTools.VecJac, we created the kwarg interface for FunctionOperator (#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 (JuliaDiff/SparseDiffTools.jl#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]

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

Let's standardize our notation.

w = L(u,p,t)*v

[ChrisRackauckas]

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]

L should have two calls:

• L(v,u,p,t) out of place
• L(w,v,u,p,t) in place

[rveltz]

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]

How is the progress on this issue?

[vpuri3]

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]

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

[ChrisRackauckas]

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

[rveltz]

bump ;)

[avik-pal]

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]

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

[rveltz]

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

[ChrisRackauckas]

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

[HumHongeKamyaab]

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]

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]

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

[HumHongeKamyaab]

@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