text:
Linear Algebra/Topic: Speed of Calculating Determinants

The permutation expansion formula for computing determinants is useful for proving theorems, but the method of using row operations is a much better for finding the determinants of a large matrix. We can make this statement precise by considering, as computer algorithm designers do, the number of arithmetic operations that each method uses.

The speed of an algorithm is measured by finding how the time taken by the computer grows as the size of its input data set grows. For instance, how much longer will the algorithm take if we increase the size of the input data by a factor of ten, from a $$1000$$ row matrix to a $$10,000$$ row matrix or from $$10,000$$ to $$100,000$$? Does the time taken grow by a factor of ten, or by a factor of a hundred, or by a factor of a thousand? That is, is the time taken by the algorithm proportional to the size of the data set, or to the square of that size, or to the cube of that size, etc.?

Recall the permutation expansion formula for determinants.

 * $$\begin{array}{rl}

\begin{vmatrix} t_{1,1} &t_{1,2}  &\ldots  &t_{1,n}  \\ t_{2,1} &t_{2,2}  &\ldots  &t_{2,n}  \\ &\vdots                     \\ t_{n,1} &t_{n,2}  &\ldots  &t_{n,n} \end{vmatrix} &=\displaystyle\sum_{\text{permutations }\phi}\!\!\!\! t_{1,\phi(1)}t_{2,\phi(2)}\cdots t_{n,\phi(n)} \left|P_{\phi}\right|  \\ &= \begin{align} &t_{1,\phi_1(1)}\cdot t_{2,\phi_1(2)}\cdots t_{n,\phi_1(n)}\left|P_{\phi_1}\right|      \\ &\quad +t_{1,\phi_2(1)}\cdot t_{2,\phi_2(2)}\cdots t_{n,\phi_2(n)}\left|P_{\phi_2}\right|      \\ &\quad\vdots            \\ &\quad +t_{1,\phi_k(1)}\cdot t_{2,\phi_k(2)}\cdots t_{n,\phi_k(n)}\left|P_{\phi_k}\right| \end{align} \end{array}$$

There are $$n!=n\cdot(n-1)\cdot(n-2)\cdots 2\cdot 1$$ different $$ n $$-permutations. For numbers $$n$$ of any size at all, this is a large value; for instance, even if $$n$$ is only $$10$$ then the expansion has $$10!=3,628,800$$ terms, all of which are obtained by multiplying $$n$$ entries together. This is a very large number of multiplications (for instance, suggests $$10!$$ steps as a rough boundary for the limit of practical calculation). The factorial function grows faster than the square function. It grows faster than the cube function, the fourth power function, or any polynomial function. (One way to see that the factorial function grows faster than the square is to note that multiplying the first two factors in $$n!$$ gives $$n\cdot(n-1)$$, which for large $$n$$ is approximately $$n^2$$, and then multiplying in more factors will make it even larger. The same argument works for the cube function, etc.) So a computer that is programmed to use the permutation expansion formula, and thus to perform a number of operations that is greater than or equal to the factorial of the number of rows, would take very long times as its input data set grows.

In contrast, the time taken by the row reduction method does not grow so fast. This fragment of row-reduction code is in the computer language FORTRAN. The matrix is stored in the $$N \! \times \! N$$ array A. For each ROW between $$1$$ and $$N$$ parts of the program not shown here have already found the pivot entry $$A(ROW,COL)$$. Now the program does a row pivot.

-PIVINV\cdot \rho_{ROW} +\rho_i $$

(This code fragment is for illustration only and is incomplete. Still, analysis of a finished version that includes all of the tests and subcases is messier but gives essentially the same conclusion.)  PIVINV=1.0/A(ROW,COL) DO 10 I=ROW+1, N DO 20 J=I, N A(I,J)=A(I,J)-PIVINV*A(ROW,J) 20 CONTINUE 10 CONTINUE  The outermost loop (not shown) runs through $$N-1$$ rows. For each row, the nested $$I$$ and $$J$$ loops shown perform arithmetic on the entries in A that are below and to the right of the pivot entry. Assume that the pivot is found in the expected place, that is, that $$COL=ROW$$. Then there are $$(N-ROW)^2$$ entries below and to the right of the pivot. On average, ROW will be $$N/2$$. Thus, we estimate that the arithmetic will be performed about $$(N/2)^2$$ times, that is, will run in a time proportional to the square of the number of equations. Taking into account the outer loop that is not shown, we get the estimate that the running time of the algorithm is proportional to the cube of the number of equations.

Finding the fastest algorithm to compute the determinant is a topic of current research. Algorithms are known that run in time between the second and third power.

Speed estimates like these help us to understand how quickly or slowly an algorithm will run. Algorithms that run in time proportional to the size of the data set are fast, algorithms that run in time proportional to the square of the size of the data set are less fast, but typically quite usable, and algorithms that run in time proportional to the cube of the size of the data set are still reasonable in speed for not-too-big input data. However, algorithms that run in time (greater than or equal to) the factorial of the size of the data set are not practical for input of any appreciable size.

There are other methods besides the two discussed here that are also used for computation of determinants. Those lie outside of our scope. Nonetheless, this contrast of the two methods for computing determinants makes the point that although in principle they give the same answer, in practice the idea is to select the one that is fast.

Exercises
Most of these problems presume access to a computer.

/Solutions/


text:
Linear Algebra/Topic: Projective Geometry

There are geometries other than the familiar Euclidean one. One such geometry arose in art, where it was observed that what a viewer sees is not necessarily what is there. This is Leonardo da Vinci's The Last Supper.

What is there in the room, for instance where the ceiling meets the left and right walls, are lines that are parallel. However, what a viewer sees is lines that, if extended, would intersect. The intersection point is called the vanishing point. This aspect of perspective is also familiar as the image of a long stretch of railroad tracks that appear to converge at the horizon.

To depict the room, da Vinci has adopted a model of how we see, of how we project the three dimensional scene to a two dimensional image. This model is only a first approximation — it does not take into account that our retina is curved and our lens bends the light, that we have binocular vision, or that our brain's processing greatly affects what we see —  but nonetheless it is interesting, both artistically and mathematically.

The projection is not orthogonal, it is a central projection from a single point, to the plane of the canvas. (It is not an orthogonal projection since the line from the viewer to $$C$$ is not orthogonal to the image plane.) As the picture suggests, the operation of central projection preserves some geometric properties — lines project to lines. However, it fails to preserve some others — equal length segments can project to segments of unequal length; the length of $$AB$$ is greater than the length of $$BC$$ because the segment projected to $$AB$$ is closer to the viewer and closer things look bigger. The study of the effects of central projections is projective geometry. We will see how linear algebra can be used in this study.

There are three cases of central projection. The first is the projection done by a movie projector. We can think that each source point is "pushed" from the domain plane outward to the image point in the codomain plane. This case of projection has a somewhat different character than the second case, that of the artist "pulling" the source back to the canvas. In the first case $$S$$ is in the middle while in the second case $$I$$ is in the middle. One more configuration is possible, with $$P$$ in the middle. An example of this is when we use a pinhole to shine the image of a solar eclipse onto a piece of paper.

We shall take each of the three to be a central projection by $$P$$ of $$S$$ to $$I$$.

Consider again the effect of railroad tracks that appear to converge to a point. We model this with parallel lines in a domain plane $$S$$ and a projection via a $$P$$ to a codomain plane $$I$$ (The gray lines are parallel to $$S,I$$). All three projection cases appear here. The first picture below shows $$P$$ acting like a movie projector by pushing points from part of $$S$$ out to image points on the lower half of $$I$$. The middle picture shows $$P$$ acting like the artist by pulling points from another part of $$S$$ back to image points in the middle of $$I$$. In the third picture, $$P$$ acts like the pinhole, projecting points from $$S$$ to the upper part of $$I$$. This picture is the trickiest — the points that are projected near to the vanishing point are the ones that are far out on the bottom left of $$S$$. Points in $$S$$ that are near to the vertical gray line are sent high up on $$I$$. There are two awkward things about this situation. The first is that neither of the two points in the domain nearest to the vertical gray line (see below) has an image because a projection from those two is along the gray line that is parallel to the codomain plane (we sometimes say that these two are projected "to infinity"). The second awkward thing is that the vanishing point in $$I$$ isn't the image of any point from $$S$$ because a projection to this point would be along the gray line that is parallel to the domain plane (we sometimes say that the vanishing point is the image of a projection "from infinity").

For a better model, put the projector $$P$$ at the origin. Imagine that $$P$$ is covered by a glass hemispheric dome. As $$P$$ looks outward, anything in the line of vision is projected to the same spot on the dome. This includes things on the line between $$P$$ and the dome, as in the case of projection by the movie projector. It includes things on the line further from $$P$$ than the dome, as in the case of projection by the painter. It also includes things on the line that lie behind $$P$$, as in the case of projection by a pinhole. From this perspective $$P$$, all of the spots on the line are seen as the same point. Accordingly, for any nonzero vector $$\vec v\in\R^3$$, we define the associated point $$v$$ in the projective plane to be the set $$\bigl\{k\vec v:k\in\R,k\ne0\bigr\}$$ of nonzero vectors lying on the same line through the origin as $$\vec v$$. To describe a projective point we can give any representative member of the line, so that the projective point shown above can be represented in any of these three ways.
 * $$\begin{pmatrix}1\\2\\3\end{pmatrix}\qquad\begin{pmatrix}\frac13\\\frac23\\1\end{pmatrix}\qquad\begin{pmatrix}-2\\-4\\-6\end{pmatrix}$$

Each of these is a homogeneous coordinate vector for $$v$$.

This picture, and the above definition that arises from it, clarifies the description of central projection but there is something awkward about the dome model: what if the viewer looks down? If we draw $$P$$'s line of sight so that the part coming toward us, out of the page, goes down below the dome then we can trace the line of sight backward, up past $$P$$ and toward the part of the hemisphere that is behind the page. So in the dome model, looking down gives a projective point that is behind the viewer. Therefore, if the viewer in the picture above drops the line of sight toward the bottom of the dome then the projective point drops also and as the line of sight continues down past the equator, the projective point suddenly shifts from the front of the dome to the back of the dome. This discontinuity in the drawing means that we often have to treat equatorial points as a separate case. That is, while the railroad track discussion of central projection has three cases, the dome model has two.

We can do better than this. Consider a sphere centered at the origin. Any line through the origin intersects the sphere in two spots, which are said to be antipodal. Because we associate each line through the origin with a point in the projective plane, we can draw such a point as a pair of antipodal spots on the sphere. Below, the two antipodal spots are shown connected by a dashed line to emphasize that they are not two different points, the pair of spots together make one projective point. While drawing a point as a pair of antipodal spots is not as natural as the one-spot-per-point dome mode, on the other hand the awkwardness of the dome model is gone, in that if as a line of view slides from north to south, no sudden changes happen on the picture. This model of central projection is uniform — the three cases are reduced to one.

So far we have described points in projective geometry. What about lines? What a viewer at the origin sees as a line is shown below as a great circle, the intersection of the model sphere with a plane through the origin. (One of the projective points on this line is shown to bring out a subtlety. Because two antipodal spots together make up a single projective point, the great circle's behind-the-paper part is the same set of projective points as its in-front-of-the-paper part.) Just as we did with each projective point, we will also describe a projective line with a triple of reals. For instance, the members of this plane through the origin in $$\R^3$$
 * $$\left\{\begin{pmatrix}x\\y\\z\end{pmatrix}:x+y-z=0\right\}$$

project to a line that we can described with the triple $$(1\ 1\ -1)$$ (we use row vectors to typographically distinguish linesfrom points). In general, for any nonzero three-wide row vector $$\vec L$$ we define the associated line in the projective plane, to be the set $$L=\bigl\{k\vec L:k\in\R,k\ne0\bigr\}$$ of nonzero multiples of $$\vec L$$.

The reason that this description of a line as a triple is convienent is that in the projective plane, a point $$v$$ and a line $$L$$ are incident — the point lies on the line, the line passes throught the point — if and only if a dot product of their representatives $$v_1L_1+v_2L_2+v_3L_3$$ is zero (Problem 4 shows that this is independent of the choice of representatives $$\vec v$$ and $$\vec L$$). For instance, the projective point described above by the column vector with components $$1,2,3$$ lies in the projective line described by $$(1\ 1\ -1)$$, simply because any vector in $$\R^3$$ whose components are in ratio $$1:2:3$$ lies in the plane through the origin whose equation is of the form $$1k\cdot x+1k\cdot y-1k\cdot z=0$$ for any nonzero $$k$$. That is, the incidence formula is inherited from the three-space lines and planes of which $$v$$ and $$L$$ are projections.

Thus, we can do analytic projective geometry. For instance, the projective line $$L=(1\ 1\ -1)$$ has the equation $$1v_1+1v_2-1v_3=0$$, because points incident on the line are characterized by having the property that their representatives satisfy this equation. One difference from familiar Euclidean analytic geometry is that in projective geometry we talk about the equation of a point. For a fixed point like
 * $$v=\begin{pmatrix}1\\2\\3\end{pmatrix}$$

the property that characterizes lines through this point (that is, lines incident on this point) is that the components of any representatives satisfy $$1L_1+2L_2+3L_3=0$$ and so this is the equation of $$v$$.

This symmetry of the statements about lines and points brings up the Duality Principle of projective geometry: in any true statement, interchanging "point" with "line" results in another true statement. For example, just as two distinct points determine one and only one line, in the projective plane, two distinct lines determine one and only one point. Here is a picture showing two lines that cross in antipodal spots and thus cross at one projective point. $$(*)$$ Contrast this with Euclidean geometry, where two distinct lines may have a unique intersection or may be parallel. In this way, projective geometry is simpler, more uniform, than Euclidean geometry.

That simplicity is relevant because there is a relationship between the two spaces: the projective plane can be viewed as an extension of the Euclidean plane. Take the sphere model of the projective plane to be the unit sphere in $$\R^3$$ and take Euclidean space to be the plane $$z=1$$. This gives us a way of viewing some points in projective space as corresponding to points in Euclidean space, because all of the points on the plane are projections of antipodal spots from the sphere. $$(**)$$ Note though that projective points on the equator don't project up to the plane. Instead, these project "out to infinity". We can thus think of projective space as consisting of the Euclidean plane with some extra points adjoined — the Euclidean plane is embedded in the projective plane. These extra points, the equatorial points, are the ideal points or points at infinity and the equator is the ideal line or line at infinity (note that it is not a Euclidean line, it is a projective line).

The advantage of the extension to the projective plane is that some of the awkwardness of Euclidean geometry disappears. For instance, the projective lines shown above in $$(*)$$ cross at antipodal spots, a single projective point, on the sphere's equator. If we put those lines into $$(**)$$ then they correspond to Euclidean lines that are parallel. That is, in moving from the Euclidean plane to the projective plane, we move from having two cases, that lines either intersect or are parallel, to having only one case, that lines intersect (possibly at a point at infinity).

The projective case is nicer in many ways than the Euclidean case but has the problem that we don't have the same experience or intuitions with it. That's one advantage of doing analytic geometry, where the equations can lead us to the right conclusions. Analytic projective geometry uses linear algebra. For instance, for three points of the projective plane $$t,u,v$$, setting up the equations for those points by fixing vectors representing each, shows that the three are collinear — incident in a single line — if and only if the resulting three-equation system has infinitely many row vector solutions representing that line. That, in turn, holds if and only if this determinant is zero.
 * $$\begin{vmatrix}t_1&u_1&v_1\\t_2&u_2&v_2\\t_3&u_3&v_3\end{vmatrix}$$

Thus, three points in the projective plane are collinear if and only if any three representative column vectors are linearly dependent. Similarly (and illustrating the Duality Principle), three lines in the projective plane are incident on a single point if and only if any three row vectors representing them are linearly dependent.

The following result is more evidence of the "niceness" of the geometry of the projective plane, compared to the Euclidean case. These two triangles are said to be in perspective from $$P$$ because their corresponding vertices are collinear. Consider the pairs of corresponding sides: the sides $$T_1U_1,T_2U_2$$, the sides $$T_1V_1,T_2V_2$$ , and the sides $$U_1V_1,U_2V_2$$. Desargue's Theorem is that when the three pairs of corresponding sides are extended to lines, they intersect (shown here as the point $$TU$$, the point $$TV$$ , and the point $$UV$$), and further, those three intersection points are collinear. We will prove this theorem, using projective geometry. (These are drawn as Euclidean figures because it is the more familiar image. To consider them as projective figures, we can imagine that, although the line segments shown are parts of great circles and so are curved, the model has such a large radius compared to the size of the figures that the sides appear in this sketch to be straight.)

For this proof, we need a preliminary lemma : if $$W,X,Y,Z$$ are four points in the projective plane (no three of which are collinear) then there are homogeneous coordinate vectors $$\vec w,\vec x,\vec y,\vec z$$ for the projective points, and a basis $$B$$ for $$\R^3$$, satisfying this.

\text{Rep}_B(\vec w)=\begin{pmatrix}1\\0\\0\end{pmatrix}\quad \text{Rep}_B(\vec x)=\begin{pmatrix}0\\1\\0\end{pmatrix}\quad \text{Rep}_B(\vec y)=\begin{pmatrix}0\\0\\1\end{pmatrix}\quad \text{Rep}_B(\vec z)=\begin{pmatrix}1\\1\\1\end{pmatrix} $$ The proof is straightforward. Because $$W,X,Y$$ are not on the same projective line, any homogeneous coordinate vectors $$\vec{w}_0,\vec{x}_0,\vec{y}_0$$ do not line on the same plane through the origin in $$\R^3$$ and so form a spanning set for $$\R^3$$. Thus any homogeneous coordinate vector for $$Z$$ can be written as a combination $$\vec{z}_0=a\vec{w}_0+b\vec{x}_0+c\vec{y}_0$$. Then, we can take
 * $$\begin{align}\vec w&=a\vec{w}_0\\\vec x&=b\vec{x}_0\\\vec y&=c\vec{y}_0\\\vec z&=\vec{z}_0\end{align}$$

where the basis is $$B=\langle\vec w,\vec x,\vec y\rangle$$.

Now, to prove of Desargue's Theorem, use the lemma to fix homogeneous coordinate vectors and a basis.

\text{Rep}_B(\vec{t}_1)=\begin{pmatrix}1\\0\\0\end{pmatrix}\quad \text{Rep}_B(\vec{u}_1)=\begin{pmatrix}0\\1\\0\end{pmatrix}\quad \text{Rep}_B(\vec{v}_1)=\begin{pmatrix}0\\0\\1\end{pmatrix}\quad \text{Rep}_B(\vec{o})=\begin{pmatrix}1\\1\\1\end{pmatrix} $$ Because the projective point $$T_2$$ is incident on the projective line $$OT_1$$, any homogeneous coordinate vector for $$T_2$$ lies in the plane through the origin in $$\R^3$$ that is spanned by homogeneous coordinate vectors of $$O$$ and $$T_1$$ :
 * $$\text{Rep}_B(\vec{t}_2)=a\begin{pmatrix}1\\1\\1\end{pmatrix}+b\begin{pmatrix}1\\0\\0\end{pmatrix}$$

for some scalars $$a$$ and $$b$$. That is, the homogenous coordinate vectors of members $$T_2$$ of the line $$OT_1$$ are of the form on the left below, and the forms for $$U_2,V_2$$ are similar.

\text{Rep}_B(\vec{t}_2)=\begin{pmatrix}t_2\\1\\1\end{pmatrix}\qquad \text{Rep}_B(\vec{u}_2)=\begin{pmatrix}1\\u_2\\1\end{pmatrix}\qquad \text{Rep}_B(\vec{v}_2)=\begin{pmatrix}1\\1\\v_2\end{pmatrix} $$ The projective line $$T_1U_1$$ is the image of a plane through the origin in $$\R^3$$. A quick way to get its equation is to note that any vector in it is linearly dependent on the vectors for $$T_1$$ and $$U_1$$ and so this determinant is zero.
 * $$\begin{vmatrix}1&0&x\\0&1&y\\0&0&z\end{vmatrix}=0\qquad\implies\qquad z=0$$

The equation of the plane in $$\R^3$$ whose image is the projective line $$T_2U_2$$ is this.
 * $$\begin{vmatrix}t_2&1&x\\1&u_2&y\\1&1&z\end{vmatrix}=0\qquad\implies\qquad(1-u_2)x+(1-t_2)y+(t_2u_2-1)z=0$$

Finding the intersection of the two is routine.
 * $$T_1U_1\cap T_2U_2=\begin{pmatrix}t_2-1\\1-u_2\\0\end{pmatrix}$$

(This is, of course, the homogeneous coordinate vector of a projective point.) The other two intersections are similar.
 * $$T_1V_1\cap T_2V_2=\begin{pmatrix}1-t_2\\0\\v_2-1\end{pmatrix}\qquad U_1V_1\cap U_2V_2=\begin{pmatrix}0\\u_2-1\\1-v_2\end{pmatrix}$$

The proof is finished by noting that these projective points are on one projective line because the sum of the three homogeneous coordinate vectors is zero.

Every projective theorem has a translation to a Euclidean version, although the Euclidean result is often messier to state and prove. Desargue's theorem illustrates this. In the translation to Euclidean space, the case where $$O$$ lies on the ideal line must be treated separately for then the lines $$T_1T_2,U_1U_2,V_1V_2$$ are parallel.

The parenthetical remark following the statement of Desargue's Theorem suggests thinking of the Euclidean pictures as figures from projective geometry for a model of very large radius. That is, just as a small area of the earth appears flat to people living there, the projective plane is also "locally Euclidean".

Although its local properties are the familiar Euclidean ones, there is a global property of the projective plane that is quite different. The picture below shows a projective point. At that point is drawn an $$xy$$-axis. There is something interesting about the way this axis appears at the antipodal ends of the sphere. In the northern hemisphere, where the axis are drawn in black, a right hand put down with fingers on the $$x$$-axis will have the thumb point along the $$y$$-axis. But the antipodal axis has just the opposite: a right hand placed with its fingers on the $$x$$-axis will have the thumb point in the wrong way, instead, it is a left hand that works. Briefly, the projective plane is not orientable: in this geometry, left and right handedness are not fixed properties of figures. The sequence of pictures below dramatizes this non-orientability. They sketch a trip around this space in the direction of the $$y$$ part of the $$xy$$-axis. (Warning: the trip shown is not halfway around, it is a full circuit. True, if we made this into a movie then we could watch the northern hemisphere spots in the drawing above gradually rotate about halfway around the sphere to the last picture below. And we could watch the southern hemisphere spots in the picture above slide through the south pole and up through the equator to the last picture. But: the spots at either end of the dashed line are the same projective point. We don't need to continue on much further; we are pretty much back to the projective point where we started by the last picture.) $$\implies$$ $$\implies$$ At the end of the circuit, the $$x$$ part of the $$xy$$-axes sticks out in the other direction. Thus, in the projective plane we cannot describe a figure as right-{} or left-handed (another way to make this point is that we cannot describe a spiral as clockwise or counterclockwise).

This exhibition of the existence of a non-orientable space raises the question of whether our universe is orientable: is it possible for an astronaut to leave right-handed and return left-handed? An excellent nontechnical reference is. A classic science fiction story about orientation reversal is.

So projective geometry is mathematically interesting, in addition to the natural way in which it arises in art. It is more than just a technical device to shorten some proofs. For an overview, see. The approach we've taken here, the analytic approach, leads to quick theorems and — most importantly for us — illustrates the power of linear algebra (see, , and ). But another approach, the synthetic approach of deriving the results from an axiom system, is both extraordinarily beautiful and is also the historical route of development. Two fine sources for this approach are or. An interesting and easy application is.

Exercises
/Solutions/

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Linear Algebra/Topic: Speed of Gauss' Method

We are using Gauss' Method to solve the linear systems in this book because it is easy to understand, easily shown to give the right answers, and fast. It is fast in that, in all the by-hand calculations we have needed, we have gotten the answers in only a few steps, taking only a few minutes. However, scientists and engineers who solve linear systems in practice must have a method that is fast enough for large systems, with 1000 equations or 10,000 equations or even 100,000 equations. These systems are solved on a computer, so the speed of the machine helps, but nonetheless the speed of the method used is a major consideration, and is sometimes the factor limiting which problems can be solved.

The speed of an algorithm is usually measured by finding how the time taken to solve problems grows as the size of the input data set grows. That is, how much longer will the algorithm take if we increase the size of the input data by a factor of ten, say from a 1000-equation system to a 10,000-equation system, or from 10,000 to 100,000? Does the time taken grow ten times, or a hundred times, or a thousand times? Is the time taken by the algorithm proportional to the size of the data set, or to the square of that size, or to the cube of that size, etc.?

Here is a fragment of Gauss' Method code, implemented in the computer language FORTRAN. The coefficients of the linear system are stored in the $$N \! \times \! N$$ array A, and the constants are stored in the $$N \! \times \! 1$$ array B. For each ROW between $$1$$ and $$N$$ this program has already found the pivot entry $$A(ROW,COL)$$. Now it will pivot.

-PIVINV\cdot \rho_{ROW} +\rho_{i} $$

(This code fragment is for illustration only, and is incomplete. For example, see the later topic on the Accuracy of Gauss' Method. Nonetheless, this fragment will do for our purposes because analysis of finished versions, including all the tests and sub-cases, is messier but gives essentially the same result.)

The outermost loop (not shown) runs through $$N-1$$ rows. For each of these rows, the shown loops perform arithmetic on the entries in A that are below and to the right of the pivot entry (and also on the entries in B, but to simplify the analysis we will not count those operations---see Exercise ). We will assume the pivot is found in the usual place, that is, that $$COL=ROW$$ (as above, analysis of the general case is messier but gives essentially the same result). That means there are $$(N-ROW)^2$$ entries to perform arithmetic on. On average, ROW will be $$N/2$$. Thus we estimate the nested loops above will run something like $$(N/2)^2$$ times, that is, will run in a time proportional to the square of the number of equations. Taking into account the outer loop that is not shown, we get the estimate that the running time of the algorithm is proportional to the cube of the number of equations.

Algorithms that run in time directly proportional to the size of the data set are fast, algorithms that run in time proportional to the square of the size of the data set are less fast, but typically quite usable, and algorithms that run in time proportional to the cube of the size of the data set are still reasonable in speed.

Speed estimates like these are a good way of understanding how quickly or slowly an algorithm can be expected to run on average. There are special cases, however, of systems on which the above Gauss' method code is especially fast, so there may be factors about a problem that make it especially suitable for this kind of solution.

In practice, the code found in computer algebra systems, or in the standard packages, implements a variant on Gauss' method, called triangular factorization. To state this method requires the language of matrix algebra, which we will not see until Chapter Three. Nonetheless, the above code is conceptually quite close to that usually used in applications.

There have been some theoretical speed-ups in the running time required to solve linear systems. Algorithms other than Gauss' method have been invented that take a time proportional not to the cube of the size of the data set, but instead to the (approximately) $$2.7$$ power (this is still under active research, so this exponent may come down somewhat over time). However, these theoretical improvements have not come into widespread use, in part because the new methods take a quite large data set before they overtake Gauss' method (although they will outperform Gauss' method on very large sets, there is some startup overhead that keeps them from being faster on the systems that have, so far, been solved in practice).

Exercises
/Solutions/


text:
Linear Algebra/Complex Vector Spaces

This chapter requires that we factor polynomials. Of course, many polynomials do not factor over the real numbers; for instance, $$ x^2+1 $$ does not factor into the product of two linear polynomials with real coefficients. For that reason, we shall from now on take our scalars from the complex numbers.

That is, we are shifting from studying vector spaces over the real numbers to vector spaces over the complex numbers— in this chapter vector and matrix entries are complex.

Any real number is a complex number and a glance through this chapter shows that most of the examples use only real numbers. Nonetheless, the critical theorems require that the scalars be complex numbers, so the first section below is a quick review of complex numbers.

In this book we are moving to the more general context of taking scalars to be complex only for the pragmatic reason that we must do so in order to develop the representation. We will not go into using other sets of scalars in more detail because it could distract from our goal. However, the idea of taking scalars from a structure other than the real numbers is an interesting one. Delightful presentations taking this approach are in and.

title: Linear Algebra/Complex Vector Spaces


text:
Linear Algebra/Factoring and Complex Numbers: A Review

''This subsection is a review only and we take the main results as known. For proofs, see or .''

Just as integers have a division operation— e.g., "$$ 4 $$ goes $$ 5 $$ times into $$ 21 $$ with remainder $$ 1 $$"— so do polynomials.

In this book constant polynomials, including the zero polynomial, are said to have degree $$ 0 $$. (This is not the standard definition, but it is convienent here.)

The point of the integer division statement "$$ 4 $$ goes $$ 5 $$ times into $$ 21 $$ with remainder $$ 1 $$" is that the remainder is less than $$ 4 $$— while $$ 4 $$ goes $$ 5 $$ times, it does not go $$ 6 $$ times. In the same way, the point of the polynomial division statement is its final clause.

If a divisor $$ m(x) $$ goes into a dividend $$ c(x) $$ evenly, meaning that $$ r(x) $$ is the zero polynomial, then $$ m(x) $$ is a factor of $$ c(x) $$. Any root of the factor (any $$ \lambda\in\mathbb{R} $$ such that $$ m(\lambda)=0 $$) is a root of $$ c(x) $$ since $$ c(\lambda)=m(\lambda)\cdot q(\lambda)=0 $$. The prior corollary immediately yields the following converse.

Finding the roots and factors of a high-degree polynomial can be hard. But for second-degree polynomials we have the quadratic formula: the roots of $$ ax^2+bx+c $$ are

\lambda_1=\frac{-b+\sqrt{b^2-4ac}}{2a} \qquad \lambda_2=\frac{-b-\sqrt{b^2-4ac}}{2a} $$

(if the discriminant $$ b^2-4ac $$ is negative then the polynomial has no real number roots). A polynomial that cannot be factored into two lower-degree polynomials with real number coefficients is irreducible over the reals.

Note the analogy with the prime factorization of integers. In both cases, the uniqueness clause is very useful.

While $$ x^2+1 $$ has no real roots and so doesn't factor over the real numbers, if we imagine a root— traditionally denoted $$ i $$ so that $$ i^2+1=0 $$— then $$ x^2+1 $$ factors into a product of linears $$ (x-i)(x+i) $$.

So we adjoin this root $$ i $$ to the reals and close the new system with respect to addition, multiplication, etc. (i.e., we also add $$ 3+i $$, and $$ 2i $$, and $$ 3+2i $$, etc., putting in all linear combinations of $$1$$ and $$i$$). We then get a new structure, the complex numbers, denoted $$ \mathbb{C} $$.

In $$\mathbb{C}$$ we can factor (obviously, at least some) quadratics that would be irreducible if we were to stick to the real numbers. Surprisingly, in $$ \mathbb{C} $$ we can not only factor $$ x^2+1 $$ and its close relatives, we can factor any quadratic.

ax^2+bx+c= a\cdot \big(x-\frac{-b+\sqrt{b^2-4ac}}{2a}\big) \cdot \big(x-\frac{-b-\sqrt{b^2-4ac}}{2a}\big) $$


text:
Linear Algebra/Complex Representations

Recall the definitions of the complex number addition

(a+bi)\,+\,(c+di)=(a+c)+(b+d)i $$

and multiplication.

 * $$\begin{array}{rl}

(a+bi)(c+di) &=ac+adi+bci+bd(-1) \\ &=(ac-bd)+(ad+bc)i \end{array}$$

Handling scalar operations with those rules, all of the operations that we've covered for real vector spaces carry over unchanged.

Everything else from prior chapters that we can, we shall also carry over unchanged. For instance, we shall call the ordered set of vectors

\left\langle \begin{pmatrix} 1+0i \\ 0+0i \\ \vdots \\ 0+0i \end{pmatrix}, \begin{pmatrix} 0+0i \\ 1+0i \\ \vdots \\ 0+0i \end{pmatrix}, \dots, \begin{pmatrix} 0+0i \\ 0+0i \\ \vdots \\ 1+0i \end{pmatrix} \right\rangle $$

the standard basis for $$ \mathbb{C}^n $$ as a vector space over $$\mathbb{C}$$ and again denote it $$\mathcal{E}_n $$.

title: Linear Algebra/Complex Representations


text:
Linear Algebra/Introduction to Similarity

While studying matrix equivalence, we have shown that for any homomorphism there are bases $$B$$ and $$D$$ such that the representation matrix has a block partial-identity form.

{\rm Rep}_{B,D}(h) = \left(\begin{array}{c|c} \textit{Identity} &\textit{Zero}   \\ \hline \textit{Zero}      &\textit{Zero} \end{array}\right) $$

This representation describes the map as sending $$ c_1\vec{\beta}_1+\dots+c_n\vec{\beta}_n $$ to $$c_1\vec{\delta}_1+\dots+c_k\vec{\delta}_k+\vec{0}+\dots+\vec{0} $$, where $$n$$ is the dimension of the domain and $$ k $$ is the dimension of the range. So, under this representation the action of the map is easy to understand because most of the matrix entries are zero.

This chapter considers the special case where the domain and the codomain are equal, that is, where the homomorphism is a transformation. In this case we naturally ask to find a single basis $$ B $$ so that $$ {\rm Rep}_{B,B}(t) $$ is as simple as possible (we will take "simple" to mean that it has many zeroes). A matrix having the above block partial-identity form is not always possible here. But we will develop a form that comes close, a representation that is nearly diagonal.


text:
Linear Algebra/Definition and Examples of Similarity

Definition and Examples
We've defined $$ H $$ and $$ \hat{H} $$ to be matrix-equivalent if there are nonsingular matrices $$ P $$ and $$ Q $$ such that $$ \hat{H}=PHQ $$. That definition is motivated by this diagram

 * [[Image:Linalg__matrix_equivalent_cd_1.png|x150px]]

showing that $$H$$ and $$\hat{H}$$ both represent $$h$$ but with respect to different pairs of bases. We now specialize that setup to the case where the codomain equals the domain, and where the codomain's basis equals the domain's basis.

 * [[Image:Linalg__matrix_equivalent_cd_2.png|x150px]]

To move from the lower left to the lower right we can either go straight over, or up, over, and then down. In matrix terms,

{\rm Rep}_{D,D}(t) ={\rm Rep}_{B,D}(\mbox{id})\;{\rm Rep}_{B,B}(t)\;\bigl({\rm Rep}_{B,D}(\mbox{id})\bigr)^{-1} $$

(recall that a representation of composition like this one reads right to left).

Since nonsingular matrices are square, the similar matrices $$T$$ and $$S$$ must be square and of the same size.

Since matrix similarity is a special case of matrix equivalence, if two matrices are similar then they are equivalent. What about the converse: must matrix equivalent square matrices be similar? The answer is no. The prior example shows that the similarity classes are different from the matrix equivalence classes, because the matrix equivalence class of the identity consists of all nonsingular matrices of that size. Thus, for instance, these two are matrix equivalent but not similar.

T= \begin{pmatrix} 1 &0  \\ 0  &1 \end{pmatrix} \qquad S= \begin{pmatrix} 1 &2  \\ 0  &3 \end{pmatrix} $$

So some matrix equivalence classes split into two or more similarity classes— similarity gives a finer partition than does equivalence. This picture shows some matrix equivalence classes subdivided into similarity classes.

To understand the similarity relation we shall study the similarity classes. We approach this question in the same way that we've studied both the row equivalence and matrix equivalence relations, by finding a canonical form for representatives of the similarity classes, called Jordan form. With this canonical form, we can decide if two matrices are similar by checking whether they reduce to the same representative. We've also seen with both row equivalence and matrix equivalence that a canonical form gives us insight into the ways in which members of the same class are alike (e.g., two identically-sized matrices are matrix equivalent if and only if they have the same rank).

Exercises
/Solutions/


text:
Linear Algebra/Diagonalizability

The prior subsection defines the relation of similarity and shows that, although similar matrices are necessarily matrix equivalent, the converse does not hold. Some matrix-equivalence classes break into two or more similarity classes (the nonsingular $$n \! \times \! n$$ matrices, for instance). This means that the canonical form for matrix equivalence, a block partial-identity, cannot be used as a canonical form for matrix similarity because the partial-identities cannot be in more than one similarity class, so there are similarity classes without one. This picture illustrates. As earlier in this book, class representatives are shown with stars. We are developing a canonical form for representatives of the similarity classes. We naturally try to build on our previous work, meaning first that the partial identity matrices should represent the similarity classes into which they fall, and beyond that, that the representatives should be as simple as possible. The simplest extension of the partial-identity form is a diagonal form.

That example shows that a diagonal form will not do for a canonical form— we cannot find a diagonal matrix in each matrix similarity class. However, the canonical form that we are developing has the property that if a matrix can be diagonalized then the diagonal matrix is the canonical representative of the similarity class. The next result characterizes which maps can be diagonalized.

In the next subsection, we will expand on that example by considering more closely the property of Corollary 2.4. This includes seeing another way, the way that we will routinely use, to find the $$\lambda$$'s.

Exercises
/Solutions/

title: Linear Algebra/Diagonalizability


text:
Linear Algebra/Eigenvalues and Eigenvectors

In this subsection we will focus on the property of Corollary 2.4.

("Eigen" is German for "characteristic of" or "peculiar to"; some authors call these characteristic values and vectors. No authors call them "peculiar".)

That example shows why the "non-$$\vec{0}$$" appears in the definition. Disallowing $$ \vec{0} $$ as an eigenvector eliminates trivial eigenvalues.

The next example illustrates the basic tool for finding eigenvectors and eigenvalues.

Problem 11 checks that the characteristic polynomial of a transformation is well-defined, that is, any choice of basis yields the same polynomial.

Notice the familiar form of the sets of eigenvectors in the above examples.

By the lemma, if two eigenvectors $$\vec{v}_1$$ and $$\vec{v}_2$$ are associated with the same eigenvalue then any linear combination of those two is also an eigenvector associated with that same eigenvalue. But, if two eigenvectors $$ \vec{v}_1 $$ and $$ \vec{v}_2 $$ are associated with different eigenvalues then the sum $$ \vec{v}_1+\vec{v}_2 $$ need not be related to the eigenvalue of either one. In fact, just the opposite. If the eigenvalues are different then the eigenvectors are not linearly related.

Exercises
/Solutions/


text: User talk:68.42.76.159

title: User talk:68.42.76.159


text:
Linear Algebra/Nilpotence

The goal of this chapter is to show that every square matrix is similar to one that is a sum of two kinds of simple matrices. The prior section focused on the first kind, diagonal matrices. We now consider the other kind.

title: Linear Algebra/Nilpotence


text:
Adventist Youth Honors Answer Book/Insect/Diptera

|


text:
Linear Algebra/Self-Composition

This subsection is optional, although it is necessary for later material in this section and in the next one.

A linear transformations $$t:V\to V$$, because it has the same domain and codomain, can be iterated. That is, compositions of $$t$$ with itself such as $$ t^2=t\circ t $$ and $$ t^3=t\circ t\circ t $$ are defined. Note that this power notation for the linear transformation functions dovetails with the notation that we've used earlier for their squared matrix representations because if $${\rm Rep}_{B,B}(t)=T$$ then $$ {\rm Rep}_{B,B}(t^j)=T^j $$.

These examples suggest that on iteration more and more zeros appear until there is a settling down. The next result makes this precise.

This graph illustrates Lemma 1.3. The horizontal axis gives the power $$j$$ of a transformation. The vertical axis gives the dimension of the rangespace of $$t^j$$ as the distance above zero— and thus also shows the dimension of the nullspace as the distance below the gray horizontal line, because the two add to the dimension $$n$$ of the domain. As sketched, on iteration the rank falls and with it the nullity grows until the two reach a steady state. This state must be reached by the $$n$$-th iterate. The steady state's distance above zero is the dimension of the generalized rangespace and its distance below $$n$$ is the dimension of the generalized nullspace.

Exercises
/Solutions/

title: Linear Algebra/Self-Composition


text:
Linear Algebra/Strings

This subsection is optional, and requires material from the optional Direct Sum subsection.

The prior subsection shows that as $$ j $$ increases, the dimensions of the $$\mathcal{R}(t^j)$$'s fall while the dimensions of the $$\mathcal{N}(t^j)$$'s rise, in such a way that this rank and nullity split the dimension of $$V$$. Can we say more; do the two split a basis— is $$ V=\mathcal{R}(t^j)\oplus\mathcal{N}(t^j) $$?

The answer is yes for the smallest power $$j=0$$ since $$ V=\mathcal{R}(t^0)\oplus\mathcal{N}(t^0)=V\oplus\{\vec{0}\} $$. The answer is also yes at the other extreme.

In contrast to the $$j=0$$ and $$j=n$$ cases, for intermediate powers the space $$V$$ might not be the direct sum of $$\mathcal{R}(t^j)$$ and $$\mathcal{N}(t^j)$$. The next example shows that the two can have a nontrivial intersection.

In those three examples all vectors are eventually transformed to zero.

Not all nilpotent matrices are all zeros except for blocks of subdiagonal ones.

The goal of this subsection is Theorem 2.13, which shows that the prior example is prototypical in that every nilpotent matrix is similar to one that is all zeros except for blocks of subdiagonal ones.

We shall show that every nilpotent map has an associated string basis. Then our goal theorem, that every nilpotent matrix is similar to one that is all zeros except for blocks of subdiagonal ones, is immediate, as in Example 2.5.

Looking for a counterexample, a nilpotent map without an associated string basis that is disjoint, will suggest the idea for the proof. Consider the map $$ t:\mathbb{C}^5\to \mathbb{C}^5 $$ with this action. $${\rm Rep}_{\mathcal{E}_5,\mathcal{E}_5}(t)= \begin{pmatrix} 0 &0  &0  &0  &0  \\ 0  &0  &0  &0  &0  \\ 1  &1  &0  &0  &0  \\ 0  &0  &0  &0  &0  \\ 0  &0  &0  &1  &0 \end{pmatrix}$$ Even after omitting the zero vector, these three strings aren't disjoint, but that doesn't end hope of finding a $$t$$-string basis. It only means that $$ \mathcal{E}_5 $$ will not do for the string basis.

To find a basis that will do, we first find the number and lengths of its strings. Since $$t$$'s index of nilpotency is two, Lemma 2.12 says that at least one string in the basis has length two. Thus the map must act on a string basis in one of these two ways.

\begin{array}{ccccc} \vec{\beta}_1 &\mapsto &\vec{\beta}_2 &\mapsto &\vec{0} \\ \vec{\beta}_3 &\mapsto &\vec{\beta}_4 &\mapsto &\vec{0} \\ \vec{\beta}_5 &\mapsto &\vec{0} \end{array}$$
 * $$\begin{array}{ccccc}

\vec{\beta}_1 &\mapsto &\vec{\beta}_2 &\mapsto &\vec{0} \\ \vec{\beta}_3 &\mapsto &\vec{0}  \\ \vec{\beta}_4 &\mapsto &\vec{0}  \\ \vec{\beta}_5 &\mapsto &\vec{0} \end{array} $$

Now, the key point. A transformation with the left-hand action has a nullspace of dimension three since that's how many basis vectors are sent to zero. A transformation with the right-hand action has a nullspace of dimension four. Using the matrix representation above, calculation of $$t$$'s nullspace

\mathcal{N}(t)= \{\begin{pmatrix} x \\ -x \\ z \\ 0 \\ r \end{pmatrix}\,\big|\, x,z,r\in\mathbb{C} \} $$

shows that it is three-dimensional, meaning that we want the left-hand action.

To produce a string basis, first pick $$ \vec{\beta}_2 $$ and $$ \vec{\beta}_4 $$ from $$ \mathcal{R}(t)\cap\mathcal{N}(t) $$

\vec{\beta}_2=\begin{pmatrix} 0 \\ 0 \\ 1 \\ 0 \\ 0 \end{pmatrix}\qquad \vec{\beta}_4=\begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 1 \end{pmatrix} $$

(other choices are possible, just be sure that $$ \{\vec{\beta}_2,\vec{\beta}_4\} $$ is linearly independent). For $$ \vec{\beta}_5 $$ pick a vector from $$ \mathcal{N}(t) $$ that is not in the span of $$ \{ \vec{\beta}_2,\vec{\beta}_4 \} $$.

\vec{\beta}_5=\begin{pmatrix} 1 \\ -1 \\ 0 \\ 0 \\ 0 \end{pmatrix} $$

Finally, take $$ \vec{\beta}_1 $$ and $$ \vec{\beta}_3 $$ such that $$ t(\vec{\beta}_1)=\vec{\beta}_2 $$ and $$ t(\vec{\beta}_3)=\vec{\beta}_4 $$.

\vec{\beta}_1=\begin{pmatrix} 0 \\ 1 \\ 0 \\ 0 \\ 0 \end{pmatrix}\qquad \vec{\beta}_3=\begin{pmatrix} 0 \\ 0 \\ 0 \\ 1 \\ 0 \end{pmatrix} $$

Now, with respect to $$ B=\langle \vec{\beta}_1,\ldots,\vec{\beta}_5 \rangle $$, the matrix of $$t$$ is as desired.

{\rm Rep}_{B,B}(t)= \left(\begin{array}{cc|cc|c} 0 &0  &0  &0  &0  \\ 1  &0  &0  &0  &0  \\  \hline 0  &0  &0  &0  &0  \\ 0  &0  &1  &0  &0  \\  \hline 0  &0  &0  &0  &0 \end{array}\right) $$

This illustrates the proof. Basis vectors are categorized into kind $$1$$, kind $$2$$, and kind $$3$$. They are also shown as squares or circles, according to whether they are in the nullspace or not.

 * [[Image:Linalg nilpotent has string basis.png|x300px]]

This form is unique in the sense that if a nilpotent matrix is similar to two such matrices then those two simply have their blocks ordered differently. Thus this is a canonical form for the similarity classes of nilpotent matrices provided that we order the blocks, say, from longest to shortest.

Exercises
/Solutions/

title: Linear Algebra/Strings


text:
Linear Algebra/Jordan Form

This section uses material from three optional subsections: Direct Sum, Determinants Exist, and Other Formulas for the Determinant.

The chapter on linear maps shows that every $$h:V\to W$$ can be represented by a partial-identity matrix with respect to some bases $$B\subset V$$ and $$D\subset W$$. This chapter revisits this issue in the special case that the map is a linear transformation $$t:V\to V$$. Of course, the general result still applies but with the codomain and domain equal we naturally ask about having the two bases also be equal. That is, we want a canonical form to represent transformations as $${\rm Rep}_{B,B}(t)$$.

After a brief review section, we began by noting that a block partial identity form matrix is not always obtainable in this $$B,B$$ case. We therefore considered the natural generalization, diagonal matrices, and showed that if its eigenvalues are distinct then a map or matrix can be diagonalized. But we also gave an example of a matrix that cannot be diagonalized and in the section prior to this one we developed that example. We showed that a linear map is nilpotent— if we take higher and higher powers of the map or matrix then we eventually get the zero map or matrix— if and only if there is a basis on which it acts via disjoint strings. That led to a canonical form for nilpotent matrices.

Now, this section concludes the chapter. We will show that the two cases we've studied are exhaustive in that for any linear transformation there is a basis such that the matrix representation $${\rm Rep}_{B,B}(t)$$ is the sum of a diagonal matrix and a nilpotent matrix in its canonical form.

title: Linear Algebra/Jordan Form


text:
Linear Algebra/Polynomials of Maps and Matrices

Recall that the set of square matrices is a vector space under entry-by-entry addition and scalar multiplication and that this space $$ \mathcal{M}_{n \! \times \! n} $$ has dimension $$ n^2 $$. Thus, for any $$ n \! \times \! n $$ matrix $$T$$ the $$ n^2+1 $$-member set $$ \{I,T,T^2,\dots,T^{n^2} \} $$ is linearly dependent and so there are scalars $$ c_0,\dots,c_{n^2} $$ such that $$c_{n^2}T^{n^2}+\dots+c_1T+c_0I$$ is the zero matrix.

Of course, if $$ T={\rm Rep}_{B,B}(t) $$ then $$ f(T)={\rm Rep}_{B,B}(f(t)) $$, which follows from the relationships $$ T^j={\rm Rep}_{B,B}(t^j) $$, and $$ cT={\rm Rep}_{B,B}(ct) $$, and $$ T_1+T_2 ={\rm Rep}_{B,B}(t_1+t_2) $$.

As Example 1.2 shows, there may be polynomials of degree smaller than $$n^2$$ that zero the map or matrix.

A minimal polynomial always exists by the observation opening this subsection. A minimal polynomial is unique by the "with leading coefficient $$ 1 $$" clause. This is because if there are two polynomials $$ m(x) $$ and $$ \hat{m}(x) $$ that are both of the minimal degree to make the map or matrix zero (and thus are of equal degree), and both have leading $$ 1 $$'s, then their difference $$ m(x)-\hat{m}(x) $$ has a smaller degree than either and still sends the map or matrix to zero. Thus $$ m(x)-\hat{m}(x) $$ is the zero polynomial and the two are equal. (The leading coefficient requirement also prevents a minimal polynomial from being the zero polynomial.)

Using the method of that example to find the minimal polynomial of a $$ 3 \! \times \! 3 $$ matrix would mean doing Gaussian reduction on a system with nine equations in ten unknowns. We shall develop an alternative. To begin, note that we can break a polynomial of a map or a matrix into its components.

In particular, if a minimial polynomial $$m(x)$$ for a transformation $$t$$ factors as $$m(x)=(x-\lambda_1)^{q_1}\cdots (x-\lambda_\ell)^{q_\ell}$$ then $$ m(t)=(t-\lambda_1)^{q_1}\circ \cdots\circ (t-\lambda_\ell)^{q_\ell} $$ is the zero map. Since $$ m(t) $$ sends every vector to zero, at least one of the maps $$ t-\lambda_i $$ sends some nonzero vectors to zero. So, too, in the matrix case— if $$m$$ is minimal for $$T$$ then $$ m(T)=(T-\lambda_1I)^{q_1}\cdots (T-\lambda_\ell I)^{q_\ell} $$ is the zero matrix and at least one of the matrices $$T-\lambda_iI$$ sends some nonzero vectors to zero. Rewording both cases: at least some of the $$ \lambda_i $$ are eigenvalues. (See Problem 17.)

Recall how we have earlier found eigenvalues. We have looked for $$\lambda$$ such that $$T\vec{v}=\lambda\vec{v}$$ by considering the equation $$\vec{0}=T\vec{v}-x\vec{v}=(T-xI)\vec{v}$$ and computing the determinant of the matrix $$T-xI$$. That determinant is a polynomial in $$x$$, the characteristic polynomial, whose roots are the eigenvalues. The major result of this subsection, the next result, is that there is a connection between this characteristic polynomial and the minimal polynomial. This results expands on the prior paragraph's insight that some roots of the minimal polynomial are eigenvalues by asserting that every root of the minimal polynomial is an eigenvalue and further that every eigenvalue is a root of the minimal polynomial (this is because it says "$$1\leq q_i$$" and not just "$$0\leq q_i$$").

The proof takes up the next three lemmas. Although they are stated only in matrix terms, they apply equally well to maps. We give the matrix version only because it is convenient for the first proof.

The first result is the key— some authors call it the Cayley-Hamilton Theorem and call Theorem 1.8 above a corollary. For the proof, observe that a matrix of polynomials can be thought of as a polynomial with matrix coefficients.

\begin{pmatrix} 2x^2+3x-1 &x^2+2    \\ 3x^2+4x+1 &4x^2+x+1 \end{pmatrix} = \begin{pmatrix} 2 &1  \\ 3  &4 \end{pmatrix}x^2 + \begin{pmatrix} 3 &0  \\ 4  &1 \end{pmatrix}x + \begin{pmatrix} -1 &2  \\ 1  &1 \end{pmatrix} $$

We sometimes refer to that lemma by saying that a matrix or map satisfies its characteristic polynomial.

Combining the prior two lemmas gives that the minimal polynomial divides the characteristic polynomial. Thus, any root of the minimal polynomial is also a root of the characteristic polynomial. That is, so far we have that if $$ m(x)=(x-\lambda_1)^{q_1}\dots(x-\lambda_i)^{q_i} $$ then $$ c(x) $$ must has the form $$ (x-\lambda_1)^{p_1}\dots(x-\lambda_i)^{p_i} (x-\lambda_{i+1})^{p_{i+1}}\dots(x-\lambda_\ell)^{p_\ell} $$ where each $$ q_j $$ is less than or equal to $$ p_j $$. The proof of the Cayley-Hamilton Theorem is finished by showing that in fact the characteristic polynomial has no extra roots $$\lambda_{i+1}$$, etc.

Exercises
/Solutions/


text:
Linear Algebra/Jordan Canonical Form

This subsection moves from the canonical form for nilpotent matrices to the one for all matrices.

We have shown that if a map is nilpotent then all of its eigenvalues are zero. We can now prove the converse.

We have a canonical form for nilpotent matrices, that is, for each matrix whose single eigenvalue is zero: each such matrix is similar to one that is all zeroes except for blocks of subdiagonal ones. (To make this representation unique we can fix some arrangement of the blocks, say, from longest to shortest.) We next extend this to all single-eigenvalue matrices.

Observe that if $$ t $$'s only eigenvalue is $$ \lambda $$ then $$ t-\lambda $$'s only eigenvalue is $$ 0 $$ because $$ t(\vec{v})=\lambda\vec{v} $$ if and only if $$ (t-\lambda)\,(\vec{v})=0\cdot\vec{v} $$. The natural way to extend the results for nilpotent matrices is to represent $$t-\lambda$$ in the canonical form $$N$$, and try to use that to get a simple representation $$T$$ for $$t$$. The next result says that this try works.

An array that is all zeroes, except for some number $$\lambda$$ down the diagonal and blocks of subdiagonal ones, is a Jordan block. We have shown that Jordan block matrices are canonical representatives of the similarity classes of single-eigenvalue matrices.

We will now finish the program of this chapter by extending this work to cover maps and matrices with multiple eigenvalues. The best possibility for general maps and matrices would be if we could break them into a part involving their first eigenvalue $$ \lambda_1 $$ (which we represent using its Jordan block), a part with $$ \lambda_2 $$, etc.

This ideal is in fact what happens. For any transformation $$ t:V\to V $$, we shall break the space $$ V $$ into the direct sum of a part on which $$ t-\lambda_1 $$ is nilpotent, plus a part on which $$ t-\lambda_2 $$ is nilpotent, etc. More precisely, we shall take three steps to get to this section's major theorem and the third step shows that $$ V=\mathcal{N}_\infty(t-\lambda_1)\oplus\cdots\oplus \mathcal{N}_\infty(t-\lambda_\ell) $$ where $$ \lambda_1,\ldots,\lambda_\ell $$ are $$ t $$'s eigenvalues.

Suppose that $$ t:V\to V $$ is a linear transformation. Note that the restriction of $$ t $$ to a subspace $$ M $$ need not be a linear transformation on $$ M $$ because there may be an $$ \vec{m}\in M $$ with $$ t(\vec{m})\not\in M $$. To ensure that the restriction of a transformation to a "part" of a space is a transformation on the partwe need the next condition.

Two examples are that the generalized null space $$\mathcal{N}_\infty(t)$$ and the generalized range space $$\mathcal{R}_\infty(t)$$ of any transformation $$t$$ are invariant. For the generalized null space, if $$\vec{v}\in\mathcal{N}_\infty(t)$$ then $$t^n(\vec{v})=\vec{0}$$ where $$n$$ is the dimension of the underlying space and so $$t(\vec{v})\in\mathcal{N}_\infty(t)$$ because $$t^n(\,t(\vec{v})\,)$$ is zero also. For the generalized range space, if $$\vec{v}\in\mathcal{R}_\infty(t)$$ then $$\vec{v}=t^n(\vec{w})$$ for some $$\vec{w}$$ and then $$t(\vec{v})=t^{n+1}(\vec{w})=t^n(\,t(\vec{w})\,)$$ shows that $$t(\vec{v})$$ is also a member of $$\mathcal{R}_\infty(t)$$.

Thus the spaces $$\mathcal{N}_\infty(t-\lambda_i)$$ and $$\mathcal{R}_\infty(t-\lambda_i)$$ are $$t-\lambda_i$$ invariant. Observe also that $$t-\lambda_i$$ is nilpotent on $$\mathcal{N}_\infty(t-\lambda_i)$$ because, simply, if $$\vec{v}$$ has the property that some power of $$t-\lambda_i$$ maps it to zero— that is, if it is in the generalized null space— then some power of $$t-\lambda_i$$ maps it to zero. The generalized null space $$\mathcal{N}_\infty(t-\lambda_i)$$ is a "part" of the space on which the action of $$t-\lambda_i$$ is easy to understand.

The next result is the first of our three steps. It establishes that $$ t-\lambda_j $$ leaves $$ t-\lambda_i $$'s part unchanged.

The second step of the three that we will take to prove this section's major result makes use of an additional property of $$ \mathcal{N}_\infty(t-\lambda_i) $$ and $$ \mathcal{R}_\infty(t-\lambda_i) $$, that they are complementary. Recall that if a space is the direct sum of two others $$ V=\mathcal{N}\oplus \mathcal{R} $$ then any vector $$ \vec{v} $$ in the space breaks into two parts $$ \vec{v}=\vec{n}+\vec{r} $$ where $$ \vec{n}\in \mathcal{N} $$ and $$ \vec{r}\in \mathcal{R} $$, and recall also that if $$ B_{\mathcal{N}} $$ and $$ B_{\mathcal{R}} $$ are bases for $$ \mathcal{N} $$ and $$ \mathcal{R} $$ then the concatenation $$ B_{\mathcal{N}}\!\mathbin{{}^\frown}\!B_{\mathcal{R}} $$ is linearly independent (and so the two parts of $$ \vec{v} $$ do not "overlap"). The next result says that for any subspaces $$ \mathcal{N} $$ and $$ \mathcal{R} $$ that are complementary as well as $$ t $$ invariant, the action of $$ t $$ on $$ \vec{v} $$ breaks into the "non-overlapping" actions of $$ t $$ on $$ \vec{n} $$ and on $$ \vec{r} $$.

{{TextBox|1= Let $$ t:V\to V $$ be a transformation and let $$ \mathcal{N} $$ and $$ \mathcal{R} $$ be $$ t $$ invariant complementary subspaces of $$ V $$. Then $$ t $$ can be represented by a matrix with blocks of square submatrices $$T_1$$ and $$T_2$$
 * Lemma 2.8{{anchor|le:InvCompSubspSplitTrans}}:

\left(\begin{array}{c|c} T_1  &Z_2  \\  \hline Z_1 &T_2 \end{array}\right) \begin{array}{ll} \} \dim(\mathcal{N})\text{-many rows} \\ \} \dim(\mathcal{R})\text{-many rows} \end{array} $$

where $$ Z_1 $$ and $$ Z_2 $$ are blocks of zeroes. }}

To see that $$ t $$ has been decomposed into its action on the parts, observe that the restrictions of $$ t $$ to the subspaces $$ \mathcal{N} $$ and $$ \mathcal{R} $$ are represented, with respect to the obvious bases, by the matrices $$ T_1 $$ and $$ T_2 $$. So, with subspaces that are invariant and complementary, we can split the problem of examining a linear transformation into two lower-dimensional subproblems. The next result illustrates this decomposition into blocks.

From Lemma 2.9 we conclude that if two subspaces are complementary and $$ t $$ invariant then $$ t $$ is nonsingular if and only if its restrictions to both subspaces are nonsingular.

Now for the promised third, final, step to the main result.

{{TextBox|1= Because $$ \dim (V) $$ is the degree $$ p_1+\cdots+p_\ell $$ of the characteristic polynomial, to establish statement (1) we need only show that statement (2) holds and that $$ \mathcal{N}_\infty(t-\lambda_i)\cap\mathcal{N}_\infty(t-\lambda_j) $$ is trivial whenever $$ i\neq j $$.
 * Proof:

For the latter, by Lemma 2.7, both $$ \mathcal{N}_\infty(t-\lambda_i) $$ and $$ \mathcal{N}_\infty(t-\lambda_j) $$ are $$ t $$ invariant. Notice that an intersection of $$ t $$ invariant subspaces is $$ t $$ invariant and so the restriction of $$ t $$ to $$ \mathcal{N}_\infty(t-\lambda_i)\cap\mathcal{N}_\infty(t-\lambda_j) $$ is a linear transformation. But both $$ t-\lambda_i $$ and $$ t-\lambda_j $$ are nilpotent on this subspace and so if $$ t $$ has any eigenvalues on the intersection then its "only" eigenvalue is both $$ \lambda_i $$ and $$ \lambda_j $$. That cannot be, so this restriction has no eigenvalues: $$ \mathcal{N}_\infty(t-\lambda_i)\cap\mathcal{N}_\infty(t-\lambda_j) $$ is trivial (Lemma V.II.3.10 shows that the only transformation without any eigenvalues is on the trivial space).

To prove statement (2), fix the index $$ i $$. Decompose $$ V $$ as $$ \mathcal{N}_\infty(t-\lambda_i)\oplus\mathcal{R}_\infty(t-\lambda_i) $$

and apply Lemma 2.8.

T=

\left(\begin{array}{c|c} T_1  &Z_2  \\  \hline Z_1   &T_2 \end{array}\right) \begin{array}{ll} \} \dim(\,\mathcal{N}_\infty(t-\lambda_i)\,)\text{-many rows} \\ \} \dim(\,\mathcal{R}_\infty(t-\lambda_i)\,)\text{-many rows} \end{array} $$

By Lemma 2.9, $$ \left|T-xI\right|=\left|T_1-xI\right|\cdot\left|T_2-xI\right| $$. By the uniqueness clause of the Fundamental Theorem of Arithmetic, the determinants of the blocks have the same factors as the characteristic polynomial $$ \left|T_1-xI\right|=(x-\lambda_1)^{q_1}\dots(x-\lambda_\ell)^{q_\ell} $$ and $$ \left|T_2-xI\right|=(x-\lambda_1)^{r_1}\dots(x-\lambda_\ell)^{r_\ell} $$, and the sum of the powers of these factors is the power of the factor in the characteristic polynomial: $$ q_1+r_1=p_1 $$, ..., $$ q_\ell+r_\ell=p_\ell $$. Statement (2) will be proved if we will show that $$q_i=p_i$$ and that $$q_j=0$$ for all $$j\neq i$$, because then the degree of the polynomial $$\left|T_1-xI\right|$$— which equals the dimension of the generalized null space— is as required.

For that, first, as the restriction of $$ t-\lambda_i $$ to $$ \mathcal{N}_\infty(t-\lambda_i) $$ is nilpotent on that space, the only eigenvalue of $$ t $$ on it is $$ \lambda_i $$. Thus the characteristic equation of $$ t $$ on $$ \mathcal{N}_\infty(t-\lambda_i) $$ is $$ \left|T_1-xI\right|=(x-\lambda_i)^{q_i} $$. And thus $$q_j=0$$ for all $$j\neq i$$.

Now consider the restriction of $$ t $$ to $$ \mathcal{R}_\infty(t-\lambda_i) $$. By Note V.III.2.2, the map $$ t-\lambda_i $$ is nonsingular on $$ \mathcal{R}_\infty(t-\lambda_i) $$ and so $$ \lambda_i $$ is not an eigenvalue of $$ t $$ on that subspace. Therefore, $$ x-\lambda_i $$ is not a factor of $$ \left|T_2-xI\right| $$, and so $$ q_i=p_i $$. }}

Our major result just translates those steps into matrix terms.

Jordan form is a canonical form for similarity classes of square matrices, provided that we make it unique by arranging the Jordan blocks from least eigenvalue to greatest and then arranging the subdiagonal $$1$$ blocks inside each Jordan block from longest to shortest.

We close with the statement that the subjects considered earlier in this Chpater are indeed, in this sense, exhaustive.

Exercises
/Solutions/

title: Linear Algebra/Jordan Canonical Form


text:
Linear Algebra/Topic: Geometry of Eigenvalues

--Refer to Topic on Geometry of Linear Transformations---

The characterization of linear transformations in terms of the elementary operations is nice in some ways (for instance, we can easily see that lines are mapped to lines because each of the operations of projection, dilation, reflection, and skew maps lines to lines), but when a map is expressed as a composition of many small operations---no matter how simple---the description is less than ideal. We finish with another way, a somewhat more holistic way, of picturing the geometric effect of transformations of $$\mathbb{R}^2$$.

The pictures in that area give the action of the map on just one or two members of the domain. Although we know that a transformation is described completely by its action on a basis, and so to describe a transformation of $$\mathbb{R}^2$$ therefore, strictly speaking, requires only a description of where it sends the two vectors from any basis, those pictures seem not to convey much geometric intuition. Can we make clear a linear map's geometry by putting in more information, but not so much information that the picture gets confused?

A transformation of $$\mathbb{R}^2$$ sends lines through the origin to lines through the origin. Thus, two points on a line $$y=k_1x$$ will both be sent to the line, say, $$y=k_2x$$. Consider two such points. One is a multiple of the other, so we can write them with the second one as $$r$$ times the first, for some scalar $$r$$.

\begin{pmatrix} x \\ k_1x \end{pmatrix}\quad\mbox{and}\quad\begin{pmatrix} (rx) \\ k_1(rx) \end{pmatrix} $$

Compare their images.

\begin{pmatrix} a &c  \\ b &d \end{pmatrix} \begin{pmatrix} x \\ k_1x \end{pmatrix} = \begin{pmatrix} ax+ck_1x \\ bx+dk_1x \end{pmatrix} \qquad \begin{pmatrix} a &c  \\ b &d \end{pmatrix} \begin{pmatrix} (rx) \\ k_1(rx) \end{pmatrix} = \begin{pmatrix} a(rx)+ck_1(rx) \\ b(rx)+dk_1(rx) \end{pmatrix} $$

The second vector is $$r$$ times the first, and the image of the second is $$r$$ times the image of the first. Not only does the transformation preserve the fact that the vectors are colinear, it also preserves the relative scale of the vectors. That is, a transformation treats the points on a line through the origin uniformily. To describe the effect of the map on the entire line, we need only describe its effect on a single non-zero point in that line.

Since every point in the space is on some line through the origin, to understand the action of a linear transformation of $$\mathbb{R}^2$$, it is sufficient to pick one point from each line through the origin (say the point that is on the upper half of the unit circle) and show how the map's effect on that set of points.

Here is such a picture for a straightforward dilation. Below, the same map is shown with the circle and its image superimposed. Certainly the geometry here is more evident. For example, we can see that some lines through the origin are actually sent to themselves: the $$x$$-axis is sent to the $$x$$-axis, and the $$y$$-axis is sent to the $$y$$-axis.

This is the flip shown earlier, here with the circle and its image superimposed. And this is the skew shown earlier. Contrast the picture of this map's effect on the unit square with this one.

Here is a somewhat more complicated map (the second coordinate function is the same as the map in the prior picture, but the first coordinate function is different). Observe that some vectors are being both dilated and rotated through some angle

\begin{pmatrix} x \\ 2x \end{pmatrix}\mapsto\begin{pmatrix} x \\ k_1x \end{pmatrix} $$

while others are just being dilated, not rotated at all.

\begin{pmatrix} x \\ 3x \end{pmatrix}\mapsto\begin{pmatrix} x  \\ 3x \end{pmatrix} $$

Exercises
/Solutions/


text:
Linear Algebra/Topic: Stable Populations

Imagine a reserve park with animals from a species that we are trying to protect. The park doesn't have a fence and so animals cross the boundary, both from the inside out and in the other direction. Every year, 10% of the animals from inside of the park leave, and 1% of the animals from the outside find their way in. We can ask if we can find a stable level of population for this park: is there a population that, once established, will stay constant over time, with the number of animals leaving equal to the number of animals entering?

To answer that question, we must first establish the equations. Let the year $$n$$ population in the park be $$p_n$$ and in the rest of the world be $$r_n$$.

 * $$\begin{array}{rl}

p_{n+1} &=.90p_n+.01r_n   \\ r_{n+1} &=.10p_n+.99r_n \end{array}$$

We can set this system up as a matrix equation (see the Markov Chain topic).

\begin{pmatrix} p_{n+1} \\ r_{n+1} \end{pmatrix} =\begin{pmatrix} .90 &.01  \\ .10  &.99 \end{pmatrix} \begin{pmatrix} p_{n} \\ r_{n} \end{pmatrix} $$

Now, "stable level" means that $$p_{n+1}=p_n$$ and $$r_{n+1}=r_n$$, so that the matrix equation $$\vec{v}_{n+1}=T\vec{v}_{n}$$ becomes $$\vec{v}=T\vec{v}$$. We are therefore looking for eigenvectors for $$T$$ that are associated with the eigenvalue $$1$$. The equation $$(I-T)\vec{v}=\vec{0}$$ is

\begin{pmatrix} .10 &-.01  \\ -.10  &.01 \end{pmatrix} \begin{pmatrix} p \\ r \end{pmatrix}= \begin{pmatrix} 0 \\ 0 \end{pmatrix} $$

which gives the eigenspace: vectors with the restriction that $$p=.1r$$. Coupled with additional information, that the total world population of this species is is $$p+r=110\,000$$, we find that the stable state is $$p=10,000$$ and $$r=100,000$$.

If we start with a park population of ten thousand animals, so that the rest of the world has one hundred thousand, then every year ten percent (a thousand animals) of those inside will leave the park, and every year one percent (a thousand) of those from the rest of the world will enter the park. It is stable, self-sustaining.

Now imagine that we are trying to gradually build up the total world population of this species. We can try, for instance, to have the world population grow at a rate of 1% per year. In this case, we can take a "stable" state for the park's population to be that it also grows at 1% per year. The equation $$\vec{v}_{n+1}=1.01\cdot\vec{v}_n=T\vec{v}_{n}$$ leads to $$((1.01\cdot I)-T)\vec{v}=\vec{0}$$, which gives this system.

\begin{pmatrix} .11 &-.01  \\ -.10  &.02 \end{pmatrix} \begin{pmatrix} p \\ r \end{pmatrix}= \begin{pmatrix} 0 \\ 0 \end{pmatrix} $$

The matrix is nonsingular, and so the only solution is $$p=0$$ and $$r=0$$. Thus, there is no (usable) initial population that we can establish at the park and expect that it will grow at the same rate as the rest of the world.

Knowing that an annual world population growth rate of 1% forces an unstable park population, we can ask which growth rates there are that would allow an initial population for the park that will be self-sustaining. We consider $$\lambda\vec{v}=T\vec{v}$$ and solve for $$\lambda$$.

0=\begin{vmatrix} \lambda-.9 &-.01  \\ -.10        &\lambda-.99 \end{vmatrix} =(\lambda-.9)(\lambda-.99)-(.10)(.01) =\lambda^2-1.89\lambda+.89 $$

A shortcut to factoring that quadratic is our knowledge that $$\lambda=1$$ is an eigenvalue of $$T$$, so the other eigenvalue is $$.89$$. Thus there are two ways to have a stable park population (a population that grows at the same rate as the population of the rest of the world, despite the leaky park boundaries): have a world population that is does not grow or shrink, and have a world population that shrinks by 11% every year.

So this is one meaning of eigenvalues and eigenvectors— they give a stable state for a system. If the eigenvalue is $$1$$ then the system is static. If the eigenvalue isn't $$1$$ then the system is either growing or shrinking, but in a dynamically-stable way.

Exercises
/Solutions/


text:
Linear Algebra/Topic: The Method of Powers

In practice, calculating eigenvalues and eigenvectors is a difficult problem. Finding, and solving, the characteristic polynomial of the large matrices often encountered in applications is too slow and too hard. Other techniques, indirect ones that avoid the characteristic polynomial, are used. Here we shall see such a method that is suitable for large matrices that are "sparse" (the great majority of the entries are zero).

Suppose that the $$n \! \times \! n$$ matrix $$T$$ has the $$n$$ distinct eigenvalues $$\lambda_1$$, $$\lambda_2$$, ..., $$\lambda_n$$. Then $$\mathbb{R}^n$$ has a basis that is composed of the associated eigenvectors $$\langle \vec{\zeta}_1,\dots,\vec{\zeta}_n \rangle $$. For any $$\vec{v}\in\mathbb{R}^n$$, where $$\vec{v}=c_1\vec{\zeta}_1+\dots+c_n\vec{\zeta}_n$$, iterating $$T$$ on $$\vec{v}$$ gives these.

 * $$\begin{array}{rl}

T\vec{v} &=c_1\lambda_1\vec{\zeta}_1+c_2\lambda_2\vec{\zeta}_2+ \dots+c_n\lambda_n\vec{\zeta}_n \\ T^2\vec{v} &=c_1\lambda_1^2\vec{\zeta}_1+c_2\lambda_2^2\vec{\zeta}_2+ \dots+c_n\lambda_n^2\vec{\zeta}_n \\ T^3\vec{v} &=c_1\lambda_1^3\vec{\zeta}_1+c_2\lambda_2^3\vec{\zeta}_2+ \dots+c_n\lambda_n^3\vec{\zeta}_n \\ &\vdots                                           \\ T^k\vec{v} &=c_1\lambda_1^k\vec{\zeta}_1+c_2\lambda_2^k\vec{\zeta}_2+ \dots+c_n\lambda_n^k\vec{\zeta}_n \end{array}$$

If one of the eigenvalues, say, $$\lambda_1$$, has a larger absolute value than any of the other eigenvalues then its term will dominate the above expression. Put another way, dividing through by $$\lambda_1^k$$ gives this,

\frac{T^k\vec{v}}{\lambda_1^k} =c_1\vec{\zeta}_1+c_2\frac{\lambda_2^k}{\lambda_1^k}\vec{\zeta}_2+ \dots+c_n\frac{\lambda_n^k}{\lambda_1^k}\vec{\zeta}_n $$

and, because $$\lambda_1$$ is assumed to have the largest absolute value, as $$k$$ gets larger the fractions go to zero. Thus, the entire expression goes to $$c_1\vec{\zeta}_1$$.

That is (as long as $$c_1$$ is not zero), as $$k$$ increases, the vectors $$T^k\vec{v}$$ will tend toward the direction of the eigenvectors associated with the dominant eigenvalue, and, consequently, the ratios of the lengths $$|\,T^{k}\vec{v}\,|/|\,T^{k-1}\vec{v}\,|$$ will tend toward that dominant eigenvalue.

For example, (sample computer code for this follows the exercises), because the matrix

T=\begin{pmatrix} 3 &0  \\ 8  &-1 \end{pmatrix} $$

is triangular, its eigenvalues are just the entries on the diagonal, $$3$$ and $$-1$$. Arbitrarily taking $$\vec{v}$$ to have the components $$1$$ and $$1$$ gives $$ \begin{array}{c|ccccc} \vec{v} &T\vec{v}  &T^2\vec{v} &\cdots &T^9\vec{v} &T^{10}\vec{v}       \\ \hline \begin{pmatrix} 1 \\ 1 \end{pmatrix} &\begin{pmatrix} 3 \\ 7 \end{pmatrix} &\begin{pmatrix} 9 \\ 17 \end{pmatrix} &\cdots &\begin{pmatrix} 19\,683 \\ 39\,367 \end{pmatrix} &\begin{pmatrix} 59\,049 \\ 118\,097 \end{pmatrix} \end{array} $$ and the ratio between the lengths of the last two is $$2.999\,9$$.

Two implementation issues must be addressed. The first issue is that, instead of finding the powers of $$T$$ and applying them to $$\vec{v}$$, we will compute $$\vec{v}_1$$ as $$T\vec{v}$$ and then compute $$\vec{v}_2$$ as $$T\vec{v}_1$$, etc. (i.e., we never separately calculate $$T^2$$, $$T^3$$, etc.). These matrix-vector products can be done quickly even if $$T$$ is large, provided that it is sparse. The second issue is that, to avoid generating numbers that are so large that they overflow our computer's capability, we can normalize the $$\vec{v}_i$$'s at each step. For instance, we can divide each $$\vec{v}_i$$ by its length (other possibilities are to divide it by its largest component, or simply by its first component). We thus implement this method by generating

 * $$\begin{array}{rl}

\vec{w}_0 &=\vec{v}_0/|\vec{v}_0| \\ \vec{v}_1 &=T\vec{w}_0                 \\ \vec{w}_1 &=\vec{v}_1/|\vec{v}_1| \\ \vec{v}_2 &=T\vec{w}_2                 \\ &\vdots   \\ \vec{w}_{k-1} &=\vec{v}_{k-1}/|\vec{v}_{k-1}| \\ \vec{v}_k &=T\vec{w}_k \end{array}$$

until we are satisfied. Then the vector $$\vec{v}_k$$ is an approximation of an eigenvector, and the approximation of the dominant eigenvalue is the ratio $$|\vec{v}_{k}|/|\vec{w}_{k-1}|=|\vec{v}_k|$$.

One way we could be "satisfied" is to iterate until our approximation of the eigenvalue settles down. We could decide, for instance, to stop the iteration process not after some fixed number of steps, but instead when $$|\vec{v}_k|$$ differs from $$|\vec{v}_{k-1}|$$ by less than one percent, or when they agree up to the second significant digit.

The rate of convergence is determined by the rate at which the powers of $$|\lambda_2/\lambda_1|$$ go to zero, where $$\lambda_2$$ is the eigenvalue of second largest norm. If that ratio is much less than one then convergence is fast, but if it is only slightly less than one then convergence can be quite slow. Consequently, the method of powers is not the most commonly used way of finding eigenvalues (although it is the simplest one, which is why it is here as the illustration of the possibility of computing eigenvalues without solving the characteristic polynomial). Instead, there are a variety of methods that generally work by first replacing the given matrix $$T$$ with another that is similar to it and so has the same eigenvalues, but is in some reduced form such as tridiagonal form: the only nonzero entries are on the diagonal, or just above or below it. Then special techniques can be used to find the eigenvalues. Once the eigenvalues are known, the eigenvectors of $$T$$ can be easily computed. These other methods are outside of our scope. A good reference is.

Exercises
/Solutions/

This is the code for the computer algebra system Octave that was used to do the calculation above. (It has been lightly edited to remove blank lines, etc.)

Computer Code  >T=[3, 0; 8, -1] T= 3  0 8 -1 >v0=[1; 2] v0= 1 1 >v1=T*v0 v1= 3 7 >v2=T*v1 v2= 9 17 >T9=T**9 T9= 19683 0 39368 -1 >T10=T**10 T10= 59049 0 118096  1 >v9=T9*v0 v9= 19683 39367 >v10=T10*v0 v10= 59049 118096 >norm(v10)/norm(v9) ans=2.9999

Remark: we are ignoring the power of Octave here; there are built-in functions to automatically apply quite sophisticated methods to find eigenvalues and eigenvectors. Instead, we are using just the system as a calculator.


text:
Linear Algebra/Topic: Linear Recurrences

In 1202 Leonardo of Pisa, also known as Fibonacci, posed this problem. A certain man put a pair of rabbits in a place surrounded on all sides by a wall. How many pairs of rabbits can be produced from that pair in a year if it is supposed that every month each pair begets a new pair which from the second month on becomes productive? This moves past an elementary exponential growth model for population increase to include the fact that there is an initial period where newborns are not fertile. However, it retains other simplyfing assumptions, such as that there is no gestation period and no mortality.

The number of newborn pairs that will appear in the upcoming month is simply the number of pairs that were alive last month, since those will all be fertile, having been alive for two months. The number of pairs alive next month is the sum of the number alive current month and the number of newborns.

f(n+1)=f(n)+f(n-1) \qquad \text{where }f(0)=1\text{, }f(1)=1 $$

The is an example of a recurrence relation (it is called that because the values of $$f$$ are calculated by looking at other, prior, values of $$f$$). From it, we can easily answer Fibonacci's twelve-month question. $\begin{array}{r

The sequence of numbers defined by the above equation (of which the first few are listed) is the Fibonacci sequence. The material of this chapter can be used to give a formula with which we can can calculate $$f(n+1)$$ without having to first find $$f(n)$$, $$f(n-1)$$, etc.

For that, observe that the recurrence is a linear relationship and so we can give a suitable matrix formulation of it.

\begin{pmatrix} 1 &1   \\ 1  &0 \end{pmatrix} \begin{pmatrix} f(n) \\ f(n-1) \end{pmatrix} = \begin{pmatrix} f(n+1) \\ f(n) \end{pmatrix} \qquad\text{where }\begin{pmatrix} f(1) \\ f(0) \end{pmatrix}=\begin{pmatrix} 1 \\ 1 \end{pmatrix} $$

Then, where we write $$T$$ for the matrix and $$\vec{v}_{n}$$ for the vector with components $$f(n+1)$$ and $$f(n)$$, we have that $$\vec{v}_n=T^n\vec{v}_0$$. The advantage of this matrix formulation is that by diagonalizing $$T$$ we get a fast way to compute its powers: where $$T=PDP^{-1}$$ we have $$T^n=PD^nP^{-1}$$, and the $$n$$-th power of the diagonal matrix $$D$$ is the diagonal matrix whose entries that are the $$n$$-th powers of the entries of $$D$$.

The characteristic equation of $$T$$ is $$\lambda^2-\lambda-1$$. The quadratic formula gives its roots as $$(1+\sqrt{5})/2$$ and $$(1-\sqrt{5})/2$$. Diagonalizing gives this.

\begin{pmatrix} 1 &1  \\ 1  &0 \end{pmatrix} =\begin{pmatrix} \frac{1+\sqrt{5}}{2} &\frac{1-\sqrt{5}}{2} \\ 1                    &1 \end{pmatrix} \begin{pmatrix} \frac{1+\sqrt{5}}{2} &0   \\ 0                    &\frac{1-\sqrt{5}}{2} \end{pmatrix} \begin{pmatrix} \frac{1}{\sqrt{5}}    &-\frac{1-\sqrt{5}}{2\sqrt{5}}  \\ \frac{-1}{\sqrt{5}}   &\frac{1+\sqrt{5}}{2\sqrt{5}} \end{pmatrix} $$

Introducing the vectors and taking the $$n$$-th power, we have

 * $$\begin{array}{rl}

\begin{pmatrix} f(n+1) \\ f(n) \end{pmatrix} &=\begin{pmatrix} 1 &1  \\ 1  &0 \end{pmatrix}^n \begin{pmatrix} f(1) \\ f(0) \end{pmatrix}                                         \\ &=\begin{pmatrix} \frac{1+\sqrt{5}}{2} &\frac{1-\sqrt{5}}{2} \\ 1                    &1 \end{pmatrix} \begin{pmatrix} \frac{1+\sqrt{5}}{2}^n &0   \\ 0                      &\frac{1-\sqrt{5}}{2}^n \end{pmatrix} \begin{pmatrix} \frac{1}{\sqrt{5}}    &-\frac{1-\sqrt{5}}{2\sqrt{5}}  \\ \frac{-1}{\sqrt{5}}   &\frac{1+\sqrt{5}}{2\sqrt{5}} \end{pmatrix} \begin{pmatrix} f(1) \\ f(0) \end{pmatrix} \end{array}$$

We can compute $$f(n)$$ from the second component of that equation.

f(n)=\frac{1}{\sqrt{5}}\left[\left(\frac{1+\sqrt{5}}{2}\right)^{n+1} -\left(\frac{1-\sqrt{5}}{2}\right)^{n+1}\right] $$

Notice that $$f$$ is dominated by its first term because $$(1-\sqrt{5})/2$$ is less than one, so its powers go to zero.

In general, a linear recurrence relation has the form

f(n+1)=a_nf(n)+a_{n-1}f(n-1)+\dots+a_{n-k}f(n-k) $$

(it is also called a difference equation). This recurrence relation is homogeneous because there is no constant term; i.e, it can be put into the form $$0=-f(n+1)+a_nf(n)+a_{n-1}f(n-1)+\dots+a_{n-k}f(n-k)$$. This is said to be a relation of order $$k$$. The relation, along with the initial conditions $$f(0)$$, ..., $$f(k)$$ completely determine a sequence. For instance, the Fibonacci relation is of order $$2$$ and it, along with the two initial conditions $$f(0)=1$$ and $$f(1)=1$$, determines the Fibonacci sequence simply because we can compute any $$f(n)$$ by first computing $$f(2)$$, $$f(3)$$, etc. In this Topic, we shall see how linear algebra can be used to solve linear recurrence relations.

First, we define the vector space in which we are working. Let $$V$$ be the set of functions $$f$$ from the natural numbers $$\mathbb{N}=\{0,1,2,\ldots\}$$ to the real numbers. (Below we shall have functions with domain $$\{1,2,\ldots\}$$, that is, without $$0$$, but it is not an important distinction.)

Putting the initial conditions aside for a moment, for any recurrence, we can consider the subset $$S$$ of $$V$$ of solutions. For example, without initial conditions, in addition to the function $$f$$ given above, the Fibonacci relation is also solved by the function $$g$$ whose first few values are $$g(0)=1$$, $$g(1)=3$$, $$g(2)=4$$, and $$g(3)=7$$.

The subset $$S$$ is a subspace of $$V$$. It is nonempty because the zero function is a solution. It is closed under addition since if $$f_1$$ and $$f_2$$ are solutions, then

a_{n+1}(f_1+f_2)(n+1)+\dots+a_{n-k}(f_1+f_2)(n-k)$$

\begin{align} &=(a_{n+1}f_1(n+1)+\dots+a_{n-k}f_1(n-k))         \\ &\quad+(a_{n+1}f_2(n+1)+\dots+a_{n-k}f_2(n-k))    \\ &=0. \end{align} $$

And, it is closed under scalar multiplication since

a_{n+1}(rf_1)(n+1)+\dots+a_{n-k}(rf_1)(n-k)$$

\begin{align} &=r(a_{n+1}f_1(n+1)+\dots+a_{n-k}f_1(n-k))  \\ &=r\cdot 0                                   \\ &=0. \end{align} $$

We can give the dimension of $$S$$. Consider this map from the set of functions $$S$$ to the set of vectors $$\mathbb{R}^k$$.

f \mapsto \begin{pmatrix} f(0) \\ f(1) \\ \vdots \\ f(k) \end{pmatrix} $$

Problem 3 shows that this map is linear. Because, as noted above, any solution of the recurrence is uniquely determined by the $$k$$ initial conditions, this map is one-to-one and onto. Thus it is an isomorphism, and thus $$S$$ has dimension $$k$$, the order of the recurrence.

So (again, without any initial conditions), we can describe the set of solutions of any linear homogeneous recurrence relation of degree $$k$$ by taking linear combinations of only $$k$$ linearly independent functions. It remains to produce those functions.

For that, we express the recurrence $$f(n+1)=a_nf(n)+\dots+a_{n-k}f(n-k)$$ with a matrix equation.

\begin{pmatrix} a_n &a_{n-1}  &a_{n-2}  &\ldots  &a_{n-k+1} &a_{n-k}  \\ 1   &0        &0        &\ldots  &0         &0        \\ 0   &1        &0                                      \\ 0    &0        &1                                      \\ \vdots &\vdots &         &\ddots  &           &\vdots  \\ 0   &0        &0        &\ldots   &1          &0 \end{pmatrix} \begin{pmatrix} f(n) \\ f(n-1) \\ \vdots \\ f(n-k) \end{pmatrix} = \begin{pmatrix} f(n+1) \\ f(n) \\ \vdots \\ f(n-k+1) \end{pmatrix} $$

In trying to find the characteristic function of the matrix, we can see the pattern in the $$2 \! \times \! 2$$ case

\begin{pmatrix} a_n-\lambda &a_{n-1} \\ 1           &-\lambda \end{pmatrix} =\lambda^2-a_n\lambda-a_{n-1} $$

and $$3 \! \times \! 3$$ case.

\begin{pmatrix} a_n-\lambda &a_{n-1}   &a_{n-2}  \\ 1           &-\lambda  &0        \\ 0           &1         &-\lambda \end{pmatrix} =-\lambda^3+a_n\lambda^2+a_{n-1}\lambda+a_{n-2} $$

Problem 4 shows that the characteristic equation is this.

\begin{vmatrix} a_n-\lambda &a_{n-1} &a_{n-2}  &\ldots  &a_{n-k+1} &a_{n-k}  \\ 1   &-\lambda &0        &\ldots  &0         &0        \\ 0   &1        &-\lambda                                      \\ 0   &0        &1                                      \\ \vdots &\vdots &         &\ddots   &           &\vdots  \\ 0   &0        &0        &\ldots   &1          &-\lambda \end{vmatrix}$$

=\pm(-\lambda^k+a_n\lambda^{k-1}+a_{n-1}\lambda^{k-2} +\dots+a_{n-k+1}\lambda+a_{n-k}) $$

We call that the polynomial "associated" with the recurrence relation. (We will be finding the roots of this polynomial and so we can drop the $$\pm$$ as irrelevant.)

If $$-\lambda^k+a_n\lambda^{k-1}+a_{n-1}\lambda^{k-2} +\dots+a_{n-k+1}\lambda+a_{n-k}$$ has no repeated roots then the matrix is diagonalizable and we can, in theory, get a formula for $$f(n)$$ as in the Fibonacci case. But, because we know that the subspace of solutions has dimension $$k$$, we do not need to do the diagonalization calculation, provided that we can exhibit $$k$$ linearly independent functions satisfying the relation.

Where $$r_1$$, $$r_2$$, ..., $$r_k$$ are the distinct roots, consider the functions $$f_{r_1}(n)=r_1^n$$ through $$f_{r_k}(n)=r_k^n$$ of powers of those roots. Problem 2 shows that each is a solution of the recurrence and that the $$k$$ of them form a linearly independent set. So, given the homogeneous linear recurrence $$f(n+1)=a_nf(n)+\dots+a_{n-k}f(n-k)$$ (that is, $$0=-f(n+1)+a_nf(n)+\dots+a_{n-k}f(n-k)$$) we consider the associated equation $$0=-\lambda^k+a_n\lambda^{k-1}+\dots+a_{n-k+1}\lambda+a_{n-k}$$. We find its roots $$r_1$$, ..., $$r_k$$, and if those roots are distinct then any solution of the relation has the form $$f(n)=c_1r_1^n+c_2r_2^n+\dots+c_kr_k^n$$ for $$c_1, \dots, c_n\in\mathbb{R}$$. (The case of repeated roots is also easily done, but we won't cover it here— see any text on Discrete Mathematics.)

Now, given some initial conditions, so that we are interested in a particular solution, we can solve for $$c_1$$, ..., $$c_n$$. For instance, the polynomial associated with the Fibonacci relation is $$-\lambda^2+\lambda+1$$, whose roots are $$(1\pm\sqrt{5})/2$$ and so any solution of the Fibonacci equation has the form $$f(n)=c_1((1+\sqrt{5})/2)^n+c_2((1-\sqrt{5})/2)^n$$. Including the initial conditions for the cases $$n=0$$ and $$n=1$$ gives

\begin{array}{*{2}{rc}r} c_1              &+  &c_2               &=  &1  \\ (1+\sqrt{5}/2)c_1 &+ &(1-\sqrt{5}/2)c_2 &=  &1 \end{array} $$

which yields $$c_1=1/\sqrt{5}$$ and $$c_2=-1/\sqrt{5}$$, as was calculated above.

We close by considering the nonhomogeneous case, where the relation has the form $$f(n+1)=a_nf(n)+a_{n-1}f(n-1)+\dots+a_{n-k}f(n-k)+b$$ for some nonzero $$b$$. As in the first chapter of this book, only a small adjustment is needed to make the transition from the homogeneous case. This classic example illustrates.

In 1883, Edouard Lucas posed the following problem. In the great temple at Benares, beneath the dome which marks the center of the world, rests a brass plate in which are fixed three diamond needles, each a cubit high and as thick as the body of a bee. On one of these needles, at the creation, God placed sixty four disks of pure gold, the largest disk resting on the brass plate, and the others getting smaller and smaller up to the top one. This is the Tower of Bramah. Day and night unceasingly the priests transfer the disks from one diamond needle to another according to the fixed and immutable laws of Bramah, which require that the priest on duty must not move more than one disk at a time and that he must place this disk on a needle so that there is no smaller disk below it. When the sixty-four disks shall have been thus transferred from the needle on which at the creation God placed them to one of the other needles, tower, temple, and Brahmins alike will crumble into dusk, and with a thunderclap the world will vanish. (Translation of from .)

How many disk moves will it take? Instead of tackling the sixty four disk problem right away, we will consider the problem for smaller numbers of disks, starting with three.

To begin, all three disks are on the same needle.

After moving the small disk to the far needle, the mid-sized disk to the middle needle, and then moving the small disk to the middle needle we have this.

Now we can move the big disk over. Then, to finish, we repeat the process of moving the smaller disks, this time so that they end up on the third needle, on top of the big disk.

So the thing to see is that to move the very largest disk, the bottom disk, at a minimum we must: first move the smaller disks to the middle needle, then move the big one, and then move all the smaller ones from the middle needle to the ending needle. Those three steps give us this recurrence.

T(n+1)=T(n)+1+T(n)=2T(n)+1 \quad \text{where }T(1)=1 $$

We can easily get the first few values of $$T$$. $\begin{array}{r

We recognize those as being simply one less than a power of two.

To derive this equation instead of just guessing at it, we write the original relation as $$-1=-T(n+1)+2T(n)$$, consider the homogeneous relation $$0=-T(n)+2T(n-1)$$, get its associated polynomial $$-\lambda+2$$, which obviously has the single, unique, root of $$r_1=2$$, and conclude that functions satisfying the homogeneous relation take the form $$T(n)=c_12^n$$.

That's the homogeneous solution. Now we need a particular solution.

Because the nonhomogeneous relation $$-1=-T(n+1)+2T(n)$$ is so simple, in a few minutes (or by remembering the table) we can spot the particular solution $$T(n)=-1$$ (there are other particular solutions, but this one is easily spotted). So we have that— without yet considering the initial condition— any solution of $$T(n+1)=2T(n)+1$$ is the sum of the homogeneous solution and this particular solution: $$T(n)=c_12^n-1$$.

The initial condition $$T(1)=1$$ now gives that $$c_1=1$$, and we've gotten the formula that generates the table: the $$n$$-disk Tower of Hanoi problem requires a minimum of $$2^n-1$$ moves.

Finding a particular solution in more complicated cases is, naturally, more complicated. A delightful and rewarding, but challenging, source on recurrence relations is ., For more on the Tower of Hanoi, or are good starting points. So is. Some computer code for trying some recurrence relations follows the exercises.

Exercises
Computer Code This code allows the generation of the first few values of a function defined by a recurrence and initial conditions. It is in the Scheme dialect of LISP (specifically, it was written for A. Jaffer's free scheme interpreter SCM, although it should run in any Scheme implementation).

First, the Tower of Hanoi code is a straightforward implementation of the recurrence.

(define (tower-of-hanoi-moves n) (if (= n 1) 1 (+ (* (tower-of-hanoi-moves (- n 1)) 2) 1) ) )

(Note for readers unused to recursive code: to compute $$T(64)$$, the computer is told to compute $$2*T(63)-1$$, which requires, of course, computing $$T(63)$$. The computer puts the "times $$2$$" and the "plus $$1$$" aside for a moment to do that. It computes $$T(63)$$ by using this same piece of code (that's what "recursive" means), and to do that is told to compute $$2*T(62)-1$$. This keeps up (the next step is to try to do $$T(62)$$ while the other arithmetic is held in waiting), until, after $$63$$ steps, the computer tries to compute $$T(1)$$. It then returns $$T(1)=1$$, which now means that the computation of $$T(2)$$ can proceed, etc., up until the original computation of $$T(64)$$ finishes.)

The next routine calculates a table of the first few values. (Some language notes:  is the empty list, that is, the empty sequence, and   pushes something onto the start of a list. Note that, in the last line, the procedure  is called on argument  .)

(define (first-few-outputs proc n) (first-few-outputs-helper proc n ') ) ; (define (first-few-outputs-aux proc n lst) (if (< n 1) lst (first-few-outputs-aux proc (- n 1) (cons (proc n) lst)) ) )

The session at the SCM prompt went like this.

>(first-few-outputs tower-of-hanoi-moves 64) Evaluation took 120 mSec (1 3 7 15 31 63 127 255 511 1023 2047 4095 8191 16383 32767 65535 131071 262143 524287 1048575 2097151 4194303 8388607 16777215 33554431 67108863 134217727 268435455 536870911 1073741823 2147483647 4294967295 8589934591 17179869183 34359738367 68719476735 137438953471 274877906943 549755813887 1099511627775 2199023255551 4398046511103 8796093022207 17592186044415 35184372088831 70368744177663 140737488355327 281474976710655 562949953421311 1125899906842623 2251799813685247 4503599627370495 9007199254740991 18014398509481983 36028797018963967 72057594037927935 144115188075855871 288230376151711743 576460752303423487 1152921504606846975 2305843009213693951 4611686018427387903 9223372036854775807 18446744073709551615)

This is a list of $$T(1)$$ through $$T(64)$$. (The $$120$$ mSec came on a 50 mHz '486 running in an XTerm of XWindow under Linux. The session was edited to put line breaks between numbers.) /Solutions/


text:
User talk:MasterTama

Welcome message
Hello, and welcome to Wikibooks!

Here are some tips to help you get started:
 * To learn the wiki-markup-language syntax, see Help:Editing.
 * Make sure to sign your posts and comments with four tildes, like this: ~~~~
 * There is a box at the top of the edit window (if javascript is enabled on your browser) that will insert it too (looks like part of a signature). This will let others know who left it, and make it easy to reply back to you.
 * Remember to conduct any editing experiments in the sandbox.
 * You can tell the community something about yourself in your userpage.
 * You can get to this page by clicking the tab at the top of the page labeled with your registered username.
 * Wikibooks is a collection of open-source textbooks. (What is Wikibooks).
 * If you are a Wikipedian, see Wikibooks for Wikipedians for a primer on how things work here (it's a little different).
 * If you want to base your work here on materials from Wikipedia, please use WB:RFI (requests for importation) to bring the material over in compliance with the GFDL.
 * If you're an instructor and plan on using Wikibooks for a class project, see Guidelines for class projects
 * Please say hello at the Reading Room with any questions or ideas.
 * Eventually, you might want to read the Manual of Style and Policies and Guidelines.
 * Help us by participating in policy and guideline creation.
 * Please take a look at Naming policy before starting a new book.
 * Remember to maintain a Neutral point of view.
 * Explore, be bold in editing pages, and have fun!

You will find more resources in Community Portal. If you want to ask a question, visit the Study help desk, the Reading Room, IRC channel or ask me personally on my talk page. For site news, see the Bulletin board. It might be a good idea to add that page to your "watchlist" so that you can see when any new information is posted there. You can do that by clicking the tab labeled "watch" at the top of the page.

Good luck! --Panic (talk) 01:30, 12 March 2008 (UTC)

title: User talk:MasterTama