feat: Add Covariant Lyapunov Vectors (CLV) via Ginelli algorithm #357
Conversation
Implement the clv function for computing Covariant Lyapunov Vectors using the forward-backward Ginelli algorithm from Ginelli et al., Phys. Rev. Lett. 99, 130601, 2007. Features: - Works with both discrete and continuous dynamical systems - Supports in-place and out-of-place system definitions - Uses TangentDynamicalSystem for automatic or hand-coded Jacobians - Returns CLVs, Lyapunov exponents, trajectory states, and timestamps Implementation details: - Forward pass: QR decomposition with positive-diagonal sign correction - Backward pass: Back-substitution to obtain coefficient matrices - Proper orthonormalization (signed-corrected QR) to ensure CLV uniqueness Tests: - Diagonal system test verifying CLVs align with coordinate axes - Lyapunov exponent consistency with lyapunovspectrum - Covariance property verification for Henon map - Neutral CLV alignment with flow direction for Lorenz system Documentation: - Comprehensive docstrings with algorithm description - Convergence troubleshooting notes - Visualized examples with Henon 2D and Lorenz 3D attractors
|
Hello @hinsley thank you very much for the contribution. I've been away from GitHub for a while but I am back now and I will review your PR soon, or as s oon as I can! May I also ask how you use this package or thi new functionality? Perhaps there is an opportunity to give a talk at an upcoming JuliaDynamics meeting. |
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #357 +/- ##
===========================================
+ Coverage 54.86% 65.35% +10.49%
===========================================
Files 22 27 +5
Lines 904 1146 +242
===========================================
+ Hits 496 749 +253
+ Misses 408 397 -11 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
|
Sorry, I was distracted for a few days.
I mainly use this algorithm in parameter sweeps of maps & flows to look for (pseudo)hyperbolic attractors. A good example of how this can be done with CLVs is found in this 2018 paper, though they employ some tricks related to those from the Kuptsov & Parlitz paper to avoid explicitly computing vectors, and they're only looking at hand-selected parameter values rather than doing sweeps. The basic idea is to consider some stable and unstable tangent subspaces spanned by a certain number* of leading and lagging CLVs (corresponding to greatest/least Lyapunov exponents) and ask whether these subspaces intersect nontrivially anywhere, or alternatively whether the minimal angle between the subspaces stays bounded away from zero. If there's a nontrivial intersection, you have a violation of pseudohyperbolicity and homoclinic tangencies indicating lack of structural stability (pseudohyperbolicity guarantees structural stability). In the additions to the docs for this PR I included a couple of visualizations of the CLVs so one can see what these subspace intersections look like in practice. *For example, in the Hénon map (2D), there's just one leading CLV, corresponding to a positive Lyapunov exponent and indicating the direction of greatest expansion, and one lagging CLV, with a negative exponent and indicating greatest contraction. In Lorenz, you've got a positive exponent, a zero exponent, and a negative exponent. In general (see the paper linked above) the "contracting subspace" needs to contain only CLVs associated with negative Lyapunov exponents, so we take 2 leading CLVs to span the "center-unstable subspace" One could think of the minimum angle between contracting and expanding subspaces as a kind of degree of hyperbolicity of an invariant set, but there are some technicalities to keep in mind about non-uniqueness of dimension of these subspaces if the codimension of the orbits is greater than 2 (maps dim > 2, flows dim > 3). I have seen these CLVs used in climate sciences in more sophisticated ways, though I believe they're working in higher-dimensional spaces where you maybe don't want to use the Ginelli algorithm. Presumably it is first and foremost used as a dimensionality reduction technique for forecasting, since you can calculate errors coming just from the center-unstable bundle instead of the full tangent bundle. I don't know anything about this area though, so this is just a guess. I am pretty interested in also understanding the Kuptsov & Parlitz approach to CLVs, which is what it seems like the applied dynamicists tend to reach for. There's an existing Python implementation at https://github.com/ThomasSavary08/Lyapynov I'd probably study alongside the original paper. If the JuliaDynamics talks are intended to be more than 5-10 minutes, I think I would feel more prepared waiting until I had also played with uses of the K&P algorithm and maybe contributed it to this repo in another PR. If there's any interest in seeing how to use the CLV algorithm from this PR to study degree of hyperbolicity in dynamical attractors, I could give a talk. But this would really feel more like a brief presentation of a special method rather than a survey of the uses of CLVs. Up to you (or anyone else interested). |
|
Hello, sorry for being late I am swamped! I think this sounds generally interesting and if you want to chat about this you are welcomed to do so. Feel free even to pick a slot and write your name directly in this sheet here: https://docs.google.com/spreadsheets/d/19MrgrHdhy6r1x8OUS-B3Gj1WEMIAToCNRabBimCzuCk/edit?gid=1016767766#gid=1016767766 ! |
|
I'll now go ahead to finally review this PR and sorry for being late. |
There was a problem hiding this comment.
This is a fantastic contribution. Thank you very much. It is rare to see this level of quality for the first contribution of someone! Congratulations! I only left minor comments on the PR. I did go through the source code, and all performance fundamentals of type stability and unecessary allocations seem to be covered. Did you by any chance also do some profiling yourself, or something similar, to optimize the performance of the code? EDIT: I read in the first post:
Type stability verified with @code_warntype.
Docs all are great.
Besides the comments below, can you also add changelog entry and bump minor version in Proejct.toml?
Lastly, I am sorry but I need to ask this in this day and age. Are parts of this code generated by GenAI? If so, this will be the very first code in DynamicalSystems.jl that has this property, and so we would need to put a disclaimer somewhere.
| Covariant Lyapunov vectors (CLVs) are tangent vectors of a dynamical system corresponding to the directions of contraction and expansion described by the Lyapunov exponents. In 2007, Ginelli et al. [Ginelli2007](@cite) introduced a numerical method for computing CLVs, based on the Benettin method used for computing Lyapunov exponents. We do not currently implement the alternative algorithm from the 2007 Wolfe and Samelson paper. | ||
|
|
||
| Note that the most non-central directions converge first, and the least expanding/contracting directions converge last. In particular, the most contracting direction converges first in the forward pass, and the most expanding direction converges first in the backward pass (both passes are computed in `clv()`). To see if you're getting the correct CLVs, you can check the Lyapunov exponents of the computed CLVs against a much longer trajectory using `lyapunovspectrum`, and in continuous-time systems you can check that the CLV corresponding to the flow direction is parallel to the flow direction at each point in the trajectory; if not, try the following: | ||
|
|
||
| - Use a different integrator | ||
| - Increase the length of forward or backward transients | ||
| - Decrease the timestep size | ||
| - Increase the orthonormalization/QR frequency by decreasing `Δt` |
There was a problem hiding this comment.
The docs of ChaosTools.jl are admittedly the olders in JuliaDynamics and do not follow many of the improvments done... For now, can you please move such discussions to the ## Description section inside the documentation string of the function?
There was a problem hiding this comment.
Additionally, can you please use some tooling that limits lines of source code to 92 when you do port this to the source code?
| For continuous systems this is approximate. | ||
| * `show_progress = false`: Display a progress bar. | ||
|
|
||
| ## Algorithm |
There was a problem hiding this comment.
| ## Algorithm | |
| ## Description |
| The CLVs have the key property of being *covariant*: if V(t) is a CLV at time t, | ||
| then after evolution by the tangent dynamics, the resulting vector is parallel to V(t+Δt). |
There was a problem hiding this comment.
I am a bit confused reading this. I am not sure where "the resulting vector" alludes to. After evolution by tangent dynamics, V(t) has become V(t+Δt). So the resulting vector is exactly V(t+Δt), not just parallel to it. I am misundersrtanding something.
| ## Returns | ||
|
|
||
| A named tuple `(V, λ, x, t)` where: | ||
| - `V::Vector{Matrix}`: CLVs at each stored time step. `V[i]` is a `D×k` matrix whose columns | ||
| are the `k` CLVs at time `t[i]`, ordered from most expanding to most contracting. | ||
| - `λ::Vector`: Lyapunov exponents (growth rates along CLVs), ordered from largest to smallest. | ||
| - `x::Vector`: State vectors at each stored time step. | ||
| - `t::Vector`: Time stamps corresponding to each stored CLV snapshot. |
There was a problem hiding this comment.
I would challenge the named tuple return type as being an uncommon practice. I would favour returning a tuple, and not saying that explicitly as that is what all functions do. So do not say anything, just start wit hthe bullet point directly.
| Use this method for looping over different initial conditions or parameters by | ||
| calling [`reinit!`](@ref) on `tands`, or when you have a hand-coded Jacobian. | ||
|
|
||
| ## Returns |
There was a problem hiding this comment.
| ## Returns | |
| ## Return |
| b = Statistics.covm(x, mx, y, my) / Statistics.varm(x, mx) | ||
| a = my - b * mx |
There was a problem hiding this comment.
| b = Statistics.covm(x, mx, y, my) / Statistics.varm(x, mx) | |
| a = my - b * mx | |
| b = Statistics.covm(x, mx, y, my)/Statistics.varm(x, mx) | |
| a = my - b*mx |
| # Creation of a trivial system with one coordinate going to infinity | ||
| # and the other to zero. Lyapunov exponents are known analytically | ||
| trivial_rule(x, p, n) = SVector{2}(p[1]*x[1], p[2]*x[2]) | ||
| trivial_rule(x, p, n) = SVector{2}(p[1] * x[1], p[2] * x[2]) |
There was a problem hiding this comment.
Please revert all automatic formatting done in this file.
|
oh additionally the |
Store snapshots of current_state with recursivecopy so x_history does not alias a mutable view and collapse to the final state.
| end | ||
| end | ||
|
|
||
| @testset "CLV Lyapunov consistency" begin |
There was a problem hiding this comment.
Hm, I took a closer look at the tests after your latest commit. I have two comments: 1) there should be a test that access the history of the CLVs, not just the final point, so that we would have caught what you just fixed. 2) can we reduce the computational cost of the tests? At the moment to test each property you are interested in, the whole CLV estimation is run. But I'd argue for most properties you can run CLVs once and then check each property in an inner testset.
3 participants
Add Covariant Lyapunov Vectors (CLV) via Ginelli algorithm
This PR implements the
clv()function for computing Covariant Lyapunov Vectors using the forward-backward Ginelli algorithm (Ginelli et al., Phys. Rev. Lett. 99, 130601, 2007; arXiv).What are CLVs?
Covariant Lyapunov Vectors are tangent vectors over a trajectory that characterize the local stability directions of a dynamical system (i.e., they are associated to the Lyapunov exponents). Unlike Gram-Schmidt vectors from the standard Lyapunov exponent computation, CLVs are covariant under the dynamics: if V(t) is a CLV at time t, then after evolution by the tangent dynamics, the resulting vector is parallel to V(t+Δt).
API
Features
DeterministicIteratedMap) and continuous (CoupledODEs) systemsTangentDynamicalSystemfor automatic or hand-coded JacobianslyapunovspectrumconventionsTests
lyapunovspectrumfor Hénon map and Lorenz systemDocumentation
Added a new "Covariant Lyapunov Vectors" section to
lyapunovs.md("Chaos detection" in the docs) with:Performance
Benchmarks for the documentation examples (on Apple M4):
clv(ds, 3000; Ttr=1000, Ttr_bkw=200)clv(ds, 500; Ttr=1000, Ttr_bkw=1000, Δt=0.05)Type stability verified with
@code_warntype.