the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 15 Nov 2023
Computation of covariant lyapunov vectors using data assimilation
Abstract. Computing Lyapunov vectors from partial and noisy observations is a challenging problem. We propose a method using data assimilation to approximate the Lyapunov vectors using the estimate of the underlying trajectory obtained from the filter mean. We then extensively study the sensitivity of these approximate Lyapunov vectors and the corresponding Oseledets' subspaces to the perturbations in the underlying true trajectory. We demonstrate that this sensitivity is consistent with and helps explain the errors in the approximate Lyapunov vectors from the estimated trajectory of the filter. Using the idea of principal angles, we demonstrate that the Oseledets' subspaces defined by the LVs computed from the approximate trajectory are less sensitive than the individual vectors.
- Preprint
(14033 KB) - Metadata XML
- BibTeX
- EndNote
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2023-2168', Anonymous Referee #1, 10 Dec 2023
The manuscript investigates the estimation of the Lyapunov exponents from noisy system trajectory. The noisy trajectory is generated in two approaches: 1) posterior analysis mean of an ensemble Kalman filter by assimilating partially observed noisy observations. 2) direct perturbation of the model trajectory at observation time without data assimilation (DA).
The manuscript presents the changes of Lyapunov vectors and Oseledets’ subspace when erroneous trajectory is used compared to a reference (true) trajectory in these two scenarios. The problem is worth investigating. However, there are major concerns about the manuscript:
1. The manuscript claims it presents a new algorithm to compute the covariance Lyapunov vectors using partial and noisy observations. I find the claim troublesome for the following reasons: 1) the method relies on an underlying model to reconstruct the full system state resulting an algorithm that is exactly the same as Ginelli et al (2013); This is equivalent to using a 'surrogate model' to compute LVs; 2) although experiments use partial and noisy observations, the error of the model itself, comes entirely from model instability arising from inaccurate initial conditions. That is, it is doubtful that the approach can behave well if the model itself is erroneous/biased. Hence, I suggest the authors remove this claim and rephrase the objectives to the differences of the LV between true trajectories and the one estimated by DA.
2. I am in general in support of using first-person pronouns for scientific writing. However, I feel using ‘we’ in Section 2 is slightly misleading. For example, Section 2.2 is describing the algorithm proposed by Ginelli et al (2013) and using ‘we’ can give readers a feeling that it is the approach proposed by the authors.
3. The results show that the estimated trajectory works better for the first and the last LVs than other LVs and the Oseledets’ subspace is aligned for most dimensions. Based on Bocquet et al, 2017, the error of EnKF should converge to the unstable-neutral subspace. Similar conclusions were found by some other studies, for example, Chen, Y., Carrassi, A., and Lucarini, V.: Inferring the instability of a dynamical system from the skill of data assimilation exercises, Nonlin. Processes Geophys., 28, 633–649, https://doi.org/10.5194/npg-28-633-2021, 2021. ; Dan Crisan, Michael Ghil; Asymptotic behavior of the forecast–assimilation process with unstable dynamics. Chaos 1 February 2023; 33 (2): 023139. https://doi.org/10.1063/5.0105590;
Based on this line of findings, can the authors give better explanations for your results? Is it because that when the Ginelli et al (2013) is used, the errors in the mean state are projected to stable components of the LVs as well? If the observations are not partial but just noisy, will the results change?
4. I feel I don’t fully understand the principal angle metric used by the authors. Consider the BLVs are orthonormal, based on Eq 14, is the principal angles just a reorder of the angle between BLVs? Can authors give a better explanation for this?
Minor comments:
1. L97, I’m not sure if we should call f a vector field as I think it is a map.
2. L147, I feel the statement that ‘the dimensions of i-th forward and backward Oseledets subspace S^+_i and S^-_i are, …, is n+d_i’ needs some revision. Based on Eq. (5), dim(S^+_i) = n, and if dim(S^+_i) = d_i, the sum of these dimensions must be > n + d_i. That said, I think the symbol d_i is not defined clearly. This is only clear once d_i is, for example, for the dimension of S^+_i\S^+_i+1, I think.
3. In Fig 4, I think there is no inverse sign in B_j+lR_j, j+1
4. Equation 11 can be better explained by saying D is local growth rate.
5. L225, ‘This leads to a our proposed modification ….’ I don’t need a is needed.
6. L228, ‘with B_j as the initial condition for (2)’; When you compute CLVs, do you use a new matrix or is B_j simply the ensemble perturbations in the EnKF?
7. L234, isn’t l2 and RMSE the same in this context?
8. Eq. 14, what is y_k here?
9. L327, if only y is observed, I believe H = [0, 1, 0]^T instead of [1, 0, 1]
10. L328, what is x_0?
11. L343, are 25 ensemble members still being used for L96? What is the ensemble size being used? The ensemble size depends on the size of the unstable-neutral space, but not necessarily the system dimension.
12. Is it possible to give the RMSEs for each variable in L63 model?
13. L390, this instead of his
Citation: https://doi.org/10.5194/egusphere-2023-2168-RC1 -
RC2: 'Comment on egusphere-2023-2168', Anonymous Referee #2, 15 Dec 2023
This manuscript investigates the sensitivity of Lyapunov vectors computation (using the QR and Ginelli algorithms) with respect to perturbations introduced along trajectories of chaotic dynamical systems. These algorithms are used respectively to compute Backward Lyapunov Vectors (BLVs) and Covariant Lyapunov Vectors (CLVs), and the present work studies how the sought Lyapunov vectors and Osedelets subspaces are impacted by the perturbations.
The perturbations considered are of two kinds:
- obtained through an ensemble Kalman filter, filtering observational errors introduced with respect to a reference trajectory.
- obtained by perturbing directly the said reference trajectory by Gaussian white noise.
This research question is interesting and has many "real world" applications, however the present analysis is rather shallow, the methodology is inadequate, and the presentation of the results suffers from many problems. Before explaining all of this in more detail, I must say directly that my recommendation is to do a very major revision of the paper, which is not suitable for publication in the present state.
The manuscript is actually at the limit of needing a complete resubmission, and it is up to the authors to reflect on whether this is needed.
At the end of this report I propose some further studies and improvements which could lead to a suitable manuscript for the next revision, at least in my opinion.
# General comments
In general the manuscript was quite difficult to read, and I had to go back and forth constantly (a bit of "back and forth" is ok of course).
## Title and abstract
- The title of the manuscript itself is misleading since it tends to indicate that this is a new method to compute LVs, while in fact state-of-the-art methods are used and simple sensitivity analysis about them are performed. Something like "Sensitivity analysis of Lyapunov Vectors computation with respect to perturbations" would be more accurate.
- Again, in the abstract, it is stated that
"We propose a method using data assimilation to approximate the Lyapunov vectors using the estimate of the underlying trajectory obtained from the filter mean."
and this is misleading since the method to compute the LVs are the well-known and unmodified QR and Ginelli algorithms, the core of the study here being the sensitivity analysis. Thus this abstract needs to be clearer about what the subject of the article really is.
## Introduction
- In the introduction, it is asserted that the used trajectories are not shadowing trajectories. Could you please cite some works exploring this aspect in data assimilation (DA)? As stated, the noise is injected at discrete observation times, therefore in between the evolution is deterministic. How can you conclude that the portion of trajectories in between observation times are not shadowing?
## Section 2
- Section 2.1 : You do not mention or explain the problem of degeneracies of the Lyapunov exponents. Also, related to that, please define clearly the dimension d_i of the Osedelets subspaces. This may help the reader understand the dimension of the intersection of the Osedelets subspaces.
- Line 138: "see the review in Kuptsov..."
- Line 146: "only the latter of the two limits above."
- Eqs. 5 and 6: it should not be the Greek letter phi but the non-empty symbol instead.
- Section 2.2: This is a rather tedious section to explain an already well-documented algorithm in the literature. Therefore I would move most of this part to an appendix, keeping only a qualitative sketch in order to explain the following experiments.
- Section 2.2: Please use capital letters for "Gram-Schmidt".
- Section 2.4: A general problem with the methodology of the study is that it is difficult to grasp what the different amplitudes of noise used to perturb the trajectories represent with respect to the typical variability of the model. Also, using the same amplitude for each variable - assuming they have the same variability - may bias the results. One way to solve this would be to express the perturbation in each variable as a percentage of its variability, and then vary this percentage, and not the absolute noise amplitude.
- Section 2.4: Also, I would move section 2.4 to after section 3.2, since this one is more about the study methodology than about the methods used. (I got lost looking for the details about the perturbations while reading the results section 4.2.).
- Section 2.4 : There is a risk of confusion between the x_k, y_k ... and the symbol used in the model's equations (x, y, z). Similarly, sigma can represent the noise amplitude or is a parameter of the L63 model. Notations need thus to be revisited over the whole manuscript.
- Eq. 13 : Please specify explicitly the symbols meaning, like the distribution you are using and I_d.
- Line 246: Why do you use bold I for the identity here, it is the only place where it is noted this way.
- Eq. 14 : y_k should be q_k
- Lines 293-297: It should be specified that subspaces hence constructed approximate the Osedelet subspaces.
## Section 3
- Section 3.2: Is that correct that you assimilate every Delta t timeunit? If so, you should specify it, by saying that you observe AND assimilate.
- Section 3.2: Again, here, mu should be expressed as percentage of the variability of the concerned variables, otherwise it is difficult to estimate what this quantity means with respect to the system at hand.
- Line 327: With respect to equation 15 order, H should be [0,1,0].
## Section 4
- Figure 2: There is a problem with the left panel. BLV 2 curve should be BLV 1 (compare with right panel CLV 1, it should be CLV 1 = BLV 1). Actual results for BLV 2 seem to be missing. Also the RMSE of the analysis displayed on top of the plots as a scale is a very bad practice, please find another way to represent it (for example you can put the numbers on top of higher 75th percentile).
- Line 370 to 372: The phrase "We first note that since BLV are orthonormal, two of these angles are necessarily equal which happen to be those between the second and third BLV and they are quite small even for the largest observational noise strength we have used." is difficult to understand if not false. Figure 2 shows the angles between reference and perturbed vectors, not the angle between vectors of different rank. This sentence needs to be rephrased. Also what does "quite small" mean here? A more general comment is that the left panel of Figure 2 shows that comparing BLVs in such a constrained 3-dimensional system does not make a lot of sense. Also, in rather specific systems like L63, it is well possible that the results you obtain here are not generic, with very low sensitivity compared to other higher-dimensional models.
- Line 372 to 373: "In addition we note that the median of angle between the first - most unstable - BLV is also within 15 degrees and does not increase rapidly with the observation noise strength µ." Again, this sentence makes sense only if BLV 2 is BLV 1 on Figure 2.
- Figure 4: It should be specified somewhere that i is a number indexing the principal angles.
- Figure 6: This shows why using L63 for sensitivity analysis of the stability directions of trajectories lead to non-generic very low sensitivity, with large contiguous regions of the attractor where the stability changes smoothly and slowly. In my view this is why studying L63 in the present framework is not very informative.
- Figure 7: Suffer from the same problem as Figure 2 with BLV1 being represented as BLV2, and actual BLV 2 is missing. Also, the x axis scale is wrong. Also, it is not specified what the insets are representing.
- Line 412-413: "The 1 st and 2 nd BLV have the same rate, whereas the rates are different for the 1st and 2nd CLVs." It is impossible to verify with the BLV 2 results missing.
- Line 413-415: The plot of the exponents seems to be missing. Therefore it is difficult to understand this statement about the order 0.2 of the exponents.
- Lines 429-432: It is concluded that the Lyapunov spectrum is quite robust to perturbations of the trajectory, but this is not surprising since there is no model error and the Lyapunov spectrum is a global property of the system. What would have been more interesting to see here is the impact of the perturbations on the local exponents.
- Line 434-435: 'instead of the individual vectors' what does it mean? Please rephrase.
- line 445: 'ture' should be 'true'
- Figure 9: Again, i should be identified as the index of the principal angle. Also, I would collapse both panels into one figure for a proper comparison.
## Conclusion
- Line 475: 'quanify' -> 'qualify'
-----------------------------------------------------------------------------
# Explicit reasons for major revision
Beside the problems of the manuscript and of the figures (but which could be improved), I recommend a complete major revision of the paper for the following reasons:
## The study is not comprehensive enough
The only message that the paper finally conveys is that the Osedelets subspaces spanned by the BLVs are less sensitive to perturbation than particular stability directions, which can intuitively be understood as “it is harder to perturb volumes than directions”. Therefore the results obtained in this paper are not surprising. More precisely it is stated starting line 395 that:
"embedded within any high dimensional Oseledets subspace, there are lower dimensional subspaces which are close to the true [BLVs] subspaces. In order to understand this behaviour more clearly, we now study the dependence of this approximation on the strength of perturbations of the trajectories."
The ensuing perturbation study does not give a clue on the characterization of the said embedded subspaces but rather confirms that perturbing (thus apparently in any way) the trajectories perturb the stability directions a lot, but less the volume spanned by the BLVs. In my opinion, the study about the principal angle is not enough. This paper is not about introducing and testing a new method, so what is left is a rather simple sensitive analysis in two very idealized models which confirms and quantifies an intuition.
## Methodology
What could have been interesting is to compare the two perturbation methods, because as acknowledged in the conclusion, investigating the sensitivity conditional on the stability directions of the perturbations is interesting. But EnKF is already doing this in part, because in perfect chaotic models the analysis error covariance matrix converges to the unstable-neutral subspace. So making the two experiments more comparable could have already led to a result about the impact of the structure of the perturbations.
# Way to improve and resubmit the manuscript
- I would regroup the results per model, and not per perturbation method. Since L63 is not very interesting for the purpose of the present study, therefore I would keep this part as small as possible.
- I would suggest to try to characterize the subspaces which remain closely aligned under perturbation. One way to do this would be to project the true BLVs or/and CLVs on the perturbed Osedelets subspaces, and compute thus the angles between the true BLVs or/and CLVs and these subspaces, but there might be other more insightful methods.
- I would investigate what happens if perturbations are done in particular stability directions/subspaces.
- I would also make the two already performed experiments more comparable (and do this comparison per model), such that a comparison between Figures 4 and 9 is possible.
- I would also investigate other models, closer to geophysical applications, to see if some genericity of the results can be found.
Because this might require a lot of additional work, I recommend at least a very major revision.
Citation: https://doi.org/10.5194/egusphere-2023-2168-RC2
Model code and software
Mr-Markovian/Reconstruct_CLV_via_enkf: Published version of the code Shashank Kumar Roy https://zenodo.org/record/8396549
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 116 | 49 | 11 | 176 | 8 | 8 |
- HTML: 116
- PDF: 49
- XML: 11
- Total: 176
- BibTeX: 8
- EndNote: 8
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1