Multi-variate factorisation of numerical simulations

This paper presents a meaningful step forward and an important clarification to the methodology of (paleo)climate factorisations. As such I would be happy to see it published following some small revisions. I found the introduction of the synergy term in the Lunt et al (2012) factorisation a little bewildering at first. The manuscript didn’t seem to contain an explanation as to why there is no need for a synergy term in the two-dimensional case. Is it because the averaging process effectively removes one degree of freedom (dimension)? This study is well written, and very interesting. It presents a new development to Lunt et al., (2012) investigation of the factorization method used in paleoclimate modeling. The goal is to achieve completeness, uniqueness and symmetry of the factorization, and eliminate the synergy term. Mathematically, it is a fine solution. But I do have a couple questions: The completeness is achieved by averaging out different additive paths of applying climate forcings (eq. 6) or sharing the synergy term proportionally with the generated warming by the individual “forcing” factors (eq. 15). I actually quite like the synergy term, because the synergy term has the physical meaning of capturing the non-linear effects of changing vegetation, ice sheet, topo/geography and CO2. Abstract. Factorisation is widely used in the analysis of numerical simulations. It allows changes in properties of a system to be attributed to changes in multiple variables associated with that system. There are many possible factorisation methods; here we discuss three previously-proposed factorisations that have been applied in the ﬁeld of climate modelling: the linear factorisation, the Stein and Alpert (1993) factorisation, and the Lunt et al. (2012) factorisation. We show that, when more than two variables are being considered, none of these three methods possess all three ::: four : properties of ‘uniqueness’, ‘symmetry’, 5 and ‘completeness’, ::: and ::::::: ‘purity’. Here, we extend each of these factorisations so that they do possess these properties for any number of variables, resulting in three factorisations – the ‘linear-sum’ factorisation, the ‘shared-interaction’ factorisation, and the ‘scaled-total’ factorisation. We show that the linear-sum factorisation and the shared-interaction factorisation reduce to be identical :: in ::: the :::: case :: of :::: four :: or ::::: fewer ::::::::: variables, ::: and ::: we ::::::::: conjecture ::: that :::: this ::::: holds ::: for ::: any ::::::: number :: of :::::::: variables. We present the results of the factorisations in the context of studies that used the previously-proposed factorisations. This reveals that only the 10 linear-sum/shared-interaction factorisation possesses a fourth ::: ﬁfth : property – ‘boundedness’, and as such we recommend the use of this factorisation in applications for which these properties are desirable.


We agree that the synergy/interaction terms can be interesting and useful in some cases. We state this in the manuscript: "These interaction terms are important and interesting, but there are cases where it can be useful or essential to only include attributable terms in the factorisation.". See lines 125-126.
Beyond what is shown in Fig. 2, additional axes are needed to capture these non-linear effects. In fact, the points seemingly overlapping at T111 could be a visual illusion, and separable with additional axes of non-linear interactions. Wouldn't it be better to leave the synergy terms alone?
In Figure 2, the lengths of the various lines do not represent the magnitude of the individual temperature changes. It would require a 4 th dimension show the temperature variable (challenging on a 2-dimensional pdf!). Similarly, in Figure 1 the temperature itself can be considered as a surface in a third dimension sitting above the 2-dimensional plane of CO2 and ice sheet. We clarified this in the caption of Figure 1. However, even if we could represent this fourth dimension, then it would still be the case that the three lines in the top-right corner of Figure 2 would converge on a single point, T111, because T111 is single-valued. Taking two example "routes" from the bottom-left (T000) to the top-right (T111), the overall change is, for one route: (T100-T000)+(T101-T100)+(T111-T101). For another route, it is (T010-T000)+(T011-T010)+(T111-T011). For both of these, terms cancel and they reduce to T111-T000. This is a mathematical property, not one associated with the climate system or the presence or absence of physical nonlinearities.
Additionally, Fig. 4 is a natural result of absorbing residuals into the calculation with corrections. I am worried that the physical meaning of factorization is contaminated through this absorption instead of being enhanced.
The differences between the left, middle, and right columns are related to the size of the synergy terms. With larger synergy terms (strong interactions between different factors), then the three columns would be more different to each other. The fact that the three columns look very similar is because actually, in this case, the synergy terms are relatively small. Added: "The first thing to note is that the three factorisations all have very similar results; visually it is difficult to tell them apart on a regional scale, and they result in global means for each factor that differ by less than 10%. This is because, in this example, the non-linearities (i.e. the interaction terms) are relatively small." See lines 289-292. For LGM, whether the symmetry is a feature of climate response to CO2 forcing is questioned (e.g., Zhu et al., 2020, Clim. Past. https://doi.org/10.5194/cp-2020. This study shows that changes in CO2 and ice sheet have different forcing efficacies under the LGM and preindustrial climate conditions. Similarly, asymmetric vegetation, ice sheet, and CO2 forcings might be prevalent for past climates. Would it be more useful to use the proposed framework to understand the asymmetry of climate forcings and responses instead of trying to force symmetry, which might not be a real feature in climate system? We agree that there may be (in fact, almost certainly is!) an asymmetry in forcing and/or response of the climate system when considering the transition from preindustrial to a warm climate (e.g. the Pliocene), compared to a transition from the preindustrial to a cold climate (e.g. the LGM). This is a physical property of the climate system, related to the nonlinearity and/or state dependence of forcings and feedbacks. However, the "symmetry" that we discuss and define in this paper is not a physical property, but a mathematical one. It is the symmetry between the transition from preindustrial to LGM versus LGM to preindustrial. This is a mathematical requirement (A-B = -(B-A) !) . Here we "force" this symmetry, which we believe is the correct thing to do from a mathematical viewpointit does not preclude analysis of the type of physical "asymmetry" that the author refers to (B-A  -(A-C)). In fact, the reviewer is basically saying that (in our 2-dimensional) case, T10-T00  T11-T01 ; which is certainly true in general, and fully consistent with our factorisation methodologies.
Lastly, this framework is described in the context of LGM, showing the LGM results would be more consistent.
Correspondence: Daniel J. Lunt (d.j.lunt@bristol.ac.uk) Abstract. Factorisation is widely used in the analysis of numerical simulations. It allows changes in properties of a system to be attributed to changes in multiple variables associated with that system. There are many possible factorisation methods; here we discuss three previously-proposed factorisations that have been applied in the field of climate modelling: the linear factorisation, the Stein and Alpert (1993) factorisation, and the Lunt et al. (2012) factorisation. We show that, when more than two variables are being considered, none of these three methods possess all three 'purity'. Here, we extend each of these factorisations so that they do possess these properties for any number of variables, resulting in three factorisations -the 'linear-sum' factorisation, the 'shared-interaction' factorisation, and the 'scaled-total' factorisation. We show that the linear-sum factorisation and the shared-interaction factorisation reduce to be property -'boundedness', and as such we recommend the use of this factorisation in applications for which these properties are desirable.
Copyright statement. TEXT

Previous factorisation methods
In order to introduce and discuss previous factorisation methods, we use an example case study from the field of climate science.

The linear factorisation
The simplest factorisation that can be carried out is a linear one. For the LGM ::::::: Pliocene : example with 2 factors, 3 simulations 60 are carried out in which variables are changed consecutively; for example, T 00 , T 10 , and T 11 . The factorisation of the total change, ∆T , between contributions due to CO 2 (∆T 1 ) and ice (∆T 2 ) would then be: This factorisation is illustrated graphically in Figure 1(a). However, an equally valid linear factorisation would be and in a non-linear system this would in general give a different answer to Equation 1. In this sense, the linear factorisation method is not 'unique'. However, it is 'complete' in the sense that the individual factors sum to the total change, ∆T exactly, i.e.
∆T 1 +∆T 2 = T 11 −T 00 . Considering the linear factorisation as a 'path' starting at T 00 and ending at T 11 , it is also 'symmetric ', 70 in that if we instead started from T 11 we would retrieve the same numerical values for the two linear factorisations (differing just by a minus sign for the numerical value of each factor). complete.

The Stein and Alpert (1993) factorisation
Stein and Alpert (1993) proposed an alternative factorisation method, illustrated in Figure 1 In contrast to the linear factorisation, the Stein and Alpert (1993) factorisation is unique. It is also complete because ∆T 1 + ∆T 2 + S = T 11 − T 00 (in fact, S is defined such that the factoristion is complete). However,  2012) proposed another factorisation, in which the factorisation for a particular variable is defined as the mean of the difference between each pair of simulations that differ by just that variable. This is illustrated in Figure 1(c); the factorisation of ice is represented by the mean of the two blue lines and the factorisation of CO 2 is represented by the mean of the two red lines: The N = 2 factorisation in Equation 4 is unique, complete, and symmetric ::::::::: symmetric, :::: and :::: pure. It is worth noting that Equation 4 can be interpreted in multiple ways -either (i) as described above, the factorisation averages all the possible pairs of simulations that differ solely by a change in that variable, i.e. for a particular variable it is the mean of either the horizontal or term, S, shared equally between the two factors.
In extending to N = 4 variables, Lunt et al. (2012) assumed that the first of these interpretations would still hold for any 100 number of variables. However, consider the N = 3 case illustrated in Figure 2, in which we have added vegetation as a third variable to contribute to LGM cooling ::::::: Pliocene :::::::: warming. Averaging the edges (interpretation (i) above) would result in a 4 factorisation: Although this is uniqueand symmetric, : , ::::::::: symmetric, ::: and ::::: interesting, but there are cases where it can be useful or essential to only include attributable terms in the factorisation.

Extensions to the previous factorisations
Here we discuss possible extensions to the three previous factorisations discussed above, that are unique, symmetric, pure. For three dimensions, this is illustrated by the dotted lines in Figure   2.
Each possible linear factorisation can be represented as a non-returning 'path' from the vertex T 000 to the opposite vertex 135 T 111 , traversing edges along the way (dotted lines in Figure 2). When considering the sum of all possible paths, some edges are traversed more than others. In general, those edges near the initial or final vertices are traversed more times than edges that are further away from these vertices. As such, when we average the possible linear factorisations, different edges (corresponding to different terms in the factorisation) will have different weightings. This is in contrast to Equation 5 where each term (i.e. edge of the cube) has the same weighting. For three dimensions, Figure 2 shows that the 6 edges adjacent to the initial and final 140 vertex are traversed twice, whereas the 6 other edges are traversed only once. Therefore, the factorisation is : This factorisation is complete (∆T 1 + ∆T 2 + ∆T 3 = T 111 − T 000 ), unique, and symmetric pure.

145
To generalise to N dimensions, consider an N -dimensional cube, which has a total of 2 N vertices and N × 2 N −1 edges.
There are 2 N −1 edges in each dimension. There are N ! paths from the initial vertex of the cube to the final opposite vertex, each of which consists of a traverse along N edges. Therefore, in each dimension there are a total of N ! edges traversed for all paths combined.
As for the 3-dimensional case above, let us label each vertex, V , of this N -dimensional cube as V a1···aN , where each a i is 150 either 0 or 1. A value a i = 0 represents the first value for variable i, and a i = 1 represents the second value for variable i. Each vertex is also associated with a system value, denoted T a1···aN (see Figure 2 for the case N = 3).
All factorisations consist of partitioning the total change, ∆T = T 1···1 − T 0···0 between N factors. Each factor is associated with a dimension, i, in the N -dimensional cube. The factorisation for dimension i is ∆T i .
For the linear-sum factorisation, all paths that we consider start at the origin vertex, 0 · · · 0, and end at the opposite vertex vertices for an edge that is oriented in the ith dimension, so that Y i is related to X i by changing the ith subscript of each 160 element from 0 to 1. For example, for N = 3 and i = 2, Y 2 = {010, 011, 110, 111}. Order X i and Y i so that their elements correspond. Then we write X j i to indicate the jth element of X i , and Y j i as indicating the jth element of Y i . For example, for the X 2 defined above, X 3 2 = 100.

6
The Lunt et al. (2012) factorisation averages along each edge oriented in dimension i: For the linear-sum factorisation, we instead carry out a weighted average, with the weight for each edge in dimension i given by how many times it is traversed in all N ! paths. Consider all the paths that traverse an edge which starts at a vertex defined by k subscripts of '1' and N − k subscripts of '0'. There are k! possible paths to the start of this edge, and (N − k − 1)! paths from the end of this edge to the final corner (defined by N subscripts of '1'). Therefore, there are k! ×(N − k − 1)! paths that use this edge. As such, we can write the linear-sum factorisation as: where k j i is the number of subscripts of '1' in X j i . For example, for N = 4 and i = 1, we have N ! = 24 edges traversed in this dimension, and 2 N −1 = 8 edges. mensions. For consistency, we use the same notation as (Stein and Alpert, 1993) ::::::::::::::::::: Stein and Alpert (1993). In their notation,f 1 represents the difference between a simulation in which only factor i is modified with a simulation in which no factors are modified, andf ijk··· represents interaction terms between the different factors. For example, for our original N = 2 example illustrated in Figure 1 and given in Equation 3, ∆T 1 ≡f 1 , ∆T 2 ≡f 3 :::::::: ∆T 2 ≡f 2 , and S ≡f 12 .

Pliocene
: example for N = 3,f 12 is the interaction term (i.e. the synergy) between factors 1 and 2 (CO 2 and ice),f 13 is the interaction term between factors 1 and 3 (CO 2 and vegetation),f 23 is the interaction term between factors 2 and 3 (ice and vegetation), andf 123 is the interaction term between all three factors. In this case, Stein and Alpert (1993) give that ∆T =f 1 +f 2 +f 3 +f 12 +f 13 +f 23 +f 123 pure (because we are just re-partitioning the interaction terms). It turns out that it is also symmetric. For example for CO 2 , This factorisation for N = 3 is represented visually in Figure 3(a). Equations 10 and 11 give that, for CO 2 , 210 This is identical to the equivalent term in Equation 6, indicating that the shared-interaction and linear-sum interpretations are identical for N = 3, and that therefore for N = 3 the shared-interaction factorisation is unique, symmetric, :::: pure, : and complete. Stein and Alpert (1993) give the generalisation of their factorisation to N factors (their Equations 11-16). For N = 4, the interaction terms are shared so that, for example for CO 2 , This factorisation for N = 4 is represented visually in Figure 3(b). Again, for N = 4 this is the same as the :::::: identical ::: to ::: the linear-sum interpretation (Equation 9). We conjecture that for any N these two interpretations will give identical results. In the scaled-total factorisation, the Lunt et al. (2012) factorisation is modified so that it is complete :::: (and ::::::: remains :::: pure). This is 220 achieved by taking the total residual term required for completeness(the 'synergy', S in the sense of (Stein and Alpert, 1993) R, is defined such that where the ∆T ′ i are defined in Equation 5. We then share the synergy ::: this ::::::: proportionally across the three factors, such :::::::::::::::: :::::::::::::::: :::::::::::::::: Equations 14 and 15 reduce to: This shows that this factorisation can also be interpreted as simply scaling the Lunt et al. (2012) factorisation so that the sum of the factors equals T 111 − T 000 . In N dimensions, this generalises to: where ∆T ′ i is defined in Equation 7. For example, for N = 4 and i = 1 we have: and similarly for ∆T ′ 2 , ∆T ′ 3 , and ∆T ′ 4 .

9
Here we discuss three examples of papers in which the Lunt et al. (2012) factorisation has been used. For each, we show how 245 using our factorisations would affect the results in that paper. ( pure, but not complete. Using their notation, their factorisation for the CO 2 variable is (equivalent to Equation 9 in their paper):

Implications for Lunt et al
The equivalent linear-sum/shared-interaction factorisation is given by Equation 9, which in the notation of Lunt et al. (2012) 255 is: and similarly for the other three variables.
The equivalent scaled-total factorisation is given by Equation 18, which in the notation of Lunt et al. (2012) is: where dT ′ CO2 is given in Equation 19; and similarly for the other three variables. In Lunt et al. (2012), although Equation 19 (Equation 9 in their paper) was presented, the four variables were actually factorised by two N = 2 factorisations for all the analysis in that paper (Equation 13 in their paper). Because for N = 2 265 dimensions the Lunt et al. (2012), linear-sum/shared-interaction, and scaled-total factorisations are identical, the actual results related to Pliocene temperature change presented in Lunt et al. (2012) would not be affected by using our proposed  by Chandan and Peltier (2017). (b-m) The top three rows present factorisations of the total anomaly into contributions arising from changes to CO2 (upper,(b,f,j)), orography (middle, (c,g,k)) and ice sheets (lower, (d,h,l)), while the bottom row shows the residual :::::: to ::::::::::::::::::::::::::::::::: The first column (b,c,d,e) shows results using the methodology of Lunt et al. (2012) and is identical to results reported in Figure 7 of Chandan and Peltier (2018). The second column (f,g,h,i) shows results from the linear-sum/shared-interaction factorisation (Eq. :::::: Equation 6) and the third column (j,k,l,m) shows results of the scaled-total factorisation (Eq. An alternative, using the linear-sum/shared-interaction factorisation that is complete, is obtained from Equation 6, which in their notation is, for CO 2 (and analogously for the other two components): Another alternative, using the scaled-total factorisation that is complete, is obtained from Equation 16, which in their notation is, for CO 2 (and analogously for the other two components): originating from a change in CO 2 , orography and ice sheets, and is reproduced here in Figure 4(a). Figure 4(b-d) shows the results of the original factorisation and is identical to those presented in Figure 7 of Chandan and Peltier (2018). Figure 4(f-h) shows the factorisation of the same anomaly using the linear-sum/shared-interaction method (Equations 22) while Figure 4(jm) shows the results of employing the scaled-total method (Equations 23). methodology.
The bottom row in Figure 4 shows, for the case of each method, the residual difference between the sum of all the factors and the total change (i.e. the synergy ::::::::::::::: interaction/synergy ::::: terms : in the sense of Stein and Alpert (1993)). The Lunt et al. (2012) 295 method yields spatially coherent structures in the residual whose magnitude can be comparable to those of the factorized components, whereas the residuals for the other two methods are zero by definition, because they are complete  18. The linear-sum/shared-interaction factorisation is by definition bounded, because it consists of a weighted average of those very terms. In contrast, the scaled-total factorisation is not bounded, and as such it should only be used with caution. It is worth noting that if absolute weightings were used in Equation 15, such that the scaled-total factorisation became (e.g. for CO 2 ): :::::::::::::::::: , were sufficiently large.

Conclusions
In this paper, we have reviewed three previously-proposed factorisations, and extended them to produce factorisations that are unique, symmetric, :::: pure, and complete. We have presented them for 3 dimensions (i.e. 3 factors), and generalised to N dimen-   sults of these extended factorisations in the context of previous work carried out by Lunt et al. (2012), Haywood et al. (2016), and Chandan and Peltier (2018) in the context of Pliocene climate change. This reveals that the scaled-total factorisation is not bounded, and therefore can lead to anomalous results that are hard to interpret. Therefore we recommend the use of the linear-sum/shared-interaction factorisation for cases where the properties of uniqueness, symmetry, and ::::: purity, completeness, and boundedness are desirable. identical.

:
The properties of all the factorisations discussed in this paper are shown in Table 1 for 2,3, :: 4, and N dimensions. The methods that we present here will be of particular use in the analysis of systems with multiple variables, and have application beyond solely climate science.
Code and data availability. The model fields underlying Figure 4 are available from the University of Toronto Dataverse in netcdf for-335 mat: https://doi.org/10.5683/SP2/QGK5B0 . The python code used to calculate the factorisations illustrated in Figure 4 is available in the Supplement.
Author contributions. DJL led the study and wrote the first draft of the paper. GAS made Figure 3 and DC made Figure 4. DJL devised the linear-sum factorisation, GAS devised the shared-interaction factorisation, and GML devised the scaled-total factorisation. JCR provided the derivation of Equation 8. HJD, US, PJV, and AMH developed the boundary conditions and early Pliocene modelling that underpin the 340 Pliocene simulations discussed. All authors contributed to writing the final version of the paper.
Competing interests. The authors declare no competing interests.