[DOC] what is the expected output type for the Jacobian?

[Datseris]

Hi there,

I'm trying to come up with a simple example for the bifurcation of 1-dimensional system to contribute to the documentation, under guidance of @gszep . I have the example running, by altering the code https://rveltz.github.io/BifurcationKit.jl/dev/iterator/# , but there are some things that are not yet clear enough for me.

What I really want to know is what is the allowed return types for both the vector field f and the Jacobian. The code has:

Jac_m = (x, p) -> diagm(0 => 1.0  .- x.^k)

but it was confusing for me because when I tried other types of Jac_m I got errors (I am not reporting errors here until I first see what is the expected form of Jacobian). I also tried to use the output of ForwardDiff.jacobian, also getting errors.

So, can you please document somewhere centrally in the documentation what is the expected forms for f, J? Not only what are the input arguments, but what is the expected output. The documentation of e.g. https://rveltz.github.io/BifurcationKit.jl/dev/library/#BifurcationKit.continuation does not discuss with clarity what J should be. There you can also provide tips and tricks, as e.g. the above usage of diagm is definitely something advanced for me, I would expected a 1-element matrix as return type in this example.

---

A comment on how the docstring is structured: it provides:

F = (x, p) -> F(x, p)

This, almost always, is invalid Julia code. If F is an existing function, here you are redefining it as an anonymous function F, and Julia will just error. It also doesn't highlight what is the return type of F, yet repeats what are the expected input arguments.

Why not consider replacing it with the more valid F(x,p)::AbstractVector, or simply explicitly write things out as e.g. "F is a function with input arguments (x, p) returning a vector r that represents the vector field of the dynamics".

---

Are static vectors supported?

[rveltz]

Hi,

For F, it is like you said. "F is a function with input arguments (x, p) returning a vector r that represents the vector field of the dynamics and for type stability, the types of x and r should match".

For the jacobian J, it should be a function J(x,p) which returns a function taking one argument dx and returns dr of the same type of dx. It just happens that when J(x,p) returns a matrix, you can use it like J(x,p)(dx) and so the above works.

When J(x,p) is a fancy function, you can pass your own linear solver if you have a better one for your problem.

---

Jac_m = (x, p) -> diagm(0 => 1.0  .- x.^k)

This iterator docs, uses the above definition for the jacobian. It returns a diagonal matrix, as explained above.

For the use or ForwardDiff, please check the docs online, it has been used many times in the tutorials either in the matrix way or functions.

---

For static vectors, I have not tried and it is unlikely to work I would guess. For example, this line would break. This could be corrected by specializing the methods BorderedArrays.jl of to static arrays.

---

I will modify the docs with your suggestions. Thank you for opening this issue.

[Datseris]

For the jacobian J, it should be a function J(x,p) which returns a function taking one argument dx and returns dr of the same type of dx. It just happens that when J(x,p) returns a matrix, you can use it like J(x,p)(dx) and so the above works.

I don't understand this at all :( What is Jacobian for you? I recall that Jacobian is the derivative of the vector field F towards each of the variables, like so:

In our discussion so far, only two things exist: x, p. The Jacobian is derivatives of F towards x. What is r? What is dx and what is dr?

There is another problem: Jac_m = (x, p) -> diagm(0 => 1.0 .- x.^k) is not a function that returns a function. It is a function that returns a matrix.

I also don't understand

It just happens that when J(x,p) returns a matrix, you can use it like J(x,p)(dx)

rand(10,10)(whatever_argument) errors because matrices are not callable in Julia...

[rveltz]

I was wrong!

The jacobian object is used for 2 things

1. solve J * dx = dy in dx given dy (basically in Newton)
2. compute the eigenvalues of J (for bifurcation detection)

For example, the first part is done in the default linear solver implemented as callable struct DefaultLS which itself calls \.

Thus, 2 cases are allowed:

• either J is a function and J(x,p) returns a ::Matrix. In this case, the default arguments of continuation will make everything work. If you provide an object with a backslash operator, it should work too (I guess).
• or J is a function and J(x,p) returns a function taking one argument dx and returns dr of the same type of dx. In your notations, dr = Jf * dx. In this case, the default parameters of continuation will not work and you have to use a Matrix Free linear solver, for example GMRESIterativeSolvers.

[Datseris]

AHAAAAA!

Alright, but then, for the first case, shouldn't in-place versions be much more performant? I.e. a function

Jacobian!(J, x, p)

that updates an existing matrix J in-place? We have such versions in both DifferentialEquations.jl as well as DynamicalSystems.jl

[rveltz]

Probably but for the matrix-free case, it should not matter much as you have to build a Krylov space anyway: this is the usual bottleneck for the continuation.

Now I can put everything inplace for the jacobian solver without breaking much in the package.

[Datseris]

I am thinking of making an interface from DynamicalSystems.jl to BifurcationKit.jl, because there for any DynamicalSystem, we anyways create a Jacobian for all systems using ForwardDiff (if the user doesn't provide a manual one, fastest version). But if you have a look at system creation

It seems that at the moment neither of the two approaches would work here. Thus, if any of these two Jacobian forms became possible here, please do tag me! Such an interface would be very helpful, given how well bifurcation analysis ties in with dynamical systems theory. Even though the main goal of this package is PDEs, nothing stops us from using it for ODEs as well!

[rveltz]

I just corrected the docs of newton and continuation and hope it makes things clearer.

What is n in jacobian(x,p,n)?

I dont understand why it wouldnot work

[Datseris]

n is just the step number, which is not used in BifurcationKit.jl, because it is assumed explicitly that no time dependence exists (which is fine).

It won't work because the out of place version returns a static matrix (performance gains) instead of Matrix. On the other hand, the in place version, that uses Matrix, doesn't have the syntax you expect.

[rveltz]

Coming back to this. Can't you do something along the lines of:

const Jnew = 1 #put something meaningfull

res = newton(F, (x,p) -> (jacobian!(Jnew, x, p, 0); Jnew), ...)

[Datseris]

Indeed! Since I anyways have to make a transformation to drop the t argument, then I can do this directly.

[rveltz]

Can we close this issue?

[Datseris]

Yeap! If I have other doc questions or suggestions I'll open new issues!

[rveltz]

👍