An Inquiry-Based Introduction to Linear Algebra and Applications

Subsection Application: Leslie Matrices and Population Modeling

The Leslie Matrix (also called the Leslie Model) is a powerful model for describing an age distributed growth of a population that is closed to migration. In a Leslie model, it is usually the case that only one gender (most often female) is considered. As an example, we will later consider a population of sheep that is being grown commercially. A natural question that we will address is how we can harvest the population to build a sustainable environment.

When working with populations, the matrices we use are often large. For large matrices, using the characteristic polynomial to calculate eigenvalues is too time and resource consuming to be practical, and we generally cannot find the exact values of the eigenvalues. As a result, approximation techniques are very important. In this section we will explore a method for approximating eigenvalues. The eigenvalues of a Leslie matrix are important because they describe the limiting or steady-state behavior of a population. The matrix and model were introduced by Patrick H. Leslie in “On the Use of Matrices in Certain Population Mathematics”, Leslie, P.H., Biometrika, Volume XXXIII, November 1945, pp. 183-212.

Subsection Introduction

We have used the characteristic polynomial to find the eigenvalues of a matrix, and for each eigenvalue row reduced a corresponding matrix to find the eigenvectors This method is only practical for small matrices — for more realistic applications approximation techniques are used. We investigate one such technique in this section, the power method.

Subsection The Power Method

While the examples we present in this text are small in order to highlight the concepts, matrices that appear in real life applications are often enormous. For example, in Google's PageRank algorithm that is used to determine relative rankings of the importance of web pages, matrices of staggering size are used (most entries in the matrices are zero, but the size of the matrices is still huge). Finding eigenvalues of such large matrices through the characteristic polynomial is impractical. In fact, finding the roots of all but the smallest degree characteristic polynomials is a very difficult problem. As a result, using the characteristic polynomial to find eigenvalues and then finding eigenvectors is not very practical in general, and it is often a better option to use a numeric approximation method. We will consider one such method in this section, the power method.

In Preview Activity 20.1, we saw an example of a matrix \(A = \left[ \begin{array}{cc} 2\amp 6 \\ 5\amp 3 \end{array} \right]\) so that the sequence \(\{\vx_k\}\text{,}\) where \(\vx_k = A \vx_{k-1}\text{,}\) converged to a dominant eigenvector of \(A\) for an initial guess vector \(\vx_0 = [1 \ 0]^{\tr}\text{.}\) The vectors \(\vx_i\) for \(i\) from \(1\) to \(6\) (with scaling) are approximately

Numerically we can see that the sequence \(\{\vx_k\}\) approaches the vector \([1 \ 1]^{\tr}\text{,}\) and Figure 20.1 illustrates this geometrically as well.

This method of successive approximations \(\vx_k = A \vx_{k-1}\) is called the power method (since we could write \(\vx_k\) as \(A^k \vx_0\)). Our task now is to show that this method works in general. In the next activity we restrict our argument to the \(2 \times 2\) case, and then discuss the general case afterwards.

Let \(A\) be an arbitrary \(2 \times 2\) matrix with two linearly independent eigenvectors \(\vv_1\) and \(\vv_2\) and corresponding eigenvalues \(\lambda_1\) and \(\lambda_2\text{,}\) respectively. We will also assume \(|\lambda_1| > |\lambda_2|\text{.}\) An eigenvalue whose absolute value is larger than that of any other eigenvalue is called a dominant eigenvalue. Any eigenvector for a dominant eigenvalue is called a dominant eigenvector. Before we show that our method can be used to approximate a dominant eigenvector, we recall that since \(\vv_1\) and \(\vv_2\) are eigenvectors corresponding to distinct eigenvalues, then \(\vv_1\) and \(\vv_2\) are linearly independent. So there exist scalars \(a_1\) and \(a_2\) such that

We have seen that for each positive integer \(k\) we can write \(\vx_n\) as

With this representation of \(\vx_0\) we can now see why the power method approximates a dominant eigenvector of \(A\text{.}\)

The power method is straightforward to implement, but it is not without its drawbacks. We began by assuming that we had a basis of eigenvectors of a matrix \(A\text{.}\) So we are also assuming that \(A\) is diagonalizable. We also assumed that \(A\) had a dominant eigenvalue \(\lambda_1\text{.}\) That is, if \(A\) is \(n \times n\) we assume that \(A\) has eigenvalues \(\lambda_1\text{,}\) \(\lambda_2\text{,}\) \(\ldots\text{,}\) \(\lambda_n\text{,}\) not necessarily distinct, with

and with \(\vv_i\) an eigenvector of \(A\) with eigenvalue \(\lambda_i\text{.}\) We could then write any initial guess \(\vx_0\) in the form

The initial guess is also called a seed.

Then

and

Notice that we are not actually calculating the vectors \(\vx_k\) here — this is a theoretical argument and we don't know \(\lambda_1\) and are not performing any scaling like we did in Preview Activity 20.1. We are assuming that \(\lambda_1\) is the dominant eigenvalue of \(A\text{,}\) though, so for each \(i\) the terms \(\left(\frac{\lambda_i}{\lambda_1}\right)^k\) converge to \(0\) as \(k\) goes to infinity. Thus,

for large values of \(k\text{,}\) which makes the sequence \(\{\vx_k\}\) converge to a vector in the direction of a dominant eigenvector \(\vv_1\) provided \(a_1 \neq 0\text{.}\) So we need to be careful enough to choose a seed that has a nonzero component in the direction of \(\vv_1\text{.}\) Of course, we generally don't know that our matrix is diagonalizable before we make these calculations, but for many matrices the sequence \(\{\vx_k\}\) will approach a dominant eigenvector.

The power method approximates a dominant eigenvector, and there are ways that we can approximate the dominant eigenvalue. Exercise 8 presents one way — by keeping track of the components of the \(\vx_k\) that have the largest absolute values, and the next activity shows another.

The result of Activity 20.3 is that, when the vectors in the sequence \(\{\vx_k\}\) approximate a dominant eigenvector of a matrix \(A\text{,}\) the quotients

approximate the dominant eigenvalue of \(A\text{.}\) The quotients in (20.4) are called Rayleigh quotients.

To summarize, the procedure for applying the power method for approximating a dominant eigenvector and dominant eigenvalue of a matrix \(A\) is as follows.

Step 1

Select an arbitrary nonzero vector \(\vx_0\) as an initial guess to a dominant eigenvector.

Step 2

Let \(\vx_1 = A\vx_0\text{.}\) Let \(k = 1\text{.}\)

Step 3

To avoid having the magnitudes of successive approximations become excessively large, scale this approximation \(\vx_k\text{.}\) That is, find the entry \(\alpha_k\) of \(\vx_k\) that is largest in absolute value. Then replace \(\vx_k\) by \(\frac{1}{\alpha_k} \vx_k\text{.}\)

Step 4

Calculate the Rayleigh quotient \(r_k = \frac{(A\vx_{k}) \cdot \vx_k}{\vx_k \cdot \vx_k}\text{.}\)

Step 5

Let let \(\vx_{k+1} = A \vx_k\text{.}\) Increase \(k\) by \(1\) and repeat Steps 3 through 5.

If the sequence \(\{\vx_k\}\) converges to a dominant eigenvector of \(A\text{,}\) then the sequence \(\{r_k\}\) converges to the dominant eigenvalue of \(A\text{.}\)

The power method can be useful for approximating a dominant eigenvector as long as the successive multiplications by \(A\) are fairly simple — for example, if many entries of \(A\) are zero.³⁸ The rate of convergence of the sequence \(\{\vx_k\}\) depends on the ratio \(\frac{\lambda_2}{\lambda_1}\text{.}\) If this ratio is close to \(1\text{,}\) then it can take many iterations before the power \(\left(\frac{\lambda_2}{\lambda_1}\right)^k\) makes the \(\vv_2\) term negligible. There are other methods for approximating eigenvalues and eigenvectors, e.g., the QR factorization, that we will not discuss at this point.

Subsection The Inverse Power Method

The power method only allows us to approximate the dominant eigenvalue and a dominant eigenvector for a matrix \(A\text{.}\) It is possible to modify this method to approximate other eigenvectors and eigenvalues under certain conditions. We consider an example in the next activity to motivate the general situation.

Activity 20.4 provides evidence that we can translate the matrix \(A\) having a dominant eigenvalue to a different matrix \(B\) with the same eigenvectors as \(A\) and with a dominant eigenvalue of our choosing. To see why, let \(A\) be an \(n \times n\) matrix with eigenvalues \(\lambda_1\text{,}\) \(\lambda_2\text{,}\) \(\ldots\text{,}\) \(\lambda_n\text{,}\) and let \(\alpha\) be any real number distinct from the eigenvalues. Let \(B = (A - \alpha I_n)^{-1}\text{.}\) In our example in Activity 20.4 the numbers

were the eigenvalues of \(B\text{,}\) and that if \(\vv_i\) is an eigenvector for \(A\) corresponding to the eigenvalue \(\lambda_i\text{,}\) then \(\vv_i\) is an eigenvector of \(B\) corresponding to the eigenvalue \(\frac{1}{\lambda_i - \alpha}\text{.}\) To see why, let \(\lambda\) be an eigenvalue of an \(n \times n\) matrix \(A\) with corresponding eigenvector \(\vv\text{.}\) Let \(\alpha\) be a scalar that is not an eigenvalue of \(A\text{,}\) and let \(B = (A - \alpha I_n)^{-1}\text{.}\) Now

So \(\frac{1}{\lambda - \alpha}\) is an eigenvalue of \(B\) with eigenvector \(\vv\text{.}\)

Now suppose that \(A\) is an \(n \times n\) matrix with eigenvalues \(\lambda_1\text{,}\) \(\lambda_2\text{,}\) \(\ldots\text{,}\) \(\lambda_n\text{,}\) and that we want to approximate an eigenvector and corresponding eigenvalue \(\lambda_i\) of \(A\text{.}\) If we can somehow find a value of \(\alpha\) so that \(|\lambda_i - \alpha| \lt |\lambda_j - \alpha|\) for all \(j \neq i\text{,}\) then \(\left| \frac{1}{\lambda_i - \alpha} \right| > \left| \frac{1}{\lambda_j - \alpha} \right|\) for any \(j \neq i\text{.}\) Thus, the matrix \(B = (A - \alpha I_n)^{-1}\) has \(\frac{1}{\lambda_i - \alpha}\) as its dominant eigenvalue and we can use the power method to approximate an eigenvector and the Rayleigh quotient to approximate the eigenvalue \(\frac{1}{\lambda_i - \alpha}\text{,}\) and hence approximate \(\lambda_i\text{.}\)

Subsection Examples

What follows are worked examples that use the concepts from this section.

Subsection Summary

• The power method is an iterative method that can be used to approximate the dominant eigenvalue of an \(n \times n\) matrix \(A\) that has \(n\) linearly independent eigenvectors and a dominant eigenvalue.

• To use the power method we start with a seed \(\vx_0\) and then calculate the sequence \(\{\vx_k\}\) of vectors, where \(\vx_k = A \vx_{k-1}\text{.}\) If \(\vx_0\) is chosen well, then the sequence \(\{\vx_k\}\) converges to a dominant eigenvector of \(A\text{.}\)

• If \(A\) is an \(n \times n\) matrix with eigenvalues \(\lambda_1\text{,}\) \(\lambda_2\text{,}\) \(\ldots\text{,}\) \(\lambda_n\text{,}\) to approximate an eigenvector of \(A\) corresponding to the eigenvalue \(\lambda_i\text{,}\) we apply the power method to the matrix \(B = (A - \alpha I_n)^{-1}\text{,}\) where \(\alpha\) is not a eigenvalue of \(A\) and \(\left| \frac{1}{\lambda_i - \alpha} \right| > \left| \frac{1}{\lambda_j - \alpha} \right|\) for any \(j \neq i\text{.}\)

Subsection Project: Managing a Sheep Herd

Sheep farming is a significant industry in New Zealand. New Zealand is reported to have the highest density of sheep in the world. Sheep can begin to reproduce after one year, and give birth only once per year. Table 20.8 gives Birth and Survival Rates for Female New Zealand Sheep (from G. Caughley, “Parameters for Seasonally Breeding Populations,” Ecology, 48, (1967), 834-839). Since sheep hardly ever live past 12 years, we will only consider the population through 12 years.

As sheep reproduce, they add to the 0-1 sheep (lamb) population. The potential to produce offspring is called fecundity (derived from the word fecund which generally refers to reproductive ability) and determines how many lamb are added to the population. Let \(F_k\) (the fecundity rate) be the rate at which females in age class \(k\) give birth to female offspring. Not all members of a given age group survive to the next age groups, so let \(s_k\) be the fraction of individuals that survives from age group \(k\) to age group \(k+1\text{.}\) With these ideas in mind, we can create a life cycle chart as in Figure 20.9 that illustrates how the population of sheep changes on a farm (for the sake of space, we illustrate with four age classes).

To model the sheep population, we need a few variables. Let \(n_1^{(0)}\) be the number of sheep in age group 0-1, \(n_2^{(0)}\) the number in age group 1-2, \(n_3\) the number in age group 2-3 and, in general, \(n_k^{(0)}\) the number of sheep in age group \((k-1)\)-\(k\) at some initial time (time \(0\)), and let

We wish to determine the populations in the different groups after one year. Let

where \(n_1^{(1)}\) denotes the number of sheep in age group 0-1, \(n_2^{(1)}\) the number of sheep in age group 1-2 and, in general, \(n_{k}^{(1)}\) the number of tilapia in age group \((k-1)\)-\(k\) after one year.

Notice that our matrix \(L\) has the form

Such a matrix is called a Leslie matrix.

Leslie matrices have certain useful properties, and one eigenvalue of a Leslie matrix can tell us a lot about the long-term behavior of the situation being modeled. You can take these properties as fact unless otherwise directed.

1. A Leslie matrix \(L\) has a unique positive eigenvalue \(\lambda_1\) with a corresponding eigenvector \(\vv_1\) whose entries are all positive.

2. If \(\lambda_i\) (\(i > 1\)) is any other eigenvalue (real or complex) of \(L\text{,}\) then \(|\lambda_i| \leq \lambda_1\text{.}\) If \(\lambda_1\) is the largest magnitude eigenvalue of a matrix \(L\text{,}\) we call \(\lambda_1\) a dominant eigenvalue of \(L\text{.}\)

3. If any two successive entries in the first row of \(L\) are both positive, then \(|\lambda_i| \lt \lambda_1\) for every \(i > 1\text{.}\) In this case we say that \(\lambda_1\) is a strictly dominant eigenvalue of \(L\text{.}\) In a Leslie model, this happens when the females in two successive age classes are fertile, which is almost always the case.

4. If \(\lambda_1\) is a strictly dominant eigenvalue, then \(\vx_k\) is approximately a scalar multiple of \(\vv_1\) for large values of \(k\text{,}\) regardless of the initial state \(\vx_0\text{.}\) In other words, large state vectors are close to eigenvectors for \(\lambda_1\text{.}\)

We can use these properties to determine the long-term behavior of the sheep herd.

Project Activity 20.7 indicates that, unchecked, the sheep population will grow without bound, roughly exponentially with ratio equal to the dominant eigenvalue of our Leslie matrix \(L\text{.}\) Of course, a sheep farmer cannot provide the physical environment or the resources to support an unlimited population of sheep. In addition, most sheep farmers cannot support themselves only by shearing sheep for the wool. Consequently, some harvesting of the sheep population each year for meat and skin is necessary. A sustainable harvesting policy allows for the regular harvesting of some sheep while maintaining the population at a stable level. It is necessary for the farmer to find an optimal harvesting rate to attain this stable population and the following activity leads us through an analysis of how such a harvesting rate can be determined.