Averaging the population projection matrices: Heuristics against uncertainty and nonexistence

Ecological Complexity 33 (2018) 66–74

Contents lists available at ScienceDirect

Ecological Complexity journal homepage: www.elsevier.com/locate/ecocom

Original Research Article

Averaging the population projection matrices: Heuristics against uncertainty and nonexistence Dmitrii O. Logofeta,b,* a Laboratory of Mathematical Ecology, A.M. Obukhov Institute of Atmospheric Physics, Russian Academy of Sciences, 3 Pyzhevskii Pereulok, Moscow, 119017, Russia b Institute of Forest Science, Russian Academy of Sciences, Sovetskaya 21, Uspenskoe, Moscow Region, 143030, Russia

A R T I C L E I N F O

Article history: Received 26 September 2017 Received in revised form 2 December 2017 Accepted 7 December 2017 Available online xxx Keywords: Matrix population models Uncertainty in data Dominant eigenvalue Nonautonomous model Matrix average Age-specific traits Absorbing Markov chains Eritrichium caucasicum Short-lived perennial

A B S T R A C T

In matrix population models, the process of population “projection” through a number of time steps is fundamentally multiplicative, hence the arithmetic mean of the consecutive matrices is of doubtful meaning, while the geometric mean quite corresponds to the multiplication principle. The geometric mean of positive numbers does not bear any problem, but that of matrices does. The “population projection matrices” (PPMs) are rather nonnegative than positive, with the allocation of non-zeros that is predetermined by the life cycle graph (LCG) reflecting the development biology of a given species, and this graph is principally incomplete. The average matrix A should logically have the same fixed pattern of zeros as those to be averaged, and this causes the averaging matrix equation to be overdetermined as a system of element-wise algebraic equations for the unknown positive elements. Therefore, the exact solution to the problem of pattern-geometric averaging does generally not exist, while the classical (leastsquares) approximation leads to a significant error. My heuristic approach to finding a better approximate A for the PPMs in the form of L = T(transition) + F(fertility) is to solve the problem in a combined way: the pattern-geometric approximation for the T part and the exact arithmetic mean for F (as population recruitment is an additive process). As a result, the approximation error decreases drastically due to matrix T being always substochastic, while the combined, TF-averaging, method turns out efficient even under ‘reproductive uncertainty’ in data, i.e., for the whole families of feasible matrices F in the sum L = T + {F}. I illustrate the method of TF-averaging with 5 matrices L(t) calibrated for each pair of consecutive years from a 6-year period of observation in a case study of Eritrichium caucasicum, a shortlived perennial herbaceous species. The approximate TF-average enables gaining the ‘age-specific traits from stage-specific models’ (Caswell, 2001, p. 116) that are characteristic of the entire period, and I discuss other motivations/advantages for/of pattern-geometric means in matrix population models. © 2017 Elsevier B.V. All rights reserved.

1. Introduction The population projection matrix (PPM) is an n
n nonnegative matrix L that governs the discrete-time dynamics in a singlespecies population structured into n stage-specific groups. The stage is understood in a generalized sense, as any discrete (or discretized) characteristic that can be used to classify the status of individuals in the population (Caswell, 2001), such as age, size, or stage of ontogeny, etc. The population structure is then represented

Abbreviations: LCG, life cycle graph; PPM, population projection matrix. * Corresponding author at: Laboratory of Mathematical Ecology, A.M. Obukhov Institute of Atmospheric Physics, Russian Academy of Sciences, 3 Pyzhevskii Pereulok, Moscow, 119017, Russia. E-mail address: [email protected] (D.O. Logofet). https://doi.org/10.1016/j.ecocom.2017.12.003 1476-945X/© 2017 Elsevier B.V. All rights reserved.

by vector x(t) 2 Rnþ whose components are the (absolute or relative) numbers of individuals in the corresponding stagespecific groups at time moment t. These moments are normally associated with consecutive (e.g., annual) censuses of the population and with the basic model equation, x(t + 1) = L(t) x(t), t = 0, 1, 2, . . . .,

(1)

where the elements of matrix L(t), or the vital rates (Caswell, 2001), depend generally on t and may also depend on the population density or the densities of some status-specific groups. However, in a growing number of practical applications (MPIDR, 2017), the PPM represents a linear transformation of the vector space. When considered as an operator in a vector space, an n
n matrix A is called a projection matrix if AA = A. The meaning of this definition is obvious: once projected to a subspace, any

D.O. Logofet / Ecological Complexity 33 (2018) 66–74

Fig. 1. An example of the LCG for a hypothetical stage-structured population with n stages, of which the (n  1)-th and the n-th stages are reproductive.

vector x 2 Rn can no longer be changed by the same operator A, but this matrix property could hardly be motivated in model (1). Therefore, the “Projection Matrices” in the title are not quite legal, though well-established among the matrix models of population dynamics (Caswell, 1989, 2001), following the idea that the matrix “projects” the current population structure for one time step further. When the knowledge of species biology is expressed as a life cycle graph (LCG, Caswell, 2001) reflecting the transitions among structural groups that may occur for one time step and the population recruitment that may appear for the same period (Fig. 1), the pattern of matrix L (i.e., the allocation of its positive and zero elements) is determined by the LCG as the directed graph associated formally to the matrix (Harary et al., 1965). For example, the LCG in Fig. 1 generates a so-called Lefkovitch (1965) matrix 2 3 0    bn1 bn 6 s1 r 2 0 07 6 7 ð2Þ L¼6 . .. 7 4 .. } } . 5 0



sn1

rn

with positive survival rates si > 0 (i = 1, 2, . . . , n  1) on the first subdiagonal, nonnegative1 stage delay rates rj  0 (j = 1, 2, . . . , n) on the principal diagonal, two positive reproduction rates bn1, bn > 0 in the first row, and zeroes elsewhere. The actual quantitative meanings of these rates follow evidently from Eq. (1) explicated for each vector component (Caswell, 2001, Ch. 4). Other LCGs may generate matrices with more complicated patterns (see, e.g., Logofet, 2008). When matrix L(t) = L remains constant over time, we obviously have x(t) = Lt x(0) 8 t = 0, 1, . . . . Perron–Frobenius theorem for nonnegative matrices, the mathematical ground of matrix population models, provides for the existence of the dominant eigenvalue coincident with the spectral radius, l1(L) = r(L) > 0, and guaranteeing, in the primitive case, the convergence x(t)/l1t ! x* as t ! 1,

Shown is the case where r1 = 0 and all the rest rj > 0.

expectancy and the mean age at first reproduction (Caswell, 2001, Ch. 5), by means of constructing some virtual absorbing Markov chains that can reach certain age-specific events in the life cycle of individuals for random times (random numbers of time steps), while the mean and variance of these random variables are amenable to calculation by known formulae (Kemeny and Snell, 1960), thus resulting in the estimates desired. However, these age-specific traits, as well as the adaptation measure l1(L), can only be relevant to quite a limited interval of real time, namely, to the period defined by the data used to calibrate the projection matrix L. This period may often reduce to just one discrete-time step of the model, e.g., the interval between two consecutive censuses at time moment t = 0 and at t = 1. The situation is typical for the type of data called ‘identified individuals’ (Caswell, 2001, p. 134), where the ‘individuals are marked and followed over time’ (Caswell, 2001, p. 134). When available is a time series of data for several time steps, it is logical to seek for the population characteristics that cover all the period of time. In terms of a matrix population model calibrated on several time steps, it means that we have a finite set of one-step matrices, L(t), each obeying Eq. (1) with the known data x(t), t = 0, 1, . . . , M  1, where M denotes the final moment in the time series of data. In general terms, we are dealing with a nonautonomous matrix model and seeking for the population characteristics averaged over the time series. The nearest way to reveal those expanded characteristics is to extract them from a matrix A that represents an average of the calibrated one-step matrices, and the question is ‘What kind of matrix averaging does correspond to the idea of time- averaged characteristics?’ Averaging issues get little attention in the literature on matrix population models, while authors use the arithmetic mean by default (Logofet, 2013a) since it poses no problems due to the linearity of matrices as operators in a vector space. In the present paper, I show, first, that the correct answer to the above question is geometric mean and, second, that it does typically not exist as an exact solution to the system of averaging equations for the PPMs of matrix population models. An approximate solution remediates the situation, and I illustrate this in Section 3, where a nonautonomous matrix model of stage-structured dynamics in a local population of a short-lived perennial plant species is reported as a case study. The field data were of the ‘identified individuals with uncertain parents’ type (Logofet, 2010, p. 33), and the uncertainty in data complicated the task to average five calibrated annual projection matrices L(t). Their own l1(L)s were localized differently with regard to 1, thus depriving the model of the ability to forecast the population dynamics in terms of asymptotic increase, steady state, or decline. Averaging has enabled the forecast to be certain as well as a certain answer to ‘How specifically short does the short-lived perennial live?’ Finally, I discuss more motivations to seek for the geometric mean of the time-specific PPMs.

(3)

where x* is a corresponding positive eigenvector whose length depends on the initial structure x(0) (Logofet, 1993; Cushing, 1998; Li and Schneider, 2002; Logofet and Belova, 2008). Thus, l1(L) represents the asymptotic growth rate of the population, hence, in applications, i.e., whenever L has been calibrated from data, it serves as a measure at which the local population is adapted to the environment, or an efficient measuring tool in comparative demography (see, e.g., Klimas et al., 2012; Logofet, 2013b, 2016). In addition, a number of ‘age-specific traits’ can be extracted ‘from stage-specific models’ (Caswell, 2001, p. 116), such as the life

1

67

2. How to average a nonautonomous matrix model Suppose we have a nonautonomous matrix model for stagestructured population dynamics in the form of Eq. (1), where vector x(t) 2 Rnþ represents the population structure at a discrete time moment t = 0, 1, . . . If the empirical data are of the ‘identified individuals with uncertain parents’ type (Logofet, 2010, p. 33), then the vectors x(t) should be known for a finite number M + 1  3 of consecutive moments t = 0, 1, . . . , M in spite of the ‘uncertain parents’ (Caswell, 2001, Ch. x6.1.1). If we have managed to calibrate the matrices L(t) at M those moments, then it follows from Eq. (1) that x(M) = L(M  1) L(M  2) . . . L(1) L(0) x(0)

(4)

68

D.O. Logofet / Ecological Complexity 33 (2018) 66–74

for the initial, x(0), and terminal, x(M), vectors served for the calibration. Note that each of the M matrices in the right-hand side of (4) are also known quantitatively from the calibration, hence so is their product. Now, the average, G, of M matrices should obviously project x(0) to the very same vector x(M), when applied M times. It should do so both for the calibrating vectors in (4) and for the model forecast over M time steps from any feasible initial structure x(0) 2 Rnþ . Therefore, we have GM – L(M  1) L(M  2) . . . L(1) L(0) = 0,

(5)

as the averaging equation, where 0 denotes the n
n matrix of zeroes. Eq. (5) represents the idea of geometric mean, but the average matrix G has also to correspond to the same LCG as each of the matrices L(t) does correspond to. Otherwise, matrix G would yield some characteristics of another population with a distinct LCG. Thus, matrix G must have the same pattern of zero/ nonzero elements, so that we are looking for a fixed-pattern mean of M given matrices. I call this kind of averaging the pattern-geometric average to contrast it with the unconstraint geometric. Remark 1. As a reflexion of species biology, the LCG is typically not complete, i.e., some transitions make no sense, hence those arcs are absent in the graph and the corresponding entries are zero in the PPMs L(t). In particular, any progressive LCG (Logofet, 2008, 2013c) implies the triangular set of zero elements between the first row and the principal diagonal (see, e.g., pattern (2)), with the number of zeros being equal to (n  1)
(n  2)/2. & Remark 2. Classical Perron–Frobenius theory for nonnegative matrices specifies an irreducible matrix A to be primitive if l1(A) = r(A) > 0 strictly dominates all the rest eigenvalues in modulus (e.g., Horn and Johnson, 1990, p. 516). Any positive element, ajj > 0, in the principal diagonal of A is a sufficient condition for A to be primitive (Horn and Johnson, 1990, Theorem 8.5.3), and this condition holds true for any stagestructured model with a PPM of the Lefkovitch pattern (2) or a more “saturated” one with positive entries (Logofet and Belova, 2008). &

So, when considered element-wise, Eq. (5) represents a system of n2 algebraic equations for a lesser number of unknown positive elements in G, i.e., an overdetermined system. Such a system is generally inconsistent (it has no solution) except for quite special cases where some equation occurs several times in the system, or if some equations are linear combinations of the others. No reason exists to assume these cases for system (5) when the calibrated matrices L(t) originate from real data. In the applications, when no exact solution exists at the background of usual data errors, we might be content with an approximate solution to the averaging Eq. (5), hence with an approximate pattern-geometric mean of given PPMs. The wellknown least-squares method gives a standard tool to search for the best approximation, i.e., a set of positive entries to matrix G that minimizes the Euclidean norm of the left-hand side in Eq. (5). Also, matrix G being a PPM imposes certain constrains on solutions of the minimization problem to be considered below. 3. How to approximate the pattern-geometric mean of population projection matrices 3.1. General consideration It is traditional to treat the projection matrix as a sum of its transition (T) and fertility (F) parts, L = T + F,

(9)

both being nonnegative nonzero n
n matrices (Cushing and Yicang, 1994; Li and Schneider, 2002). Each column of T represents transition rates for all the transitions outgoing from the corresponding node of the LCG (see Fig. 1 and matrix (2)). The rates are essentially some fractions by their meanings and, when summed up in a column j, they indicate a fraction of those individuals at stage j that survive for the next time step and can be found in any of those stage-specific groups that are prescribed by the LCG. When all the individuals at stage j do survive, the sum of column j elements equals 1. So, for any PPM (9) and its transition part T = [tij], we have tij  0 8 i; j ¼ 1; 2; . . . ; n;

n X

tij  1 8j ¼ 1; 2; . . . ; n

ð10Þ

i¼1

Remark 3. If A is a given primitive matrix, then the least k such that Ak > 0 is called the index of primitivity and is denoted by g (A) (Horn and Johnson, 1990, p. 519). Any positive element ajj > 0 in the principal diagonal of A provides for the estimation

g (A)  2n  2

(6)

(Horn and Johnson, 1990, Theorem 8.5.8). & Below, I explain why the averaging Eq. (5) does normally not have any exact solution when the product of M matrices is positive, for instance, when M  2n  2.

(7) M

Due to estimates (6) and (7), we have L(0) > 0, and, in terms of the corresponding LCG, it means that any vertex of the graph can be reached from any other one for M successive steps following the graph arrows (Harary et al., 1965). Since this observation is predetermined exclusively by the LCG, which is the same for any other L(t), t = 1, . . . , M  1, it follows that L(M  1) L(M  2) . . . L(1) L(0) > 0, 2

(8)

meaning that each of n elements in the product is positive. Meanwhile, the number of positive elements in matrix G is less than n2 by Remark 1.

Matrices of this kind are called substochastic (Katz, 1972). As regards matrix F, the second summand in sum (9), it contains only one nonzero row (the 1st one for matrix (2)) when the population recruitment occurs in a single stage only (Fig. 1). This is typical for biological populations, yet not obligatory (Logofet, 2016). If, in a nonautonomous matrix model, we have M matrices L(t) calibrated in a unique way from data, then finding the unique approximate solution to the averaging problem (5) with the pattern constraint poses no problem at all – neither theoretical due to the foundations of the well-known least squares method, nor practical one due to the corresponding routines available in the modern software systems. If however the calibrating data are of the ‘identified individuals with uncertain parents’ type (Logofet, 2010, p. 33), the multivariate calibration at each step t may result in either a unique matrix L(t) under the hypothesis of ‘maximal adaptation’(Logofet et al., 2016) or, without any additional hypotheses, in a unique transition part T(t) plus a family of fertility parts {F(t)} bearing the ‘reproductive uncertainty’ (Logofet et al., 2012, p. 99). Whether the maximaladaptation hypothesis can or can not be accepted is a matter of theoretical discourse (Metz et al., 2008a,b; Gyllenberg and Service, 2011) and that of expert judgement in each practical case. However, without any additional assumptions, seeking for the

D.O. Logofet / Ecological Complexity 33 (2018) 66–74

69

pattern-geometric mean faces a nontrivial mathematical problem presented below. The family of fertility parts {F(t)} is determined by the ‘recruitment balance equation’ (Logofet, 2008, p. 220), an obvious consequence of Eqs. (1) and (9): x(t + 1)  T(t) x(t) = F(t) x(t).

(11)

Here the left-hand side is calculable from the data and calibration outcome for T(t), while the uncertain, stage-specific reproduction rates take part in the linear expression(s) of the right-hand side with the coefficients known by the data as the ‘identified’ components of x(t). The constraint maximization problem w.r.t. the uncertain reproduction rates is generally not convex (Logofet, 2008), while avoiding this problem leads just to certain ranges for the adaptation measure l1(L) (segments [l1min, l1max]), rather than a certain number (Logofet, 2016), with the corresponding certain ranges for the uncertain rates within {F(t)}. This circumstance has prompted the following heuristic way to solve the averaging problem for a nonautonomous matrix population model under ‘reproductive uncertainty’ in data: we represent the average of L(t) as LTF = Tpg + {Fav},

(12)

where Tpg is the approximate pattern-geometric mean of M substochastic matrices T(t) (t = 0, 1, . . . , M  1) to be found by the least-squares method, and {Fav} denotes the averaged family of matrices {F(t)} to be found by interval computations (e.g., Neumaier, 1990). Since this method relies principally on the decomposition (12), it can be referred to as TF-averaging, and there are at least three arguments in favour of this approach: 1) the approximate geometric mean of M matrices L(t), each representing a family, {L(t)} = T(t) + {F(t)}, because of uncertain reproduction rates, does invoke theoretical and technical complications, while the approximation for the fixed feasible L(t)s yields a great approximation error (Logofet et al., 2017c, Appendix A); 2) the error in approximating the pattern-geometric mean of transition parts, T(t), is rather small, perhaps, due to their spectral radii being bounded, r(T(t))  1 (Horn and Johnson, 1990, Theorem 8.1.18, Corollary 8.1.30; Protasov and Logofet, 2014); 3) the technique I use further to gain ‘age-specific traits from stage-specific models’ (Caswell, 2001, p. 116) exploits only the transition parts of the PPMs. The next Subsection illustrates how the TF-averaging method applies in a case study. 3.2. Case study with reproductive uncertainty Fig. 2 presents the LCG for a local population of Eritrichium caucasicum, a short-lived perennial plant inhabiting alpine lichen heaths. Severe wintering conditions cause some of generative plants to not flower in the upcoming season, and this phenomenon causes the backward transition va g to appear in the LCG, hence a positive entry to matrix L between its first row and principal diagonal: 2 3 2 3 0 0 0 0 0 0 a b 6 c d e 07 60 0 0 07 7 6 7 L¼T þF ¼6 ð13Þ 4 0 f h 0 5 þ 4 0 0 0 0 5: 0 k l 0 0 0 0 0

Fig. 2. LCG for a local population of Eritrichium caucasicum: j denotes young plants (seedlings and juveniles); va, virignal plants and the adult non-flowering ones; g, generative plants; gt, terminally generative plants. Solid arrows indicate ontogenetic transitions occurring for one year (the lack of transition, in particular); dashed arrows correspond to annual recruitment (Logofet et al., 2016).

Correspondingly, x(t) = [j(t), va(t), g(t), gt(t)]T 2R4þ , and the data to calibrate matrix L are of the ‘identified individuals with uncertain parents’ type, mined from M + 1 = 6 annual censuses on permanent sample plots in the period of 2009–2014 (Logofet et al., 2016). After matrix (13) has been substituted for L(t) into Eq. (1), the ‘recruitment balance Eq.’ (11) takes on the following form: 3 2 2 3 2 3 32 jðtÞ jðt þ 1Þ ajðtÞ þ bgðtÞ 0 0 a b 6 7 6 7 6 7 6 7 0 0 6 7 ¼ 6 0 0 0 0 76 vaðtÞ 7 ¼ 6 7 ð14Þ 4 5 4 0 0 0 0 54 gðtÞ 5 4 5 0 0 gtðtÞ 0 0 0 0 0 0 thus restricting the uncertain reproduction rates a, b just to one linear constraint of the equality type. Therefore, the nonautonomous matrix population model represents a set of M = 5 calibrated annual matrices L(t) = T (t) + {F(t)} (Table 1; transition rates are given in their original rational forms as frequencies). We see quite a diverse allocation of 5 adaptation ranges relative to l1 = 1: two to the left of 1, two to the right, and one covering the l1 = 1. This gives no chance for any certain judgement on the future fate of the population in terms of asymptotic growth, decline, or steady state to be based upon the 6year observations. Now, the averaging Eq. (5) for the transition part Tpg, with due regard to pattern (13), takes on the following form: 2 35 0 0 0 0 6 c d e 07 7 0 ¼ T 5pg  Tð4ÞTð3ÞTð2ÞTð1ÞTð0Þ ¼ 6 40 f h 05 0 k l 0  Tð4ÞTð3ÞTð2ÞTð1ÞTð0Þ ð15Þ where matrices T(0), . . . , T(4) are numerically presented in the second column of Table 1, inside of the proper L(t)s. The product T (4)T(3)T(2)T(1)T(0) has obviously the same first row and last column of zeros as each of the multipliers does have, yet the other nine elements are all positive (Appendix A). Thus, in the element-wise form, the matrix Eq. (15) contains 9 nontrivial equations w.r.t. 7 unknowns c, d, . . . , k (omitting i, j, g), hence it is an overdetermined system, indeed. The traditional leastsquares method to find the approximate solution ought to be supplemented with certain constraints, namely, with those of matrix Tpg being substochastic and those that ensue from the sense of averaging. The latter suggests, besides the pattern conservation, that the averaging elements should occur within the same bounds as those to be averaged. In the formal, element-wise terms, we have mint{T(t)}  Tpg  maxt{T(t)}.

(16)

70

D.O. Logofet / Ecological Complexity 33 (2018) 66–74

Table 1 PPMs L(t) calibrated on the time step from t to t + 1, the bounds of reproductive uncertainty, and the ensuing ranges of l1(t). Current year

2009

2010

2011

2012

2013

Matrix L(t)

Eq. (14); the range of a, the range of b

2

0 6 68 6 149 Lð0Þ ¼ 6 4 0 0 2 0  6 17 6 31 Lð1Þ ¼ 6 4 0 0 2 0  6 76 6 150 Lð2Þ ¼ 6 4 0 0 2 0  6 137 6 211 Lð3Þ¼6 0 4 0 2 0 6 23 6 119 Lð4Þ¼6 4 0 0

0   80 0

a 10 10 1

0  106 136  9 136 2

a  6 9  1 9 1

63

5

6

3

80

10

136

0  101 129  7 129 4

9

a 4 10 2 10 3

129

0  153  181 6 181

0 0  139  296 9 296 4 296

3

10

a  6 9 1

0  9

a  4 6 2 6

0

b 07 7 07 5 0 3 b 07 7 07 5 0 3 b 07 7 7 05 0 3 b 07 7 07 5 0 3 b 07 7 07 5 0

The range of l1(t)

l1min

l1max

10a + 4b = 31; [0,3.1], [0,7.75]

0.9035

0.9949

9a + b = 150; [0,150/9], [0,150]

1.2460

1.5201

10a + 3b = 211; [0,21.1], [0,211/3]

1.2476

1.4948

9a + 7b = 119; [0,119/9], [0,17]

0.9213

1.1004

6a + b = 99; [0,99/6], [0,99]

0.7864

0.8588

adapted from Logofet et al., 2017a,b,c.

Table 2 The outcome of averaging the annual matrices in the form 2 0 0 Fav 60 0 6 40 0 0 0 a, b

of LTF = Tpg + {Fav}. 3 a b 0 07 7 0 05 0 0

a/14 2 + b/69 = 1 0 0 6 0:408294 0:708987 6 4 0:062115 0 0:015220 0

Tpg

Approx. error

l1(LTF)min xmin ; %

l1(LTF)max ;% xmax 

0 0:648420 0:044791 0:066129

3 0 07 7 05 0

1.1942
105 1.071169 [42.31 53.49 3.24 0.96]T 1.149345 [47.98 48.50 2.73 0.80]T Fig. 3. Digraph of the virtual absorbing Markov chain as a complement to the ontogenetic part of the LCG in Fig. 2.

So, to find the best approximations for matrix Tpg means to solve the constraint minimization problem for the Euclidean norm of the right-hand side of (15) (approximation error) with constraints (10) and (16). The numerical solution can be found by means of a constraint minimization routine in a modern software system, such as fmincon in MATLAB (Appendix A). The results are summarised in Table 2. The error of approximation has turned out to be surprisingly small, while the range of l1(Lav) to locate entirely to the right of l1 = 1. It means that the population increases in the long term in spite of declining for some years. The stable population structures x* for both minimal and maximal values of l1(Lav) feature strong domination by the juvenile and vegetative stages, while the share of reproductive stages being quite small. 3.3. How short does the short-lived perennial live? When asked on the average, the above question is nontrivial as far as the field data for a 6-year period of observations reveal

individuals of quite different ages (Kazantseva et al., 2016), while some of them are still alive and supposed to survive even longer. Cochran and Ellner (1992) proposed an efficient method to estimate the mean age-based life history parameters for stagestructured populations. The core idea was to treat a substochastic matrix T, the transition part of the projection matrix L = T + F for a stage-structured population, as a submatrix of the stochastic transition matrix P of a virtual absorbing Markov chain that represents the same transitions as those in the LCG plus those to an absorbing state meaning death (Fig. 3). Thereafter, known formulae from the theory of finite Markov chains (Kemeny and Snell, 1960) enable calculating the average number of time steps that the chain sojourns in nonabsorbing states before absorption, hence the parameters desired. Again, the result is predetermined by matrix T, hence it varies with t in a nonautonomous model T(t). In the E. caucasicum case, the life expectancy at the juvenile stage does vary roughly from 1 to 6 years (Logofet et al., 2017a, Table 4), thus giving another motivation to find the correct average of matrices T(t).

D.O. Logofet / Ecological Complexity 33 (2018) 66–74

With the digraph shown in Fig. 3, the stochastic transition matrix P has the following pattern:

(17)

for which the fundamental matrix is N = (I  Q)1,

(18)

where Q, the submatrix responsible for the transitions among nonabsorbing states, coincides with T. With the time-dependent matrices P(t) that are caused by T(t), the transitions for 5 time steps are described by the product P(4)P (3)P(2)P(1)P(0), which reduces to that in Eq. (15) due to the special pattern (17). This gives another motivation to average the transition parts geometrically. With the average matrix Tpg = Q shown in Table 2 and with an accuracy to 104, matrix N = [nij] (18) is presented in Table 3. Each of its elements nij represents the mean number of time steps the Markov chain sojourns in state i after starting from state j (Kemeny and Snell, 1960); element n44 is equal non-stochastically to 1 in strict accordance with the LCG (Fig. 2). In terms of demography, therefore, each column sum represents the number of time steps that individuals of the corresponding stage are expected to survive for, or the life expectancy, ES, at the stage s = j, va, g, gt. So, this short-lived perennial lives about 3 years on the average. At the first glance, it is paradoxical that Ej is markedly less than Eva, i.e., the plants of the next stage are expected to live longer than those of the current one. The reason is however in greater mortality at the juvenile stage (roughly 0.6) versus that at the vegetative stage (0.2, see Table 2). 4. Discussion Although nonlinearity is a symbol of theoretical ecology since the last century, hundreds of practical ecologists do essentially linear modelling when calibrating PPMs from data and calculating their characteristics to be used in comparative demography of the plant or animal species. The proper case studies have been contributing to the COMPADRE and COMADRE databases (MPIDR, 2017) of the corresponding publications on plant and animal species, respectively. The linear-vs.-nonlinear controversy is however not destructive as it is never destructive that we use the differential of a nonlinear function (which is a linear operation) to calculate the major part of the function increment at a finite increment of its argument. In our case, the argument increment coincides with the time step, while the PPMs L(t), when being calibrated according to Eq. (1) at each time step, give the exact fit for the nonlinear “curve” x(t) at the discrete moments. This is certainly the best fit ever reachable. Cited in the present study, the data of ‘identified individuals’ type are advantageous as enabling the calibration of the PPM to be done just for one time step, i.e., from two consecutive censuses of Table 3 Stage-specific life expectancies in the local population of E. caucasicum as a result of averaging the transition matrices T(t) over 6 years of observation. Fundamental matrix 2 1 0 6 1:6407 400184 6 4 0:1067 0:2613 0:0320 0:0784 2.7794 Ej

N = (I  Q)1 3 0 0 2:7278 0 7 7 1:2243 0 5 0:1225 1 4.3581 Eva

4.0745 Eg

1 Egt

71

Table 4 Reproduction rates a(t) and b(t) as the interval numbers with certain ranges (according to Table 1) and the outcome of averaging. t

Eq. (14)

Range of a(t)

Range of b(t)

0 (2009) 1 2 3 4 (2013) Mean ranges Approx. average

10a + 4b = 31 9a + b = 150 10a + 3b = 211 9a + 7b = 119 6a + b = 99

[0, 31/10] [0, 150/9] [0, 211/10] [0, 119/9] [0, 33/2] [0, 6353/450] [0, 14]

[0, 31/4] [0, 150] [0, 211/3] [0, 17] [0, 99] [0, 4129/60] [0, 69]

69a + 14b = 1469

the local population. Implicitly, each of those year-specific PPMs, L (t), reflects a particular environment that was effecting (or affecting) survival and reproduction during the time interval from t to t + 1. But this advantage turns into a serious problem when the environment varies in time, hence so do the population characteristics gained from the year-specific PPMs. To characterize the population over the whole period of observation, we have to average all the calibrated PPMs. Whether the geometric average exists for nonnegative matrices of a fixed pattern depends on whether the associated digraph is complete or incomplete. The latter case is typical for projection matrices of matrix population models as the biology excludes certain transitions from the LCG. The issue is therefore open in matrix theory (Logofet, 2015), where a progress has only been achieved for those positive definite matrices that commute in the product (Ando et al., 2004) and a few results found for those that do not (Bhatia and Holbrook, 2006). But the PPMs can neither be symmetric, nor commutative because of, e.g., decomposition (9). In mathematical terms, to find a solution, G, for the averaging Eq. (5) is equivalent to finding the M-th root for the product of M given matrices, each being a primitive nonnegative matrix; the product, B, is therefore a primitive nonnegative (or even a positive) matrix, too. The problem has been permanently attracting certain attention in the literature on matrix theory (Jain et al., 1979; Denman, 1981; Reams, 1997; Higham and Lin, 2009; McDonald et al., 2014; Politi and Popolizio, 2015; Tam and Huang, 2016), treated mostly as a particular case of finding a matrix function of a given matrix (Horn and Johnson, 1994, Ch. 9; Higham, 2008, Ch. 7). In particular, it was noted that many questions about nonnegative matrices with a positive eigenvector could be reduced to the questions about stochastic matrices (Horn and Johnson, 1990, x8.7. P1), linking to applications with discrete-time Markov chain models. Therefore, the published results cover a variety of problem formulations and solution methods, but none of them imposes the constraint of a fixed associated digraph for the M-th root of a given matrix B, which means, in terms of matrix population models, the same LCG for the geometric average as that for the individual PPMs to be averaged. Moreover, numerical examples demonstrate typically the same (or even greater) number of nonzero elements in the matrix roots as (than) that in matrix B (Appendix C; Denman, 1981; Higham and Lin, 2009; McDonald et al., 2014; Politi and Popolizio, 2015), whereas these numbers ought to differ essentially in the present case study (Table 2). Quantitative constraints such as bounds (16), which are sound from the practical viewpoint, seem to have never been imposed in the mathematical formulations. This justifies my introducing the pattern-geometric averaging as a new mode of averaging motivated in matrix population models for the case of several time-specific PPMs. The approximate solution I suggest as TF-averaging meets all the constraints mentioned above. Besides providing for a certain answer to the question in the title of Section 3.3, the patterngeometric mean of substochastic transition matrices T(t) is also

72

D.O. Logofet / Ecological Complexity 33 (2018) 66–74

crucial for other ‘age-specific traits from stage-specific models’ (Caswell, 2001, p. 116), such as the mean age at first reproduction, in nonautonomous cases. The technique has revealed it to be 12 years in the E. caucasicum case study (Logofet et al., 2017b, Table 2) in spite of 6 years only available in the observation data. In the theory of population growth in stochastic environments (Caswell, 2001, Ch. 14), the time series of population data is supposed to implicitly represent the range of temporal random variations in the environmental conditions. Once calibrated on a time step, each time-specific projection matrix L(t) is considered as an explicit expression of the environment that has realised during the time step, while predicting the population fate under random fluctuations of the environment, a concept of stochastic growth rate, lS (Tuljapurkar, 1986, 1990; Caswell, 2001), is applied elsewhere. The direct calculation of lS suggests a long series of Monte Carlo simulations, while theoretical approximations of lS (Tuljapurkar, 1986, 1990; Caswell, 2001) just assume the mean matrix to exist rather than propose any practical recipes to average the timespecific matrices in hand. The concept of pattern-geometric mean might fill this gap in theory and practice. When the population experiences stochastic environments, the order at which time-specific PPMs occur in their product is no longer fixed (as it was in the definition of the pattern-geometric mean), but obeys a random choice. However, the number of different products (of a finite number of PPMs) is still finite, hence there exist a product yielding the least approximation error in the proper pattern-geometric mean or a mean matrix that is the closest to the centre of the ensemble in a metric of the matrix space. To study the ensuing issues of existence and uniqueness is a challenge to the theory of matrix population models motivated in the modelling practice. Meanwhile, the accuracy of the best approximation in the problem of pattern-geometric average has turned out surprisingly high in the E. caucasicum case (Table 2). Whether it is a consequence of matrices T(t) being substochastic, with their r(T) < 1, is a particular question to be answered. No question is there about whether the other heuristic modes of averaging, such as arithmetic or element-wise geometric, might achieve any better approximation accuracy when they substitute for Tpg in the overdetermined system (15) (Appendix B). The arithmetic mean has evidently turned out the worst both in the approximation error and in the overestimation of the range of l1(LTF), with LTF = Tpg + {Fav}. This standard decomposition prompts a heuristic explanation for the success of Tpg: while the transition probabilities for several steps are essentially multiplicative, the population recruitments are additive in essence. Further studies with the direct and approximate calculations of stochastic growth rate lS by means of the approximate geometric mean might reveal further advantages of this practical approach. It was a school exercise to prove that geometric mean of positive numbers could never be greater than their arithmetic mean. Table 5 illustrates this rule with the elements of two

alternative average matrices, while the ensuing l1s illustrate the monotonicity of l1(A), the dominant eigenvalue of a nonnegative matrix, as a function of matrix elements (Horn and Johnson, 1990, Ch. 8). To prove that l1 of the real geometric mean (rather than the element-wise) of nonnegative matrices can never be greater than l1 of their arithmetic mean is therefore another challenge to matrix theory. 5. Conclusion A simple task to average several (M) given nonnegative matrices turns out complicated in matrix population models, both in theoretical and practical aspects. While the correct mode of averaging is theoretically shown to be pattern-geometric, its mathematical expression, evolving from the M-th root of a given nonnegative matrix, encounters certain obstacles in the practical application. In general, the exact mathematical solution to the M-th root problem cannot guarantee the same pattern of the average matrix as the pattern of those to be averaged, which is predetermined by a given LCG, the (incomplete) directed graph associated with each of the year-specific PPMs. Adopting this constraint leads to the concept of pattern-geometric averaging with the nonexistence of exact solution to the averaging problem. The approximate solution ensues from a constraint minimization problem whose formulation adopts also the observable bounds for the matrix elements. The TF-averaging, a heuristic method proposed here to solve the approximation problem, makes use of the standard PPM decomposition into the transition (T) and fertility (F) parts and combines the pattern-geometric mean of the transition parts with the arithmetic mean of the fertility ones. The resulting error of approximation has turned out to be of the same order of magnitude as the direct approximation error in the Eritrichium caucasicum case study. The approximate patterngeometric mean has enabled gaining the ‘age-specific traits from stage-specific models’ that are characteristic of the entire observation period, and there are other motivations/ advantages for/of pattern-geometric means in matrix population models. Acknowledgments The study is supported by the Russian Foundation for Basic Research, project No. 16-04-00832. Author is grateful to three anonymous reviewers, whose critical remarks have markedly improved the manuscript. Appendix A. Approximate solution to the averaging Eq. (15) Let a vector of unknown parameters be denoted as p = [c, d, . . . , k] (omitting i, j, g). Then p 2 B R7þ where B is defined by the substochasticity constraints: 0  c  1, 0  d + f + k  1, 0  e + h + l  1,

(A1)

Table 5 Alternative modes of averaging (aa) and the corresponding ranges of l1(Laa). Mode of averaging

Approx. error [l1min, l1max]

Arithmetic 2 0 0 6 0:4708 0:7330 6 4 0:0518 0 0:0118 0 0.0648 [1.1038, 1.1860]

0 0:5800 0:1889 0:1244

3 0 07 7 05 0

Element-wise geometric 2 0 0 0 6 0:4369 0:7180 0:5683 6 4 0:0486 0:1246 0 0:0162 0:0871 0 0.0100 [1.0523, 1.1812]

3 0 07 7 05 0

D.O. Logofet / Ecological Complexity 33 (2018) 66–74

73

and by the lower and upper bounds (16), which take on the following numeric form (see Table 1): 3 2 3 2 0 0 0 0 0 0 0 0   6 23 119 139 296 2 5 0 7 6 c d e 0 7 7 6 76  9 4 0 0 05 40 f h 05 296 0 0 0 0 20 k l 0 3 0 0 0 0    153 2 6 137 07 6 7 211  181 1 3 6 ðA2Þ 3 0 07 4 5 40 3   3 4 0 0 129 10

the errors of approximations when those averages substitute for Tpg in the overdetermined system (15). Because of occasional zeros in the entries (4, 2), (4, 3), and (3, 3) to matrices T(t), their actual Hadamard product does have zeros at these places, too, which contradict the idea of averaging. Therefore, the minimal positive entry 4/296 (Table 1) is substituted for those zeros in the multipliers T(0), T(3), and T(4).

Now, the problem is to find 2 35 0 0 0 0 6 c d e 07 7 minB k 6 4 0 f h 0 5  Tð4ÞTð3ÞTð2ÞTð1ÞTð0Þ k 0 k l 0

Extreme values of l1s, such as l1max(0), l1min(1), etc. (Table 1), require reproduction rate b = 0, hence the 4-th column of zeros and the corresponding matrix L(t) being singular. To avoid manipulation with singular matrices, we select those L(t)s of the feasible ones that produce a non-singular product, e.g., those with a = b.

ðA3Þ

Appendix C. Calculating the 5th root of the product L(4)L(3)L(2) L(1)L(0)

>> Pr aa ¼ ½42874249=243530070;

759338657=1270185600;

17695857883=7779886800;

1256827=2288202

½ 606733743973=1606463501760;

589066799 . . . =586520902 . . . ;

1391122005 . . . =513205789 . . . ;

329677277=580549536

½3081824857=133871958480;

330710802887=4887674188800;

860405230819=4276714915200;

3510007=96758256

½176799191=20080793772;

90184185389=3665755641600;

216247822177=3207536186400;

831181=72568692;

under constraints (A1) and (A2), where || . . . || denotes the Euclidean norm, while the product of 5 matrices takes on the following numeric form: 2

0 6 2003449367  6 15446764440 6 6 6 57631247  6 5148921480 4  2739937 772338222

137029965913

0 

3940949473



37482403



5141191577

0 

147748897



563962406400

23498433600

187987468800 7032559

5639624064



7832811200 1174921680

0

3

7 07 7 7; 7 07 5

>> evalðPr aa^ =5Þ ans ¼ 0:0352 0:0755 0:0046 0:0018

0:1196 0:2009 0:0135 0:0049

0:4549 0:5421 0:0402 0:0135

0:1099 0:1136 0:0073 0:0023;

0

ðA4Þ

with the 9 > 7 positives entries, indeed. To solve the problem numerically, MATLAB offers a routine called fmincon (MathWorks, 2016) that minimizes a given nonlinear function under given constraints. The function is written down explicily in (A3) and (A4), while the constraints in (A1) and (A2). The function in (A3) is readily implemented with standard MATLAB tools, say, as a user-defined function normT5_ProT. However, the syntax of fmincon supposes the function argument to be a vector, rather than a matrix, and its upper and lower bounds to be a special kind of constraint. This required a pair of mutually reciprocal functions that convert a matrix T of the fixed pattern (13) into the corresponding column vector p 2 R7þ and vice versa. Eventually, the proper call looks like

[p, FVAL] = fmincon(@normT5_ProT, p0, A, b, [], [], minTvec, maxTvec),

Calculated in symbolic MATLAB, the product then takes on the following rational form: where from the 5-th root is calculated as

(A5)

where the algorithm starts at p0, inequality Ap  b implements the constraints (A1), two empty input arguments [] designate the absence of equality-type constraints, minTvec and maxTvec realise the bounds (A2). The output arguments p and FVAL represent the minimizing solution and the minimal value (approximation error) respectively, see Table 2. The outcome of interval computations for matrices {F(t)}is shown in Table 4 Appendix B. Arithmetic and element-wise geometric averaging Table 5 contains the arithmetic and element-wise (Hadamard) geometric averages of matrices T(t) presented in Table 1, as well as

and it is positive in contrast to the nonnegative pattern-geometric mean to be found. The symbolic rational form of matrix elements is crucial for the accuracy of calculations as far as the condition number of the product,

eval(cond(Pr_aa)) ans = 3.9165e + 06  3.7709e  13i, is fairly great, hence the product is an ill-conditioned matrix. References Ando, T., Li, C.-K., Mathias, R., 2004. Geometric means. Linear Algebra Appl. 385, 305–334. Bhatia, R., Holbrook, J., 2006. Noncommutative geometric means. Math. Intell. 28, 32–39. Caswell, H., 1989. Matrix Population Models: Construction, Analysis, and Interpretation. Sinauer, Sunderland, MA. Caswell, H., 2001. Matrix Population Models: Construction, Analysis and Interpretation, 2nd ed. Sinauer, Sunderland, MA. Cochran, M.E., Ellner, S., 1992. Simple methods for calculating age-based life history parameters for stage-structured populations. Ecol. Monogr. 62 (3), 3455–3464. Cushing, J.M., 1998. An Introduction to Structured Population Dynamics (Conference Series in Applied Mathematics, vol. 71. SIAM, Philadelphia, PA. Cushing, J.M., Yicang, Z., 1994. The net reproductive value and stability in matrix population models. Nat. Resour. Model. 8, 297–333. Denman, E.D., 1981. Roots of real matriuces. Linear Algebra Appl. 36, 133–139. Gyllenberg, M., Service, R., 2011. Necessary and sufficient conditions for the existence of an optimisation principle in evolution. J. Math. Biol. 62, 359–369. Harary, F., Norman, R.Z., Cartwright, D., 1965. Structural Models: An Introduction to the Theory of Directed Graphs. John Wiley, New York. Higham, N.J., 2008. Functions of Matrices: Theory and Computation. SIAM, Philadelphia, PA. Higham, N.J., Lin, L., 2009. On pth Roots of Stochastic Matrices. MIMS EPrint, pp. 21. http://eprints.maths.manchester.ac.uk/1241/.

74

D.O. Logofet / Ecological Complexity 33 (2018) 66–74

Horn, R.A., Johnson, C.R., 1990. Matrix Analysis. Cambridge University Press, Cambridge, UK. Horn, R.A., Johnson, C.R., 1994. Topics in Matrix Theory, 2nd ed. Cambridge University Press, Cambridge, UK. Jain, S.K., Goel, V.K., Kwak, E.K., 1979. Nonnegative mth roots of nonnegative 0symmetric idempotent matrices. Linear Algebra Appl. 37–51. Katz, M., 1972. On the extreme points of the set of substochastic and symmetric matrices. J. Math. Anal. Appl. 37, 576–579. Kazantseva, E.S., Onipchenko, V.G., Kipkeev, A.M., 2016. Age at the first flowering in alpine herbaceous short-lived perennials of Northwest Caucasus. Bull. MOIP Seriya Biologiya 121 (2), 73–80 (in Russian). Kemeny, J.G., Snell, J.L., 1960. Finite Markov Chains. Van Nostrand, Princeton, NJ. Klimas, C.A., Cropper Jr., W.P., Kainer, K.A., Wadt de Oliveira, L.H., 2012. Viability of combined timber and non-timber harvests for one species: a Carapa guianensis case study. Ecol. Modell. 246, 147–156. Lefkovitch, L.P., 1965. The study of population growth in organisms grouped by stages. Biometrics 21, 1–18. Li, C.-K., Schneider, H., 2002. Application of Perron–Frobenius theory to population dynamics. J. Math. Biol. 44, 450–462. Logofet, D.O., 1993. Matrices and Graphs: Stability Problems in Mathematical Ecology. CRC Press, Boca Raton, FL p. 308. Logofet, D.O., 2008. Convexity in projection matrices: projection to a calibration problem. Ecol. Modell. 216, 217–228. Logofet, D.O., 2010. Svirezhev’s substitution principle and matrix models for dynamics of populations with complex structures. Zhurnal Obschei Biologii (J. Gen. Biol.) 71 (1), 30–40 (in Russian, with English summary). Logofet, D.O., 2013a. Projection matrices in variable environments: l1 in theory and practice. Ecol. Modell. 251, 307–311. Logofet, D.O., 2013b. Complexity in matrix population models: polyvariant ontogeny and reproductive uncertainty. Ecol. Complex. 15, 43–51. Logofet, D.O., 2013c. Calamagrostis model revisited: matrix calibration as a constraint maximization problem. Ecol. Modell. 254, 71–79. Logofet, D.O., 2015. Current problems in matrix population models: expanding classics and new discoveries. 4th International Conference on Matrix Methods in Mathematics and Applications (MMMA-2015). Program and Abstracts 39– 40. http://matrix.inm.ras.ru/program_and_abstracts.pdf. Logofet, D.O., 2016. Estimating the fitness of a local discrete-structured population: from uncertainty to an exact number. Ecol. Modell. 329 (10), 112–120. Logofet, D.O., Belova, I.N., 2008. Nonnegative matrices as a tool to model population dynamics: classical models and contemporary expansions. J. Math. Sci. 155 (6), 894–907.

Logofet, D.O., Ulanova, N.G., Belova, I.N., 2012. Two paradigms in mathematical population biology: an attempt at synthesis. Biol. Bull. Rev. 2 (1), 89–104. Logofet, D.O., Ulanova, N.G., Belova, I.N., 2016. Polyvariant ontogeny in woodreeds: novel models and new discoveries. Biol. Bull. Rev. 6 (5), 365–385. Logofet, D.O., Belova, I.N., Kazantseva, E.S., Onipchenko, V.G., 2017a. Local population of Eritrichium caucasicum as an object of mathematical modelling. I. Life cycle graph and a nonautonomous matrix model. Biol. Bull. Rev. 77 (2), 106– 121 (in press). Logofet, D.O., Kazantseva, E.S., Belova, I.N., Onipchenko, V.G., 2017b. Local population of Eritrichium caucasicum as an object of mathematical modelling. II. How short does the short-lived perennial live? Zhurnal Obschei Biologii 78 (1), 56–66 (in Russian, with English summary). Logofet, D.O., Kazantseva, E.S., Belova, I.N., Onipchenko, V.G., 2017c. How long does a short-lived perennial live? A modeling approach. Zhurnal Obschei Biologii 78 (in press). MathWorks, 2016. http://www.mathworks.com/help/optim/ug/fmincon.html. McDonald, J.J., Paparella, P., Michael, J., Tsatsomeros, M.J., 2014. Matrix roots of eventually positive matrices. Linear Algebra Appl. 456, 122–137. Metz, J.A.J., Mylius, S.D., Diekmann, O., 2008a. When does evolution optimize? Evol. Ecol. Res. 10, 629–654. Metz, J.A.J., Mylius, S.D., Diekmann, O., 2008b. Even in the odd cases when evolution optimizes, unrelated population dynamical details may shine through in the ESS. Evol. Ecol. Res. 10, 655–666. MPIDR, 2017. http://www.demogr.mpg.de/en/media/1867_media_2293.htm. Neumaier, A., 1990. Interval Methods for Systems of Equations (Encyclopedia of Mathematics and its Applications), vol. 37. Cambridge Univ. Press, Cambridge, UK. Politi, T., Popolizio, M., 2015. On stochasticity preserving methods for the computation of the matrix pth root. Math. Comp. Simul. 110, 53–68. Protasov, V.Yu., Logofet, D.O., 2014. Rank-one corrections of nonnegative matrices, with an application to matrix population models. SIAM J. Matrix Anal. Appl. 35 (2), 749–764. Reams, R., 1997. A Galois approach to mth roots of matrices with rational entries. Linear Algebra Appl. 258, 1287–1294. Tam, B.-S., Huang, P.-R., 2016. Nonnegative square roots of matrices. Linear Algebra Appl. 498, 404–440. Tuljapurkar, S.D., 1986. Demography in stochastic environments. II. Growth and convergence rates. J. Math. Biol. 24, 569–581. Tuljapurkar, S.D., 1990. Population Dynamics in Variable Environments. Springer, New York.