Pattern-Multiplicative Average of Nonnegative Matrices Revisited: Eigenvalue Approximation Is the Best of Versatile Optimization Tools

1. Introduction

Pattern-multiplicative average (PMA) of nonnegative matrices is a concept motivated in the field of Matrix Population Models by the need to assess the viability of a local population from the outcome of its long-term monitoring. Matrix Population Models represent a wide and rapidly growing area of applications [CA23, CP23], where the theory of nonnegative matrices applies to population dynamics [BP94]. When biologists study a population of a plant or animal species and observe the population at discrete time moments, they classify individuals according to an observable trait such as age, body size, ontogenetic stage, etc., and divide the entire range of trait values into certain discrete classes, thus dealing with a discrete population structure [Cas01]. It is expressed mathematically as a vector x(t) ϵ ${\mathbb{R}}_{+}^n$ (n ≥ 2) representing the (absolute or relative) abundances of n class-specific groups of individuals, whereafter the Matrix Population Model is a system of first-order difference equations,

x(t + 1) = L(t) x(t), t = 0, 1, 2, …,

(1)

linking the population vector at the next moment, t + 1, to that at the current one, t.

Matrix L(t) is therefore called the population projection matrix (PPM) [Cas01], despite the projection term has quite another meaning in matrix algebra [HoJ90]. All the matrix entries are nonnegative and called vital rates as they have certain demographic sense [Cas01]. The pattern of how zero/nonzero elements are arranged within the matrix is associated, in a unique way, to what is called the life cycle graph (ibidem) as it bears the biological knowledge of how the individuals grow/develop for one time step and provide for the population recruitment.

Being nonnegative, the PPM becomes a legal target for the Perron–Frobenius Theorem [Gan59, HoJ90], whereby a rich repertoire of model properties and population characteristics can be obtained after the PPM having been calibrated from data [Cas01, LSG21]. In particular, the dominant eigenvalue, λ₁ > 0, of matrix L(t) serves as a “measure of how the local population is adapted to the environment where and when the data have been mined to calibrate the PPM” [L19] (p.1). When the data are of the ‘identified individuals’ type [Cas01], a great advantage of MPMs is that the calibration can be done just from two successive observations, thus meeting Eq.(1) and yielding a time-dependent L(t). Correspondingly, the adaptation measure, λ₁, inherits the time dependence, and the great advantage turns dialectically into a great trouble or challenge when the task is to assess the adaptation from a time series of observation data, hence from a finite set of one-step PPMs.

The historically first response to the challenge was grounded on the theory of random sequences of vectors L(t)x(t), t → ∞ [Co79, TuO80, Tu86, Tu90], and it has led to what is now called the stochastic growth rate of a population in a randomly changing environment (see [Cas01] and refs therein). Another, more recent approach reverts the task just to averaging a given number, say M (M ≥ 2), of one-step PPMs and calculating the dominant eigenvalue, λ₁(G), of the average matrix G. Letter “G” here bears a hint to the geometric (or multiplicative) nature of averaging, and the hint will be revealed further in Section 2.2. This approach has led to a concept that was proposed 5 years ago with regard to matrix population models [L18], and the task of averaging was reduced to a constrained nonlinear minimization problem for the matrix norm [L19, L18] (see Section 2.3). The norm has however not appeared quite fit for global optimization, but required versatile optimization tools [PZL22] to solve the problem.

This paper aims at expanding those “versatile optimization tools” with another target of optimization, namely, the deviation of λ₁(G_opt) from the ideal calculable λ₁(G) (Section 2.4). The approach originates from the fact the ideal λ₁(G) is known:

λ₁(G) = (Prod)^1/M,

(2)

where Prod denotes the product of M one-step PPMs. This leads to a problem that appears to be formalizable as a standard problem of linear programming (LP), hence solvable by standard routines. The solution is exemplified by a case study of monitoring a local population of an alpine short-lived perennial for 13 years [Let23] (Section 2.1). Sections 2.3 and 2.4 contain respectively a short presentation of the published “versatile optimization tools” and an in-depth one of the original LP approach. The outcomes of the both are presented in Section 3, enabling the comparison of results, which reveals a surprising contrast and induces some comments (Section 4), in particular, on how general the proposed LP technique might be. $\sout{\mathrm{Section}} \sout{5}$ contains some ensuing conclusions.

2. Materials and Methods

2.1. Case Study

A local population of Androsace albana, an alpine short-lived perennial plant species, was monitored on permanent sample plots once a year [Let23] during 14 years (2009–2022). Although A. albana reproduces actually by seeds, the stage of dormant seeds was deliberately excluded from the life cycle graph by ontogenetic stages as the seed-related vital rates cannot be reliably estimated in the field [LKBO18]. Later on, we proved such an exclusion to be correct within the integer-valued formalism [LKO20], which can always be applied for the ‘identified individuals’ data. The life cycle graph looks correspondingly as shown in Figure 1, with the “virtual” rather than real recruitment arcs. These three arcs ingo respectively to three early stages since the newly recruited plants can be observed in any of them at the moment of observation.

Now, we have x(t) ϵ ${\mathbb{R}}_{+}^5$, and the natural order of components in x(t) leads to the following PPM pattern:

$$\mathit{L}=\left[\begin{array}{ccc}\begin{array}{cc}0& 0\\ d& 0\end{array}& \begin{array}{cc}0& 0\\ 0& 0\end{array}& \begin{array}{c}a\\ b\end{array}\\ \begin{array}{cc}e& f\\ 0& 0\end{array}& \begin{array}{cc}h& 0\\ k& l\end{array}& \begin{array}{c}c\\ 0\end{array}\\ \begin{array}{cc}0& 0\end{array}& \begin{array}{cc}0& m\end{array}& 0\end{array}\right]; a, b, \dots , l, m\ge 0,$$

(3)

with 10 vital rates to be calibrated. The calibration for 13 successive pairs of observation years results in 13 annual PPMs, L(t) = [l_ij(t)], presented in Table 1. Each of them obeys Eq. (1). There are 6 annual PPMs with λ₁(L(t)) > 1 and 7 of those with λ₁(L(t)) < 1, so that even an “educated guess” of what should be λ₁(G), the dominant eigenvalue of the average matrix, remains uncertain.

2.2. Pattern-Multiplicative Average of Annual PPMs

The logic of PMA ensues from the following apparent observation: given the 13 annual PPMs calibrated according to Eq. (1) (Table 1), we have the initial population structure projected to the terminal one,

x(2022) = L(2021) L(2020) _{. . .}L(2009) x(2009),

by

Prod = L(2021) L(2020) _{. . .} L(2009),

(4)

the product of 13 annual PPMS in the chronological order.

x(2022) = Prod x(2009).

(5)

Since the average matrix ought to do absolutely the same when raised to the 13th power, we have the following basic equation of averaging:

G¹³ = Prod,

(6)

which justifies Eq. (2) of the Introduction. As is noted ibidem, a routine extraction of the 13th root from a given, even positive, Prod can hardly return a nonnegative matrix, nor guarantee the fixed pattern (3) for the root.

A workaround leads to solving matrix equation (6) as a system of scalar polynomial equations for unknown matrix elements. In fact, there are 5×5 – 10 = 15 elements that are prescribed to be zero by pattern (3), while the corresponding 15 scalar equations bind the unknown elements too, in addition to the 10 equations for the 10 unknowns. It means that Eq. (6) is overdetermined as a system of scalar equations, hence it may have an exact solution only under special relations between the equations depriving those extra 15 ones of their independence. No reason however to seek this kind of relation among the PPMs, so that Eq. (6), the basic equation of averaging has no solution over the set of pattern (3) PPMs.

The simple and clear logic of PPM averaging thus faces a principal obstacle in matrix algebra, and a workaround comes to an approximate solution.

2.3. Approximate PMA as a Nonlinear Constraned Minimization Problem

Any approximation invokes the question of how close is the approximate solution to the exact one. Although the latter is unknown, the equivalent form of Eq, (6):

G¹³– Prod = 0,

(7)

just illustrates the ideal, yet unreachable, closeness as a 5×5 matrix of zeros. Retaining the same notation G = [g_ij] for the approximate average matrix, we can see the elements of G¹³ as certain polynomials of the 5×13 = 65th degree. If we consider the quality of approximation to be a scalar value, then it should be the error of approximation. The error can be measured as the (squared) matrix norm, whereby we come to the following nonlinear minimization problem:

$$\mathsf{\Phi }\left(\mathit{G}\right)\to \underset{\mathbb{B}}{\min }{\mathsf{\Phi }\left(\mathit{G}\right)=\left|\right|{\mathit{G}}^{13}–\mathit{P}\mathit{r}\mathit{o}\mathit{d}\left|\right|}^2.$$

(8)

Here ${\mathbb{B}⊏\mathbb{R}}_{+}^{10}$ defines the problem to be a constrained one as it represents the constraints that ensue from the logics of averaging: the average vital rates should not be less, nor greater than those fixed in observations. In matrix terms, it means that each entry of the annual PPM has to lay within the bounds defined by the minimal value and the maximal one among the 13 annual values of that entry. In formal terms,

$$\underset{t}{\min }{l}_{ij}\left(t\right)\le g_{ij}\le \underset{t}{\max }{l}_{ij}\left(t\right), i, j=1,\dots , 5$$

(9)

(see Table 2).

So, we come to the constrained nonlinear minimization problem (8–9), where the bounds

of B are defined in Table 2. Both standard routines for global optimization [MW23G] and more sophisticated ones [PZL22] can be applied for solving problem (8–9).

2.4. Approximate PMA as an Eigenvalue Constrained Optimization Problem

Noted in the Introduction, expression (2) is actually a simple algebraic consequence of (6), the basic equation of averaging, that would be true if the exact solution existed. Since it does not exist, another measure of approximation quality can be found by the following speculations.

Suppose v to be a dominant eigenvector of the average matrix G, i.e., we have

G v = λ₁(G) v.

(10)

Acting by operator G successively on the both sides of Eq. (10), we get

(G)² v = λ₁(G)² v,

(G)³ v = λ₁(G)³ v,

…,

(G)¹³ v = λ₁(G)¹³ v.

(11)

In combination with the averaging Eq.(6), Eq.(11) signifies that v is a dominant eigenvector of the matrix product, Prod, too. Hereby, the other measure of approximation error is

|λ₁(G) – ρ₀|, where ρ₀ = ρ(Prod)

(12)

can be calculated from the data presented in Table 1. The optimal solution should be searched for among those matrices G that has the same eigenvector v. To exemplify how matrices of a single pattern with different eigenvalues may have the same dominant eigenvector is a student level exercise (A hint to solution in Appendix A).

To be unique, this eigenvector, v*, can, for instance, be expressed in percent and calculated like the vectors x* presented in Table 1. It can thereafter serve as a constraint in the search of approximation among various λ₁ values associated with a single v*. (To give an example of matrix eigenvalue being nonlinear as a function of matrix elements is another exercise.) Hence, we have a nonlinear minimization problem again. However, the change of variables, a beloved mathematical trick,

y = |λ₁(G) – ρ₀|,

(13)

turns the problem into a linear one. More specifically, it turns the task into the following two linear problems: to maximize the λ₁(G) when ρ₀ ≥ λ₁(G) and to minimize it when λ₁(G) ≥ ρ₀.

Moreover, vector v* = [v₁, v₂, …, v₅], the known positive eigenvector of G, and its dominant eigenvalue λ₁(G) are bound by definition to the matrix equation

Gv* = λ₁(G)v*.

With due regard to the fixed pattern (3) and with omitting subscript 1 at λ₁, this matrix equation reduces to the five scalar linear equations for a, b, …, l, m, λ, the ten unknown entries of G and the 11th formal variable λ:

av₅ – λv₁ = 0,
bv₅ + dv₁ – λv₂ = 0,
cv₅ + ev₁+ fv₂+ hv₃ – λv₃ = 0,
kv₃ + lv₄ – λv₄ = 0,
mv₄ – λv₅ = 0.

(14)

The 11th variable has its own bounds: as the dominant eigenvalue of the average matrix, the unknown λ has to locate within the following bounds:

$$\lambda \in \left[\underset{t}{\mathit{min}}{\lambda }_1\left(\mathit{L}\left(t\right)\right), \underset{t}{\mathit{max}}{\lambda }_1\left(\mathit{L}\left(t\right)\right)\right], t=2009, \dots , 2021,$$

to be determined from Table 1. The two specific problems specify these bounds further:

$$\lambda \in \left\{\begin{array}{c}\left[\underset{t}{\mathit{min}}{\lambda }_1\left(\mathit{L}\left(t\right)\right), {\rho }_0\right] if \lambda \le {\rho }_0;\\ \left[\underset{t}{{\rho }_0,\mathit{min}}{\lambda }_1\left(\mathit{L}\left(t\right)\right)\right] if \lambda \ge {\rho }_0.\end{array}\right.$$

(15)

Finally, we have the following problem formulation:

$$y\to \underset{p}{\min }{y}_p \left(p=1, 2\right);y_p\to \underset{\mathit{x}\in {\mathbb{B}}_p}{\min }{y}_p\left(\mathit{x}\right),$$

(16)

where x = [a, b, …, l, m, λ]^T, while the polygons ${\mathbb{B}}_p\subset {\mathbb{R}}_{+}^{11}$(p = 1, 2) are defined in (15) and the equality-type constraints (14). The two inner problems (16) represent the standard LP ones to be solved by standard software routines, such as function linprog in MATLAB^® [MW23L] (technical details in Appendix B).

3. Results

3.1. Case Study Outcome

Calibrated with the absolute accuracy, 13 annual PPMs for A. albana are presented in Table 1 with their dominant eigenvalues and corresponding eigenvectors.

3.2. Pattern-Multiplicative Averaging as an Approxiation Problem

The logics of approximate PMA have been implemented as problem settings in Subsections 2.3 and 2.4. Matrix Prod (4), the product of 13 annual PPMs, is too cumbersome to be presented here in its original rational form. Its numerical form up to the 15th significant digit is shown in Table 3 together with the dominant eigenvalue and the corresponding eigenvector.

3.3. Minimizing the Approximation Error as a Matrix Norm.

Solution to the constraint minimization problem (8) has been obtained in [PZL22] by the “versatile optimization tools” (Table 4) for the case of 12 annual PPMs. The quality of approximation ranges from 0.002374 to 0.002379 as a function of the optimization method (ibidem). In the format of Table 1, the obtained average matrices are presented in Table 4.

3.4. Eigenvalue Constrained Optimization Problem

Solutions to the problem (14–16) are shown in Table 5.

The solutions to the two optimization problems look indistinguishable within the publication format allowed, and they are really different at the 16th decimal digit alone, which may well be attributed to the computer round-off errors. It means that both problems have the same unique optimal solution, i.e., the corresponding vertex is in common of the polygons B₁ and B ₂.

4. Discussion and Conclusions

The case study practice has not motivated any logics of PMA beyond the population structure projection from the initial vector to the terminal one, i.e., Eqs. (4-5). Therefore, we deal here with single-objective problems alone (in the terminology of [DRH21]). Once having been formulated in mathematical terms [L18], this objective comes immediately to the difference norm as the distance between two matrices following the fundamental of real analysis [K22]. However, even the squared matrix norm has required “versatile optimization tools” [PZL22] to achieve even moderate approximation quality in the corresponding nonlinear constrained minimization problem (Table 4).

Beside of the moderate approximation quality, the basin-hopping global search, as an essential part of the “versatile optimization tools”, leaves no chance for the results being quantitatively reproducible as it is essentially stochastic [RC04]. In contrast, the LP approach is essentially deterministic [Lu73], hence the results of Section 3.4 can be reproduced any number of times.

What does not contrast the LP approach to the PMA task is the eigenvalue optimization problem being nonlinear, too. Since the matrix eigenvalue in general and the dominant eigenvalue in particular is not linear as a function of matrix elements (see Appendix A), the use of linear programming seems paradoxical at the first glance. However, beside of the trick with change of variables (11), the very nature of eigenvalue–eigenvector relations is linear, resulting in the linear relationships among matrix elements and the linear equality-type constraints in the LP setting. The trick has nevertheless failed to cheat a software routine formality, namely, to express the loss function as a linear combination of formal variables at fixed coefficients (see Appendix B). The reason is quite evident: function (11) is not a linear function, albeit consisting of two linear parts.

These two parts have logically resulted in the two adjacent LP problems differing in the range (15) of the lambda alone. The both yield a single solution surprisingly coincident with the ideal λ₁(G) (Table 5). On the one hand, the theory of linear programming is known to admit nonunique solutions [Lu73]. On the other, different matrices might well have the same dominant eigenvector (see Appendix A) in our practical case, too. Fortunately, the solution returned by the routine has turned out unique, hence the unique best solution to the averaging problem.

Moreover, this best solution is absolutely accurate (Table 5), thus turning a solution to the approximate PMA problem into the exact one and suggesting the best alternative to our former ecological practice [L19, Let23, LKBO18]. Therefore, it makes no sense to quantitatively compare the outcomes of the “versatile optimization tools” and the LP approach – the comparison should rather be qualitative, with an obvious conclusion.

this best solution is absolutely accurate (Table 5), thus turning a solution to the approximate PMA problem into the exact one and suggesting the best alternative to our former ecological practice [L19, Let23, LKBO18]. Therefore, it makes no sense to quantitatively compare the outcomes of the “versatile optimization tools” and the LP approach – the comparison should rather be qualitative, with an obvious conclusion, thus denying any possible reproaches to comparing the case of 12 annual PPMs [PZL22] with the current case of 13 PPMs.

The only “spoon of tar” in this “barrel of honey” is due to certain zeros in the optimal solution, i.e., zero entries to the average matrix G. While there are totally 4 vital rates that have zero minimal values among the 13 annual ones (Table 2), three of them, namely, for c, d, e, (Figure 1), are reached by the “ideal” solution (Table 5). The shares of c = 0, d = 0, and d = 0, in the set of 13 annual PPMs (Table 1) are equal respectively to 9/13, 7/14 and 1/13 (Appendix C). Whether the average matrix should take into account the most frequent zeros or all of them is a dialectical issue. The paradigm of biological diversity in general and the concept of polyvariant ontogeny [Zh86,95] in particular dictate retaining all the links in the life cycle graph (Figure1), hence the absence of occasional zeros in the average matrix.

Therefore, an optimistic mathematical hypothesis is that the same absolute result can be obtained in the eigenvalue constrained optimization problem when we substitute the minimal positive bounds for the zero ones. Testing this hypothesis might be an object of further studies.

Note also that the LP technique presented in Section 2.4 for a specific case study is general enough to be applied for any finite set of irreducible nonnegative fixed-pattern matrices, the specificity being only confined to the bounds and equality-type constraints. While the task is always to find the best approximation, a supertask is to motivate younger mathematicians for giving a theoretical explanation of why the best approximation turns into the exact solution.