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

Abstract

Given several nonnegative matrices with a single pattern of allocation among their zero/nonzero elements, the average matrix should have the same pattern, too. This is the first tenet of the pattern-multiplicative average (PMA) concept, while the second one suggests the multiplicative (or geometric) nature of averaging. The original concept of PMA was motivated by the practice of matrix population models as a tool to assess the population viability from long-term monitoring data. The task has reduced to searching for an approximate solution to an overdetermined system of polynomial equations for unknown elements of the average matrix (G), and hence to a nonlinear constrained minimization problem for the matrix norm. Former practical solutions faced certain technical problems, which required sophisticated algorithms but returned acceptable estimates. Now, we formulate (for the first time in ecological modeling and nonnegative matrix theory) the PMA problem as an eigenvalue approximation one and reduce it to a standard problem of linear programing (LP). The basic equation of averaging also determines the exact value of λ₁(G), the dominant eigenvalue of matrix G, and the corresponding eigenvector. These are bound by the well-known linear equations, which enable an LP formulation of the former nonlinear problem. The LP approach is realized for 13 fixed-pattern matrices gained in a case study of Androsace albana, an alpine short-lived perennial, monitored on permanent plots over 14 years. A standard software routine reveals the unique exact solution, rather than an approximate one, to the PMA problem, which turns the LP approach into ‘’the best of versatile optimization tools”. The exact solution turns out to be peculiar in reaching zero bounds for certain nonnegative entries of G, which deserves modified problem formulation separating the lower bounds from zero.

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 [1,2], where the theory of nonnegative matrices applies to population dynamics [3]. 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 [4]. 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) [4], despite the projection term having quite another meaning in matrix algebra [5]. All the matrix entries are nonnegative and called vital rates as they have a certain demographic sense [4]. 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 [5,6], whereby a rich repertoire of model properties and population characteristics can be obtained after the PPM has been calibrated from data [4,7]. 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” [8] (p. 1). When the data are of the “identified individuals” type [4], a great advantage of MPMs is that the calibration can be performed just from two successive observations, thus meeting Equation (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 problem 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 → ∞ [9,10,11,12], and it has led to what is now called the stochastic growth rate of a population in a randomly changing environment (see [4] 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 [13], and the task of averaging was reduced to a constrained nonlinear minimization problem for the matrix norm [8,13] (see Section 2.3). The norm has however not appeared quite fit for global optimization but required versatile optimization tools [14] 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), and 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 [15] (Section 2.1). Section 2.3 and Section 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 both are presented in Section 3, enabling the comparison of results, which reveals a surprising contrast and induces certain comments (Section 4), in particular, on how general the proposed LP technique might be.

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 [15] for 14 years (2009–2022). Although A. albana actually reproduces 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 [16]. Later, we proved such an exclusion to be correct within the integer-valued formalism [17], 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 because 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];\hspace{1em}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. Outcome of PPM calibration for A. albana based on 2009–2022 data (updating Table 2 from [15]).

Census Year, t	Matrix L(t): t → t + 1	λ1(L(t))	Vector x *, %
2009 j = 0	L0$=\left[\begin{array}{ccccc}0& 0& 0& 0& \raisebox{1ex}{$30$}\!\left/ \!\raisebox{-1ex}{$13$}\right.\\ \raisebox{1ex}{$8$}\!\left/ \!\raisebox{-1ex}{$37$}\right.& 0& 0& 0& \raisebox{1ex}{$40$}\!\left/ \!\raisebox{-1ex}{$13$}\right.\\ \raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$37$}\right.& \raisebox{1ex}{$22$}\!\left/ \!\raisebox{-1ex}{$110$}\right.& \raisebox{1ex}{$28$}\!\left/ \!\raisebox{-1ex}{$99$}\right.& 0& \raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$13$}\right.\\ 0& 0& \raisebox{1ex}{$7$}\!\left/ \!\raisebox{-1ex}{$99$}\right.& \raisebox{1ex}{$19$}\!\left/ \!\raisebox{-1ex}{$35$}\right.& 0\\ 0& 0& 0& \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$35$}\right.& 0\end{array}\right]$	0.5661	$\left[\begin{array}{c}10.61\\ 18.20\\ 16.99\\ 51.59\\ 2.60\end{array}\right]$
2010 j = 1	L1 = $\left[\begin{array}{ccccc}0& 0& 0& 0& \raisebox{1ex}{$19$}\!\left/ \!\raisebox{-1ex}{$1$}\right.\\ \raisebox{1ex}{$14$}\!\left/ \!\raisebox{-1ex}{$30$}\right.& 0& 0& 0& \raisebox{1ex}{$31$}\!\left/ \!\raisebox{-1ex}{$1$}\right.\\ \raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$30$}\right.& \raisebox{1ex}{$22$}\!\left/ \!\raisebox{-1ex}{$48$}\right.& \raisebox{1ex}{$17$}\!\left/ \!\raisebox{-1ex}{$55$}\right.& 0& \raisebox{1ex}{$0$}\!\left/ \!\raisebox{-1ex}{$1$}\right.\\ 0& 0& \raisebox{1ex}{$34$}\!\left/ \!\raisebox{-1ex}{$55$}\right.& \raisebox{1ex}{$23$}\!\left/ \!\raisebox{-1ex}{$26$}\right.& 0\\ 0& 0& 0& \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$26$}\right.& 0\end{array}\right]$	1.2283	$\left[\begin{array}{c}15.90\\ 31.99\\ 18.25\\ 32.83\\ 1.03\end{array}\right]$
2011 j = 2	L2 = $\left[\begin{array}{ccccc}0& 0& 0& 0& \raisebox{1ex}{$49$}\!\left/ \!\raisebox{-1ex}{$1$}\right.\\ \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$19$}\right.& 0& 0& 0& \raisebox{1ex}{$85$}\!\left/ \!\raisebox{-1ex}{$1$}\right.\\ \raisebox{1ex}{$6$}\!\left/ \!\raisebox{-1ex}{$19$}\right.& \raisebox{1ex}{$35$}\!\left/ \!\raisebox{-1ex}{$45$}\right.& \raisebox{1ex}{$21$}\!\left/ \!\raisebox{-1ex}{$43$}\right.& 0& \raisebox{1ex}{$25$}\!\left/ \!\raisebox{-1ex}{$1$}\right.\\ 0& 0& \raisebox{1ex}{$10$}\!\left/ \!\raisebox{-1ex}{$43$}\right.& \raisebox{1ex}{$48$}\!\left/ \!\raisebox{-1ex}{$57$}\right.& 0\\ 0& 0& 0& \raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$57$}\right.& 0\end{array}\right]$	1.5779	$\left[\begin{array}{c}17.20\\ 30.40\\ 39.39\\ 12.45\\ 0.55\end{array}\right]$
2012 j = 3	L3 = $\left[\begin{array}{ccccc}0& 0& 0& 0& \raisebox{1ex}{$19$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\\ \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$49$}\right.& 0& 0& 0& \raisebox{1ex}{$136$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\\ \raisebox{1ex}{$10$}\!\left/ \!\raisebox{-1ex}{$49$}\right.& \raisebox{1ex}{$45$}\!\left/ \!\raisebox{-1ex}{$86$}\right.& \raisebox{1ex}{$39$}\!\left/ \!\raisebox{-1ex}{$87$}\right.& 0& \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\\ 0& 0& \raisebox{1ex}{$28$}\!\left/ \!\raisebox{-1ex}{$87$}\right.& \raisebox{1ex}{$45$}\!\left/ \!\raisebox{-1ex}{$58$}\right.& 0\\ 0& 0& 0& \raisebox{1ex}{$6$}\!\left/ \!\raisebox{-1ex}{$58$}\right.& 0\end{array}\right]$	1.2641	$\left[\begin{array}{c}6.01\\ 43.15\\ 29.67\\ 19.56\\ 1.60\end{array}\right]$
2014 j = 5	L5 = $\left[\begin{array}{ccccc}0& 0& 0& 0& \raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$3$}\right.\\ 0& 0& 0& 0& \raisebox{1ex}{$19$}\!\left/ \!\raisebox{-1ex}{$3$}\right.\\ \raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$16$}\right.& \raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$98$}\right.& \raisebox{1ex}{$6$}\!\left/ \!\raisebox{-1ex}{$34$}\right.& 0& \raisebox{1ex}{$0$}\!\left/ \!\raisebox{-1ex}{$3$}\right.\\ 0& 0& \raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$34$}\right.& \raisebox{1ex}{$16$}\!\left/ \!\raisebox{-1ex}{$50$}\right.& 0\\ 0& 0& 0& \raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$50$}\right.& 0\end{array}\right]$	0.3988	$\left[\begin{array}{c}11.71\\ 56.64\\ 11.69\\ 17.46\\ 3.50\end{array}\right]$
2015 j = 6	L6 = $\left[\begin{array}{ccccc}0& 0& 0& 0& \raisebox{1ex}{$10$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\\ 0& 0& 0& 0& \raisebox{1ex}{$29$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\\ 0& \raisebox{1ex}{$10$}\!\left/ \!\raisebox{-1ex}{$19$}\right.& \raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$10$}\right.& 0& \raisebox{1ex}{$0$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\\ 0& 0& \raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$10$}\right.& \raisebox{1ex}{$17$}\!\left/ \!\raisebox{-1ex}{$20$}\right.& 0\\ 0& 0& 0& \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$10$}\right.& 0\end{array}\right]$	1.0679	$\left[\begin{array}{c}9.19\\ 26.66\\ 18.28\\ 41.94\\ 3.93\end{array}\right]$
2016 j = 7	L7 = $\left[\begin{array}{ccccc}0& 0& 0& 0& \raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\\ 0& 0& 0& 0& \raisebox{1ex}{$8$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\\ \raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$10$}\right.& \raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$29$}\right.& \raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$13$}\right.& 0& 0\\ 0& 0& \raisebox{1ex}{$8$}\!\left/ \!\raisebox{-1ex}{$13$}\right.& \raisebox{1ex}{$20$}\!\left/ \!\raisebox{-1ex}{$22$}\right.& 0\\ 0& 0& 0& \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$22$}\right.& 0\end{array}\right]$	0.9611	$\left[\begin{array}{c}5.26\\ 14.04\\ 6.02\\ 71.30\\ 3.37\end{array}\right]$
2017 j = 8	L8 =$\left[\begin{array}{ccccc}0& 0& 0& 0& \raisebox{1ex}{$12$}\!\left/ \!\raisebox{-1ex}{$1$}\right.\\ 0& 0& 0& 0& \raisebox{1ex}{$23$}\!\left/ \!\raisebox{-1ex}{$1$}\right.\\ \raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$3$}\right.& \raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$8$}\right.& \raisebox{1ex}{$8$}\!\left/ \!\raisebox{-1ex}{$12$}\right.& 0& 0\\ 0& 0& \raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$12$}\right.& \raisebox{1ex}{$21$}\!\left/ \!\raisebox{-1ex}{$28$}\right.& 0\\ 0& 0& 0& \raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$28$}\right.& 0\end{array}\right]$	1.1206	$\left[\begin{array}{c}12.93\\ 24.78\\ 42.13\\ 18.95\\ 1.21\end{array}\right]$
2018 j = 9	L9 = $\left[\begin{array}{ccccc}0& 0& 0& 0& \raisebox{1ex}{$13$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\\ 0& 0& 0& 0& \raisebox{1ex}{$38$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\\ \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$12$}\right.& \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$23$}\right.& \raisebox{1ex}{$0$}\!\left/ \!\raisebox{-1ex}{$13$}\right.& 0& 0\\ 0& 0& \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$13$}\right.& \raisebox{1ex}{$22$}\!\left/ \!\raisebox{-1ex}{$23$}\right.& 0\\ 0& 0& 0& \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$23$}\right.& 0\end{array}\right]$	0.9617	$\left[\begin{array}{c}13.23\\ 38.65\\ 2.89\\ 43.27\\ 1.96\end{array}\right]$
2019 j = 10	L10 = $\left[\begin{array}{ccccc}0& 0& 0& 0& \raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$1$}\right.\\ 0& 0& 0& 0& \raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$1$}\right.\\ \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$13$}\right.& \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$19$}\right.& \raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$4$}\right.& 0& 0\\ 0& 0& \raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$4$}\right.& \raisebox{1ex}{$18$}\!\left/ \!\raisebox{-1ex}{$23$}\right.& 0\\ 0& 0& \raisebox{1ex}{$0$}\!\left/ \!\raisebox{-1ex}{$4$}\right.& \raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$23$}\right.& 0\end{array}\right]$	0.8496	$\left[\begin{array}{c}18.45\\ \begin{array}{c}18.45\\ \begin{array}{c}6.84\\ 51.04\end{array}\end{array}\\ 5.22\end{array}\right]$
2020 j = 11	L11 = $\left[\begin{array}{ccccc}0& 0& 0& 0& \raisebox{1ex}{$8$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\\ \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.& 0& 0& 0& \raisebox{1ex}{$9$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\\ \raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.& \raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.& \raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$4$}\right.& 0& 0\\ 0& 0& \raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$4$}\right.& \raisebox{1ex}{$13$}\!\left/ \!\raisebox{-1ex}{$19$}\right.& 0\\ 0& 0& \raisebox{1ex}{$0$}\!\left/ \!\raisebox{-1ex}{$4$}\right.& \raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$19$}\right.& 0\end{array}\right]$	1.3008	$\left[\begin{array}{c}15.88\\ \begin{array}{c}21.94\\ \begin{array}{c}31.49\\ 25.53\end{array}\end{array}\\ 5.16\end{array}\right]$
2021 j = 12	L12 = $\left[\begin{array}{ccccc}0& 0& 0& 0& \raisebox{1ex}{$29$}\!\left/ \!\raisebox{-1ex}{$5$}\right.\\ \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$8$}\right.& 0& 0& 0& \raisebox{1ex}{$44$}\!\left/ \!\raisebox{-1ex}{$5$}\right.\\ \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$8$}\right.& \raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$5$}\right.& 0& 0& 0\\ 0& 0& \raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$6$}\right.& \raisebox{1ex}{$14$}\!\left/ \!\raisebox{-1ex}{$15$}\right.& 0\\ 0& 0& \raisebox{1ex}{$0$}\!\left/ \!\raisebox{-1ex}{$6$}\right.& \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$15$}\right.& 0\end{array}\right]$	1.1143	$\left[\begin{array}{c}14.86\\ \begin{array}{c}24.21\\ \begin{array}{c}10.36\\ 47.72\end{array}\end{array}\\ 2.85\end{array}\right]$

* To be explained in Section 2.4.

. Each of them obeys Equation (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 Equation (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 chronological order.

x(2022) = Prod x(2009).

(5)

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

G¹³ = Prod,

(6)

which justifies Equation (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 Equation (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. There is no reason however to seek this kind of relation among the PPMs, so that Equation (6), the basic equation of averaging, has no solution over the set of the 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 Constrained Minimization Problem

Any approximation invokes the question of how close the approximate solution is to the exact one. Although the latter is unknown, the equivalent form of Equation (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 }(\mathit{G})\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 lie 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}(t)\le g_{ij}\le \underset{t}{\max }{l}_{ij}\left(t\right),i,j=1,\dots ,5$$

(9)

(see Table 2).

Table 2. Bounds for the 10 unknown variables ensuing from 13 annual PPMs (Table 1).

Minimal	Vital Rate	Maximal
4/3	a	49
3	b	85
0	c	25
0	d	7/15
0	e	1
1/49	f	7/9
0	h	2/3
6/95	k	5/6
8/25	l	22/23
1/35	m	5/15

Table 2. Bounds for the 10 unknown variables ensuing from 13 annual PPMs (Table 1).

Minimal	Vital Rate	Maximal
4/3	a	49
3	b	85
0	c	25
0	d	7/15
0	e	1
1/49	f	7/9
0	h	2/3
6/95	k	5/6
8/25	l	22/23
1/35	m	5/15

So, we come to the constrained nonlinear minimization problem (8 and 9), where the bounds of B are defined in Table 2. Both standard routines for global optimization [18] and more sophisticated ones [14] can be applied for solving the problem (8 and 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 is a dominant eigenvector of the average matrix G, i.e., we have

G v = λ₁(G) v.

(10)

Acting by operator G successively on both sides of Equation (10), we obtain

(G)² v = λ₁(G)² v,
(G)³ v = λ₁(G)³ v,
…,
(G)¹³ v = λ₁(G)¹³ v.

(11)

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

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

(12)

which can be calculated from the data presented in Table 1. The optimal solution should be searched for among those matrices G that have 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 the 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 a 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, λ:

av₅ − λv₁ = 0,
bv₅ + dv₁ − λv₂ = 0,
cv₅ + ev₁+ fv₂+ hv₃ − λv₃ = 0,
kv₃ + lv₄ − λv₄ = 0,
mv₄ − λv₅ = 0,
i.e., for the ten unknown entries of G and the 11th formal variable λ.

(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[{}_0,\underset{t}{\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 the function linprog in MATLAB^® [19] (technical details in Appendix B).

3. Results

3.1. Case Study Outcome

Calibrated with 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 Section 2.3 and Section 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. Numerical inputs for the optimization problems.

Matrix Prod					λ1(G13), ρ0	Vector v*, %
0.021185585295608	0.039538528446472	0.086369212318887	0.321576284325812	0.312397941844407	0.31893645391 0.91584799085	29.2923
0.032875920640909	0.061354070510449	0.134019914355007	0.498983075380778	0.484789540051083		45.4532
0.003661007803335	0.006824023057585	0.014887199373520	0.055368303504470	0.054013601100565		5.0480
0.013845526439184	0.025919668563603	0.056669501627721	0.210813605518226	0.203676548113200		19.1962
0.000729124010650	0.001363464687049	0.002980796589266	0.011096313418442	0.010736266465211		1.0103

together with the dominant eigenvalue and the corresponding eigenvector.

3.3. Minimizing the Approximation Error as a Matrix Norm

The solution to the constraint minimization problem (8) has been obtained in [14] by the “versatile optimization tools” (

Table 4. Solutions to the PMA problem by various optimization techniques (adapted from [14]).

Optimization Method, Loss Function	Matrix G					λ1(G)	Approximation Error
Basin hopping, ${\mathsf{\Phi }\left(\mathit{G}\right)=\left|\right|{\mathit{G}}^{12}–\mathit{P}\mathit{r}\mathit{o}\mathit{d}\left|\right|}^2$	0	0	0	0	3.3309	0.8585	0.002374
	0.4530	0	0	0	7.8767
	0.0288	0.2936	0.1474	0	0
	0	0	0.1726	0.7589	0
	0	0	0	0.1034	0
Basin hopping, ${\mathsf{\Phi }\left(\mathit{G}\right)=\left|\right|{\mathit{G}}^{12}–\mathit{P}\mathit{r}\mathit{o}\mathit{d}\left|\right|}^2$ and penalty for constraint violations	0	0	0	0	3.3309	0.8585	0.002374
	0.4533	0	0	0	7.8757
	0.0287	0.22936	0.1474	0	0
	0	0	0.1726	0.7589	0
	0	0	0	0.1034	0
Basin hopping, S(G) = Φ(G)/Ψ(G)′, where $\Psi \left(\mathit{G}\right)={‖{\sum }_{j=0}^{11}{\mathit{L}}_{12}\dots {\mathit{L}}_{j+1}‖}^2$	0	0	0	0	3.3348	0.8584	0.002379
	0.4322	0	0	0	7.9666
	0.0363	0.28976	0.1485	0	0.0022
	0	0	0.1728	0.7587	0
	0	0	0	0.1034	0

) 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. Solutions to the PMA problem by the LP technique.

Inner Problem (16)	Matrix G(x)					λ1(G)	Approximation Error
$y_1\to \underset{\mathit{x}\in {\mathbb{B}}_1}{\min }{y}_1\left(\mathit{x}\right)$, $y_1\left(\mathit{x}\right)={\rho }_0-{\lambda }_1\left(\mathit{G}\left(\mathit{x}\right)\right)$	0	0	0	0	26.5530	0.915847990853247	0
	0	0	0	0	41.2023
	0	0.2777	0.6667	0	0
	0	0	0.0632	0.8992	0
	0	0	0	0.0482	0
$y_2\to \underset{\mathit{x}\in {\mathbb{B}}_2}{\min }{y}_2\left(\mathit{x}\right)$, $y_2\left(\mathit{x}\right)={\lambda }_1\left(\mathit{G}\left(\mathit{x}\right)\right)-{\rho }_0$	0	0	0	0	26.5530	0.915847990853247	1.1102 × 10−16
	0	0	0	0	41.2023
	0	0.2777	0.6667	0	0
	0	0	0.0632	0.8992	0
	0	0	0	0.0482	0

.

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., Equations (4) and (5). Therefore, we deal here with single-objective problems alone (in the terminology of [20]). Once having been formulated in mathematical terms [13], this objective comes immediately to the difference norm as the distance between two matrices following the fundamental of real analysis [21]. However, even the squared matrix norm has required “versatile optimization tools” [14] to achieve even moderate approximation quality in the corresponding nonlinear constrained minimization problem (Table 4).

Besides 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 [22]. In contrast, the LP approach is essentially deterministic [23]; hence, the results in 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. Because the matrix eigenvalue in general and the dominant eigenvalue in particular are not linear as a function of matrix elements (see Appendix A), the use of linear programming seems paradoxical at first glance. However, besides the trick with a change of variables (13), the very nature of eigenvalue–eigenvector relations is linear, resulting in the linear relationships among the 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 (12) 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. They 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 [23]. 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 to be 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 [8,15,16]. 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 [14] 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 in total 4 vital rates that have zero minimal values among the 13 annual ones (Table 2), 3 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 [24,25] in particular dictate retaining all the links in the life cycle graph (Figure 1), 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, the supertask is to motivate younger mathematicians to give a theoretical explanation of why the best approximation turns into the exact solution.