In the documentation it is specified the following about the Jacobian:
Jacobian
The keyword J provides the Jacobian function. It must be a Julia function in the same form as f, the dynamic_rule. Specifically, J(u, p, n) -> M::SMatrix for the out-of-place version or J(M, u, p, n) for the in-place version acting in-place on M. in both cases M is a matrix whose columns are the deviation vectors.
The bold part does not look consistent with what I can see in the source code and in the examples. Function J seems to return the Jacobian matrix J_f(x). The time evolution of the deviation vectors is then obtained multiplying the previous vectors by this matrix:
\frac{dY}{dt} = J_f(x)\cdot Y
If J was what the documentation says, then it would be possible to provide the evolution in tangent space as any type of function. I was planning to exploit this to speed up my computations. However, I don’t think that this tangent evolution works as it is explained in the docs.
Could anyone clarify whether my interpretation is right? If so, the docs should be modified…
What do you mean? What do you see in the docs that is inconsistent?
I don’t understand what you are proposing here. Can you give an example? What does “any type of function” mean? You hope to provide a function f(Y, u, p, t) that returns directly dY/dt ? If so, why?
Evolution in tangent space works as explained in the docs, which is the equation you state in your post. I am very confident on that. Why do you think it is wrong?
Oh sh*t you are right, the docstring doesn’t make sense. Of course the Jacobian function should return the Jacobian matrix, not a matrix with columns the deviation vectors… The jacobian matrix multiplies the latter matrix. that docstirng statement is indeed dead wrong Do you mind doing a quick PR to fix it?
coincidentally, two weeks ago, I was nearly about to re-write the continuous time version for Tangent systems but I decided I better focus on my job instead… But a rework may also allow more flexibility in how to specify the Jacobian.
Just did it. I hope I did it well, it is my first PR ever…
I don’t understand what you are proposing here. Can you give an example? What does “any type of function” mean? You hope to provide a function f(Y, u, p, t) that returns directly dY/dt ? If so, why?
Exactly, I am thinking on providing a function that gives the evolution of Y directly, like what you said. The idea would be to speed up the computations. I am computing the Lyapunov spectra of a very large system in which the Jacobian is very sparse, but sparse does not improve the performance. I was hoping that something hard coded would help, as I can exploit some properties of the system and avoid multiplying by so many zeros.
I understand that this is a very particular case, and the performance might end up not improving. Anyhow, I will try to do this by redefining TangentDynamicalSystem so that tangentrule is my new function.
That’s a good point, I did not benchmark in detail before. I just ran a simulation using Julia’s profiler and it does look like step! is the bottleneck (I attach a screenshot of ProfView).
Ok, for the moment I will try to speed things up by tweaking the implementation in tangent space, but in the future I might want to have a deeper look.
Do you mind doing a quick PR to fix it?
