List of Figures
1.1 We think that gold is buried under the sand so we make measurements of gravity at various locations on the surface. . . . . . . . . . . . . . . . 1.2 Inverse problems usually start with some procedure for predicting the response of a physical system with known parameters. Then we ask: how can we determine the unknown parameters from observed data? . 1.3 An idealized view of the beach. The surface is ﬂat and the subsurface consists of little blocks containing either sand or gold. . . . . . . . . . . 1.4 Our preconceptions as to the number of bricks buried in the sand. There is a possibility that someone has already dug up the gold, in which case the number of gold blocks is zero. But we thing it’s most likely that there are 6 gold blocks. Possibly 7, but deﬁnitely not 3, for example. Since this preconception represents information we have independent of the gravity data, or prior to the measurements, it’s an example of what is called a priori information. . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Pirate chests were well to move around much. are clustered together. dispersed, but it seems made. And gold, being rather heavy, is unlikely So we think it’s mostly likely that the gold bars It’s not impossible that the bars have become unlikely. . . . . . . . . . . . . . . . . . . . . . .

2.6 PT |O , the probability that the true density is x given some observed value. 18 2.7 A priori we know that the density of kryptonite cannot be less than 5.1 or greater than 5.6. If we’re sure of this than we can reject any observed density outside of this region. . . . . . . . . . . . . . . . . . . . . . . .

20

3.1 Simple model of a vertical seismic proﬁle (VSP). An acoustic source is at the surface of the Earth near a vertical bore-hole (left side). A receiver is lowered into the bore-hole, recording the pulses of down-going sound at various depths below the surface. From these recorded pulses (right) we can extract the travel time of the ﬁrst-arriving energy. These travel times are used to construct a best-ﬁtting model of the subsurface wavespeed (velocity). Here vi refers to the velocity in discrete layers, assumed to be constant. How we discretize a continuous velocity function into a ﬁnite numer of discrete values is tricky. But for now we will ignore this issue and just assume that it can be done. . . . . . . . . . . . . . . . . . . . 3.2 Noise is just that portion of the data we have no interest in explaining. The x’s indicate hypothetical measurements. If the measurements are very noisy, then a model whose response is a straight line might ﬁt the data (curve 1). The more precisely the data are known, the more structure is required to ﬁt them. . . . . . . . . . . . . . . . . . . . . . . 3.3 Observed data (solid curve) and predicted data for two diﬀerent assumed levels of noise. In the optimistic case (dashed curve) we assume the data are accurate to 0.3 ms. In the more pessimistic case (dotted curve), we assume the data are accurate to only 1.0 ms. In both cases the predicted travel times are computed for a model that just ﬁts the data. In other words we perturb the model until the RMS misﬁt between the observed and predicted data is about N times 0.3 or 1.0, where N is the number of observations. Here N = 78. I.e., N χ2 = 78 × 1.0 for the pessimistic case, and N χ2 = 78 × .3 for the optimistic case. . . . . . . . . . . . . .

26

27

30

LIST OF FIGURES 3.4 The true model (solid curve) and the models obtained by a truncated SVD expansion for the two levels of noise, optimistic (0.3 ms, dashed curve) and pessimistic (1.0 ms, dotted curve). Both of these models just ﬁt the data in the sense that we eliminate as many singular vectors as possible and still ﬁt the data to within 1 standard deviation (normalized χ2 = 1). An upper bound of 4 has also been imposed on the velocity. The data ﬁt is calculated for the constrained model. . . . . . . . . . . . 4.1 Family of p norm solutions to the optimization problem for various values of the parameter λ. In accordance with the uniqueness theorem, we can see that the solutions are indeed unique for all values of p > 1, but that for p = 1 this breaks down at the point λ = 1. For λ = 1 there is a cusp in the curve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Shape of the generalized Gaussian distribution for several values of p. . 4.3 Let a and b be any two vectors. We can always represent one, say b, in terms of its components parallel and perpendicular to the other. The length of the component of b along a is b cos θ which is also bT a/ a . 6.1 Examples of the intersection, union, and complement of sets. . . . . . . 6.2 The title of Bayes’ article, published posthumously in the Philosophical Transactions of the Royal Society, Volume 53, pages 370–418, 1763 . . . 6.3 Bayes’ statement of the problem. . . . . . . . . . . . . . . . . . . . . . 6.4 A normal distribution of zero mean and unit variance. Almost all the area under this curve is contained within 3 standard deviations of the mean. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Ouput from the coin-ﬂipping program. The histograms show the outcomes of a calculation simulating the repeated ﬂipping of a fair coin. The histograms have been normalized by the number of trials, so what we are actually plotting is the relative probability of of ﬂipping k heads out of 100. The central limit theorem guarantees that this curve has a Gaussian shape, even though the underlying probability of the random variable is not Gaussian. . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Two Gaussian sequences (top) with approximately the same mean, standard deviation and 1D distributions, but which look very diﬀerent. In the middle of this ﬁgure are shown the autocorrelations of these two sequences. Question: suppose we took the samples in one of these time series and sorted them in order of size. Would this preserve the nice bell-shaped curve? . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

LIST OF FIGURES 8.10 The distribution of singular values (top). A well resolved model singular vector (middle) and a poorly resolved singular vector (bottom). In this cross well experiment, the rays travel from left to right across the ﬁgure. Thus, features which vary with depth are well resolved, while features which vary with the horizontal distance are poorly resolved. . . . . . . 10.1 For square error loss, the Bayes risk associated with a uniform prior is shown along with the upper and lower bounds on the minimax risk as a function of the size of the bounding interval [−β, β]. When β is comparable to or less than the variance (1 in this case), the risk associated with a uniform prior is optimistic . . . . . . . . . . . . . . . . . . . . . 11.1 Contours of the quadratic form associated with the linear system Ax = h where A = diag(10, 1) and h = (1, −1). Superposed on top of the contours are the solution vectors for the ﬁrst few iterations. . . . . . . .

xi

129

143

159

xii

LIST OF FIGURES

Chapter 1 What Is Inverse Theory
This course is an introduction to some of the balkanized family of techniques and philosophies that reside under the umbrella of inverse theory. In this section we present the central threads that bind all of these singular items together in a harmonious whole. That’s impossible of course, but what we will do is provide a point of view that, while it will break from time-to-time, is good enough to proceed with. The goal of this chapter is to introduce a real inverse problem and explore some of the issues that arise in a non-technical way. Later, we explore the resulting complications in greater depth. Suppose that we ﬁnd ourselves on a gleaming white beach somewhere in the Caribbean with • time on our hands, • a gravimeter (a little device that measures changes in gravitational acceleration), and • the certain conviction that a golden blob of pirate booty lies somewhere beneath us. In pursuit of wealth we make a series of measurements of gravity at several points along the surface. Our mental picture looks like Figure 1.1. And although we don’t know where the gold actually is, or what amount is present, we’re pretty sure something is there. How can we use these observations to decide where the pirate gold lies and how much is present? It’s not enough to know that gold (ρ = 19.3gm/cm 3 ) is denser than sand (ρ = 2.2gm/cm3 ) and that the observed gravity should be greater above our future wealth. Suppose that we observe relative gravity values of (from left to right) 22, 34, 30, 24, and 55µ gals
1

2

What Is Inverse Theory

Measurements

Surface

x

x x

x

Sand

Figure 1.1: We think that gold is buried under the sand so we make measurements of gravity at various locations on the surface. respectively.a There’s no simple formula, (at least not that we know) into which we can plug ﬁve observed gravity observations and receive in return the depth and size of our target. So what shall we do? One thing we do know is φ(r) = Gρ(r ) dV |r − r | (1.1)

that is, Newtonian gravitation. (If you didn’t know it before, you know it now.) Equation 1.1 relates the gravitational potential, φ, to density, ρ. Equation 1.1 has two interesting properties: • it expresses something we think is true about the physics of a continuum, and • it can be turned into an algorithm which we can apply to a given density ﬁeld So although we don’t know how to turn our gravity measurements into direct information about the density in the earth beneath us, we do know how to go in the other direction: given the density in the earth beneath us, we know how to predict the gravity ﬁeld we should observe. Inverse theory begins here, as in Figure 1.2. For openers, we might write a computer program that accepts densities as inputs and produces predicted gravity values as outputs. Once we have such a tool we can play with diﬀerent density values to see what kind of gravity observations we would get. We might assume that the gold is a rectangular block of the same dimensions as a standard
a A gal is a unit of acceleration equal to one centimeter per second per second. It is named after Galileo but was ﬁrst used in this century.

1

3
observed gravity

want:

density

have:

predicted gravity

density

Figure 1.2: Inverse problems usually start with some procedure for predicting the response of a physical system with known parameters. Then we ask: how can we determine the unknown parameters from observed data?
surface x x x x predictions x unknown

sand

Figure 1.3: An idealized view of the beach. The surface is ﬂat and the subsurface consists of little blocks containing either sand or gold. pirate’s chest and we could move the block to diﬀerent locations, varying both depth and horizontal location, to see if we can match our gravity observations. Part of writing the gravity program is deﬁning the types of density models we’re going to use. We’ll use a simpliﬁed model of the beach that has a perfectly ﬂat surface, and has a subsurface that consists of a cluster of little rectangles of variable density surrounded by sand with a constant density. We’ve chosen the cluster of little rectangles to include all of the likely locations of the buried treasure. (Did we mention we have a manuscript fragment which appears to be part of a pirate’s diary?) In order to model having the buried treasure at a particular spot in the model we’ll set the density in those rectangles to be equal to the density of gold and we’ll set the density in the rest of the little rectangles to the density of sand. Here’s what the model looks like: The x’s are the locations for which we’ll compute the gravitational ﬁeld. Notice that the values produced by our program are referred to as predictions, rather than observations. Now we have to get down to business and use our program to ﬁgure out where the treasure is located. Suppose we embed our gravity program into a larger program which will
1

4

What Is Inverse Theory • generate all possible models by trying all combinations of sand and gold densities in our little rectangles, and • compare the predicted gravity values to the observed gravity values and tell us which models, if any, agreed well with the observations.

Model space and data space In the beach example a model consists of 45 parameters, namely the content (sand or gold) of each block. We could represent this mathematically as a 45-tuple containing the densities of each block. For example, (2.2, 2.2, 2.2, 19.3, 2, 2 . . .) is an example of a model. Moreover, since we’re only allowing those densities to be that of gold and sand, we might as well consider the 45-tuple as consisting of zeros and ones. Therefore all possible models of the subsurface are elements of the set of 45-tuples whose elements are 0 or 1. There are 245 such models. We call this the model space for our problem. On the other hand, the data space consists of all possible data predictions. For this example there are 5 gravity measurements, so the data space consists of all possible 5-tuples whose elements vary continuously between 0 and some upper limit; i.e., a subset of R5 , the 5-dimensional Euclidean space.

1.1

Too many models

The ﬁrst problem is that there are forty-ﬁve little rectangles under our model beach and so there are 245 ≈ 3 × 1013 (1.2)

models to inspect. If we can evaluate a thousand models per second, it will still take us about 1100 years to complete the search. It is almost always impossible to examine more than the tiniest fraction of the possible answers (models) in any interesting inverse calculation.

1.2

No unique answer

We have forty-ﬁve knobs to play with in our model (one for each little rectangle) and only ﬁve observations to match. It is very likely that there will be more than one bestﬁtting model. This likelihood increases to near certainty once we admit the possibility of noise in the observations. There are almost always many possible answers to an inverse problem which cannot be distinguished by the available observations.
1

1.3 Implausible models
relative likelihood it’s all here

5

someone else found it it’s a little bigger than we thought

0

1 2 3 4 5 6 7 number of gold rectangles in model

8

Figure 1.4: Our preconceptions as to the number of bricks buried in the sand. There is a possibility that someone has already dug up the gold, in which case the number of gold blocks is zero. But we thing it’s most likely that there are 6 gold blocks. Possibly 7, but deﬁnitely not 3, for example. Since this preconception represents information we have independent of the gravity data, or prior to the measurements, it’s an example of what is called a priori information.

1.3

Implausible models

On the basis of outside information (which we can’t reproduce here because we unfortunately left it back at the hotel), we think that the total treasure was about the equivalent of six little rectangles worth of gold. We also think that it was buried in a chest which is probably still intact (they really knew how to make pirate’s chests back then). We can’t, however, be absolutely certain of either belief because storms could have rearranged the beach or broken the chest and scattered the gold about. It’s also possible that someone else has already found it. Based on this information we think that some models are more likely to be correct than others. If we attach a relative likelihood to diﬀerent number of gold rectangles, our prejudices might look like Figure 1.4. You can imagine a single Olympic judge holding up a card as each model is displayed. Similarly, since we think the chest is probably still intact we favor models which have all of the gold rectangles in the two-by-three arrangement typical of pirate chests, and we will regard models with the gold spread widely as less likely. Qualitatively, our thoughts tend towards some speciﬁcation of the relative likelihood of models, even before we’re made any observations, as illustrated in Figure 1.5. This distinction is hard to capture in a quasi-quantitative way.
1

6

What Is Inverse Theory

x x

x x

x x

plausible

x x

x x

x

x

possible

x x

x

x x

x

unlikely

Figure 1.5: Pirate chests were well made. And gold, being rather heavy, is unlikely to move around much. So we think it’s mostly likely that the gold bars are clustered together. It’s not impossible that the bars have become dispersed, but it seems unlikely.

A priori information Information which is independent of the observations, such as that models with the gold bars clustered are more likely than those in which the bars are dispersed, is called a priori information. We will continually make the distinction between a priori (or simply prior, meaning before) and a posteriori (or simply posterior, meaning after) information. Posterior information is the result of the inferences we make from data and the prior information. What we’ve called plausibility really amounts to information about the subsurface that is independent of the gravity observations. Here the information was historic and took the form of prejudices about how likely certain model conﬁgurations were with respect to one another. This information is independent of, and should be used in addition to, the gravity observations we have.

1.4

Observations are noisy

Most observations are subject to noise and gravity observations are particularly delicate. If we have two models that produce predicted values that lie within reasonable errors of the observed values, we probably don’t want to put much emphasis on the possibility that one of the models may ﬁt slightly better than the other. Clearly learning what the observations have to tell us requires that we take account of noise in the observations.
1

1.5 The beach is not a model

7

Nature (real beach) real physics real gravity

transducer observed gravity (gravimeter)

corrections for reality

corrected observed gravity

Figure 1.6: The path connecting nature and the corrected observations is long and diﬃcult.

1.5

The beach is not a model

A stickier issue is that the real beach is deﬁnitely not one of the possible models we consider. The real beach • is three-dimensional, • has an irregular surface, • has objects in addition to sand and gold within it (bones and rum bottles, for example) • has an ocean nearby, and is embedded in a planet that has lots of mass of its own and which is subject to perceptible gravitational attraction by the Moon and Sun, • etc Some of these eﬀects, such as the beach’s irregular surface and the gravitational eﬀects due to things other than the beach (the ocean, earth, Moon, Sun), we might try to eliminate by correcting the observations (it would probably be more accurate to call it erroring the observations). We would change the values we are trying to ﬁt and, likely, increasing their error estimates. The observational process looks more or less like Figure 1.6 The wonder of it is that it works at all.
1

8

What Is Inverse Theory

Other eﬀects, such as the three-dimensionality of reality, we might handle by altering the model to make each rectangle three-dimensional or by attaching modeling errors to the predicted values.

1.6

Summary

Inverse theory is concerned with the problem of making inferences about physical systems from data (usually remotely sensed). Since nearly all data are subject to some uncertainty, these inferences are usually statistical. Further, since one can only record ﬁnitely many (noisy) data and since physical systems are usually modeled by continuum equations, if there is a single model that ﬁts the data there will be an inﬁnity of them. To make these inferences quantitative one must answer three fundamental questions. How accurately are the data known? I.e., what does it mean to “ﬁt” the data. How accurately can we model the response of the system? In other words, have we included all the physics in the model that contribute signiﬁcantly to the data. Finally, what is known about the system independent of the data? Because for any suﬃciently ﬁne parameterization of a system there will be unreasonable models that ﬁt the data too, there must be a systematic procedure for rejecting these unreasonable models.

1.7

Beach Example

Here we show an example of the beach calculation. With the graphical user interface shown in Figure 1.7 we can ﬁddle with the locations of the gold/sand rectangles and visually try to match the “observed” data. For this particular calculation, the true model has 6 buried gold bricks as shown in Figure 1.7. In Figure 1.8 we show but one example of a model that predicts the data essentially as well. The diﬀerence between the observed and predicted data is not exactly zero, but given the noise that would be present in our measurements, it’s almost certainly good enough. So we see that two fundamentally diﬀerent models predict the data about equally well.

1

1.7 Beach Example

9

Figure 1.7: The true distribution of gold bricks.

1

10

What Is Inverse Theory

Figure 1.8: An unreasonable model that predicts the data.

1

Chapter 2 A Simple Inverse Problem that Isn’t
Now we’re going to take a look at another inverse problem: estimating the density of the material in a body from information about the body’s weight and volume. Although this sounds like a problem that is too simple to be of any interest to real inverters, we are going to show you that it is prey to exactly the same theoretical problems as an attempt to model the three-dimensional elastic structure of the earth from seismic observations. Here’s a piece of something (Figure 2.1): It’s green, moderately heavy, and it appears to glow slightly (as indicated by the tastefully drawn rays in the ﬁgure). The chunk is actually a piece of kryptonite, one of the few materials for which physical properties are not available in handbooks. Our goal is to estimate the chunk’s density (which is just the mass per unit volume). Density is just a scalar, such as 7.34, and we’ll use ρ to denote various estimates of its value. Let’s use K to denote the chunk (so we don’t have to say chunk again and again).

Figure 2.1: A chunk of kryptonite. Unfortunately, kryptonite’s properties do not appear to be in the handbooks.
1

Figure 2.3: A scale may or may not measure mass directly. In this case, it actually measures the force of gravity on the mass. This is then used to infer mass via Hooke’s law.

2.1

A First Stab at ρ

In order to estimate the chunk’s density we need to learn its volume and its mass.

2.1.1

Measuring Volume

We measure volume with an instrument called a pycnometer. Our pycnometer consists of a calibrated beaker partially ﬁlled with water. If we put K in the beaker, it sinks (which tells us right away that K is denser than water). If the ﬂuid level in the beaker is high enough to completely cover K, and if we record the volume of ﬂuid in the beaker with and without K in it, then the diﬀerence in apparent ﬂuid volume is equal to the volume of K. Figure 2.2 shows a picture of everyman’s pycnometer. V denotes the change in volume due to adding K to the beaker.

2.1.2

Measuring Mass

We seldom actually measure mass. What we usually measure is the force exerted on an object by the local gravitational ﬁeld, that is, we put it on a scale and record the resultant force on the scale (Figure 2.3). In this instance, we measure the force by measuring the compression of the spring holding K up. We then convert that to mass by knowing (1) the local value of the Earth’s gravitational ﬁeld, and (2) the (presumed linear) relation between spring extension and
1

For many purposes, this story could end now. We have found an answer to our original problem (measuring the density of K). We don’t know anything (yet) about the shortcomings of our answer, but we haven’t had to do much work to get this point. However, we, being scientists, are perforce driven to consider this issue at a more fundamental level.

2.2.1

Errors in Mass Measurement

For simplicity, let’s stipulate that the volume measurement is essentially error-free, and let’s focus on errors in the measurement of mass. To estimate errors due to the scale, we can take an object that we knowa and measure its mass a large number of times. We then plot the distribution (relative frequency) of the measured masses when we had a ﬁxed standard mass. The results looks like Figure 2.4.
a An object with known properties is a standard. Roughly speaking, an object functions as a standard if the uncertainty in knowledge of the object’s properties is at least ten times smaller than the uncertainty in the current measurement. Clearly, a given object can be a standard in some circumstances and the object of investigation in others.

1

14
p(x) probability of measuring x when the correct value is 5.2

A Simple Inverse Problem that Isn’t

probability of measuring 5.2

probability of measuring 5.4

x 5.2 5.4

Figure 2.4: Pay careful attention to the content of this ﬁgure: It tells us the distribution of measurement outcomes for a particular true value.

Physics News Number 183 by Phillip F. Schewe Improved mass values for nine elements and for the neutron have been published by an MIT research team, opening possibilities for a truly fundamental deﬁnition of the kilogram as well as the most precise direct test yet of Einstein’s equation E = mc2 . The new mass values, for elements such as hydrogen, deuterium, and oxygen-16, are 20-1000 times more accurate than previous ones, with uncertainties in the range of 100 parts per trillion. To determine the masses, the MIT team, led by David Pritchard, traps single ions in electric and magnetic ﬁelds and obtains each ion’s mass-to-charge ratio by measuring its cyclotron frequency, the rate at which it circles about in the magnetic ﬁeld. The trapped ions, in general, are charged molecules containing the atoms of interest, and from their measurements the researchers can extract values for individual atomic masses. One important atom in the MIT mass table is silicon-28. With the new mass value and comparably accurate measurements of the density and the lattice spacing of ultrapure Si-28, a new fundamental deﬁnition of the kilogram (replacing the kilogram artifact in Paris) could be possible. The MIT team also plans to participate in a test of E = mc2 by using its mass values of nitrogen-14, nitrogen-15, and a neutron. When N-14 and a neutron combine, the resulting N- 15 atom is not as heavy as the sum of its parts, because it converts some of its mass into energy by releasing gamma rays. In an upcoming experiment in Grenoble, France there are plans to measure the ”E” side of the equation by making highly accurate measurements of these gamma rays. (F. DeFilippo et al, Physical Review Letters, 12 September.)
1

2.3 What is an Answer?

15

2.3

What is an Answer?

Let’s consider how we can use this information to reﬁne the results of our experiment. Since we have an observation (namely 5.2) we’d like to know the probability that the true density has a particular value, say 5.4. This is going to be a little tricky, and it’s going to lead us into some unusual topics. We need to proceed with caution, and for that we need to sort out some notation.

2.3.1

Conditional Probabilities

Let ρO be the value of density we compute after measuring the volume and mass of K; we will refer to ρO as the observed density. Let ρT be the actual value of K’s density; we will refer to ρT as the true density.b Let PO|T (ρO , ρT ) denote the conditional probability that we would measure ρO if the true density was ρT . The quantity plotted above is PO|T (ρO , 5.2), the probability that we would observe ρO if the true density was 5.2.

A few observations First, keep in mind that in general we don’t know what the true value of the density is. But if we nonetheless made repeated measurements we would still be mapping out PO|T , only this time it would be PO|T (ρO , ρT ). And secondly, you’ll notice in the ﬁgure above that the true value of the density does not lie exactly at the peak of our distribution of observations. This must be the result of some kind of systematic error in the experiment. Perhaps the scale is biased; perhaps we’ve got a bad A/D converter; perhaps there was a steady breeze blowing in the window of the lab that day. A distinction is usually made between modeling or theoretical errors and random errors. A good example of a modeling error, would be assuming that K were pure kryptonite, when in fact it is an alloy of kryptonite and titanium. So in this case our theory is slightly wrong. In fact, we normally think of random noise as being the small scale ﬂuctuations which occur when a measurement is repeated. Unfortunately this distinction is hard to maintain in practice. Few experiments are truly repeatable. So when we try to repeat it, we’re actually introducing small changes into the assumptions; as we repeatedly pick up K and put it back down on the scale, perhaps little bits ﬂeck oﬀ, or some perspiration from our hands sticks to the sample, or we disturb the balance of the scale slightly by touching it. An even better example would be the positions of the gravimeters in the buried treasure example. We need to know these to do the modeling.
b We will later consider whether this deﬁnition must be made more precise, but for now we will avoid the issue.

1

16

A Simple Inverse Problem that Isn’t

But every time we pick up the gravimeter and put it back to repeat the observation, we misposition it slightly. Do we regard these mispositionings as noise or do we regard them as actual model parameters that we wish to infer? Do we regard the wind blowing near the trees during our seismic experiment as noise, or could we actually infer the speed of the wind from the seismic data? In fact, recent work in meterology has shown how microseismic noise (caused by waves at sea) can be used to make inferences about climate. As far as we can tell, the distinction between random errors and theoretical errors is somewhat arbitrary and up to us to decide on a case by case. What it boils down to are: what features are we really interested in? Noise consists of those features of the data we have no intest in explaining. For more details see the commentary: What is Noise? [SS98].

2.3.2

What We’re Really (Really) After

What we want is PT |O (ρT , ρO ), the probability that ρT has a particular value given that we have the observed value ρO . Because PT |O and PO|T appear to be relations between the same quantities, and because they look symmetric, it’s tempting to make the connection PT |O (ρT , ρO ) = PO|T (ρO , ρT ) ? but unfortunately it’s not true. What is the correct expression for PT |O ? More important, how can we think our way through issues like this? We’ll start with the last question. One fruitful way to think about these issues is in terms of a simple, repeated experiment. Consider the quantity we already have: PO|T , which we plotted earlier. It’s easy to imagine the process of repeatedly weighing a mass and recording the results. If we did this, we could directly construct tables of PO|T .

2.3.3

A (Short) Tale of Two Experiments

Now consider repeatedly estimating density. There are two ways we might think of this. In one experiment we repeatedly estimate the density of a particular, given chunk of kryptonite. In the second experiment we repeatedly draw a chunk of kryptonite from some source and estimate its density. These experiments appear to be quite diﬀerent. The ﬁrst experiment sounds just like the measurements we (or someone) made to estimate errors in the scale, except in this case we don’t know the object’s mass to begin with. The second experiment has an
1

2.3 What is an Answer?

17

Experiment 1 given a chunk: 1. estimate its density. 2. go to 1.

Experiment 2 1. get a chunk. 2. estimate its density. 3. go to 1.

one chunk

many chunks

Figure 2.5: Two apparently diﬀerent experiments. entirely new aspect: selecting a chunk from a pool or source of chunks.c Now we’re going to do two things: • We’re going to persuade you (we hope) that both experiments are in fact the same, and they both involve acquiring (in principle) multiple chunks from some source. • We’re going to show you how to compute PT |O when the nature of the source of chunks is known and its character understood. After that we’ll tackle (and never fully resolve) the thorny but very interesting issue of dealing with sources that are not well-understood.

2.3.4

The Experiments Are Identical

Repetition Doesn’t Aﬀect Logical Structure In the ﬁrst experiment we accepted a particular K and measured its density repeatedly by conducting repeated weighings. The number of times we weigh a given chunk aﬀects the precision of the measurement but it does not aﬀect the logical structure of the experiment. If we weigh each chunk (whether we use one chunk or many) one hundred times and average the results, the mass estimate for each chunk will be more precise, because we have reduced uncorrelated errors through averaging; we could achieve the
c

The Edmund Scientiﬁc catalog might be a good bet, although we didn’t ﬁnd kryptonite in it.
1

18
p(x) probability of a true density of x when the observed value is 5.2

A Simple Inverse Problem that Isn’t

probability of true density of 5.2

probability of true density of 5.4

x 5.2 5.4

Figure 2.6: PT |O , the probability that the true density is x given some observed value. same eﬀect by using a correspondingly better scale. This issue is experimentally significant but it is irrelevant to understanding the probabilistic structure of the experiment. For simplicity, then, we will assume that in both experiments, a particular chunk is measured only once.

Answer is Always a Distribution In the (now slightly modiﬁed) ﬁrst experiment, we are given a particular chunk, K, and we make a single estimate of its mass, namely ρO . Since the scale is noisy, we have to express our knowledge of ρT , the true density, as a distribution showing the probability that the true density has some value given that the observed density has some other value. Our ﬁrst guess is that it might have the gaussianish form that we had for PO|T inFigure 2.4. So Figure 2.6 shows the suggested form for PT |O constructed by cloning the earlier ﬁgure.

A Priori Pops Up This looks pretty good until we consider whether or not we know anything about the density of kryptonite outside of the measurements we have made.
1

2.3 What is an Answer? Suppose ρT is Known Suppose that we know that the density of kryptonite is exactly ρT = 1.7π In that case, we must have PT |O (ρT , ρO ) = δ(ρT − 1.7π) (where δ(x) is the Dirac delta-function) no matter what the observed value ρ O is.

19

We are not asserting that the observed densities are all equal to 1.7π: the observations are still subject to measurement noise. We do claim that the observations must always be consistent with the required value of ρT (or that some element of this theory is wrong). This shows clearly that PT |O = PO|T since one is a delta function, while the other must show the eﬀects of experimental errors. Suppose ρT is Constrained Suppose that we don’t know the true density of K exactly, but we’re sure it lies within some range of values: CK if 5.6 > ρT > 5.1 P (ρT ) = 0 otherwise where CK is a constant and P refers to the probability distribution of possible values of the density. In that case, we’d expect PT |O must be zero for impossible values of ρT but should have the same shape everywhere else since the density distribution of chunks taken from the pool is ﬂat for those values. (The distribution does have to be renormalized, so that the probability of getting some value is one, but we can ignore this for now.) So we’d expect something like Figure 2.7. What Are We Supposed to Learn from All This? We hope it’s clear from these examples that the ﬁnal value of PT |O depends upon both the errors in the measurement process and the distribution of possible true values determined by the source from which we acquired our sample(s). This is clearly the case for the second type of experiment (in which we draw multiple samples from a pool), but we have just shown above that it is also true when we have but a single sample and a single measurement. One of the reasons we aﬀord so much attention to the simple one-sample experiment is that in geophysics we typically have only one sample, namely Earth. What we’re supposed to learn from all this, then, is
1

20
P(x)

A Simple Inverse Problem that Isn’t

Probability of a true density of x when the observed value is 5.2 probability of a true density of 5.2

Figure 2.7: A priori we know that the density of kryptonite cannot be less than 5.1 or greater than 5.6. If we’re sure of this than we can reject any observed density outside of this region.

Conclusion 1: The correct a posteriori conditional distribution of density, PT |O , depends in part upon the a priori distribution of true densities. Conclusion 2: This connection holds even if the experiment consists of a single measurement on a single sample.

2.4

What does it mean to condition on the truth?

The kryptonite example hinges on a very subtle idea: when we make repeated measurements of the density of the sample, we are mapping out the probability PO|T even though we don’t know the true density. How can this be? We have a state of knowledge about the kryptonite density that depends on measurements and prior information. If we treat the prior information as a probability, then we are considering a hypothetical range of kryptonite densities any one of which, according to the prior probability, could be the true value. So the variability in our knowledge of the density is partly due to the range of possible a priori true density values, and partly due to the experimental variation in the measurements. However, when we make repeated measurements of a single chunk of kryptonite, we are not considering the universe of possible kryptonites, but just the one we are measuring. And so this repeated measurement is in fact conditioned on the true value of the density even though we don’t know it.
1

2.4 What does it mean to condition on the truth?

21

Let us consider the simplest possible case, one observation, one parameter connected by the forward problem: d=m+ . Assume that the prior distribution for m is N (0, β 2 ) (the normal or Gaussian probability with 0 mean and variance β 2 ). Assume that the experimental error is N (0, σ 2 ). If we make repeated measurement of d on the same physical system (ﬁxed m), then the measurements will be centered about m (assuming no systematic errors) with variance just due to the experimental errors, σ 2 . So we conclude that the probability (which we will call f ) of d given m is f (d|m) = N (m, σ 2 ). (2.3) The deﬁnition of conditional probability is that f (d, m) = f (d|m)f (m) (2.4)

So, if measuring the density repeatedly maps out f (d|m), then what is f (d)? We can get f (d) formally by just integrating f (d, m) over all m:
∞ −∞

f (d) ≡

f (d, m)dm =

exp −

1 1 (d − m)2 × exp − 2 m2 dm. 2σ 2 2β

This is the deﬁnition of a marginal probability. But now you can see that the variations in f (d) depend on the a priori variations in m—we’re integrating over the universe of possible m values. This is deﬁnitely not what we do when we make a measurement.

2.4.1

Another example

Here is a more complicated example of the same idea, which we extend to the solution of a toy “inverse” problem. It involves using n measurements and a normal prior to estimate a normal mean. Assume that there are n observations d = (d1 , d2 , ...dn ) which are iidd N (a, σ 2 ) and that we want to estimate the mean a given that the prior on a f (a) is N (µ, β 2 ). Up to a constant factor, the joint distribution for a and d is:
d The term iid is used to denote independent, identically distributed random variables. This means that the random variables are statistically independent of one another and they all have the same probability law.

1

22

A Simple Inverse Problem that Isn’t

f (d, m) = exp −

1 2σ 2

n i=1

(di − m)2 exp −

1 (m − µ)2 , 2β 2

(2.6)

As we saw above, the ﬁrst term on the right is the probability f (d|m) Now the following result, known as Bayes theorem, is treated in detail later in book, but it is easy to derive from the deﬁnition of conditional probability, so we’ll give it here too. In a joint probability distribution (i.e., a probability involving more than one random variable), the order of the random variables doesn’t matter, so f (d, m) is the same as f (m, d). Using the deﬁnition of conditional probability twice we have f (d, m) = f (d|m)f (m) and f (m, d) = f (m|d)f (d). So, since f (d, m) = f (m, d), it is clear that f (d|m)f (m) = f (m|d)f (d) from which it follows that f (m|d) = f (d|m)f (m) . f (d) Bayes Theorem (2.7)

The term f (m|d) is traditionally called the posterior (or a posteriori) probability since it is conditioned on the data. Later we will see another interpretation of Bayesian inversion in which f (m|d) is not the posterior. But for now we’ll assume that’s what we’re after, as in the kryptonite study where we called it PT |O . We have everything we need to evaluate f (m|d) except the marginal f (d). So here are the steps in the calculation: • compute f (d) by integrating the joint distribution f (d, m) with respect to m. • form f (m|d) =
f (d|m)f (m) . f (d)

• from f (m|d) compute a “best” estimated value of m by computing the mean of f (m|d). We will discuss later why the posterior mean is what you want to have. If you do this correctly you should get the following for the posterior mean: ¯ nd/σ 2 + µ/β 2 , n/σ 2 + 1/β 2 (2.8)

¯ where d is the mean of the data. By a similar calculation the posterior variance is n/σ 2 1 . + 1/β 2
1

(2.9)

BIBLIOGRAPHY

23

Notice that the posterior variance is always reduced by the presence of a nonzero β. The posterior mean can also be written as 1/β 2 n/σ 2 ¯+ d µ. n/σ 2 + 1/β 2 n/σ 2 + 1/β 2 Later we will see that the posterior mean has a special signiﬁcance in that it minimizes a certain average error (called the risk). Because of this, the posterior mean has its own name: it is called the Bayes estimator. In this example the Bayes estimator is a weighted average of the mean of the data and the mean of the Bayesian prior distribution; the latter is the Bayes estimator before any data have been recorded. Note also that as β → 0, increasingly strong prior information, the estimate converges to the prior mean. As β → ∞, increasingly weak prior information, the Bayes estimate converges to the mean of the data.

Chapter 3 Example: A Vertical Seismic Proﬁle
Here we will look at another simple example of a geophysical inverse calculation. We will cover the technical issues in due course. The goal here is simply to illustrate the fundamental role of data uncertainties in any inverse calculation. In this example we will see that a certain model feature is near the limit of the resolution of the data. Depending on whether we are bold or conservative in assessing the errors of our data, this feature will or will not be required to ﬁt the data. We use a vertical seismic proﬁle (VSP–used in exploration seismology to image the Earth’s near surface) experiment to illustrate how a ﬁtted response depends on the assumed noise level in the data. Figure 3.1 shows the geometry of a VSP. A source of acoustic energy is at the surface near a vertical bore-hole (left side). A receiver is lowered into a bore-hole, recording the travel time of the down-going acoustic pulse. These times are used to construct a “best-ﬁtting” model of the wavespeed as a function of depth v(z). Of course the real velocity is a function of x, y, and z, but since in this example the rays propagate almost vertically, there will be no point in trying ot resolve lateral variations in v. If the Earth is not laterally invariant, this assumption introduces a systematic error into the calculation. For each observation (and hence each ray) the problem of data prediction boils down to computing the following integral: 1 t= d . (3.1) ray v(z) We can simplify the analysis somewhat by introducing the reciprocal velocity (called slowness): s = 1/v. Now the travel time integral is linear in slowness: t=
ray

s(z)d .

(3.2)

If the velocity model v(z) (or slowness s(z)) and the ray paths are known, then the
1

26

Example: A Vertical Seismic Proﬁle

source r1 r2 depth r3 r 4 v1 v2 v3 v4
...

downgoing pulse

vn
Figure 3.1: Simple model of a vertical seismic proﬁle (VSP). An acoustic source is at the surface of the Earth near a vertical bore-hole (left side). A receiver is lowered into the bore-hole, recording the pulses of down-going sound at various depths below the surface. From these recorded pulses (right) we can extract the travel time of the ﬁrst-arriving energy. These travel times are used to construct a best-ﬁtting model of the subsurface wavespeed (velocity). Here vi refers to the velocity in discrete layers, assumed to be constant. How we discretize a continuous velocity function into a ﬁnite numer of discrete values is tricky. But for now we will ignore this issue and just assume that it can be done.

1

27

x x

data

x x x x x

x

x x 1 x x 3 2

measurement
Figure 3.2: Noise is just that portion of the data we have no interest in explaining. The x’s indicate hypothetical measurements. If the measurements are very noisy, then a model whose response is a straight line might ﬁt the data (curve 1). The more precisely the data are known, the more structure is required to ﬁt them.

travel time can be computed by integrating the velocity along the ray path. The goal is to somehow estimate v(z) (or some function of v(z), such as the average velocity in a region), or to estimate ranges of plausible values of v(z). How well a particular v(z) model ﬁts the data depends on how accurately the data are known. Roughly speaking, if the data are known very precisely we will have to work hard to come up with a model that ﬁts them to a reasonable degree. If the data are known only imprecisely, then we can ﬁt them more easily. For example, in the extreme case of only noise, the mean of the noise ﬁts the data. separating signal from noise Consider the hypothetical measurements labeled with x’s in Figure 3.2. Suppose that we construct three diﬀerent models whose predicted data are labeled 1, 2 and 3 in the ﬁgure. If we consider the uncertainty of the measurements to be large, we might might argue that a straight line ﬁts the data (curve 1). If the uncertainties are smaller, them perhaps structure on the order of that shown in the quadratic curve is required (curve 2). If the data are even more precisely known, then more structure (such as shown in curve 3) is required. Unless we know the noise level in the data, to perform a quantitative inverse calculation we have to decide in advance which features we want to try to explain and which we do not. Just as in the gravity problem we ignored all sorts of complicating factors, such as the eﬀects of tides. Here we will ignore the fact that unless v is constant, the rays will bend (refract); this means that the domain of integration in the travel time formula (equation 3.2) depends on the velocity, which we don’t know. We will neglect this issue
1

28

Example: A Vertical Seismic Proﬁle

for now by simply asserting that the rays are straight lines. This would be a reasonable approximation for x-ray, but likely not for sound.

an example As a simple synthetic example we constructed a piecewise constant v(z) using 40 unknown layers. We computed 78 synthetic travel times and contaminated them with Gaussian noise. (The numbers 40 and 78 have no signiﬁcance whatsoever; they’re just pulled from a hat.) The level of the noise doesn’t matter for the present purposes; the point is that given an unknown level of noise in the data, diﬀerent assumptions about this noise will lead to diﬀerent kinds of reconstructions. With the constant velocity layers, the system of forward problems for all 78 rays (Equation 3.2) reduces to t=J ·s (3.3)

where s is the 40-dimensional vector of layer slownesses and J is a matrix whose (i, j) entry is the distance the i-th ray travels in the j-th layer. The details are given Bording et al. [BGL+ 87] or later in Chapter 8. For now, the main point is that Equation 3.3 is simply a numerical approximation of the continuous Equation 3.2. The data mapping, the function that maps models into data, is the inner product of the matrix J and the slowness vector s. The vector s, is another example of a model vector. It results from discretizing a function (slowness as a function of space). The ﬁrst element of s, s1 , is the slowness in the ﬁrst layer, s2 is the slowness in the second layer, and so on. Let to be the i−th observed travel time (which we get by examinging the raw data i shown in Figure 3.1. Let tc (s) be the i-th travel time calculated through an arbitrary i slowness model s (by computing J for the given geometry and taking the dot product in Equation 3.3. Finally, let σi is the uncertainty (standard deviation) of the i-th datum. If the true slowness is st , then the following model of the observed travel times is assumed to hold: to = tc (st ) + i , (3.4) i i where i is a noise term (whose standard deviation is σi ). For this example, our goal is to estimate st . A standard approach to solve this problem is to determine slowness vectors s that make a misﬁt function such as 1 χ (s) = N
2 N i=1

tc (s) − to i i σi

2

,

(3.5)

smaller than some tolerance. Here N is the number of observations. The symbol χ2 is often used to denote this sum because the sum of uncorrelated Gaussian random variables has a distribution known as χ2 by statisticians. Any statistics will have the details, for example the informative and highly entertaining [GS94]. We will come back to this idea later in the course.
1

29 We have assumed that the number of layers is known, 40 in this example, but this is usually not the case. Choosing too many layers may lead to an over-ﬁtting of the data. In other words we may end up ﬁtting noise induced structures. Using an insuﬃcient number of layers will not capture important features in the data. There are tricks and methods to try to avoid over- and under-ﬁtting. In the present example we do not have to worry since we will be using simulated data. To determine the slowness values through (3.5) we have used a truncated SVDa reconstruction, throwing away all the eigenvectors in the generalized inverse approximation of s that are not required to ﬁt the data at the χ2 = 1 level. Fitting the data this level means that, on average, all the predicted data agree with the measurements to within one σ. The resulting model is not unique, but it is representative of models that do not over-ﬁt the data (to the assumed noise level).

3.0.2

Travel time ﬁtting

We will consider the problem of ﬁtting the data under two diﬀerent assumptions about the noise. Figure 3.3 shows the observed and predicted data for models that ﬁt the travel times on average to within 0.3 ms and 1.0 ms. Remember, the actual pseudorandom noise in the data is ﬁxed throughout, all we are changing is our assumption about the noise, which is reﬂected in the data misﬁt criterion. We refer to these as the optimistic (low noise) and pessimistic (high noise) scenarios. You can clearly see that the smaller the assumed noise level in the data, the more the predicted data must follow the pattern of the observed data. It takes a complicated model to predict complicated data! Therefore, we should expect the best ﬁtting model that produced the low noise response to be more complicated than the model that produced the high noise response. If the error bars are large, then a simple model will explain the data. Now let us look at the models that actually ﬁt the data to these diﬀerent noise levels; these are shown in Figure 3.4. It is clear that if the data uncertainty is only 0.3 ms, then the model predicts (or requires) a low velocity zone. However, if the data errors are as much as 1 ms, then a very smooth response is enough to ﬁt the data, in which case a low velocity zone is not required. In fact, for the high noise case essentially a linear v(z) increase will ﬁt the data, while for the low noise case a rather complicated model is required. (In both cases, because of the singularity of J, the variances of the estimated parameters become very large near the bottom of the borehole.) Hopefully this example illustrates the importance of understanding the noise distribua We will study the singular value decomposition (SVD) in great detail later. For now just consider it to be something like a Fourier decomposition of a matrix. From it we can get an approximate inverse of the matrix, which we use to solve Equation3.3. Truncating the SVD is somewhat akin to low-pass ﬁltering a time series in the frequency domain. The more you truncate the simpler the signal.

1

30

Example: A Vertical Seismic Proﬁle

20

15 time (ms)

10 data high noise response low noise response

5

0

0

10

20 depth (m)

30

40

Figure 3.3: Observed data (solid curve) and predicted data for two diﬀerent assumed levels of noise. In the optimistic case (dashed curve) we assume the data are accurate to 0.3 ms. In the more pessimistic case (dotted curve), we assume the data are accurate to only 1.0 ms. In both cases the predicted travel times are computed for a model that just ﬁts the data. In other words we perturb the model until the RMS misﬁt between the observed and predicted data is about N times 0.3 or 1.0, where N is the number of observations. Here N = 78. I.e., N χ2 = 78 × 1.0 for the pessimistic case, and N χ2 = 78 × .3 for the optimistic case.

tion to properly interpret inversion estimates. In this particular case, we didn’t simply pull these standard deviations out of hat. The low value (0.3 ms) is what you happen to get if you assume that the only uncertainties in the data are normally distributed ﬂuctuations about the running mean of the travel times. However, keep in mind that nature doesn’t really know about travel times. Travel times are approximations to the true properties (i.e., ﬁnite bandwidth) of waveforms. Further, the travel times themselves are usually assigned by a human interpreter looking at the waveforms. Based on these considerations, one might be led to conclude that a more reasonable estimate of the uncertainties for real data would be closer to 1 ms than 0.3 ms. In any event, the interpretation of the presence of a low velocity zone should be viewed with some scepticism unless the smaller uncertainty level can be justiﬁed.
1

31

5 4
true model high noise low noise

wave speed (m/s)

3 2 1 0

0

10

20 depth (m)

30

40

Figure 3.4: The true model (solid curve) and the models obtained by a truncated SVD expansion for the two levels of noise, optimistic (0.3 ms, dashed curve) and pessimistic (1.0 ms, dotted curve). Both of these models just ﬁt the data in the sense that we eliminate as many singular vectors as possible and still ﬁt the data to within 1 standard deviation (normalized χ2 = 1). An upper bound of 4 has also been imposed on the velocity. The data ﬁt is calculated for the constrained model.

Chapter 4 A Little Linear Algebra
Linear algebra background The parts of this chapter dealing with linear algebra follow the outstanding book by Strang [Str88] closely. If this summary is too condensed, you would be well advised to spend some time working your way through Strang’s book. One diﬀerence to note however is that Strang’s matrices are m × n, whereas ours are n × m. This is not a big deal, but it can be confusing. We’ll stick with n × m because that is common in geophysics and later we will see that m is the number of model parameters in an inverse calculation.

4.1

Linear Vector Spaces

The only kind of mathematical spaces we will deal with in this course are linear vector spaces. You are already well familiar with concrete examples of such spaces, at least in the geometrical setting of vectors in three-dimensional space. We can add any two, say, force vectors and get another force vector. We can scale any such vector by a numerical quantity and still have a legitimate vector. However, in this course we will use vectors to encapsulate discrete information about models and data. If we record one seismic trace, one second in length at a sample rate of 1000 samples per second, and let each sample be deﬁned by one byte, then we can put these 1000 bytes of information in a 1000-tuple (s1 , s2 , s3 , · · · , s1000 ) (4.1)

Now, the physical vectors have a life independent of the particular 3-tuple we use to represent them. We will get a diﬀerent 3-tuple depending on whether we use cartesian or spherical coordinates, for example; but the force vector itself is independent of these considerations. On the other hand, our use of vector spaces is purely abstract. There is no physical seismogram vector; all we have is the n-tuple sampled from the recorded seismic trace. Further, the mathematical deﬁnition of a vector space is suﬃciently general to incorporate objects that you might not consider as vectors at ﬁrst glance–such as functions and matrices. The deﬁnition of such a space actually requires two sets of objects: a set of vectors V and a one of scalars F . For our purposes the scalars will always be either the real numbers R or the complex numbers C. For this deﬁnition we need the idea of a Cartesian product of two sets. Deﬁnition 1 Cartesian product The Cartesian product A × B of two sets A and B is the set of all ordered pairs (a, b) where a ∈ A and b ∈ B. Deﬁnition 2 Linear Vector Space A linear vector space over a set F of scalars is a set of elements V together with a function called addition from V × V into V and a function called scalar multiplication from F × V into V satisfying the following conditions for all x, y, z ∈ V and all α, β ∈ F : V1: (x + y) + z = x + (y + z) V2: x + y = y + x V3: There is an element 0 in V such that x + 0 = x for all x ∈ V . V4: For each x ∈ V there is an element −x ∈ V such that x + (−x) = 0. V5: α(x + y) = αx + αy V6: (α + β)x = αx + βx V7: α(βx) = (αβ)x V8: 1 · x = x The simplest example of a vector space is Rn , whose vectors are n-tuples of real numbers. Addition and scalar multiplication are deﬁned component-wise: (x1 , x2 , · · · , xn ) + (y1 , y2 , · · · , yn ) = (x1 + y1 , x2 + y2 , · · · , xn + yn ) and α(x1 , x2 , · · · , xn ) = (αx1 , αx2 , · · · , αxn ).
0

(4.3) (4.4)

4.1 Linear Vector Spaces

35

In the case of n = 1 the vector space V and the scalars F are the same. So trivially, F is a vector space over F . A few observations: ﬁrst, by adding −x to both sides of x + y = x, you can show that x + y = x if and only if y = 0. This implies the uniqueness of the zero element and also that α · 0 = 0 for all scalars α. Functions themselves can be vectors. Consider the space of functions mapping some nonempty set onto the scalars, with addition and multiplication deﬁned by: [f + g](t) = f (t) + g(t) and [αf ](t) = αf (t). (4.6) We use the square brackets to separate the function from its arguments. In this case, the zero element is the function whose value is zero everywhere. And the minus element is inherited from the scalars: [−f ](t) = −f (t). (4.5)

4.1.1

Matrices

The set of all n × m matrices with scalar entries is a linear vector space with addition and scalar multiplication deﬁned component-wise. We denote this space by Rn×m . Two matrices have the same dimensions if they have the same number of rows and columns. We use upper case roman letters to denote matrices, lower case romana to denote ordinary vectors and greek letters to denote scalars. For example, let 2 5   A =  3 8 . 1 0
 

(4.7)

Then the components of A are denoted by Aij . The transpose of a matrix, denoted by AT , is achieved by exchanging the columns and rows. In this example AT = Thus A21 = 3 = AT . 12 You can prove for yourself that (AB)T = B T AT . (4.9) 2 3 1 5 8 0 . (4.8)

a For emphasis, and to avoid any possible confusion, we will henceforth also use bold type for ordinary vectors.

0

36

A Little Linear Algebra

A matrix which equals its transpose (AT = A) is said to be symmetric. If AT = −A the matrix is said to be skew-symmetric. We can split any square matrix A into a sum of a symmetric and a skew-symmetric part via 1 1 A = (A + AT ) + (A − AT ). 2 2 4 − i 8 12 + i −12 −8 −4 − i
 

(4.10)

The Hermitian transpose of a matrix is the complex conjugate of its transpose. Thus if A= then (4.11)

Sometimes it is useful to have a special notation for the columns of a matrix. So if (4.13)

then we write where

(4.14)

2   a1 =  3  . 1

(4.15)

Addition of two matrices A and B only makes sense if they have the same number of rows and columns, in which case we can add them component-wise (A + B)ij = [Aij + Bij ] . For example if A= and 1 2 3 −3 −2 −1 0 6 2 1 1 1 1 8 5 −2 −1 0 . (4.17) (4.16)

So both matrices and vectors can be thought of as vectors in the abstract sense. Matrices can also be thought of as operators acting on vectors in Rn via the matrix-vector inner (or “dot”) product. If A ∈ Rn×m and x ∈ Rm , then A · x = y ∈ Rn is deﬁned by
m

By default, a vector x is regarded as a column vector. So this vector-vector inner product is also written as zT x or as (z, x). Similarly if A ∈ Rn×m and B ∈ Rm×p , then the matrix-matrix AB product is deﬁned to be a matrix in Rn×p with components
m

(AB)ij =
k=1

aik bkj .

(4.26)

For example, AB = 1 2 3 4 0 1 2 3 = 4 7 8 15 . (4.27)

On the other hand, note well that BA = 0 1 2 3 1 2 3 4
0

=

3 4 11 16

= AB.

(4.28)

38

A Little Linear Algebra

This deﬁnition of matrix-matrix product even extends to the case in which both matrices are vectors. If x ∈ Rm and y ∈ Rn , then xy (called the “outer” product and usually written as xyT ) is (xy)ij = xi yj . (4.29) So if x= and


−1 1


(4.30)

then xyT =

1   y= 3  0 −1 −3 0 1 3 0 .

(4.31)

(4.32)

4.1.2

Matrices With Special Structure

The identity element in the space of square n × n matrices is a matrix with ones on the main diagonal and zeros everywhere else
        

In =

1 0 0 . . .

0 1 0

0 0 1

0 ... 0

0 ...  0 ...   0 ... .   ..  .  0 1



(4.33)

Even if the matrix is not square, there is still a main diagonal of elements given by Aii where i runs from 1 to the smaller of the number of rows and columns. We can take any vector in Rn and make a diagonal matrix out of it just by putting it on the main diagonal and ﬁlling in the rest of the elements of the matrix with zeros. There is a special notation for this:
        

A matrix Q ∈ Rn×n is said to be orthogonal if QT Q = In . In this case, each column of Q is an orthornormal vector: qi · qi = 1. So why are these matrices called orthogonal? No good reason. As an example 1 Q= √ 2
0

1 −1 1 1

.

(4.35)

4.2 Matrix and Vector Norms

39

Now convince yourself that QT Q = In implies that QQT = In as well. In this case the rows of Q must be orthonormal vectors too. Another interpretation of the matrix-vector inner product is as a mapping from one vector space to another. Suppose A ∈ Rn×m , then A maps vectors in Rm into vectors in Rn . An orthogonal matrix has an especially nice geometrical interpretation. To see this ﬁrst notice that for any matrix A, the inner product (A · x) · y, which we write as (Ax, y), is equal to (x, AT y), as you will verify in one of the exercises at the end of the chapter. Similarly (AT x, y) = (x, Ay). (4.36) As a result, for an orthogonal matrix Q (Qx, Qx) = (QT Qx, x) = (x, x). (4.37)

Now, as you already know, and we will discuss shortly, the inner product of a vector with itself is related to the length, or norm, of that vector. Therefore an orthogonal matrix maps a vector into another vector of the same norm. In other words it does a rotation.

4.2

Matrix and Vector Norms

We need some way of comparing the relative “size” of vectors and matrices. For scalars, the obvious answer is the absolute value. The absolute value of a scalar has the property that it is never negative and it is zero if and only if the scalar itself is zero. For vectors and matrices both we can deﬁne a generalization of this concept of length called a norm. A norm is a function from the space of vectors onto the scalars, denoted by · satisfying the following properties for any two vectors v and u and any scalar α: Deﬁnition 3 Norms N1: v > 0 for any v = 0 and v = 0 ⇔ v = 0 N2: αv = |α| v N3: v + u ≤ v + u Here we use the symbol ⇔ to mean if and only if. Property N 3 is called the triangle inequality. The most useful class of norms for vectors in Rn is the
n 1/p p

norm deﬁned for p ≥ 1 by (4.38)

x

p

=
i=1

|xi |

p

.

0

40 For p = 2 this is just the ordinary euclidean norm: x p norm exists as p → ∞ called the ∞ norm: x
∞

A Little Linear Algebra
2

=

√

xT x. A ﬁnite limit of the

= max |xi |
1≤i≤n

(4.39)

Any norm on vectors in Rn induces a norm on matrices via A = max
x=0

Ax . x

(4.40)

A matrix norm that is not induced by any vector norm is the Frobenius norm deﬁned for all A ∈ Rn×m as A
F

=



m

n

i=1 j=1

A2  ij

1/2

.

(4.41)

Some examples (see [GvL83]): A 1 = maxj aj 1 where aj is the j-th column of A. Similarly A ∞ is the maximum 1-norm of the rows of A. For the euclidean norm we have ( A 2 )2 = maximum eigenvalue of AT A. The ﬁrst two of these examples are reasonably obvious. The third is far from so, but is the reason the 2 norm of a matrix is called the spectral norm. We will prove this latter result shortly after we’ve reviewed the properties of eigenvalues and eigenvectors.

Minor digression: breakdown of the

p

norm

Since we have alluded in the previous footnote to some diﬃculty with the p norm for p < 1 it might be worth a brief digression on this point in order to emphasize that this diﬃculty is not merely of academic interest. Rather, it has important consequences for the algorithms that we will develop in the chapter on “robust estimation” methods. For the rectangular (and invariably singular) linear systems we will need to solve in inverse calculations, it is useful to pose the problem as one of optimization; to wit, min Ax − y .
x

(4.42)

It can be shown that for the p family of norms, if this optimization problem has a solution, then it is unique: provided the matrix has full column rank and p > 1. (By full column rank we mean that all the columns are linearly independent.) For p = 1 the norm loses, in the technical jargon, strict convexity. A proof of this result can be found in [SG88]. It is easy to illustrate. Suppose we consider the one parameter linear system: 1 1 x= . (4.43) 0 λ
0

4.2 Matrix and Vector Norms
1 0.8 0.6

41

p-norm error
0.4 0.2

p=1.01
0 0.25 0.5 0.75 1 1.25

p=1.1
1.5 1.75

p=2.0 p=1.5
2

λ

Figure 4.1: Family of p norm solutions to the optimization problem for various values of the parameter λ. In accordance with the uniqueness theorem, we can see that the solutions are indeed unique for all values of p > 1, but that for p = 1 this breaks down at the point λ = 1. For λ = 1 there is a cusp in the curve. For simplicity, let us assume that λ ≥ 0 and let us solve the problem on the open interval x ∈ (0, 1). The p error function is just Ep (x) ≡ [|x − 1|p + λp |x|p ]1/p . (4.44) Restricting x ∈ (0, 1) means that we don’t have to deal with the fact that the absolute value function is not diﬀerentiable at the origin. Further, the overall exponent doesn’t aﬀect the critical points (points where the derivative vanishes) of Ep . So we ﬁnd that ∂x EP (x) = 0 if and only if 1 − x p−1 = λp (4.45) x from which we deduce that the p norm solution of the optimization problem is xp= 1 . 1 + λp/(p−1) (4.46)

But remember, λ is just a parameter. The theorem just alluded to guarantees that this problem has a unique solution for any λ provided p > 1. A plot of these solutions as a function of λ is given in Figure (4.1). This family of solutions is obviously converging to a step function as p → 1. And since this function is not single-valued at λ = 1, you can see why the uniqueness theorem is only valid for p > 1 Interpretation of the norms

p

When we are faced with optimization problems of the form min Ax − y
x
p

(4.47)

0

42

A Little Linear Algebra

the question naturally arises: which p is best? There are two aspects of this question. The ﬁrst is purely numerical. It turns out that some of the p norms have more stable numerical properties than others. In particular, as we will see, p values near 1 are more stable than p values near 2. On the other hand, there is an important statistical aspect of this question. When we are doing inverse calculations, the vector y is associated with our data. If our data have, say, a Gaussian distribution, then 2 is optimal in a certain sense to be described shortly. On the other hand, if our data have the double-exponential distribution, then 1 is optimal. This optimality can be quantiﬁed in terms of the entropy or information content of the distribution. For the Gaussian distribution we are used to thinking of this in terms of the variance or standard deviation. More generally, we can deﬁne the p norm dispersion of a given probability density ρ(x) as (σp )p ≡
∞ −∞

|x − x0 |p ρ(x) dx

(4.48)

where x0 is the center of the distribution. (The deﬁnition of the center need not concern us here. The point is simply that the dispersion is a measure of how spread out a probability distribution is.) One can show (cf. [Tar87], Chapter 1) that for a ﬁxed p norm dispersion, the probability density with the minimum information content is given by the generalized gaussian ρp (x) = −1 |x − x0 |p p1−1/p exp 2σp Γ(1/p) p (σp )p (4.49)

where Γ is the Gamma function [MF53]. These distributions are shown in Figure 4.2 for four diﬀerent values of p, 1, 2, 10, and ∞. The reason information content is so important is that being naturally conservative, we want to avoid jumping to any unduly risky conclusions about our data. One way to quantify simplicity is in terms of information content, or entropy: given two (or more) models which ﬁt the data to the same degree, we may want to choose the one with the least information content in order to avoid over-interpreting the data. This is an important caveat for all of inverse theory.b Later in the course we will come back to what it means to be “conservative” and see that the matter is more complicated than it might ﬁrst appear.

4.3

Projecting Vectors Onto Other Vectors

Figure 4.3 illustrates the basic idea of projecting one vector onto another. We can always represent one, say b, in terms of its components parallel and perpendicular to the other. The length of the component of b along a is b cos θ which is also bT a/ a
b This is a caveat for all of life too. It is digniﬁed with the title Occam’s razor after William of Occam, an English philosopher of the early 14th century. What Occam actually wrote was: “Entia non sunt multiplicanda praeter necessitatem” (things should not be presumed to exist, or multiplied, beyond necessity).

0

4.3 Projecting Vectors Onto Other Vectors

43

p=1

p=2

p=10

p=infinity

Figure 4.2: Shape of the generalized Gaussian distribution for several values of p.

0

44
y

A Little Linear Algebra

b a-b

θ

a b cos θ
x

Figure 4.3: Let a and b be any two vectors. We can always represent one, say b, in terms of its components parallel and perpendicular to the other. The length of the component of b along a is b cos θ which is also bT a/ a . Now suppose we want to construct a vector in the direction of a but whose length is the component of b along b . We did this, in eﬀect, when we computed the tangential force of gravity on a simple pendulum. What we need to do is multiply b cos θ by a unit vector in the a direction. Obviously a convenient unit vector in the a direction is a/ a , which equals a √ . aT a So a vector in the a with length b cos θ is given by b cos θˆ = a aT b a a a a aT b aaT b aaT = = T = T b a a a a a a (4.50) (4.51)

As an exercise verify that in general a(aT b) = (aaT )b. This is not completely obvious since in one expression there is an inner product in the parenthesis and in the other there is an outer product. What we’ve managed to show is that the projection of the vector b into the direction of a can be achieved with the following matrix (operator) aaT . aT a This is our ﬁrst example of a projection operator.
0

4.4 Linear Dependence and Independence

45

4.4

Linear Dependence and Independence
{x1 , x2 , · · · , xn } (4.52)

Suppose we have n vectors of the same dimension. The question is, under what circumstances can the linear combination of these vectors be zero: α1 x1 + α2 x2 + · · · αn xn = 0. (4.53)

If this is true with at least one of the coeﬃcients αi nonzero, then we could isolate a particular vector on the right hand side, expressing it as a linear combination of the other vectors. In this case the original set of n vectors are said to be linearly dependent. On the other hand, if the only way for this sum of vectors to be zero is for all the coeﬃcients themselves to be zero, then we say that the vectors are linearly independent. Now, this linear combination of vectors can also be written as a matrix-vector inner product. With a = (α1 , α2 , · · · , αn ), and X = (x1 , x2 , · · · , xn ) we have the condition for linear dependence being Xa = 0 (4.54) for some nonzero vector a, and the condition for linear independence being Xa = 0 ⇒ a = 0. As a result, if we are faced with a linear system of equations to solve Ax = b (4.56) (4.55)

we can think in two diﬀerent ways. On the one hand, we can investigate the equation in terms of the existence of a vector x satisfying the equation. On the other hand, we can think in terms of the compatibility of the right hand side with the columns of the matrix. Linear independence is also central to the notion of how big a vector space is–its dimension. It’s intuitively clear that no two linearly independent vectors are adequate to represent an arbitrary vector in R3 . For example, (1, 0, 0) and (0, 1, 0) are linearly independent, but there are no scalar coeﬃcients that will let us write (1, 1, 1) as a linear combination of the ﬁrst two. Conversely, since any vector in R3 can be written as a combination of the three vectors (1, 0, 0), (0, 1, 0), and (0, 0, 1), it is impossible to have more than three linearly independent vectors in R3 . The dimension of a space is the number of linearly independent vectors required to represent an arbitrary element.
0

46

A Little Linear Algebra

4.5

The Four Fundamental Spaces

Suppose you have an n−dimensional space and n linearly independent vectors in that space. These vectors are said to be basis vectors since any element of the space can be written as a linear combination of the basis vectors. For instance, two basis vectors for R2 are (1, 0) and (0, 1). Any element of R2 can be written as some constant times (1, 0) plus another constant times (0, 1). Any other pair of linearly independent vectors would also work, such as (2, 0) and (1, 15). OK, so take two basis vectors for R2 and consider all possible linear combinations of them. This the set of all vectors α(1, 0) + β(0, 1), where α and β are arbitrary scalars. This is called the span of the two vectors and in this case it obviously consists of all of R2 . The span of (1, 0) is just the x-axis in R2 . On the other hand, if we consider these two vectors as being in R3 , so that we write them as (1, 0, 0) and (0, 1, 0), then their span clearly doesn’t ﬁll up all of R3 . It does, however, ﬁll up a subspace of R3 , the x−y plane. The technical deﬁnition of a subspace is that it is a subset closed under addition and scalar multiplication: Deﬁnition 4 Subspaces: A subspace of a vector space is a nonempty subset S that satisﬁes S1: The sum of any two elements from S is in S, and S2: The scalar multiple of any element from S is in S. If we take a general matrix A ∈ Rn×m , then the span of the columns must be a subspace of Rn . Whether this subspace amounts to the whole of Rn obviously depends on whether the columns are linearly independent or not. This subspace is called the column space of the matrix and is usually denoted by R(A), for “range”. The dimension of the column space is called the rank of the matrix. Another fundamental subspace associated with any matrix A is associated with the solutions of the homogeneous equation Ax = 0. Why is this a subspace? Take any two such solutions, say x and y and we have A(x + y) = Ax + Ay = 0. Hence Similarly, A(αx) = αAx. (4.58) This subspace is called the nullspace or kernel and is extremely important from the point of view of inverse theory. As we shall see, in an inverse calculation the right
0

(4.57)

4.5 The Four Fundamental Spaces

47

hand side of a matrix equations is usually associated with perturbations to the data. Vectors in the nullspace have no eﬀect on the data and are therefore unresolved in an experiment. Figuring out what features of a model are unresolved is a major goal of inversion.

4.5.1

Spaces associated with a linear system Ax = y

The span of the columns is a subset of Rn and the span of the rows is a subset of Rm . In other words the rows of A have m components while the columns of A have n components. Now the column space and the nullspace are generated by A. What about the column space and the null space of AT ? These are, respectively, the row space and the left nullspace of A. The nullspace and row space are subspaces of Rm , while the column space and the left nullspace are subspaces of Rn . Here is probably the most important result in linear algebra: For any matrix whatsoever, the number of linearly independent rows equals the number of linearly independent columns. We summarize this by saying that row rank = column rank. For a generic n × m matrix, this is not an obvious result. If you haven’t encountered this before, it would be a good idea to review a good linear algebra book, such as [Str88]. We can summarize these spaces as follows: Theorem 1 Fundamental Theorem of Linear Algebra Let A ∈ Rn×m . Then 1: Dimension of column space equals r, the rank. 2: Dimension of nullspace equals m − r. 3: Dimension of row space equals r. 4: Dimension of left nullspace equals n − r. A Geometrical Picture Any vector in the null space of a matrix, must be orthogonal to all the rows (since each component of the matrix dotted into the vector is zero). Therefore all the elements in the null space are orthogonal to all the elements in the row space. In mathematical terminology, the null space and the row space are orthogonal complements of one another. Or, to say the same thing, they are orthogonal subspaces of Rm . Similarly, vectors in the left null space of a matrix are orthogonal to all the columns of this matrix. This means that the left null space of a matrix is the orthogonal complement of the column space; they are orthogonal subspaces of Rn . In other words, orthogonal complement of a subspace S consists of all the vectors x such that (x, y) = 0 for y ∈ S.
0

48

A Little Linear Algebra

4.6

Matrix Inverses

A left inverse of a matrix A ∈ Rn×m is deﬁned to be a matrix B such that BA = I. A right inverse C therefore must satisfy AC = I. (4.60) (4.59)

If there exists a left and a right inverse of A then they must be equal since matrix multiplication is associative: AC = I ⇒ B(AC) = B ⇒ (BA)C = B ⇒ C = B. (4.61)

Now if we have more equations than unknowns then the columns cannot possibly span all of Rn . Certainly the rank r must be less than or equal to n, but it can only equal n if we have at least as many unknowns as equations. The basic existence result is then [Str88]: Theorem 2 Existence of solutions to Ax = y The system Ax = y has at least one solution x for every y (there might be inﬁnitely many solutions) if and only if the columns span Rn (r = n), in which case there exists an m × n right inverse C such that AC = In . This is only possible if n ≤ m.

=

R

n xm

R

m

R

n

Don’t be mislead by the picture above into neglecting the important special case when m = n. The point is that the basic issues of existence and, next, uniqueness, depend on whether there are more or fewer rows than equations. The statement of uniqueness is [Str88]: Theorem 3 Uniqueness of solutions to Ax = y There is at most one solution to Ax = y (there might be none) if and only if the columns of A are linearly independent (r = m), in which case there exists an m × n left inverse B such that BA = Im . This is only possible if n ≥ m.

=

R

n xm

R

m

R

n

Clearly then, in order to have both existence and uniqueness, we must have that r = m = n. This precludes having existence and uniqueness for rectangular matrices. For square matrices m = n, so existence implies uniqueness and uniqueness implies existence. Using the left and right inverses we can ﬁnd solutions to Ax = y: if they exist. For example, given a right inverse A, then since AC = I, we have ACy = y. But since
0

4.7 Eigenvalues and Eigenvectors

49

Ax = y it follows that x = Cy. But C is not necessarily unique. On the other hand, if there exists a left inverse BA = I, then BAx = By, which implies that x = By. Some examples. Consider ﬁrst the case of more equations than unknowns n > m. Let −1 0  A= 0 3   0 0
 

(4.62)

Since the columns are linearly independent and there are more rows than columns, there can be at most one solution. You can readily verify that any matrix of the form −1 0 γ 0 1/3 ι (4.63)

is a left inverse. The particular left inverse given by the formula (AT A)−1 AT (cf. the exercise at the end of this chapter) is the one for which γ and ι are zero. But there are inﬁnitely many other left inverses. As for solutions of Ax = y, if we take the inner product of A with the vector (x1 , x2 )T we get


So, clearly, we must have x1 = −y1 and x2 = 1/3y2 . But, there will not be any solution unless y3 = 0. Next, let’s consider the case of more columns (unknowns) than rows (equations) n < m. Let −1 0 0 (4.65) A= 0 3 0 Here you can readily verify that any matrix of the form


y1 −x1      3x2  =  y2  y3 0







(4.64)

is a right inverse. The particular right inverse (shown in the exercise at the end of this chapter) AT (AAT )−1 corresponds to γ = ι = 0. Now if we look at solutions of the linear system Ax = y with x ∈ R3 and y ∈ R2 we ﬁnd that x1 = −y1 , x2 = 1/3y2 , and that x3 is completely undetermined. So there is an inﬁnite set of solutions corresponding to the diﬀerent values of x3 .

−1 0    0 1/3  γ ι



(4.66)

4.7

Eigenvalues and Eigenvectors

Usually when a matrix operates on a vector, it changes the direction of the vector. But for a special class of vectors, eigenvectors, the action of the matrix is to simply scale
0

50 the vector: Ax = λx.

A Little Linear Algebra

(4.67)

If this is true, then x is an eigenvector of the matrix A associated with the eigenvalue λ. Now, λx equals λIx so we can rearrange this equation and write (A − λI)x = 0. (4.68)

Clearly in order that x be an eigenvector we must choose λ so that (A − λI) has a nullspace and we must choose x so that it lies in that nullspace. That means we must choose λ so that Det(A − λI) = 0. This determinant is a polynomial in λ, called the characteristic polynomial. For example if A= then the characteristic polynomial is λ2 − 10λ + 13 whose roots are √ √ λ = 5 + 2 3, and λ = 5 − 2 3. (4.70) 5 3 4 5 (4.69)

But note well, that these eigenvectors are not unique. Because they solve a homogeneous system, we can multiply them by any scalar we like and not change the fact that they are eigenvectors. This exercise was straightforward. But imagine what would have happened if we had needed to compute the eigenvectors/eigenvalues of a 10 × 10 matrix. Can you imagine having to compute the roots of a 10-th order polynomial? In fact, once you get past order 4, there is no algebraic formula for the roots of a polynomial. The eigenvalue problem is much harder than solving Ax = y. The following theorem gives us the essential computational tool for using eigenvectors.
0

4.7 Eigenvalues and Eigenvectors

51

Theorem 4 Matrix diagonalization Let A be an n × n matrix with n linearly independent eigenvectors. Let S be a matrix whose columns are these eigenvectors. Then S −1 AS is a diagonal matrix Λ whose elements are the eigenvalues of A. The proof is easy. The elements in the ﬁrst column of the product matrix AS are precisely the elements of the vector which is the inner product of A with the ﬁrst column of S. The ﬁrst column of S, say s1 , is, by deﬁnition, an eigenvector of A. Therefore the ﬁrst column of AS is λ1 s1 . Since this is true for all the columns, it follows that AS is a matrix whose columns are λi si . But now we’re in business since [λ1 s1 λ2 s2 · · · λn sn ] = [s1 s2 · · · sn ] diag(λ1 , λ2 , · · · , λn ) ≡ SΛ. (4.75)

Therefore AS = SΛ which means that S −1 AS = Λ. S must be invertible since we’ve assumed that all it’s columns are linearly independent. Some points to keep in mind: • Any matrix in Rn×n with n distinct eigenvalues can be diagonalized. • Because the eigenvectors themselves are not unique, the diagonalizing matrix S is not unique. • Not all square matrices possess n linearly independent eigenvectors. • A matrix can be invertible without being diagonalizable. We can summarize these ideas with a theorem whose proof can be found in linear algebra books. Theorem 5 Linear independence of eigenvectors If n eigenvectors of an n × n matrix correspond to n diﬀerent eigenvalues, then the eigenvectors are linearly independent. An important class of matrices for inverse theory are the real symmetric matrices. The reason is that since we have to deal with rectangular matrices, we often end up treating the matrices AT A and AAT instead. And these two matrices are manifestly symmetric. In the case of real symmetric matrices, the eigenvector/eigenvalue decomposition is especially nice, since in this case the diagonalizing matrix S can be chosen to be an orthogonal matrix Q. Theorem 6 Orthogonal decomposition of a real symmetric matrix A real symmetric matrix A can be factored into A = QΛQT with orthonormal eigenvectors in Q and real eigenvalues in Λ.
0

(4.76)

52

A Little Linear Algebra

4.8

Orthogonal decomposition of rectangular matrices

For dimensional reasons there is clearly no hope of the kind of eigenvector decomposition discussed above being applied to rectangular matrices. However, there is an amazingly useful generalization that pertains if we allow a diﬀerent orthogonal matrix on each side of A. It is called the Singular Value Decomposition (SVD) and works for any matrix whatsoever. Essentially the singular value decomposition generates orthogonal bases of Rm and Rn simultaneously. Theorem 7 Singular value decomposition Any matrix A ∈ Rn×m can be factored as A = U ΛV T (4.77) where the columns of U ∈ Rn×n are eigenvectors of AAT and the columns of V ∈ Rm×m are the eigenvectors of AT A. Λ ∈ Rn×m is a rectangular matrix with the singular values on its main diagonal and zero elsewhere. The singular values are the square roots of the eigenvalues of AT A, which are the same as the nonzero eigenvalues of AAT . Further, there are exactly r nonzero singular values, where r is the rank of A. The columns of U and V span the four fundamental subspaces. The column space of A is spanned by the ﬁrst r columns of U . The row space is spanned by the ﬁrst r columns of V . The left nullspace of A is spanned by the last n − r columns of U . And the nullspace of A is spanned by the last m − r columns of V . A direct approach to the SVD, due to the physicist Lanczos[Lan61], is to make a symmetric matrix out of the rectangular matrix A as follows: Let S= 0 A AT 0 . (4.78)

Since A is in Rn×m , S must be in R(n+m)×(n+m) . m by n or n by m? For the rest of this book we will interpret the matrix A as mapping from the space of model parameters into the space of data–the forward problem. So there are m parameters and n data. But, obviously this is unnecessary for the interpretation of the results. Model space is simply Rm and data space is Rn . And since S is symmetric it has orthogonal eigenvectors wi with real eigenvalues λi Swi = λi wi . (4.79) If we split up the eigenvector wi , which is in Rn+m , into an n-dimensional data part and an m-dimensional model part wi =
0

ui vi

(4.80)

4.8 Orthogonal decomposition of rectangular matrices

53

then the eigenvalue problem for S reduces to two coupled eigenvalue problems, one for A and one for AT A T u i = λ i vi (4.81) Avi = λi ui . We can multiply the ﬁrst of these equations by A and the second by AT to get AT Avi = λi 2 vi AAT ui = λi 2 ui . (4.83) (4.84) (4.82)

So we see, once again, that the model eigenvectors ui are eigenvectors of AAT and the data eigenvectors vi are eigenvectors of AT A. Also note that if we change sign of the eigenvalue we see that (−ui , vi ) is an eigenvector too. So if there are r pairs of nonzero eigenvalues ±λi then there are r eigenvectors of the form (ui , vi ) for the positive λi and r of the form (−ui , vi ) for the negative λi . Keep in mind that the matrices U and V whose columns are the date and model eigenvectors are square (respectively n × n and m × m) and orthogonal. Therefore we have U T U = U U T = In and V T V = V V T = Im . But it is important to distinguish between the eigenvectors associated with zero and nonzero eigenvalues. Let U r and Vr be the matrices whose columns are the r model and data eigenvectors associated with the r nonzero eigenvalues and U0 and V0 be the matrices whose columns are the eigenvectors associated with the zero eigenvalues, and let Λr be the r × r square, diagonal matrix containing the r nonzero eigenvalues. Then we have by 4.81 and 4.82 the following eigenvalue problem AVr = Ur Λr A T Ur = V r Λr AV0 = 0 AT U0 = 0. (4.85) (4.86) (4.87) (4.88)

This is the singular value decomposition. Notice that 0 represent rectangular matrices of zeros. Since Λr is r × r and Λ is n × m then the lower left block of zeros must be n − r × r, the upper right must be r × m − r and the lower right must be n − r × m − r. It is important to keep the subscript r in mind since the fact that A can be reconstructed from the eigenvectors associated with the nonzero eigenvalues means that the
0

54

A Little Linear Algebra

experiment is unable to see the contribution due to the eigenvectors associated with zero eigenvalues. Cornelius Lanczos was born in Hungary in 1893. His family name was L¨wy, but this was changed to avoid the prevailing sentiments in o Hungary against German names. Lanczos did his university work at Budapest where he studied mathematics and physics. He did work in general relativity throughout his life but made many important contributions to numerical analysis, including the development of the Fast Fourier Transform (25 years before Tukey). Lanczos’ books are marvels of clarity. After ﬂeeing Nazi Germany in the 1930s, Lanczos took up residence ﬁrst in the US and then in Dublin, Ireland, where Schr¨dinger had built up a school of theoretical physics. He died on o a trip to his native land in 1974.

4.9

Orthogonal projections

Above we said that the matrices V and U were orthogonal so that V T V = V V T = Im and U T U = U U T = In . There is a nice geometrical picture we can draw for these equations having to do with projections onto lines or subspaces. Let vi denote the ith column of the matrix V . (The same argument applies to U of course.) The outer T product vi vi is an m × m matrix. It is easy to see that the action of this matrix on a vector is to project that vector onto the one-dimensional subspace spanned by vi :
T T vi vi x = (vi x)vi .

A “projection” operator is deﬁned by the property that once you’ve applied it to a vector, applying it again doesn’t change the result: P (P x) = P x, in other words. For T T the operator vi vi this is obviously true since vi vi = 1.
T T Now suppose we consider the sum of two of these projection operators: vi vi + vj vj . m This will project any vector in R onto the plane spanned by vi and vj . We can continue this procedure and deﬁne a projection operator onto the subspace spanned by any number p of the model eigenvectors: p T vi vi . i=1

If we let p = m then we get a projection onto all of Rm . But this must be the identity operator. In eﬀect we’ve just proved the following identity:
m T vi vi = V V T = I. i=1

0

4.10 A few examples

55

On the other hand, if we only include the terms in the sum associated with the r nonzero singular values, then we get a projection operator onto the non-null space (which is the row space). So
r T vi vi = Vr VrT i=1

is a projection operator onto the row space. By the same reasoning
m i=r+1 T vi vi = V0 V0T

is a projection operator onto the null space. Putting this all together we can say that Vr VrT + V0 V0T = I. This says that any vector in Rm can be written in terms of its component in the null space and its component in the row space of A. Let x ∈ Rm , then x = Ix = Vr VrT + V0 V0T x = (x)row + (x)null . (4.90)

4.10

A few examples

This example shows that often matrices with repeated eigenvalues cannot be diagonalized. But symmetric matrices can always be diagonalized. 3 1 0 3

A=

(4.91)

The eigenvalues of this matrix are obviously 3 and 3. This matrix has a one-dimensional family of eigenvectors; any vector of the form (x, 0)T will do. So it cannot be diagonalized, it doesn’t have enough eigenvectors. Now consider A= 3 0 0 3

(4.92)

The eigenvalues of this matrix are still 3 and 3. But it will be diagonalized by any invertible matrix! So, of course, to make our lives simple we will choose an orthogonal matrix. How about 0 1 ? 1 0
0

(4.93)

56 That will do. But so will

A Little Linear Algebra

1 √ 2

−1 1 1 1

.

(4.94)

So, as you can see, repeated eigenvalues give us choice. And for symmetric matrices we nearly always choose to diagonalize with orthogonal matrices.

Exercises 1. Give speciﬁc (nonzero) examples of 2 by 2 matrices satisfying the following properties: A2 = 0, A2 = −I2 , and AB = −BA (4.95) 2. Let A be an upper triangular matrix. Suppose that all the diagonal elements are nonzero. Show that the columns must be linearly independent and that the null-space contains only the zero vector. 3. Figure out the column space and null space of the following two matrices: 1 −1 0 0 and 0 0 0 0 0 0 (4.96)

4. Which of the following two are subspaces of Rn : the plane of all vectors whose ﬁrst component is zero; the plane of all vectors whose ﬁrst component is 1. 5. Let x= Compute x 1 , x 2 , and x 6. Deﬁne the unit
p -ball ∞.

9 −12

.

(4.97)

in the plane R2 as the set of points satisfying x
p

≤ 1.

(4.98)

Draw a picture of this ball for p = 1, 2, 3 and ∞. 7. Show that B = (AT A)−1 AT is a left inverse and C = AT (AAT )−1 is a right inverse of a matrix A, provided that AAT and AT A are invertible. It turns out that AT A is invertible if the rank of A is equal to n, the number of columns; and AAT is invertible if the rank is equal to m, the number of rows. 8. Consider the matrix a b c d

(4.99)

The trace of this matrix is a + d and the determinant is ad − cb. Show by direct calculation that the product of the eigenvalues is equal to the determinant and the sum of the eigenvalues is equal to the trace.
0

BIBLIOGRAPHY

57

9. As we have seen, an orthogonal matrix corresponds to a rotation. Consider the eigenvalue problem for a simple orthogonal matrix such as Q= 0 −1 1 0 (4.100)

How can a rotation map a vector into a multiple of itself? 10. Show that the eigenvalues of Aj are the j-th powers of the eigenvalues of A. 11. Using the SVD show that AAT = U ΛΛT U and AT A = V ΛT ΛV. (4.102) The diagonal matrices ΛΛT ∈ Rm×m and ΛT Λ ∈ Rn×n have diﬀerent dimensions, but they have the same r nonzero elements: σ1 , σ2 , · · · , σr . 12. Compute the SVD of the matrix 1 1 0  A= 0 0 1   0 0 −1
 

(4.101)

(4.103)

directly by computing the eigenvectors of AT A and AAT . Show that the pseudoinverse solution to the linear system Ax = y where y = (1, 2, 1)T is given by averaging the equations. 13. Prove that (Ax, y) = (x, AT y). 14. Prove that if Q is an orthogonal matrix, that Qx is a rotation of x. 15. What happens to the
p

Chapter 5 SVD and Resolution in Least Squares
In section 4.8 we introduced the singular value decomposition (SVD). The SVD is a natural generalization of the eigenvector decomposition to arbitrary (even rectangular) matrices. It plays a fundamental role in linear inverse problems.

The eigenvalue problem for AAT is easy; since it is diagonal, its diagonal entries are the eigenvalues. To ﬁnd the eigenvalues of AT A we need to ﬁnd the roots of the characteristic polynomial 1−λ 1 0 1 1−λ 0 = (1 − λ) (1 − λ)2 − 1 = 0 Det 0 0 1−λ which are 2, 1 and 0. Now we can compute the data eigenvectors ui by solving the eigenvalue problem
1

60

SVD and Resolution in Least Squares

AAT ui = λ2 ui i for λ2 equal to 2 and 1. So i 2 0 0 1 The only way this can be true is if u1 = Similarly, for λ2 = 1 we have i u2 = In this example, there is no data null space: Ur = U = 1 0 0 1 0 1 . 1 0 . u11 u21 =2 u11 u21 .

We could also solve the eigenvalue problem for AT A to get the model eigenvectors vi , but a shortcut is to take advantage of the coupling of the model and data eigenvectors, namely that AT Ur = Vr Λr , so all we have to do is take the inner product of AT with the data eigenvectors and divide by the corresponding singular value. But remember, the singular value is the square root of λ2 , so 1 0 1   v1 = √  1 0  2 0 1 and 1 0   v2 =  1 0  0 1 Vr =  
     

Remember, the only way that there can be no null space at all (no U0 or V0 ) is if n = m = r.

5.0.2

The Generalized Inverse

Recall the SVD of A is A = U ΛV T . U is n × n, V is m × m and Λ is n × m. If there are no zero singular values the following matrix provides a one-sided inverse of A: A† = V Λ−1 U T where Λ−1 refers to the m × n matrix with 1/λi on its main diagonal. The matrix A† is called the generalized inverse of A, or the pseudo-inverse. Be careful to keep the dimensions straight; in the SVD A = U ΛV T we know that V must be m × m (its columns span model space) and U must be n × n (its columns span data space). Therefore Λ must be n × m. Similarly if we write V Λ−1 U T it is clear that Λ−1 must refer to an m × n matrix.
a

Whether A† will be a left inverse or a right inverse depends on whether there are more equations than unknowns (n ≥ m) or fewer (m ≥ n). There is a two-sided (ordinary) inverse if and only if m = n = r, where r is the rank. To see how this goes consider a concrete case, m = 3 and n = r = 2 So Λ= λ1 0 0 0 λ2 0 n×m

a For this reason perhaps it is an abuse of notation to write Λ−1 . Perhaps we should write Λ† instead. The main danger of the current notation is that one is tempted to assume that Λ−1 Λ = ΛΛ−1 = I, which, as we have seen is not true in general.

V Λ−1 ΛV T = I. On the other hand if we multiply A on the right by A† we get AA† = U ΛΛ−1 U T . And ΛΛ−1 = λ1 0 0 0 λ2 0


So in this case we can see that A† is a right inverse but not a left inverse. You can verify for yourself that if there were more unknowns than data (n ≥ m), A† would be a left inverse of A. If there are zero singular values, then the only thing diﬀerent we must do is project out those components. The SVD then becomes: A = Ur Λr VrT . The generalized inverse is then deﬁned to be
T A† ≡ Vr Λ−1 Ur . r

1/λ1 0  0 1/λ2  =   0 0



1 0 0 1

.

Note that in this case Λr is an r × r matrix so and

A† A = Vr VrT

T AA† = Ur Ur .

The ﬁrst of these is an identity matrix only if r = m and the second only if r = n. You will show in an exercise however that in any case A† AA† = A† AA† A = A Let us explore the signiﬁcance of the generalized inverse bit by bit. This discussion is patterned on that in Chapter 12 of [AR80].
1

63 No Null Space First consider the case in which there is no data or model null space. This can only happen when r = m = n, in which case the generalized inverse is the ordinary inverse. A Data Null Space Next consider the case in which there is a data null space U0 but no model null space T (n > m). Since AT U0 = 0, it follows that U0 A = 0. And hence, the forward operator A always maps models into vectors that have no component in U0 . That means that if there is a data null space, and if the data have a component in this null space, then it will be impossible to ﬁt them exactly. That being the case, it would seem reasonable to try to minimize the misﬁt between observed and predicted data, say, min Am − d 2 , (5.1)

where m is an element of the model space. I.e., least-squares. A least-squares minimizing model must be associated with a critical point of this mis-ﬁt function. Diﬀerentiating Equation 5.1 with respect to m and setting the result equal to zero results in the normal equations: AT Am = AT d. (5.2) There are many ways to derive the normal equations. In the next section we will derive them without using any calculus. But it is not too hard to do the diﬀerentiation in Equation 5.1. First, write the norm-squared as an inner product: Am − d
2

= (Am − d, Am − d).

Expand this. You’ll get a sum of 4 inner products, such as (Am, Am). You can diﬀerentiate these with respect to each of the components of m if you like, but you can do this in vector notation with a little practice. For instance, the derivative of (m, m) with respect to m is 2m. The derivative of (m, a) (which equals (a, m)) with respect to m is a. Further, since (Am, Am) = (AT Am, m) = (m, AT Am), the derivative of (Am, Am) with respect to m is 2AT Am. You can move A back and forth across the inner product just by taking the transpose. Now, by Equation 4.89
T AT A = (Ur Λr VrT )T Ur Λr VrT = Vr ΛT Ur Ur Λr VrT . r

At this point we have to be a bit careful. We can be sure that U U T = U T U = In , an n-dimensional identity. And that V V T = V T V = Im . But this is not true of Vr if there
1

64

SVD and Resolution in Least Squares

is a V0 space, or Ur if there is a U0 space. All we can be certain of is that VrT Vr and T Ur Ur will be r−dimensional identity matrices, So we do know that AT A = Vr Λ2 VrT . r AT A is certainly invertible (since in this case there is assumed to be no model null space) so the least squares solution is
T mls = (Vr Λ2 VrT )−1 (Ur Λr VrT )T d = Vr Λ−1 Ur d. r r

But this is precisely A† d. Let us denote the generalized inverse solution by m† = A† d. In the special case that there is no model null space V0 , mls = m† .b Now we saw above that A maps arbitrary model vectors m into vectors that have no component in U0 . On the other hand it is easy to show (using the SVD) that
T T T T Ur (d − Am† ) = Ur d − Ur Ur Ur d = 0.

This means that Am† (since it lies in Ur ) must be perpendicular to d − Am† (since it lies in U0 ). A Geometrical Interpretation of Least Squares [Str88] If d were in the column space of A, then there would exist a vector m such that Am = d. On the other hand, if d is not in the column space of A a reasonable strategy is to try to ﬁnd an approximate solution from within the column space. In other words, ﬁnd a linear combination of the columns of A that is as close as possible in a least squares sense to the data. Let’s call this approximate solution mls . Since Amls is, by deﬁnition, conﬁned to the column space of A then Amls − d (the error in ﬁtting the data) must be in the orthogonal complement of the column space. (The orthogonal complement was deﬁned on page 47.) The orthogonal complement of the column space is the left null space, so Amls − d must get mapped into zero by AT : AT Amls − d = 0 or AT Amls = AT d which is just the normal equation again. Now we saw in the last chapter that the outer product of a vector or matrix with itself deﬁned a projection operator onto the subspace spanned by the vector (or columns of the matrix). If we look again at the normal equations and assume for the moment that the matrix AT A is invertible, then the least squares solution is: mls = (AT A)−1 AT d
b † m is the generalized inverse solution, A† d. It turns out this is unique, as we will prove shortly. mls is any solution of the normal equations. The complete connection between these two concepts will be made shortly when we treat the case of a model null space.

1

65 Now A applied to the least squares solution is the approximation to the data from within the column space. So Amls is precisely the projection of the data d onto the column space: Amls = A(AT A)−1 AT d. Before when we did orthogonal projections, the projecting vectors/matrices were orthogonal, so the AT A term would have been the identity, but the outer product structure in Amls is evident. The generalized inverse projects the data onto the column space of A. A few observations: • When A is invertible (square, full rank) A(AT A)−1 AT = AA−1 (AT )−1 AT = I, so every vector projects onto itself. • AT A has the same null space as A. Proof: clearly if Am = 0, then AT Am = 0. Going the other way, suppose AT Am = 0. Then mT AT Am = 0. But this can also be written as (Am, Am) = Am 2 = 0. By the properties of the norm, Am 2 = 0 ⇒ Am = 0. • As a corollary of this, if A has linearly independent columns (i.e., the rank r = m) then AT A is invertible.

A Model Null Space Now let us consider the existence of a model null space V0 (but no data null space U0 ), so m > n ≥ r. Once again, using the SVD, we can show that (since m† = A† d)
T Am† = AA† d = Ur Λr VrT Vr Λ−1 Ur d = d r T since VrT Vr = Ir and Ur Ur = Ir = In . But since m† is expressible in terms of the Vr vectors (and not the V0 vectors), it is clear that the generalized inverse solution is a model that satisﬁes Am† = d but is entirely conﬁned to Vr .

A consequence of this is that an arbitrary least squares solution (i.e., any solution of the normal equations) can be represented as the sum of the generalized solution with some component in the model null space:
M

mls = m† +
1

α i vi
i=r+1

(5.3)

66

SVD and Resolution in Least Squares

where by mls we mean any solution of the normal equations. An immediate consequence of this is that the length of mls must be at least as great as the length of m† since mls
2

= m†

M 2

+
i=r+1

2 αi .

(5.4)

To prove this just remember that mls 2 is the dot product of mls with itself. Take the dot product of the right-hand-side of Equation 5.3 with itself. Not only are the vectors vi mutually orthonormal, but they are orthogonal to m† since m† lives in Vr and Vr is orthogonal to V0 . This is referred to the minimum norm property of the generalized inverse. Of all the inﬁnity of solutions of the normal equations (assuming there is a model null space), the generalized inverse solution is the one of smallest length. Both a Model and a Data Null Space In the case of a data null space, we saw that the generalized inverse solution minimized the least squares mis-ﬁt of data and model response. While in the case of a model null space, the generalized inverse solution minimized the length of the solution itself. If there are both model and data null spaces, then the generalized inverse simultaneously optimizes these goals. As an exercise, set the derivative of Am − d
2

+ m

2

with respect to m equal to zero. The calculation is sketched on page 63. You should get the following generalization of the normal equations: AT A + I m = AT d. You can show that the matrix AT A + I is invertible for any A. How?

5.0.3

Examples

Consider the linear system


1 1 0 0 0 1 From the SVD we have

m1    m2  = m3

1 2 1 2



1 1

.

0  0 . A† =   0 1
1



67 It is obvious that m3 = 1 and that there are not enough equations to specify m1 or m2 . All we can say at this point is that m1 + m2 = 1. Some possible solutions then are: m1 = 0, m2 = 1, m1 = 1, m2 = 0, m1 = .5, m2 = .5, and so on. All of these choices explain the “data”. The generalized inverse solution is A† d = (1/2, 1/2, 1)T . Here we see the key feature of least squares (or generalized inverses): when faced with uncertainty least squares splits the diﬀerence.

5.0.4

Resolution

Resolution is all about how precisely one can infer model parameters from data. The issue is complicated by all of the uncertainties that exist in any inverse problem: uncertainties in the forward modeling, the discretization of the model itself (i.e., replacing continuous functions by ﬁnite-dimensional vectors), noise in the data, and uncertainties in the constraints or a priori information we have. This is why we need a fairly elaborate statistical machinery to tackle such problems. However, there are situations in which resolution becomes relatively straightforward–whether these situations pertain in practice is another matter. One of these occurs when the problem is linear and the only uncertainties arise from random noise in the data. In this case the true Earth model is linearly related to the observed data by d = Am+e where e is an n-dimensional vector of random errors. The meaning of this equation is as follows: if there were no random noise in the problem, e would be zero and the true Earth model would predict the data exactly (d = Am). We could then estimate the true model by applying the pseudo-inverse of A to the measurements. On the other hand, if e is nonzero, d = Am + e, we still get the generalized inverse solution by applying the pseudo-inverse to the data: m† = A† d. It follows that m† = A† (Am + e) . (5.5) Later on we will discuss the error term explicitly. For now we can ﬁnesse the issue by assuming that the errors have zero mean, in which case if we simply take the average of Equation 5.5 the error term goes away.c For now let’s simply assume that the errors are zero m† = A† Am. (5.6) This result can be interpreted as saying that the matrix A† A acts as a kind of ﬁlter relating the true Earth model to the computed Earth model. Thus, if A† A were equal
c

to the identity matrix, we would have perfect resolution. Using the SVD we have
T m† = Vr Λ−1 Ur Ur Λr VrT m = Vr VrT m. r T We can use Ur Ur = I whether there is a data null space or not. So in any case the matrix Vr VrT is the “ﬁlter” relating the computed Earth parameters to the true ones. In the example above, with 1 1 0 A= 0 0 1

This says that the model parameter m3 is perfectly well resolved, but that we can only resolve the average of the ﬁrst two parameters m1 and m2 . The more nonzero terms that appear in the rows of the resolution matrix, the more broadly averaged our inferences of the model parameters. Data resolution is connected to the fact that the observed data may be diﬀerent than the data predicted by the generalized inverse. The latter is just Am† . But this is AA† d. So if we call this d† , then we have a relation very similar to that given by the resolution matrix: T T d† = AA† d = Ur Λr VrT Vr Λ−1 Ur d = Ur Ur d r
T so we can think of the matrix Ur Ur as telling us about how well the data are predicted by the computed model. In our example above, there is no data null space, so the data T are predicted perfectly. But if there is a data null space then the row vectors of Ur Ur will represent averages of the data. Exercises

1. Verify the following two “Penrose conditions”: A† AA† = A† AA† A = A 2. Show that minimizing Am − d
2

+λ m

2

with respect to m leads to the following generalized “normal equations” AT A + λI m = AT d. 3. Show that AT A + λI is always an invertible matrix.
1

Chapter 6 A Summary of Probability and Statistics
Collected here are a few basic deﬁnitions and ideas. Fore more details consult a textbook on probability or mathematical statistics, for instance [Sin91], [Par60], and [Bru65]. We begin with a discussion of discrete probabilites, which involves counting sets. Then we introduce the notion of a numerical valued random variable and random physical phenomena. The main goal of this chapter is to develope the tools we need to characterize the uncertainties in both geophysical data sets and in the description of Earth models—at least when we have our Bayesian hats on. So the chapter culminates with a discussion of various descriptive statistics numerical data (means, variances, etc). Most of the problems that we face in geophysics involves spatially or temporally varying random phenomena, also known as stochastic processes; e.g., velocity as a function of space. For everything we will do in this course, however, it suﬃces to consider only ﬁnite dimensional vector-valued random variables, such as we would measure by sampling a random function at discrete times or spatial locations.

6.1

Sets

Probability is fundamentally about measuring sets. The sets can be ﬁnite, as in the possible outcomes of a toss of a coin, or inﬁnite, as in the possible values of a measured P-wave impedance. The space of all possible outcomes of a given experiment is called the sample space. We will usually denote the sample space by Ω. If the problem is simple enough that we can enumerate all possible outcomes, then assigning probabilities is easy.
0

where we use N (A) to denote the number of possible ways of achieving event A and N (Ω) is the size of the sample space.

Figure 6.1: Examples of the intersection, union, and complement of sets.

More on Sets

P (A) =

A Summary of Probability and Statistics

N (A) 1 = N (Ω) 6

0

B

C

(6.1)

6.1 Sets

73

The set with no elements in it is called the empty set and is denoted ∅. Its probability is always 0 P (∅) = 0. (6.2) Since, by deﬁnition, the sample space contains all possible outcomes, its probability must always be 1 P (Ω) = 1. (6.3) The other thing we need to be able to do is combine probabilities:a P (A ∪ B) = P (A) + P (B) − P (AB). (6.4)

P (AB) is the probability of the event A intersect B, which means the event A and B. In particular, if the two events are exclusive, i.e., if AB = 0 then P (A ∪ B) = P (A) + P (B). This result extends to an arbitrary number of exclusive events Ai
n

(6.5)

P (A1 ∪ A2 ∪ · · · ∪ An ) =

P (Ai ).
i=1

(6.6)

This property is called additivity. Events A and B are said to be independent if P (AB) = P (A)P (B). Example 2 Toss the fair die twice. The sample space for this experiment consists of {{1, 1}, {1, 2}, . . . {6, 5}, {6, 6}}. (6.7)

Let A be the event that the ﬁrst number is a 1. Let B be the event that the second number is a 2. The the probability of both A and B occurring is the probability of the intersection of these two sets. So P (AB) = N (AB) 1 = N (Ω) 36 (6.8)

Example 3 A certain roulette wheel has 4 numbers on it: 1, 2, 3, and 4. The even numbers are white and the odd numbers are black. The sample space associated with spinning the wheel twice is {{1, 1}, {1, 2}, {1, 3}, {1, 4}}
a

Now, in terms of black and white, the diﬀerent outcomes are {{black, black}, {black, white}, {black, black}, {black, white}} {{white, black}, {white, white}, {white, black}, {white, white}} {{white, black}, {white, white}, {white, black}, {white, white}} Let A be the event that the ﬁrst number is white, and B the event that the second number is white. Then N (A) = 8 and N (B) = 8. So P (A) = 8/16 and P (B) = 8/16. The event that both numbers are white is the intersection of A and B and P (AB) = 4/16. Suppose we want to know the probability of the second number being white given that the ﬁrst number is white. We denote this conditional probability by P (B|A). The only way for this conditional event to be true if both B and A are true. Therefore, P (B|A) is going to have to be equal to N (AB) divided by something. That something cannot be N (Ω) since only half of these have a white number in the ﬁrst slot, so we must divide by N (A) since these are the only events for which the event B given A could possibly be true. Therefore we have P (B|A) = N (AB) P (AB) = N (A) P (A) (6.9) {{black, black}, {black, white}, {black, black}, {black, white}}

assuming P (A) is not zero, of course. The latter equality holds because we can divide the top and the bottom of N (AB) by N (Ω). N (A) As we saw above, for independent events P (AB) = P (A)P (B). Therefore it follows that for independent events P (B|A) = P (B).

6.2

Random Variables

If we use a variable to denote the outcome of a random trial, then we call this a random variable. For example, let d denote the outcome of a ﬂip of a fair coin. Then d is a random variable with two possible values, heads and tails. A given outcome of a random trial is called a realization. Thus if we ﬂip the coin 100 times, the result is 100 realizations of the random variable d. Later in this book we will ﬁnd it necessary to invent a new notation so as to distinguish a realization of a random process from the random process itself, the later being usually unknown.
0

6.2 Random Variables

75

6.2.1

A Deﬁnition of Random

It turns out to be diﬃcult to give a precise mathematical deﬁnition of randomness, so we won’t try. (A brief perusal of randomness in Volume 2 of Knuth’s great The Art of Computer Programming is edifying and frustrating in equal measures.) In any case it is undoubtedly more satisfying to think in terms of observations of physical experiments. Here is Parzen’s (1960) deﬁnition, which is as good as any: A random (or chance) phenomenon is an empirical phenomenon characterized by the property that its observation under a given set of circumstances does not always lead to the same observed outcomes (so that there is no deterministic regularity) but rather to diﬀerent outcomes in such a way that there is statistical regularity. By this is meant that numbers exist between 0 and 1 that represent the relative frequency with which the diﬀerent possible outcomes may be observed in a series of observations of independent occurrences of the phenomenon. ... A random event is one whose relative frequency of occurrence, in a very long sequence of observations of randomly selected situations in which the event may occur, approaches a stable limit value as the number of observations is increased to inﬁnity; the limit value of the relative frequency is called the probability of the random event It is precisely this lack of deterministic reproducibility that allows us to reduce random noise by averaging over many repetitions of the experiment.

where a and b are constants and m is called the modulus. (E.g., 24 = 12 (mod 12).) The value at the step n is determined by the value and step n − 1. The modulus deﬁnes the maximum period of the sequence; but the multiplier a and the shift b must be properly chosen in order that the sequence generate all possible integers between 0 and m − 1. For badly chosen values of these constants there will be hidden periodicities which show up when plotting groups of k of these numbers as points in k-dimensional space. To implement Equation 6.10 We need four magic numbers: • m, the modulus m > 0
0

76 • a, the multiplier 0 ≤ a < m • c, the increment 0 ≤ c < m

A Summary of Probability and Statistics

• X(0), the starting value (seed) 0 ≤ X(0) < m. Here’s an example. Take X(0) = a = c = 7 and take m = 10. Then the sequence of numbers generated by our recursion is: 7, 6, 9, 0, 7, 6, 9, 0, ... Woops. The sequence begins repeating on the fourth iteration. This is called a periodicity. Fortunately, for better choices of the magic numbers, we can generate sequences that do not repeat until n is quite large. Perhaps as large as 264 or larger. But in all cases, such linear congruential “random number” generators are periodic. A large portion of Volume II of Knuth’s treatise on computer programming [Knu81] is devoted to computer tests of randomness and to theoretical deﬁntions. We will not discuss here how good choices of the magic numbers are made, except to quote Theorem A, from section 3.2.1.2, volume 2 of [Knu81]. Theorem 8 The linear congruential sequence deﬁned by m, a, c, and X(0) has period length m if and only if • c is relatively prime to m [i.e., if the greatest common divisor of c and m is 1] • b = a − 1 is a multiple of p, for every prime p dividing m • b is a multiple of 4, if m is a multiple of 4. In any case, the key point is that when you use a typical random number generator, you are tapping into a ﬁnite sequence of numbers. The place where you jump into the queue is called the seed. The sequence is purely deterministic (being generated by recursion), and we must rely on some other analysis to see whether or not the numbers really do look “random.” This led to a famous quote by one of the fathers of modern computing: Anyone who considers arithmetical methods of producing random digits is, of course, in a state of sin. – John von Neumann (1951) Any deterministic number generator, such as the recursion above, which is designed to mimic a random process (i.e., to produce a sequence of numbers with no apparent pattern) is called a pseudo-random number generator (PRNG). Virtually all of the
0

6.2 Random Variables

77

so-called random number generators used today are in fact pseudo-random number generators. This may seem like an minor point until you consider that many important numerical algorithms (so-called Monte Carlo methods) rely fundamentally on being able to make prodigious use of random numbers. Unsuspected periodicities in pseudorandom number generators have led to several of important physics papers being called into question. Can we do better than PRNG? Yes and no. Yes we can, but not practically at the scale required by modern computing. To see a pretty good random process, watch the lottery drawing on TV some time. State lotteries usually make use of the classic “well stirred urn”. How about aiming an arrow at a distant target. Surely the spread of the arrows is random. In fact, one of the words we use for random phenomena (stochastic) comes from a Greek word which means to aim. Or tune your radio to an unallocated channel and listen to the static. Amplify this static, feed it through an A/D board into your computer and voila: random numbers. Suspend a small particle (a grain of pollen for instance) in a dish of water and watch the particle under a microscope. (This is called Brownian motion, after Robert Brown, the 19th century Scottish botanist who discovered it.) Turn on a Geiger counter (after Hans Geiger, the 20th century German physicist) near some source of radioactivity and measure the time between clicks. Put a small marker in a turbulent ﬂuid ﬂow, then measure the position of the marker. The ﬂuctuations in sea height obtained from satellite altimetry. There are countless similarly unpredictable phenomena, but the question is: could we turn any of these into useful RNG? It’s been tried (see Knuth, 1973). But the appetite of modern physical means of computing random numbers may not be able to keep up with the voracious needs of Monte Carlo computer codes. Further, we would have to store all the numbers in a big table so that we could have them to reuse, else we couldn’t debug our codes (they wouldn’t be repeatable). Under Linux, true random numbers are available by reading /dev/random. This generator gathers environmenal noise (from hardware devices, such as mouse movements, keystrokes, etc.) and keeps track of how much disorder (entropy) is available. When the entropy pool is empty, further reads to /dev/random will be blocked. (For more information see the man page for random.) This means that the number of strong random numbers is limited; it may be inadequate for numerical simulations, such as Monte Carlo calculuations, that require vast quantities of such numbers. For a more extensive discussion of “true radom numbers” see the web page www.random.org, from which you can download true random numers or surf to other sites that provide them.

6.4 Probability Functions and Densities These are known as Bayes’ theorem.

79

Suppose we have n mutually exclusive and exhaustive events Ci . By mutually exclusive we meant that the intersection of any two of the Ci is the empty set (the set with no elements) Ci Cj = ∅. (6.16)

By exhaustive we meant that the union of all the Ci ﬁlls up the entire sample space (i.e., the certain event) C1 ∪ C2 ∪ · · · ∪ Cn = Ω. (6.17) It is not diﬃcult to see that for any event B, we have P (B) = P (BC1 ) + P (BC2 ) + · · · P (BCn ). (6.18)

You can think of this as being akin to writing a vector as the sum of its projections onto orthogonal (independent) directions (sets). Since the Ci are independent and exhaustive, every element in B must be in one of the intersections BCi ; and no element can appear in more than one. Therefore B = BC1 ∪ · · · BCn , and the result follows from the additivity of probabilities. Finally, since we know that for any Ci P (BCi ) = P (B|Ci )P (Ci ) it follows that P (B) = P (B|C1 )P (C1 ) + P (B|C2 )P (C2 ) + · · · P (B|Cn )P (Cn ). This gives us the following generalization of Bayes’ Theorem P (Ci |B) = P (BCi ) = P (B)
n j=1

(6.19)

P (B|Ci )P (Ci ) . P (B|Cj )P (Cj )

(6.20)

Thomas Bayes (1702-1761) is best known for his theory of probability outlined in his Essays towards solving a problem in the doctrine of chances published in the Philosophical transactions of the Royal Society (1763). He wrote a number of other mathematical essays but none were published during his lifetime. Bayes was a nonconcormist minister who preached at the Presbyterian Chapel in Turbridge Wells (south of London) for over 30 years. He was elected a fellow of the Royal Society in 1742.

6.4

Probability Functions and Densities

So far in this chapter we have dealt only with discrete probabilities. The sample space Ω has consisted of individual events to which we can assign probabilities. We can assign probabilities to collections of events by using the rules for the union, intersection and complement of events. So the probability is a kind of measure on sets. 1) It’s
0

80

A Summary of Probability and Statistics

Figure 6.2: The title of Bayes’ article, published posthumously in the Philosophical Transactions of the Royal Society, Volume 53, pages 370–418, 1763

These deﬁning properties of a probability function are already well known to us in other contexts. For example, consider a function M which measures the length of a subinterval of the unit interval I = [0, 1]. If 0 ≤ x1 ≤ x2 ≤ 1, then A = [x1 , x2 ] is a subinterval of I. Then M (A) = x2 − x1 is always positive unless the interval is empty, x1 = x2 , in which case it’s zero. If A = I, then M (A) = 1. And if two intervals are disjoint, the measure (length) of the union of the two intervals is the sum of the length of the individual intervals. So it looks like a probability function is just a special kind of measure, a measure normalized to one. Now let’s get fancy and write the length of an interval as the integral of some function over the interval. x2 M (A) = µ(x) dx ≡ µ(x) dx. (6.22)
A x1

In this simple example using cartesian coordinates, the density function µ is equal to a constant one. But it suggests that more generally we can deﬁne a probability density such that the probability of a given set is the integral of the probability density over that set P (A) = ρ(x) dx (6.23)
A

or, more generally, P (A) =
A

ρ(x1 , x2 , . . . , xn ) dx1 dx2 · · · dxn .

(6.24)

Of course there is no reason to restrict ourselves to cartesian coordinates. The set itself is independent of the coordinates used and we can transform from one coordinate system to another via the usual rules for a change of variables in deﬁnite integrals. Yet another representation of the probability law of a numerical valued random phenomenon is in terms of the distribution function F (x). F (x) is deﬁned as the probability that the observed value of the random variable will be less than x: F (x) = P (X < x) =
x −∞

ρ(x ) dx .

(6.25)

Clearly, F must go to zero as x goes to −∞ and it must go to one as x goes to +∞. Further, F (x) = ρ(x). Example
∞ −∞

Clearly, this probability is positive; it is normalized to one 1 P (−∞ ≤ x ≤ ∞) = √ π the probability of an empty interval is zero.
∞ −∞

e−x dx = 1;

2

(6.29)

6.4.1

Expectation of a Function With Respect to a Probability Law

Henceforth, we shall be interested primarily in numerical valued random phenomena; phenomena whose outcomes are real numbers. A probability law for such a phenomena P , can be thought of as determining a (in general non-uniform) distribution of a unit mass along the real line. This extends immediately to vector ﬁelds of numerical valued random phenomena, or even functions. Let ρ(x) be the probability density associated with P , then we deﬁne the expectation of a function f (x) with respect to P as E[f (x)] =
∞ −∞

f (x)ρ(x) dx.

(6.30)

Obviously this expectation exists if and only if the improper integral converges. The mean of the probability P is the expectation of x E[x] =
∞ −∞

xρ(x) dx.

(6.31)

For any real number ξ, we deﬁne the n-th moment of P about ξ as E[(x − ξ)n ]. The most common moments are the central moments, which correspond to E[(x − x)n ]. The ¯ second central moment is called the variance of the probability law. Keep in mind the connection between the ordinary variable x and the random variable itself; let us call the latter X. Then the probability law P and the probability density p are related by P (X < x) =
x −∞

ρ(x ) dx .

(6.32)

We will summarize the basic results on expectations and variances later in this chapter.
0

1 where I2 is a subset of the real plane R2 . So ρ(x, y) = π e−(x +y ) is an example of a joint probability density on two variables. We can extend this deﬁnition to any number of variables. In general, we denote by ρ(x1 , x2 , ...xN ) a joint density for the N -dimensional random variable. Sometimes we will write this as ρ(x). By deﬁnition, the probability that the N −vector x lies in some subset A of RN is given by:

P [x ∈ A] =

A

ρ(x) dx

(6.35)

where dx refers to some N −dimensional volume element. independence We saw above that conditional probabilities were related to joint probabilities by P (AB) = P (B|A)P (A) = P (A|B)P (B) from which result Bayes theorem follows. The same result holds for random variables. If a random variable X is independent of event Y , the probability of Y does not depend on the probability of X. That is, P (x|y) = P (x) and P (y|x) = P (y). Hence for independent events, P (x, y) = P (x)P y). Once we have two random variables X and Y , with a joint probability P (x, y), we can think of their moments. The joint n − m moment of X and Y about 0 is just E[xn y m ] = xn y m ρ(x, y)dxdy. (6.36)

The 1-1 moment about zero is called the correlation of the two random variables: ΓXY = E[xy] = xyρ(x, y)dxdy. (6.37)

On the other hand, the 1-1 moment about the means is called the covariance: CXY = E[(x − E(x))(y − E(y)) = ΓXY − E[x]E[y].
0

(6.38)

84

A Summary of Probability and Statistics

The covariance and the correlation measure how similar the two random variables are. This similarity is distilled into a dimensionless number called the correlation coeﬃcient: r= CXY σX σY (6.39)

where σX is the variance of X and σY is the variance of Y . Using Schwarz’s inequality, namely that
2

This proves that 0 ≤ r ≤ 1. A correlation coeﬃcient of 1 means that the ﬂuctuations in X and Y are essentially identical. This is perfect correlation. A correlation coeﬃcient of -1 means that the ﬂuctuations in X and Y are essentially identical but with the opposite sign. This is perfect anticorrelation. A correlation coeﬃcient of 0 means X and Y are uncorrelated. Two independent random variables are always uncorrelated. But dependent random variables can be uncorrelated too. Here is an example from [Goo00]. Let Θ be uniformly distributed on [−π/2, π/2]. Let X = cos Θ and Y = sin Θ. Since knowledge of Y completely determines X, these two random variables are clearly dependent. But CXY = . marginal probabilities From an n−dimensional joint distribution, we often wish to know the probability that some subset of the variables take on certain values. These are called marginal probabilities. For example, from ρ(x, y), we might wish to know the probability P (x ∈ I1 ). To ﬁnd this all we have to do is integrate out the contribution from y. In other words P (x ∈ I1 ) = 1 π
∞ I1 −∞

From this deﬁnition it is obvious that C is a symmetric matrix. The diagonal elements of the covariance matrix are just the ordinary variances (squares of the standard deviations) of the components: Cii (m) = (σi )2 . (6.46) The oﬀ-diagonal elements describe the dependence of pairs of components. As a concrete example, the n-dimensional normalized gaussian distribution with mean m and covariance C is given by ρ(x) = 1 1 exp − (x − m)T C −1 (x − m) . N/2 (2π det C) 2 (6.47)

This result and many other analytic calculations involving multi-dimensional Gaussian distributions can be found in [MGB74] and [Goo00]. An aside, diagonalizing the covariance Since the covariance matrix is symmetric, we can always diagonalize it with an orthogonal transformation involving real eigenvalues. It we transform to principal coordinates (i.e., rotate the coordinates using the diagonalizing orthogonal transformation) then the covariance matrix becomes diagonal. So in these coordinates correlations vanish since they are governed by the oﬀ-diagonal elements of the covariance matrix. But suppose one or more of the eigenvalues is zero. This means that the standard deviation of that parameter is zero; i.e., our knowledge of this parameter is certain. Another way to say this is that one or more of the parameters is deterministically related to the others. This is not a problem since we can always eliminate such parameters from the probabilistic description of the problems. Finally, after diagonalizing C we can scale the parameters by their respective standard deviations. In this new rotated, scaled coordinate system the covariance matrix is just the identity. In this sense, we can assume in a theoretical analysis that the covariance is the identity since in priciple we can arrange so that it is.
0

86

A Summary of Probability and Statistics

6.5

Random Sequences

Often we are faced with a number of measurements {xi } that we want to use to estimate the quantity being measured x. A seismometer recording ambient noise, for example, is sampling the velocity or displacement as a function of time associated with some piece of the earth. We don’t necessarily know the probability law associated with the underlying random process, we only know its sampled values. Fortunately, measures such as the mean and standard deviation computed from the sampled values converge to the true mean and standard deviation of the random process. The sample average or sample mean is deﬁned to be 1 x≡ ¯ N
N

xi
i=1

sample moments: Let x1 , x2 , ...xN be a random random sample from the probability density ρ. Then the r-th sample moment about 0, is given by 1 N
N

xr . i
i=1

If r = 1 this is the sample mean, x. Further, the r-th sample moment about x, is ¯ ¯ given by 1 N Mr ≡ (xi − x)r . ¯ N i=1 How is the sample mean x related to the mean of the underlying random variable ¯ (what we will shortly call the expectation, E[X])? This is the content of the law of large numbers; here is one form due to Khintchine (see [Bru65] or [Par60]): Theorem 9 Khintchine’s Theorem: If x is the sample mean of a random sample of ¯ size n from the population induced by a random variable x with mean µ, and if > 0 then: P [ x − µ ≥ ] → 0 as n → ∞. ¯ In the technical language of probability the sample mean x is said to converge in prob¯ ability to the population mean µ. The sample mean is said to be an “estimator” of true mean. A related result is Theorem 10 Chebyshev’s inequality: If a random variable X has ﬁnite mean x and ¯
0

6.5 Random Sequences

87

-5

-4

-3

-2

-1

0

1

2

3

4

5

Figure 6.4: A normal distribution of zero mean and unit variance. Almost all the area under this curve is contained within 3 standard deviations of the mean. variance σ 2 , then [Par60]: P[ X − x ≤ ] ≥ 1 − ¯ for any > 0. σ2
2

This says that in order to make the probability of X being within of the mean greater σ than some value p, we must choose at least as large as √1−p . Another way to say this would be: let p be the probability that x lies within a distance of the mean x. Then ¯ σ Chebyshev’s inequality says that we must choose to be at least as large as √1−p . For example, if p = .95, then ≥ 4.47σ, while for p = .99, then ≥ 10σ. For the normal probability, this inequality can be sharpened considerably: the 99% conﬁdence interval is = 2.58σ. But you can see this in the plot of the normal probability in Figure 6.4. This is the standard normal probability (zero mean, unit variance). Clearly nearly all the probability is within 3 standard deviations.

6.5.1

The Central Limit Theorem

The other basic theorem of probability which we need for interpreting real data is this: the sum of a large number of independent, identically distributed random variables (deﬁned on page 21), all with ﬁnite means and variances, is approximately normally distributed. This is called the central limit theorem, and has been known, more or less, since the time of De Moivre in the early 18-th century. The term “central limit theorem” was coined by George Polya in the 1920s. There are many forms of this result, for proofs you should consult more advanced texts such as [Sin91] and [Bru65].
0

Figure 6.5: Ouput from the coin-ﬂipping program. The histograms show the outcomes of a calculation simulating the repeated ﬂipping of a fair coin. The histograms have been normalized by the number of trials, so what we are actually plotting is the relative probability of of ﬂipping k heads out of 100. The central limit theorem guarantees that this curve has a Gaussian shape, even though the underlying probability of the random variable is not Gaussian. Theorem 11 Central Limit Theorem: If x is the sample mean of a sample of size n ¯ from a population with mean µ and standard deviation σ, then for any real numbers a and b with a < b 1 bσ aσ ¯ P µ+ √ <x<µ+ √ → √ n n 2π
b a

e−z

2 /2

dz.

This just says that the sample mean is approximately normally distributed. Since the central limit theorem says nothing about the particular distribution involved, it must apply to even something as apparently non-Gaussian as ﬂipping a coin. Suppose we ﬂip a fair coin 100 times and record the number of heads which appear. Now, repeat the experiment a large number of times, keeping track of how many times there were 0 heads, 1, 2, and so on up to 100 heads. Obviously if the coin is fair, we expect 50 heads to be the peak of the resulting histogram. But what the central limit theorem says is that the curve will be a Gaussian centered on 50. This is illustrated in Figure 6.5 via a little code that ﬂips coins for us. For comparison, the exact probability of ﬂipping precisely 50 heads is 100! 1 50!50! 2
100

≈ .076.

(6.48)

What is the relevance of the Central Limit Theorem to real data? Here are three conﬂicting views quoted in [Bra90]. From Scarborough (1966): “The truth is that, for the kinds of errors considered in this book (errors of measurement and observation), the Normal Law is proved by experience. Several substitutes for this law have been proposed, but none ﬁts the facts as well as it does.” From Press et al. (1986): “This infatuation [of statisticians with Gaussian statistics] tended to focus interest away from the fact that, for real data, the normal distribution
0

6.6 Expectations and Variances is often rather poorly realized, if it is realized at all.”

89

And perhaps the best summary of all, Gabriel Lippmann speaking to Henri Poincar´: e “Everybody believes the [normal law] of errors: the experimenters because they believe that it can be proved by mathematics, and the mathematicians because they believe it has been established by observation.”

6.6

Expectations and Variances

Notation: we use E[x] to denote the expectation of a random variable with respect to its probability law f (x). Sometimes it is useful to write this as Ef [x] if we are dealing with several probability laws at the same time. If the probability is discrete then E[x] =
i

xi f (xi ).

If the probability is continuous then E[x] =
∞ −∞

xf (x) dx.

Mixed probabilities (partly discrete, partly continuous) can be handled in a similar way using Stieltjes integrals [Bar76]. We can also compute the expectation of functions of random variables: E[φ(x)] =
∞ −∞

φ(x)f (x) dx.

It will be left as an exercise to show that the expectation of a constant a is a (E[a] = a) and the expectation of a constant a times a random variable x is a times the expectation of x (E[ax] = aE[x]). Recall that the variance of x is deﬁned to be V (x) = E[(x − E(x))2 ] = E[(x − µ)2 ] where µ = E[x]. Here is an important result for expectations: E[(x − µ)2 ] = E[x2 ] − µ2 . The proof is easy. E[(x − µ)2 ] = = = = E[x2 − 2xµ + µ2 ] E[x2 ] − 2µE[x] + µ2 E[x2 ] − 2µ2 + µ2 E[x2 ] − µ2
0

(6.49) (6.50) (6.51) (6.52)

90

A Summary of Probability and Statistics

An important result that we need is the variance of a sample mean. For this we use the following lemma, the proof of which will be left as an exercise: Lemma 1 If a is a real number and x a random variable, then V (ax) = a2 V (x). From this it follows immediately that V (¯) = x 1 n2
n

V (xi ).
i=1

In particular, if the random variables are identically distributed, with mean µ then V (¯) = σ 2 /n. x

6.7

Bias

In statistics, the bias of an estimator of some parameter is deﬁned to be the expectation of the diﬀerence between the parameter and the estimator: ˆ ˆ B[θ] ≡ E[θ − θ] (6.53)

ˆ where θ is the estimator of θ. In a sense, we want the bias to be small so that we have a faithful estimate of the quantity of interest. ˆ ˆ An estimator θ of θ is unbiased if E[θ] = θ For instance, it follows from the law of large numbers that the sample mean is an unbiased estimator of the population mean. In symbols, E[¯] = µ. x However, the sample variance s2 ≡ 1 N
N i=1

(6.54)

(xi − x)2 ¯

(6.55)

turns out not to be unbiased (except asymptotically) since E[s2 ] = n−1 σ 2 . To get an n n unbiased estimator of the variance we use E[ n−1 s2 ]. To see this note that s2 = 1 N
n i=1

ˆ Finally, there is the notion of the consistency of an estimator. An estimator θ of θ is consistent if for every > 0 ˆ P [|θ − θ| < ] → 1 as n → ∞. Consistency just means that if the sample size is large enough, the estimator will be close to the thing being estimated. Later on, when we talk about inverse problems we will see that bias represents a potentially signiﬁcant component of the uncertainty in the results of the calculations. Since the bias depends on something we do not know, the true value of the unknown parameter, it will be necessary to use a priori information in order to estimate it. Mean-squared error, bias and variance The mean-squared error (MSE) for an estimator m of mT is deﬁned to be ¯ MSE(m) ≡ E[(m − mT )2 ] = E[m2 − 2mmT + m2 ] = m2 − 2mmT + m2 . ¯ T T By doing a similar analysis of the variance and bias we have: Bias(m) ≡ E[m − mT ] = m − mT ¯ and ¯ Var(m) ≡ E[(m − m)2 ] = E[m2 − 2mm + m2 ] = m2 − m2 . ¯ ¯ ¯ ¯
2

(6.56)

(6.57) (6.58)

So you can see that we have: MSE = Var + Bias . As you can see, for a given meansquared error, there is a trade-oﬀ between variance and bias. The following example illustrates this trade-oﬀ. Example: Estimating the derivative of a smooth function We start with a simple example to illustrate the eﬀects of noise and prior information in the performance of an estimator. Later we will introduce tools from statistical decision theory to study the performance of estimators given diﬀerent types of prior information.
0

It is clear that choosing the smallest h possible does not lead to the best estimate; the noise has to be taken into account. The lowest upper bound is obtained with h = 21/4 σ/M . The larger the variance of the noise, the wider the spacing between the points. We have used a bound on the second derivative to bound the MSE. It is a fact that some type of prior information, in addition to model (6.59), is required to bound derivative uncertainties. Take any smooth function, g, which vanishes at the points x1 , ..., xn . ˜ Then, the function f = f + g satisﬁes the same model as f , yet their derivatives could be very diﬀerent. For example, choose an integer, m, and deﬁne g(x) = sin 2πm(x − x1 ) . h

˜ By choosing m large enough, we can make the diﬀerence, f (xmi ) − f (xmi ), as large as we want; without prior information the derivative can not be estimated with ﬁnite uncertainty.

6.8

Correlation of Sequences

Many people think that “random” and uncorrelated are the same thing. Random sequences need not be uncorrelated. Correlation of sequences is measured by looking at the correlation of the sequence with itself, the autocorrelation.c If this is approximately a δ-function, then the sequence is uncorrelated. In a sense, this means that the sequence does not resemble itself for any lag other than zero. But suppose we took a deterministic function, such as sin(x), and added small (compared to 1) random perturbations to it. The result would have the large-scale structure of sin(x) but with a lot of random junk superimposed. The result is surely still random, even though it will not be uncorrelated. If the autocorrelation is not a δ-function, then the sequence is correlated. Figure 6.6 shows two pseudo-random Gaussian sequences with approximately the same mean, standard deviation and 1D distributions: they look rather diﬀerent. In the middle of this ﬁgure are shown the autocorrelations of these two sequences. Since the autocorrelation of the right-hand sequence drops oﬀ to approximately zero in 10 samples, we say the correlation length of this sequence is 10. In the special case that the autocorrelation of a sequence is an exponential function, the the correlation length is deﬁned as the (reciprocal) exponent of the best-ﬁtting exponential curve. In other words, if the autocorrelation can be ﬁt with an exponential e−z/ , then the best-ﬁtting value of is the correlation length. If the autocorrelation is not an exponential, then the correlation length is more diﬃcult to deﬁne. We could say that it is the number of lags of the autocorrelation within which the autocorrelation has most of its energy. It is often impossible to deﬁne meaningful correlation lengths from real data. A simple way to generate a correlated sequence is to take an uncorrelated one (this is what pseudo-random number generators produce) and apply some operator that correlates the samples. We could, for example, run a length- smoothing ﬁlter over the uncorrelated samples. The result would be a series with a correlation length approximately equal to . A fancier approach would be to build an analytic covariance matrix and impose it on an uncorrelated pseudo-random sample.
c From the convolution theorem, it follows that the autocorrelation is just the inverse Fourier transform of the periodogram (absolute value squared of the Fourier transform).

Figure 6.6: Two Gaussian sequences (top) with approximately the same mean, standard deviation and 1D distributions, but which look very diﬀerent. In the middle of this ﬁgure are shown the autocorrelations of these two sequences. Question: suppose we took the samples in one of these time series and sorted them in order of size. Would this preserve the nice bell-shaped curve?

0

6.8 Correlation of Sequences

95

For example, an exponential covariance matrix could be deﬁned by Ci,j = σ 2 e− i−j / where σ 2 is the (in this example) constant variance and is the correlation length. To impose this correlation on an uncorrelated, Gaussian sequence, we do a Cholesky decomposition of the covariance matrix and dot the lower triangular part into the uncorrelated sequence [Par94]. If A is a symmetric matrix, then we can always write A = LLT , where L is lower triangular [GvL83]. This is called the Cholesky decomposition of the matrix A. You can think of the Cholesky decomposition as being somewhat like the square root of the matrix. Now suppose we apply this to the covariance matrix C = LLT . Let x be a mean zero pseudo-random vector whose covariance is the identity. We will use L to transform x: z ≡ Lx. The covariance of z is given by Cov(z) = E[zz T ] = E[(Lx)(Lx)T ] = LE[xxT ]LT = LILT = C which is what we wanted. Here is a simple Scilab code that builds an exponential covariance matrix Ci,j = σ 2 e− i−j /l and then returns n pseudo-random samples drawn from a Gaussian process with this covariance (and mean zero).

The vector z should now have a standard deviation of approximately 1 and a correlation length of 10. To see the autocorrelation of the vector you could do: plot(corr(z,200)) the autocorrelation function will be approximately exponential for short lags. For longer lags you will see oscillatory departures from exponential behavior—even though the samples are drawn from an analytic exponential covariance. This is due to the small number of realizations. From the exponential part of the autocorrelation you can estimate the correlation length by simply looking for the point at which the autocorrelation has decayed by 1/e. Or you can ﬁt the log of the autocorrelation with a straight line.

6.9

Random Fields

If we were to make N measurements of, say, the density of a substance, we could plot these N data as a function of the sample number. This plot might look something like the top left curve in Figure 6.6. But these are N measurements of the same thing, unless we belive that the density of the sample is changing with time. So a histogram of the measurements would approximate the probability density function of the parameter and that would be the end of the story. On the other hand, curves such as those in Figure 6.6 might result from measuring a random, time-varying process. For instance, these might be measurements of a noisy accelerometer, in which case the plot would be N samples of voltage versus time. But then these would not be N measurements of the same parameter, rather they would constitute a single measurement of a random function of time, sampled at N times. The distinction we are making here is the distinction between a scalar-valued random process and a random function. Now when we measure a function we always measure it at a ﬁnite number of locations (in space or time). So our measurements of random function result in ﬁnite-dimensional random vectors. This is why the sampling of a time series of voltage, say, at N times, is really the realization of an N -dimensional random process. For such a process it is not suﬃcient to simply make a histogram of the samples. We need higher order characterizations of the probability law behind the time-series in order to account for the correlations of the measured values. This is the study of random ﬁelds or stochastic processes. A real-data example of measurements of a random ﬁeld is shown in Figure 6.7. These traces represent 38 realizations of a time series recorded in a random medium. In this case the randomness is spatial: each realization is made at a diﬀerent location in the medium. But you could easily imagine the same situation arising with a temporal random process. For example, these could be measurements of wave propagation in a medium undergoing random ﬂuctuations in time, such as the ocean or the atmosphere.
0

6.9 Random Fields
realization 20

97

0 0.05 time (ms) 0.10 0.15 0.20

10

30

Figure 6.7: 38 realizations of an ultrasonic wave propagation experiment in a spatially random medium. Each trace is one realization of an unknown random process U (t). The study of random processes is both conceptually and notationally diﬃcult. There are advanced mathematical books such as Pugachev’s Theory of Random Functions [Pug65] which explain the situation, but from a physical point of view, one of the clearest explanations is in the book Statistical Optics by Goodman [Goo00]. We will follow his lead. So let us assume that there is an underlying random process U (t) (or U (r) if we want to think about spatial randomness, it doesn’t matter for our purposes), the probability law of which we do not know. The notation U (t) is used to represent the ensemble of all possible outcomes of the random process along with their probabilities. Each outcome would be a function of time, say u(t). It is as if t is an index and u labels the sample description space of the random process. If we don’t know the probability law of the random process, perhaps it is possible to nevertheless characterize U by means, covariances, that sort of thing. Here is the theorem, which you can ﬁnd in [Pug65]. To completely characterize the probability law of a stochastic process we must know the joint probability distribution ρU (u1 , u2 , . . . , un , . . . ; t1 , t2 , . . . , tn , . . .) for all n, where u1 = u(t1 ), etc. Of course, in practice, we only make measurements of U at a ﬁnite numer of samples, so in practice we must make due with the n − th order ρU .d Now, the ﬁrst-order PDF (probability density function) is ρU (u; t). Knowing this
It is sometimes convenient to label the joint distribution by the random process, and sometimes by the order. So, ρU (u1 , u2 , . . . , un ; t1 , t2 , . . . , tn ) ≡ ρn (u1 , u2 , . . . , un ; t1 , t2 , . . . , tn )
d

It is relatively uncommon to go beyond second order statistical characterizations (means and covariances) and we will not do so in this class. Time versus statistical autocorrelation Given a known function u(t) (or this could be discrete samples of a known function u(ti )), the time autocorrelation function of u is deﬁned as: 1 ˜ Γu (τ ) = lim T ←∞ T
T /2 −T /2

In geophysics there is a large amount of a priori information that could be used to inﬂuence inverse calculations. Here, a priori refers to the assumption that this information is known independently of the data. Plausible geologic models can be based on rock outcrops, models of sedimentological deposition, subsidence, etc. There are also often in situ and laboratory measurements of rock properties that have a direct bearing on macroscopic seismic observations, such as porosity, permeability, crack orientation, etc. There are other, less quantitative, forms of information as well, the knowledge of experts for instance. This prior information can be deterministic or probabilistic. Examples of deterministic information include: density is positive, wave velocity is positive (and less than the
0

speed of light!), or that the mass of the planet is bounded above. Prior information can also be probabilistic, as we have discussed from a Bayesian perspective since the beginning of this course. In a later chapter we will treat problems with deterministic information using non-Bayesian methods, but for now we will consider probabilistic prior information. There is a simple conceptual model that can be used to visualize the application of this diverse a priori information. Suppose we have a black box into which we put all sorts of information about our problem. We can turn a crank or push a button and out pops a model that is consistent with whatever information that is put in, as illustrated in Figure 6.8. If this information consists of an expert’s belief that, say, a certain sand/shale sequence is 90% likely to appear in a given location, then we must ensure that 90% of the models generated by the black box have this particular sequence. One may repeat the process indeﬁnitely, producing a collection of models that have one thing in common: they are all consistent with the information put into the black box. Let us assume, for example, that a particular layered Earth description consists of the normal-incidence P-wave reﬂection coeﬃcient r at 1000 diﬀerent depth locations in some well-studied sedimentary basin. Suppose, further, that we know from in situ measurements that r in this particular sedimentary basin almost never exceeds .1. What does it mean for a model to be consistent with this information? We can push the button on the black box and generate models which satisfy this requirement. Figure 6.9 shows
0

Figure 6.9: Three models of reﬂectivity as a function of depth which are consistent with the information that the absolute value of the reﬂection coeﬃcient must be less than .1. On the right is shown the histogram of values for each model. The top two models are uncorrelated, while the bottom model has a correlation length of 15 samples. some examples. The three models shown satisfy the hard constraint that |r| ≤ .1 but they look completely diﬀerent. In the case of the two Gaussian models, we know this is because they have diﬀerent correlation length. The real question is, which is most consistent with our assumed prior information? What do we know about the correlation length in the Earth? And how do we measure this consistency? If we make a histogram of the in situ observations of r and it shows a nice bell-shaped curve, are we justiﬁed in assuming a Gaussian prior distribution? On the other hand, if we do not have histograms of r but only extreme values, so that all we really know is that |r| ≤ .1, are we justiﬁed in thinking of this information probabilistically? If we accept for the moment that our information is best described probabilistically, then a plausible strategy for solving the inverse problem would be to generate a sequence of models according to the prior information and see which ones ﬁt the data. Assuming, of course, that we know the probability function that governs the variations of the Earth’s properties. In the case of the reﬂectivity sequence, imagine that we have surface seismic data to be inverted. For each model generated by the black box, compute synthetic seismograms, compare them to the data and decide whether they ﬁt the data well enough to be acceptable. If so, the models are saved; if not, they are rejected. Repeating this procedure many times results in a collection of models that are, by deﬁnition, a priori reasonable and ﬁt the data. If the models in this collection all look alike, then the features the models show are well-constrained by the combination of data ﬁt and a priori information. If, on the other hand, the models show a diversity of features, then
0

6.11 Other Common Analytic Distributions

101

Figure 6.10: The lognormal is a prototype for asymmetrical distributions. It arises naturally when considering the product of a number of iid random variables. This ﬁgure was generated from Equation 6.62 for s = 2. these features cannot be well-resolved.

6.11

Other Common Analytic Distributions

Finally, we close this chapter by mentioning a few other commonly used analytic distributions. Nearly as important theoretically as the Gaussian, is the lognormal distribution. A variable X is lognormally distributed if Y = ln(X) is normally distributed. The central limit theorem treats the problem of sums of random variables; but the product of exponentials is the exponential of the sum of the exponents. Therefore we should expect that the lognormal would be important when dealing the a product of iid e random variables. One of these distributions, sketched in Figure 6.10, is the prototype of asymmetrical distributions. It also will play an important role later, when we talk about so-called non-informative distributions. In fact, there is a whole family of lognormal distributions given by ρ(x) = sx 2π 1 √ exp − 1 x log 2s2 x0
2

(6.62)

where x0 plays the analog of the mean of the distribution and s governs the shape. For small values of s (less than 1), the lognormal distribution is approximately gaussian. While for large values of s (greater than about 2.5), the lognormal approches 1/x. Figure 6.10 was computed for an s value of 2. The Gaussian distribution is a member of a family of exponential distributions referred to as generalized Gaussian distributions. Four of these distributions are shown in Fige

where Γ is the Gamma function [MF53] and σp is a generalized measure of variance known in the general case as the dispersion of the distribution: (σp )p ≡
∞ −∞

|x − x0 |p ρ(x) dx

(6.64)

where x0 is the center of the distribution. See [Tar87] for more details.

0

6.11 Other Common Analytic Distributions

103

Exercises
1. Show that for any two events A and B P (AB c ) = P (A) − P (BA) 2. Show that for any event A, P (Ac ) = 1 − P (A). 3. Show that 1/x is a measure, but not a probability density. 4. Show that the truth of the following formula for any two sets A and B P (A ∪ B) = P (A) + P (B) − P (A ∪ B) follows from the fact that for independent sets A and B P (A ∪ B ) = P (A ) + P (B ). (6.67) (6.66) (6.65)

Hint. The union of any two sets A and B can be written as the sum of three independent sets of elements: the elements in A but not in B; the elements in B but not in A; and the elements in both A and B. 5. Show that all the central moments of the normal distribution beyond the second are either zero or can be written in terms of the mean and variance. 6. You have made n diﬀerent measurements of the mass of an object. You want to ﬁnd the mass that best “ﬁts” the data. Show that the mass estimator which minimizes the sum of squared errors is given by the mean of the data, while the mass estimator which minimizes the sum of the absolute values of the errors is given by the median of the data. Feel free to assume that you have an odd number of data. 7. Show that Equation 6.47 is normalized. 8. Take the n data you recorded above and put them in numerical order: x1 ≤ x2 ≤ ... ≤ xn . Compute the sensitivity of the two diﬀerent estimators, average and median, to perturbations in xn . What does this say about how least squares and least absolute values treat “outliers” in the data? 9. Find the normalization constant that will make p(x) = e−(x
2 −x 2 0 x+x0 )

(6.68)

a probability density on the real line. x0 is a constant. What are the mean and variance?
0

104

A Summary of Probability and Statistics Answer: The exponential integral is ubiquitous. You should remember the following trick. H= H2 = Therefore
∞ −∞ ∞ −∞
2

Write a program that computes the sample covariance matrix of repeated recordings of a time series. To test your code, generate 25 correlated time series of length 100 and use these as data. In other words the sample size will be 25 and the covariance matrix will be of order 100.
0

Chapter 7 Linear Inverse Problems With Uncertain Data
In Chapter 5 we showed that the SVD could be used to solve a linear inverse problem in which the only uncertainties were associated with the data. The canonical formulation of such a problem is d = Am + (7.1)

where it is assumed that the forward operator A is linear and exactly known and that the uncertainties arise from additive noise in the data. In most geophysical inverse problems, the vector m is properly deﬁned in an inﬁnite dimensional space of functions; for example, the elastic tensor as a function of space. As a practical matter the model space is usually discretized so that the problem is numerically ﬁnite dimensional. This is a potential source of error (bias, discretization error), but for now we will ignore this and assume that the discretization is very ﬁne, but nevertheless ﬁnite. This is equivalent to assuming a priori that the true Earth model is conﬁned to a ﬁnite dimensional subspace of the model space. If there are no discretization errors, and if the forward model is linear and known, then the observations d are the response of the true model mT under the action of A, provided there are no measurement or other systematic errors. As deﬁned in Section 6.5 m† is a pseudo-inverse estimator of the true model: the generalized solution of Equation 7.1 is given by m† = A† d. Since d is the response of the true model, mtrue , it follows that m† ≡ A† d = A† (Amtrue + ) = A† Amtrue + A† . In terms of the SVD, the resolution matrix A† A can be written Vr VrT and so represents a projection operator onto the non-null space of the forward problem (i.e., the row space). Since none of the columns of Vr lie in the null space of A, the net result of this
1

108

Linear Inverse Problems With Uncertain Data

is that A† A can have no component in the null space. So, apart from the noise, the matrix A† A acts as a ﬁlter through which we see the Earth. We proved in Section 4.9 that a projection operator onto the null space is V0 V0T = I − VrT VrT and therefore (A† A − I)mT = − [mT ]null . (7.3) (7.2)

The null space components of the true model have a special statistical signiﬁcance. In statistics, the bias of an estimator of some parameter is deﬁned to be the expectation of the diﬀerence between the parameter and the estimator (see Section 6.7): ˆ ˆ B[θ] ≡ E[θ − θ] (7.4)

ˆ where θ is the estimator of θ. In a sense, we want the bias to be small so that we have a faithful estimate of the quantity of interest. For instance, it follows from the linearity ¯ of the expectation that the sample mean x is an unbiased estimator of the population mean µ: E[¯ ] = µ x (7.5) and hence E[¯ − µ] = 0. x Using the previous result we can see that the bias of the generalized inverse solution as an estimator of the true earth model is just (minus) the projection of the true model onto the null space of the forward problem: B(m† ) ≡ E[m† − mtrue ] = E[A† d − mtrue ] = E[A† Amtrue + A† − mtrue ] =

A† A − I E[mtrue ] + A† E[ ]

(7.6)

and so, assuming that the noise is zero mean (E[ ] = 0), we can see that the bias is simply the projection of the true model’s expected value onto the null-space. If we assume that the true model is non-random, then E[mtrue ] = mtrue . The net result is that the bias associated with the generalized inverse solution is the component of the true model in the null space of the forward problem. Inverse problems with no null-space are automatically unbiased. But the existence of a null-space does not automatically lead to bias since the true model could be orthogonal to the null-space. If the expected value of the true model is a constant, then this orthogonality is equivalent to having the row sums of the matrix A† A − I be zero. In fact, the requirement that the row sums of this matrix be zero is sometimes stated as the deﬁnition of unbiasedness [OP95], but as we have just seen, such a deﬁnition would, in general, be inconsistent with the standard statistical use of this term.
1

7.1 The World’s Second Smallest Inverse Problem

109

7.0.1

Model Covariances

Estimators are functions of the data and therefore random variables. The covariance of a random variable x is the second central moment: C = E[(x − E[x])(x − E[x])T ]. (7.7)

The covariance of the generalized inverse estimate m† ≡ A† d is easy to compute. First realize that if d has zero mean than so does m† a and therefore assuming zero mean errors T T ˆ Cov(m) = E[m† m† ] = A† Cov(d)A† (7.8) If the data are uncorrelated, then Cov(d) is a diagonal matrix whose elements are the standard deviations of the data. If we go one step further and assume that all these 2 standard deviations are the same, σd ,b then the covariance of the generalized inverse estimate takes on an especially simple form:
2 2 Cov(m† ) = σd A† A† = σd Vr Λ−2 VrT . r

We can see that the uncertainties in the estimated model parameters (expressed as Cov(δm† )) are proportional to the data uncertainties and inversely proportional to the squared singular values. This is as one would expect: as the noise increases, the uncertainty in our parameter estimates increases; and further, the parameters associated with the smallest singular values will be less well resolved than those associated with the largest.

7.1

The World’s Second Smallest Inverse Problem

Suppose we wanted to use sound to discover the depth to bedrock below our feet. We could set oﬀ a loud bang at the surface and wait to see how long it would take for the echo from the top of the bedrock to return to the surface. Then, assuming that the geologic layers are horizontal, can we compute the depth to bedrock z from the travel time of the reﬂected bang t? Suppose we do not know the speed with which sounds propagates beneath us, so that all we can say is that the travel time must depend both on this speed and on the unknown depth t = 2z/c. Since this toy problem involves many of the complications of more realistic inverse calculations, it will be useful to go through the steps of setting up and solving the calculation. We can absorb the factor of two into a new sound speed c and write t = z/c.
a b

So the model vector m is (z, c) since c and z are both unknown, and the data vector d is simply t. The forward problem is g(m) = z/c. Notice that g is linear in depth, but nonlinear in sound speed c. We can linearize the forward problem by doing a Taylor series expansion about some model (z0 , c0 ) and retaining only the ﬁrst order term: t = t0 + 1 z0 1 ,− c0 z0 c0 δz δc (7.10)

where t0 = z0 /c0 . Pulling the t0 over to the left side and diving by t0 we have δt = [1, −1] t0
δz z0 δc c0

(7.11)

In this particular case the linearization is independent of the starting model (z0 , c0 ) since by computing the total derivative of of t we get δt δz δc = − . t z c (7.12)

In other words, by deﬁning new parameters to be the logarithms of the old parameters, or the dimensionless perturbations, (but keeping the same symbols for convenience) we have t = z − c. (7.13) In any case, the linear(-ized) forward operator is the 1 × 2 matrix A = (1, −1) and AT A = 1 −1 −1 1 . (7.14)

Let’s work out the SVD of A by hand. First, let us make the convention that model vectors and data vectors are column vectors. We could make them row vectors too, but we must keep to some convention in order to avoid getting confused. So The forward operator matrix A must be a 1 by 2 matrix A = [1 − 1] since d ∈ R1 and m ∈ R2 . Therefore AT A = and AAT = [1 − 1] 1 −1 [1 − 1] = 1 −1 1 −1 −1 1 = 2.

So the eigenvalue of a 1 × 1 matrix (a scalar) is just this number. The eigenvalues of AAT are the squares of the singular values, so the one and only non-zero singular value √ is 2.
1

7.1 The World’s Second Smallest Inverse Problem

111

Now since the data space is one-dimensional, the data-space eigenvector is just a normalized vector in R1 –which is just 1. So Ur = 1. To get Vr we don’t even need AT A since we know that A T Ur = V r Λr . So 1 −1 The complete SVD then is
T √ √ 1 1 1 A=1· 2· √ = 1 · 2 · √ [1 − 1] . 2 −1 2 √ √ where Ur = 1, Λr = 2, and Vr = 1/ 2[1 − 1]T .

·1=

√ 1 2Vr ⇒ Vr = √ 2

1 −1

.

So, the eigenvalues of AT A are 2 and 0. 2 is the eigenvalue of the (unnormalized) eigenvector (1, −1)T , while 0 is the eigenvalue of (1, 1)T . The latter follows from the fact that AV0 = 0 so v0 1 [1 − 1] = 0 ⇒ V0 = v1 1 This has a simple physical interpretation. An out-of-phase perturbation of velocity and depth (increase one and decrease the other) changes the travel time, while an in phase perturbation (increase both) does not. Since an in phase perturbation must be proportional to (1, 1)T , it stands to reason that this vector would be in the null space of A. But notice that we have made this physical argument without reference to the linearized (log parameter) problem. However, since we spoke in terms of perturbations to the model, the assumption of a linear problem was implicit. In other words, by thinking of the physics of the problem we were able to guess the singular vectors of the linearized problem without even considering the linearization explicitly. In the notation we developed for the SVD, we can say that Vr , the matrix of non-nullspace model singular vectors is (1, −1)T , while V0 , the matrix of null-space singular vectors is (1, 1)T . And hence, using the normalized singular vectors, the resolution operator is 1 1 −1 Vr VrT = √ . (7.15) 2 −1 1 The covariance matrix of the depth/velocity model is A† Cov(d)A†T = σ 2 A† A†T (7.16)

The important thing to notice here is that this says the velocity and depth are completely correlated (oﬀ diagonal entries magnitude equal to 1), and that the correlation is negative. This means that increasing one is the same as decreasing the other. The covariance matrix itself has the following eigenvalue/eigenvector decomposition. Cov(m) = σ 2
2

1 1 −1 1

1 0 0 0

1 −1 1 1

.

(7.18)

These orthogonal matrices correspond to rotations of the velocity/depth axes. These axes are associated with the line z/c = t. So for a given travel time t, we can be anywhere on the line z = tc: there is complete uncertainty in model space along this line and this uncertainty is reﬂected in the zero eigenvalue of the covariance matrix. For a two-dimensional problem such as this the correlation coeﬃcient measures the similarity in the ﬂuctuations in the two random variables. Here the two random variables are our estimates of z and c. Formally the correlation coeﬃcient is deﬁned to be: r= Czc . σz σc

It is not hard to show that for a two-dimensional Gaussian probability density, the level surfaces (contours of constant probability) are ellipses (circles and lines being special cases of ellipses). If the two random variables are zero-mean and have the same variances, then the level surfaces fall into one of three classes depending on the size of the correlation coeﬃcient. First note that the correlation coeﬃcient is always less than or equal to one in absolute value. If r = 0 then the level surfaces are circles. If 0 < |r| < 1, then the level surfaces are true ellipses. Finally, if |r| = 1 the level surfaces are lines, as in the example above.

As we have seen before, least squares tends to want to average over ignorance. Since we cannot determine velocity and depth individually, but only their ratio, least squares puts half the data into each. Damping does not change this, it is still least squares, but it does change the magnitude of the computed solution. Since damping penalizes the
1

7.1 The World’s Second Smallest Inverse Problem

113

norm of the solution, we can take an educated guess that the damped solution should, for large values of the damping parameter λ, tend to t λ 1 −1 . (7.20)

For small values of λ, the damped solution must tend to the already computed generalized inverse solution. It will be shown shortly that the damped generalized inverse solution is t 1 m† = . (7.21) λ −1 λ+2 Exact Solution The damped least squares estimator satisﬁes (AT A + λI)mλ = AT d. Since the matrix on the left is by construction invertible, we have m† λ = (AT A + λI)−1 AT d. If AT A = then (AT A + λI)−1 =

Damping changes the covariance structure of the problem too. We will not bother deriving a analytic expression for the damped covariance matrix, but a few cases will serve to illustrate the main idea. The damped problem [AT A + λI]m = AT d, is equivalent to the ordinary normal equations for the augmented matrix Aλ ≡ A √ λI (7.22)

where A is the original matrix and I is an identity matrix of dimension equal to the number of columns of A. In our toy problem this is 1 −1  √ 0  Aλ ≡  λ √  . 0 λ
1

Right away we can see that since the eigenvalues of this matrix are 1 and 3, instead of a degenerate ellipsoid (inﬁnite aspect ratio), the error ellipsoid of the damped problem has an aspect ratio of 3. As the damping increases, the covariance matrix becomes increasingly diagonal, resulting in a circular error ellipsoid. You will calculate the analytic result as an exercise. Your result should become degenerate as λ → 0.

Exercises

• Extend the two-parameter travel time inversion problem to the case in which the ray reﬂects from the ﬂat interface at an angle of θ, measured relative to the vertical. I.e, θ = 0 would correspond to a ray that goes straight up and down. Assume that the travel time can be measured with an uncertainty of σ second. • Compute the pseudoinverse and resolution matrix of 1 −1 2 0 4 −4 0 0 Assuming the right hand side is (0, 1)T , what is the least squares estimator of the 4-dimensional model vector. • Compute the pseudoinverse and resolution matrix of
    

1. −1 4 −4    0 1  0 −1



Assuming the right hand side is (0, 1, −1, 0)T , what is the least squares estimator of the 4-dimensional model vector. • Assuming the data covariance matrix is
    

1 0 0 0

0 .5 0 0

0 0 0 0 .1 0 0 .0001

    

compute the covariance of matrix of the least squares estimator.
1

BIBLIOGRAPHY

115

• Find the inverse of the covariance matrix in Equation 7.25. Now see if you can ﬁnd the square root of this matrix. I.e., a matrix such that when you square it, you get the inverse of the covariance matrix. • Compute the exact covariance and resolution for the damped two-parameter problem.

Chapter 8 Examples: Absorption and Travel Time Tomography
Before we go any further, it will be useful to motivate all the work we’re doing with an example that is suﬃciently simple that we can do all the calculations without too much stress. We will consider and example of “tomography”. The word tomography comes from the Greek tomos meaning section or slice. The idea is to use observed values of some quantity which is related via a line integral to the physical parameter we wish to infer. Here we will study seismic travel time tomography, is is a widely used method of imaging the earth’s interior. Mathematically this problem is identical to other types of tomography such as used in X-ray CAT scans.

8.1

The X-ray Absorber

Most of the examples shown here are based on an idealized two-dimensional x-ray absorption experiment. This experiment consists of the measurement of the power loss of x-ray beams passing through a domain ﬁlled with an x-ray absorber. We suppose that the absorber is conﬁned to the unit square, (x, y) ∈ [0, 1] ⊗ [0, 1]; we will represent this domain by DX . We also suppose that the sensing beam follows a perfectly straight path from transmitter to receiver and that the transmitter and receiver are located on the perimeter of the unit square. The geometry of a single x-ray absorption measurement looks like Figure 8.1.
1

118

Tomography

transmitter

receiver

Figure 8.1: An x-ray source shoots x-rays across a target to a detector where the intensity (energy) of the beam is measured.

8.1.1

The Forward Problem

Let c(x, y) be the absorption coeﬃcient in DX ; we assume that c(x, y) is non-negative everywhere. Let IT be the emitted beam intensity from a transmitter, T ; and let IR be the received intensity at some receiver, R. Then the absorption law is exactly IR = I T e−
T R

c(x,y)dλ

(8.1)

where the integral is along the (perfectly straight) path from T to R and dλ is arc-length along the path. (Note that c(x, y) = 0 in a vacuum and the exponent in equation (8.1) vanishes.) It is convenient to replace intensities with ρ= IT − I R , IT (8.2)

which is just the fractional intensity drop. ρ has the virtues that • ρ is independent of transmitter strength, IT , • ρ = 0 for a beam which passes only through a vacuum, • ρ ≥ 0 for all reasonable mediaa and, in fact, 0 ≤ ρ < 1, if c(x, y) is everywhere non-negative and ﬁnite.
a A “reasonable” medium is one which does not add energy to beams passing through. A laser is not a reasonable medium.

1

8.1 The X-ray Absorber For our uses, we will need two types of absorption calculations:

119

exact We will want an exact calculation which we can use to generate synthetic data to test our inverse algorithms. This calculation mimics the data generation process in nature. This calculation should either be exact or at least much more accurate than the associated linearized calculation (if we wish to study uncertainties and ambiguities in the inverse process, rather than errors in the synthetic data generator). linear We will also want a calculation in which the relation between a model’s parameters and the observations predicted for that model is linear. The precise linear relationship will form the heart of a linear inverse calculation. The diﬀerence between these two calculations is a measure of the problem’s nonlinearity. The next few subsections describe these two calculations in more detail. Exact Absorption Let ρexact be the exact absorption calculation. We can think of ρexact (c; T, R) as a computer program to which is given the transmitter and receiver locations and the function c(x, y) which deﬁnes the absorption coeﬃcient everywhere in DX . This program then returns the fractional intensity drop for the path T R through the medium c(x, y). The calculation itself is quite straightforward. In an actual application we would have to specify the accuracy with which the quadrature along the ray path is performed. In the calculations discussed here, we performed the quadrature by dividing the ray path into segments of a ﬁxed length and then summing the contribution to the integral from each tiny segment. We took the segment length to be about 10−3 ; recall that the sides of the model are of unit length.

8.1.2

Linear Absorption

A simple way to linearize the exact calculation, equation (8.1), is to assume that the path integral,
T R

c(x, y)dλ

is small. Since ex ≈ 1 + x for small x, we have IR ≈ I T 1 − or ρlinear ≈
T R T R

c(x, y)dλ

c(x, y)dλ

(8.3)

1

120

Tomography

0 0.0

10

20

% Error 30

40

50

0.1

0.2

0.3 Exact

0.4

0.5

0.6

Figure 8.2: The fractional error of the linearized absorption as a function of ρexact . This result is a good approximation to the extent that the total relative absorption is much less than one. An exact linearization can be had by changing the observable quantity from ρ to log(IR /IT ). Simply involves taking the logarithm of equation 8.1 leads to log(IR /IT ) = −
T R

c(x, y)dλ

which assures us that the logarithms of the observed intensities are exactly linear in the absorption coeﬃcient distribution (c(x, y)). We chose the approximate form, (8.3), when we developed the examples in order to induce some non-linearity into the calculation.

Linearization Errors It is very easy to compute the errors due to linearization in this case, since ρexact can be easily related to ρlinear as ρexact = 1 − e−ρlinear . A plot of the fractional error of the linearized absorption, ρlinear − ρexact ρexact as a function of ρexact is shown in Figure 8.2
1

8.1 The X-ray Absorber

121

Figure 8.3: Geometry of the tomography problem. The model is speciﬁed by blocks of constant absorption. Notice that the error expressed as a fraction is of the same order as the fractional absorption (ρexact ). (For example, when ρexact = 0.5, error ≈ 40%.) In a rough sense, if we think of linearization as neglecting a quadratic term which has about the same coeﬃcient as the linear term we are retaining, then we should expect an error of the same order as the quantity which has been linearized. Although this property is so simple as to be self-evident, it almost always comes as an ugly surprise in any particular application.

8.1.3

Model Representation

All of our inverse calculations use a model consisting of a regular array of homogeneous rectangular blocks. We completely specify a model’s geometry by specifying the number of blocks along the x-axis (Nx ) and the number of blocks along the y-axis (Ny ). We completely specify a model by specifying its geometry and by specifying the Nx Ny constant absorptivities of the blocks. The geometry of a model with Nx = 4, Ny = 5 is shown in Figure 8.3. We will need to map the cells in a model onto the set of integers {1, . . . Nx Ny }. Let C11 be the upper-left corner, CNx 1 be the upper-right corner, and let CNx Ny be the lower-right corner. The matrix {Cij } is mapped onto the vector {mk } a row at a time, starting with the ﬁrst (lowermost) row: {mi } = {C11 , . . . CNx 1 , C1,2 . . . CNx Ny } We chose this representation because it is very simple (possibly too simple for some applications) and it is strongly local. The latter property simply means that a perturbation in a model parameter only changes the values of c(x, y) in a limited neighborhood. Strong locality makes some results quite a bit easier to interpret; the trade-oﬀ is that locality is always associated with discontinuities in the representation’s derivatives.
1

122

Tomography

Figure 8.4: A perspective view of the model. Model Sensitivity Vector Let T R be the path joining a given transmitter-receiver pair, and let m be an absorption model, a vector of length Nx Ny . We want to ﬁnd a vector, q(T R), a function of the path T R and of dimension Nx Ny , such that ρlinear (m; T, R) = q(T R) · m. (8.4)

It is easy to see, by inspection of (8.3), that qi is simply the length of the portion of the path T R that passes through the ith block. Notice that the components of q depend only upon the model representation and the path T R. In particular, they are independent of the absorptivities.

8.1.4

Some Numerical Results

A perspective view of a model consisting of a centered disc of radius 0.25 is shown in Figure 8.4. inside the disc the absorption coeﬃcient is 0.1 and outside of the disc it vanishes. We sent nine shots through this structure. All of the shots came from a common transmitter in the upper-left corner and went to receivers spread along the right-hand side. The model and shot geometry looks like Figure 8.5. Figure 8.6 shows the computed values of ρexact as a function of receiver elevation. The numerical value of the extinction for the lowermost ray was 0.048757. This ray traveled from the point (0, 0.9), the transmitter, to (1, 0.1), the receiver. The value of the integrated absorptivities along the path, the path integral in equation (8.1), should have been exactly 0.05 (as a little contemplation should show). From this we compute
1

8.1 The X-ray Absorber

123

•

•

Figure 8.5: The model and shot geometry.

0.05

•

• • •

0.04

0.02

0.03

•

• 0.01 0.0

• 0.2 0.4 Receiver 0.6

• 0.8

•

Figure 8.6:

1

124

Tomography

s r r r r r

Figure 8.7: Plan view of the model showing one source and ﬁve receivers. ρtrue = 0.048771, which is in satisfactory agreement. (Note that the linearized estimate is exactly 0.05 and it is about 1% high.)

8.2

Travel Time Tomography

Along the same lines as the x-ray absorption problem, the time-of-ﬂight of a wave propagating through a medium with wavespeed v(x, y, z) is given by the line integral along the ray of the reciprocal wavespeed (slowness or index of refraction) dλ . v(x,y,z) v(x, y, z) The problem is to infer the unknown wavespeed of the medium from repeated observations of the time-of-ﬂight for various transmitter/detector locations. For the sake of deﬁniteness, let’s suppose that the source emits pressure pulses, the receiver is a hydrophone, and the medium is a ﬂuid. Figure (8.7) shows a 2D model of an anomaly embedded in a homogeneous medium. Also shown are 5 hypothetical rays between a source and 5 detectors. This is an idealized view on two counts. The ﬁrst is that not only is the raypath unknown–rays refract–but the raypath depends on the unknown wavespeed. This is what makes the travel time inversion problem nonlinear. On the other hand, if we can neglect the refraction of the ray, then the problem of determining the wavespeed from travel time observations is completely linear.
1

8.3 Computer Example: Cross-well tomography

125

The second complicating factor is that a “travel time” is only unambiguously deﬁned at asymptotically high frequencies. In general, we could deﬁne the travel time in various ways: ﬁrst recorded energy above a threshold, ﬁrst peak after the threshold value is surpassed, and others. Further, the travel times themselves must be inferred from the recorded data, although it is possible in some cases that this can be done automatically; perhaps by using a triggering mechanism which records a time whenever some threshold of activity is crossed. For purposes of this example, we will neglect both of these diﬃculties. We will compute the travel times as if we were dealing with an inﬁnite frequency (perfectly localized) pulse, and we will assume straight ray propagation. The ﬁrst assumption is made in most travel time inversion calculations since there is no easy way around the diﬃculty without invoking a more elaborate theory of wave propagation. The second assumption, that of linearity, is easily avoided in practice by numerically tracing rays through the approximate medium. But we won’t worry about this now.

8.3

Computer Example: Cross-well tomography

In the code directory you will ﬁnd various implementations of straight-ray tomography. These are extensive codes and will not be described in detail here. They begin by setting up the source/receiver geometry of the problem, computing a Jacobian matrix and fake travel times, adding noise to these and doing the least squares problem via SVD. Here we just show some of the results that you will be able to get. In Figure (8.8) you see a plot of the Jacobian matrix itself. The i − j element of this matrix is the length the i-th ray travels in the j-th cell. This comes from discretizing the travel time integral (t = s(r)d ) along each ray (one for each travel time). Black indicates zero elements and white nonzero. This particular matrix is about 95% sparse, so until we take advantage of this fact, we’ll be doing a lot of redundant operations, e.g., 0 × 0 = 0. Below this we show the “hit count”. This is the summation of the ray segments within each cell of the model and represents the total “illumination” of each cell. Below this we show the exact model whose features we will attempt to reconstruct via a linear inversion. Finally, before we can do an inversion, we need some data to invert. First we’ll compute the travel times in the true model shown above, then we’ll compute the travel times through a background model which is presumed to be correct except for the absence of the anomaly. It’s the diﬀerence between these two that we take to be the right hand side of the linear system Jδm = δd relating model perturbations to data perturbations. The computed solutions are shown
1

126

Tomography

600

500

400

300

Jacobian matrix

200

100

0 0 100 200 300 400

20

15

10

Illumination per cell

5

0 0 5 10 15 20

2 1.75 1.5 1.25 1 10 5 10 15 20 5

20 15

Exact model

Figure 8.8: Jacobian matrix for a cross hole tomography experiment involving 25 × 25 rays and 20×20 cells (top). Black indicates zeros in the matrix and white nonzeros. Cell hit count (middle). White indicates a high total ray length per cell. The exact model used in the calculation (bottom). Starting with a model having a constant wavespeed of 1, the task is to image the perturbation in the center.

1

8.3 Computer Example: Cross-well tomography in Figure 8.9.

127

Finally, in Figure (8.10), we show the spectrum of singular values present in the jacobian matrix, and one well resolved and one poorly resolved model singular vectors. Note well that in the cross hole situation, vertically stratiﬁed features are well resolved while horizontally stratiﬁed features are poorly resolved. Imagine the limiting case of purely horizontal rays. A v(z) model would be perfectly well resolved, but a v(x) model would be completely ambiguous.

Figure 8.10: The distribution of singular values (top). A well resolved model singular vector (middle) and a poorly resolved singular vector (bottom). In this cross well experiment, the rays travel from left to right across the ﬁgure. Thus, features which vary with depth are well resolved, while features which vary with the horizontal distance are poorly resolved.

1

130

Tomography

1

Chapter 9 From Bayes to Weighted Least Squares
In the last chapters we have developed the theory of least squares estimators for linear inverse problems in which the only uncertainty was the random errors in the data. Now we return to our earlier discussion of Bayes theorem and show how within the Bayesian strategy we can incorporate prior information on model parameters and still get away with solving weighted least squares calculations. Denote by f (m, d) the joint distribution on models and data. Recall that from Bayes’ theorem, the conditional probability on m given d is p(m|d) = f (d|m)ρ(m) , h(d)

where f (d|m) measures how well a model ﬁts the data, ρ(m) is the prior model distribution, and h(d) is the marginal density of d. The conditional probability p(m|d) is the so-called Bayesian posterior probability, expressing the idea that p(m|d) assimilates the data and prior information. For now we will assume that all uncertainties (model and data) can be described by Gaussian distributions. Since any Gaussian distribution can be characterized by its mean and covariance, this means that we must specify a mean and covariance for both the a priori distribution and the data uncertainties. In this case the Bayesian posterior probability is the normalized product of the following two functions: (2π)−n 1 −1 exp − (g(m) − dobs )T CD (g(m) − dobs ) , det CD 2 (9.1) where dobs is the vector of observed data which dimension is n, CD is the data covari1

132

From Bayes to Weighted Least Squares

ance matrix and g(m) is the forward operator; and 1 (2π)−m −1 exp − (m − mprior )T CM (m − mprior ) , det CM 2 (9.2) where m is the number of model parameters and CM is the covariance matrix describing the distribution of models about the prior model mprior . If the forward operator is linear, then the posterior distribution is itself a Gaussian. If the forward operator is nonlinear, then the posterior is non-Gaussian. The physical interpretation of Equation 9.1 is that it represents the probability that a given model predicts the data. Remember that d = g(mtrue ) + e where e is the noise. If we take expectations of both sides then E[d] = g(mtrue ) + E[e]. So if the errors are zero mean then the true model predicts the mean of the data. Of course, we don’t know the true model, but if we have an estimate of it, say m then g(m) is an estimate of the mean of the data and d − g(m) is an estimate of e. If we want to estimate the true model we still have the problem of deﬁning what sort of estimator we want to use. Maybe this is not what we want. It may suﬃce to ﬁnd regions in model space which have a high probability, as measured by the posterior. But for now let’s consider the problem of estimating the true model. A reasonable choice turns out to be: look for the mean of the posterior.a If the forward operator is linear (so that g(m) = Gm for some matrix G), then Tarantola [Tar87] shows that the normalized product σ(m) ∝ exp − 1 −1 (Gm − dobs )T CD (Gm − dobs ) 2 (9.3)

−1 + (m − mprior )T CM (m − mprior ) .

can be written as
−1 σ(m) ∝ exp (m − mmap )T CM (m − mmap ) ,

(9.4)

where
−1 −1 CM = G T CD G + C M

−1

,

(9.5)

is the covariance matrix of the posterior probability. This is approximately true even when g is nonlinear, provided it’s not too nonlinear.
a

We will take up the reasonableness of this choice later.
1

BIBLIOGRAPHY

133

Here mmap is the maximum of the posterior distribution, which for a Gaussian is also the mean. So to ﬁnd our estimator we need to optimize Equation 9.3. But that is equivalent to minimizing the exponent: min
m −1 −1 (Gm − dobs )T CD (Gm − dobs ) + (m − mprior )T CM (m − mprior ) .

(9.6)

But this is nothing but a weighted least squares problem. This is even easier to see if we introduce the square roots of the covariance matrices. Then the ﬁrst term above is (CD
−1/2 T

) (Gm − dobs ), CD

−1/2

(Gm − dobs ) = CD (m − mprior ) 2 .

−1/2

(Gm − dobs )

2

while the second term is

CM

−1/2

Here we have used two important facts. This ﬁrst is that the inverse of a symmetric matrix is symmetric. The second is that every symmetric matrix has a square root. To see this consider the diagonalization of such a matrix via an orthogonal transformation: A = QΛQT . So it is not too hard to see that A = QΛ1/2 QT QΛ1/2 QT = QΛ1/2 Λ1/2 QT = QΛQT

where the meaning of Λ1/2 is clear since it is diagonal with real elements. So QΛ1/2 QT is the square root of A.

Chapter 10 Bayesian versus Frequentist Methods of Inference
If we are Bayesians, then all prior information about models must be cast in the form of probabilities. However, in some cases, prior information is purely deterministic. For example, we know based on deﬁnition, that mass density and wavespeed are positive. There is nothing uncertain about this information. On the other hand, we know from observation that the average mass density of the Earth is less than 7g/cm2 . In priciple we can convert these pieces of deterministic information into probabilities. In this chapter we will compare and contrast the Bayesian and frequentist views of inference. This chapter is adapted from [SS97b] and [ST01]. There are two fundamentally diﬀerent meanings of the term ‘probability’ in common usage [SS97a]. If we toss a coin N times, where N is large, and see roughly N /2 heads, then we say the probability of getting a head in a given toss is about 50%. This interpretation of probability, based on the frequency of outcomes of random trails, is therefore called ‘frequentist’. On the other hand it is common to hear statements such as: ‘the probability of rain tomorrow is 50%’. Since this statement does not refer to the repeated outcome of a random trial, it is not a frequentist use of the term probability. Rather, it conveys a statement of information (or lack thereof). This is the Bayesian use of ‘probability’. Both ideas seem natural to some degree, so it is perhaps unfortunate that the same term is used to describe them. Bayesian inversion has gained considerable popularity in its application to geophysical inverse problems. The philosophy of this procedure is as follows. Suppose one knows something about a model before observing the data. This knowledge is cast in a probabilistic form and is called the prior probability model (prior means before the data have been observed.) Bayesian inversion then provides a framework for combining the probabilistic prior information with the information contained in the observed data in order to update the prior information. The updated distribution is the posterior conditional model distribution given the data; it is what we know about the model after
1

136

Bayesian versus Frequentist Methods of Inference

we have assimilated the data and the prior information. The point of using the data is that the posterior information hopefully constrains the model more tightly than the prior model distribution. However, the selection of a prior statistical model can in practice be somewhat shaky. For example, in a seismic survey we may have a fairly accurate idea of the realistic ranges of seismic velocity and density, and perhaps even of the vertical correlation length (if bore-hole measurements are available). However, the horizontal length scale of the velocity and density variation is to a large extent unknown. Given this, how can Bayesian inversion be so popular when our prior knowledge is often so poor? The reason for this is that in practice the prior model is used to regularize the posterior solution. Via a succession of diﬀerent calculations, the characteristics of the prior model are often tuned in such a way that the retrieved model has subjectively agreeable features. But logically, the prior distribution must be ﬁxed before hand. The features used to tune the prior should in fact be included as part of the prior information [GS97]. So, the practice of using the data to tune the prior suggests that the reason for the popularity of Bayesian inversion within the Earth sciences is inconsistent with the underlying philosophy. A common attitude seems to be: ‘If I hadn’t believed it, I wouldn’t have seen it.’ Since Bayesian statistics relies completely on the speciﬁcation of a prior statistical model, the ﬂexibility taken in using the prior model as a knob to tune properties of the retrieved model is completely at odds with the philosophy of Bayesian inversion. One can, however, use an empirical Bayes approach to use data to help determine a prior distribution. But having used the data to select a prior, one has to correct the uncertainty estimates so as not to be overconﬁdent [see Carlin and Louis [CL96]]. This correction is not usually done in geophysical Bayesian inversion.

10.0.1

Bayesian Inversion in Practice

There are two important questions that have to be addressed in any Bayesian inversion: • How do we represent the prior information? This applies both to the prior model information and to the description of the data statistics. • How do we summarize the posterior information? The second question is the easiest one to answer, at least in principle. It is just a matter of applying Bayes’ theorem to compute the posterior distribution. We then use this distribution to study the statistics of diﬀerent parameter estimates. For example, we can ﬁnd credible regions for the model parameters given the data, or simply use posterior means as estimates and posterior standard deviations as ‘error bars’. However, very seldom will we be able to compute all the posterior estimates analytically; we often
1

137 have to use computer-intensive approximations based on Markov Chain Monte Carlo methods [see for example, Tanner [Tan93]]. But still, a complete Bayesian analysis may be computationally intractable. The ﬁrst question is a lot more diﬃcult to answer. A ﬁrst strategy is a subjective Bayesian one: prior probabilities are designed to represent states of mind, prejudices or prior experience. But, depending on the amount and type of prior information, the choice of prior may or may not be clear. For example, if an unknown parameter, µ, must lie between a and b, is it justiﬁed to assume that µ has a uniform prior distribution on the interval [a, b]? We will address this question in an example below, but for now simply observe that there are inﬁnitely many probability distributions consistent with µ being in the interval [a, b]. To pick one may be an over-speciﬁcation of the available information. Even an apparently conservative approach, such as taking the probability distribution that maximizes the entropy subject to the constraint that µ lies in the interval, may lead to pathologies in high-dimensional problems. This shows how diﬃcult it may be to unambiguously select a prior statistical model. One way out of this dilemma is to sacriﬁce objectivity and presume that ‘probability lies in the eye of the beholder’. Of course, this means that our posterior probability will be diﬀerent from yours. A second approach attempts to make a somewhat more objective choice of priors by relying on theoretical considerations such as maximum entropy a [Jay82] , transformation invariance [Tar87], or by somehow using observations to estimate a prior. This latter approach is the empirical Bayes mentioned in a previous section. For example, suppose one is doing a gravity inversion to estimate mass density in some reservoir. Suppose further that there are available a large number of independent, identically distributed laboratory measurements of density for rocks taken from this reservoir (a big if!). Then one could use the measurements to estimate a probability distribution for mass density that could be used as a prior for the gravity inversion. This is the approach taken in [GS98], where in-situ (bore-hole) measurements are used to derive an empirical prior for surface seismic data. An empirical Bayes analysis is basically an approximation to a full hierarchical Bayes analysis based on the joint probability distribution of all parameters and available data. In other words, in a full Bayesian analysis the prior distribution may depend on some paramaters which in turn follow a second-stage prior. This latter prior can also depend on some third-stage prior, etc. This hierachical model ends when all the remaining parameters are assumed known. We can use the empirical Bayes approach when the last parameters can not be assumed to be known. Instead, we use the data to estimate the remaining parameters and stop the sequence of priors. We then proceed as in the standard Bayesian procedure. For an introduction to empirical and hierarchical Bayes methods see Casella [Cas85], Gelman et al. [GCSR97] and references therein. For a review on the development of objective priors see Kass and Wasserman [KW96].
a

See the appendix on entropy.
1

138

Bayesian versus Frequentist Methods of Inference

A third strategy is to abandon Bayes altogether and use only deterministic prior information about models: wave-speed is positive (a matter of deﬁnition); velocity is less than the speed of light (a theoretical prediction); the Earth’s average mass density is less than 7 g/cm3 (a combination of observation and theory that is highly certain). The inference problem is still statistical since random data uncertainties are taken into account. Essentially the idea is to look at the set of all models that ﬁt the data. Then perform surgery on this set, cutting away those models that violate the deterministic criteria, e.g., have negative density. The result will be a (presumably smaller and more realistic) set of models that ﬁt the data and satisfy the prior considerations. We choose any model that ﬁts the data to a desired level and satisﬁes the prior model constraints. Tikhonov’s regularization is one way of obtaining an inversion algorithm by restricting the family of models that ﬁt the data; for example, among all the models that ﬁt the data, we choose one that has particular features, the smoothest, the shortest, etc.

10.0.2

Bayes vs Frequentist

In the Bayesian approach, probability distributions are the fundamental tools. Bayesians can speak of the probability of a hypothesis given some evidence, and are able to conduct pre-data and post-data inferences. Frequentists, on the other hand, are more concerned with pre-data inference and run into diﬃculties when trying to give post-data interpretations to their pre-data formulation. In other words, uncertainty estimates, such as conﬁdence sets, are based on the error distribution, which is assumed to be known a priori, and on hypothetical repetitions of the data gathering process. We have seen that the choice of prior distributions is not always well deﬁned. In this case it would seem more reasonable to follow a frequentist approach. But it may also be the case that the determinism that frequentists rely on in the deﬁnition of parameters may be ill-deﬁned. For instance, if we are trying to estimate the mass of the earth, is this a precisly deﬁned, non-random quantity? Perhaps, but does the deﬁnition include the atmosphere? If so, how much of the atmosphere? If not, does it take into account that the mass is constantly changing (slightly) from, for example, micrometeorites? Even if you make the ‘true mass of the Earth’ well-deﬁned (it will still be arbitrary to some extent), it can never be precisely known. So, which approach is better? Bayesians are happy to point to some well known inconsistencies in the frequentist methodology and to diﬃculties frequentists face to use available prior information. Some Bayesians even go as far as claiming that anyone in her/his right frame of mind should be a Bayesian. Frequentists, on the other hand, complain about the sometimes subjective choice of priors and about the computational complexity of the Bayesian approach. In real down-to-earth data analysis we prefer to keep an open mind. Diﬀerent methods may be better than others depending on the problem. Both schools of inference have something to oﬀer. For colorful discussions on the comparison of the two approaches see Efron [Efr86] and Lindley [Lin75].
1

10.1 What Diﬀerence Does the Prior Make?

139

10.1

What Diﬀerence Does the Prior Make?

In a Bayesian calculation, whatever estimator we use depends on the prior and conditional distributions given the data. There is no clear established procedure to check how much information a prior injects into the posterior estimates. [This is one of the open problems mentioned in Kass and Wasserman [KW96].] In this example we will compare the risks of the estimators. ˆ To measure the performance of an estimator, m, of m we deﬁne the loss function, ˆ L(m, m); L is a non-negative function which is zero for the true model. That is, for any other model m1 , L(m, m1 ) ≥ 0 and L(m, m) ≡ 0. The loss is a measure of the cost ˆ of estimating the true model with m when it is actually m. For example, a common ˆ ˆ loss function is the square error: L(m, m) = (m − m)2 . But there are other choices p ˆ ˆ like p -norm error: L(m, m) = m − m . ˆ ˆ The loss, L(m, m), is a random variable since m depends on the data. We average over ˆ the data to obtain an average loss. This is called the risk of m given the model m ˆ ˆ R(m, m) = EP L(m, m), (10.1)

where P is the error probability distribution and EP the expectation with respect to this distribution. For square error loss the risk is the usual mean square error.

10.1.1

Bayes Risk

The expected loss depends on the chosen model. Some estimators may have small risks for some models but not for others. To compare estimators we need a global measure that takes all plausible models into account. A natural choice is to take the expected value of the loss with respect to the posterior distribution, p(m|d), of the model given the data. This is called the posterior risk ˆ rm|d = Em|d L[m, m(d)]. Alternatively we can take a weighted average of the risk (10.1) using the prior model distribution as weight function. This is the Bayes risk ˆ rρ = Eρ R(m, m), where ρ is the prior model distribution. An estimator with the smallest Bayes risk is called a Bayes estimator. Note that we have used a frequentist approach to deﬁne the Bayes risk, since we have not conditioned on the observed data. It does make sense, however, to expect good frequentist behavior if the Bayesian approach is to be used repeatedly with diﬀerent data sets. In addition, it can be shown that, under very general conditions, minimizing the Bayes risk is equivalent to minimizing the posterior risk [Ber85].
1

140

Bayesian versus Frequentist Methods of Inference

Let f denote the joint distribution of models and data. The distribution (marginal) of the data is obtained by integrating f over the models h(d) =
M

f (m, d)dm,

where M is the space of models. From Bayes’ theorem, the conditional distribution of m given d is f (d|m)ρ(m) p(m|d) = , h(d) where ρ(m), the prior distribution, is the marginal of f with respect to m. The conditional distribution, p(m|d), is the so-called Bayesian posterior distribution, which updates the prior information in view of the data. One can deﬁne a number of reasonable estimators of m based on p(m|d). For example, ˆ the m that maximizes p(m|d) (or that is close, in probability, to m.) Or one could compute the estimator that gives the smallest Bayes risk for a given prior and loss function. It can be shown [Lehmann [Leh83], p.239] that, for square error loss function, the Bayes estimator is the posterior mean. Here is a simple example of using a normal prior to estimate a normal mean. Assume that there are n observations, (d1 , d2 , ..., dn ) = d, which are iid N (η, σ 2 ) and that we want to estimate the mean, η, given that the prior, ρ, is N (µ, β 2 ). Up to a constant factor, the joint distribution of η and d is [Lehmann [Leh83], p.243] 1 f (d, η) = exp − 2σ The posterior mean is η = E(η|d) = ˆ
n i=1

(di − η)2 exp −

1 (η − µ)2 , 2β

¯ nd/σ 2 + µ/β 2 , n/σ 2 + 1/β 2

¯ where d is the arithmetic mean of the data. The posterior variance is 1 Var(η|d) = . 2 + 1/β 2 n/σ Notice that the posterior variance is always reduced by the presence of a nonzero β. The posterior mean, which is the Bayes estimator for square error loss, can be written as n/σ 2 1/β 2 ¯ η (d) = ˆ µ. d+ n/σ 2 + 1/β 2 n/σ 2 + 1/β 2 We see that the Bayes estimator is a weighted average of the mean of the data and the mean of the Bayesian prior distribution; the latter is the Bayes estimator before any data have been observed. The Bayes risk is the integral, over the data, of the posterior variance of η. Since the posterior variance does not depend of d, the Bayes risk is just the posterior variance. Note also that as β → 0, increasingly strong prior information, the estimate converges to the prior mean. As β → ∞, increasingly weak prior information, the Bayes estimate converges to the mean of the data. Also note that as β → ∞ the prior becomes improper (not normalizable).
1

10.2 Example: A Toy Inverse Problem

141

10.1.2

What is the Most Conservative Prior?

It often happens that there is not enough information to choose a prior density for the unknown parameters, or that the information available is not easily translated into a probabilistic statement; yet we need a prior to be able to apply Bayes’ theorem. In this case we try to ﬁnd a ‘noninformative’, or ‘conservative’, prior that will allow us to conduct the Bayesian inference while injecting a minimum of artiﬁcial information; that is, information which is not justiﬁed by the physical process. We have deﬁned the Bayes risk, rρ , and the Bayes estimator for a given prior density. It stands to reason that the more informative the prior the smaller its associated risk; we therefore say that a prior, ρ, is least favorable if rρ ≥ rρ for any other prior, ρ . A least favorable prior is associated with the greatest unavoidable loss. In the frequentist approach the greatest unavoidable loss is associated with the maximum of the risk (10.1) over all the possible models. An estimator that minimizes this maximum risk is called a minimax estimator. Under certain conditions the Bayes estimator corresponding to a least favorable prior actually minimizes the maximum risk [see Lehmann [Leh83]]. This is true, for example, when the Bayes estimator has a constant risk. In this sense we can think of a least favorable prior as being a route to the most conservative Bayesian estimator. How does one ﬁnd a conservative (noninformative) prior? There is no easy answer, even the terms ‘conservative’ and ‘noninformative’ are not well deﬁned. One possibility is to deﬁne a measure of information (e.g., entropy) and determine a prior which minimizes/maximizes this measure (e.g., maximum entropy). We could also look for priors which are invariant under some family of transformations.

10.2

Example: A Toy Inverse Problem

We consider a simple example of estimating the mean, η, of a unit variance normal distribution, N (η, 1), with an observation, d, from N (η, 1) given that |η| is known to be bounded by β. Following Stark [Sta97], we will use this as a model of an inverse problem with a prior constraint. Without the prior bound, d is an estimator of η but we hope to do better (obtain a smaller risk) by including the bound information. How can we include this information in the estimation procedure? One possibility is to use a Bayesian approach and assign a prior distribution to η which is uniform on [−β, β]. We will show that this distribution injects stronger information than might be evident.
1

142

Bayesian versus Frequentist Methods of Inference

10.2.1

Bayes Risk

Start with an observation, d, from N (η, 1) and suppose we know a priori that |η| is bounded by β. We incorporate the bound by assigning to η a prior uniformly distributed on [−β, β]. The joint distribution of η and d is then f (d, η) = 1 1 1 I[−β,β] √ exp − (d − η)2 , 2β 2 2π

where I[−β,β] (x) = 1 for x ∈ [−β, β] and zero otherwise. We reproduce Stark’s Monte Carlo calculation of the Bayes risk for this problem. Figure 10.1 shows the Bayes risk, using a uniform prior on [−β, β], and the minimax risk to be described next. As the constraint weakens (β increases) the Bayes risk gets closer to 1. (The dashed and dotted curves in this ﬁgure will be explained in the next section.)

10.2.2

The Flat Prior is Informative

We have used the uniform distribution to ‘soften’ (i.e., convert to a probabilistic statement) the constraint |η| ≤ β. Now we want to measure the eﬀect of this constraint softening. Have we included more information than we really had? Given the observation, d, from N (η, 1) and knowing that |η| ≤ β, what is the worst risk (mean square error) we may hope to achieve with the best possible estimator without imposing a prior distribution on η? In other words we want to compute the minimax risk, R(β), given the bound β R(β) = minδ maxη∈[−β,β] EP [η − δ(d)]2 . R(β) is a lower bound for the maximum risk of any other estimator. Although it is diﬃcult to compute its exact value, it is easy to see that R(β) ≤ min{β 2 , 1}. In addition, Donoho et al. [DLM90] show that 4 β2 ≤ R(β). 5 β2 + 1 Figure 1 shows upper and lower bounds for for the minimax risk as a function of β. Note that for β ≤ 3 the Bayes risk is outside the minimax bounds. This is an artifact of the way we have ‘softened’ the bound. In other words, the uniform prior distribution injects more information than the hard bound on η, as judged by comparing the most pessimistic frequentist risk with that of the Bayesian estimator. It can also be shown that R(b) → 1 as b → ∞. So, as the bound weakens the Bayes and minimax risk both approach 1.
1

10.2 Example: A Toy Inverse Problem

143

1.0

0.8

risk

0.6

0.4 Bayes MC simulation Minimax lower bound Minimax upper bound

0.2

0.0

0

2

4

6 8 interval half−width

10

12

Figure 10.1: For square error loss, the Bayes risk associated with a uniform prior is shown along with the upper and lower bounds on the minimax risk as a function of the size of the bounding interval [−β, β]. When β is comparable to or less than the variance (1 in this case), the risk associated with a uniform prior is optimistic

1

144

Bayesian versus Frequentist Methods of Inference

10.3

Priors in High Dimensional Spaces: The Curse of Dimensionality

As we have just seen, most probability distributions usually have more information than implied by a hard constraint. To say, for instance, that any model with m ≤ 1 is feasible is certainly not the same thing as saying that all models with m ≤ 1 are equally likely. And while we could look for the most conservative or least favorable such probabilistic assignment, Backus [Bac88] makes an interesting argument against any such probabilistic replacement in high- or inﬁnite-dimensional model spaces. His point can be illustrated with a simple example. Suppose that all we know about an ndimensional model vector, m, is that its length, m ≡ m , is less than some particular value–unity for the sake of deﬁniteness. In other words, suppose we know a priori that m is constrained to be within the n-dimensional unit ball, Bn . Backus considers various probabilistic replacements of this hard constraint; this is called ‘softening’ the constraint. We could for example choose a prior probability on m which is uniform on Bn . Namely, the probability that m will lie in some small volume, δV ∈ Bn , shall be equal to δV divided by the volume of Bn . Choosing this uniform prior on the ball, it is not diﬃcult to show that the expectation of m2 for an n-dimensional m is E(m2 ) = n n+2

which converges to 1 as n increases. Unfortunately, the variance of m2 goes as 1/n for large n, and thus we seem to have introduced a piece of information that was not implied by the original constraint; namely that for large n, the only likely vectors, m, will have length equal to one. The reason for this apparently strange behavior has to do with the way volumes behave in high dimensional spaces. The volume, Vn (R), of the R − diameter ball in n dimensional space is Vn (R) = Cn Rn , where Cn is a constant that depends only on the dimension n, not on the radius. [This is a standard result in statistical mechanics; e.g., Becker [Bec67].] If we compute the volume, V ,n , of an n-dimensional shell of thickness just inside an R-diameter ball we can see that V ,n ≡ Vn (R) − Vn (R − ) = Cn (Rn − (R − )n ) = Vn (R) 1 − 1 − Now, for /R 1 and n 1 we have V ,n ≈ Vn (R) 1 − e−n /R . This says that as n gets large, nearly all of the volume of the ball is compressed into a thin shell just inside the radius.
1
n

R

.

(10.2)

10.3 Priors in High Dimensional Spaces: The Curse of Dimensionality 145 But even this objection can be overcome with a diﬀerent choice of probability distribution to soften the constraint. For example, choose m to be uniformly distributed on [0, 1] and choose the n − 1 spherical polar angles uniformly on their respective domains. This probability is uniform on m , but non-uniform on the ball. However it is consistent with the constraint and has the property that the mean and variance of m2 is independent of the dimension of the space. So, as Backus has said, we must be very careful in replacing a hard constraint with a probability distribution, especially in a high-dimensional model space. Apparently innocent choices may lead to unexpected behavior.

Appendix: Entropy
This appendix is a brief introduction to entropy as it relates to inversion. For more details see [GMS01]. In Bayesian inversion we use probabilities to represent states of information. But just how does one quantify such a state? Is it possible to say that one probability has more information than another? Consider an experiment with N possible outcomes each occuring with a probability pi . In analogy with the statistical mechanical deﬁnition of entropy, [Sha48] introduced the following deﬁnition of the entropy for such discrete probabilities: H(p) = − pi log pi .
i

(10.3)

Following Shannon, three postulates should be satisﬁed by H(p) or any other measure of information. Those are: 1. Continuity; 2. Monotonicity, and 3. Composition Law. These postulates are discussed in detail in [GMS01], here a qualitative understanding is suﬃcient. The ﬁrst postulate requires that we should not gain or loose a large amount of information by making a small change to the probabilities. The second postulate, monotonicity, refers to the information associated with a collection of independent, equally likely events. It is clear that in such a case the uncertainty must increase monotonically with the number of possible outcomes. The third postulate requires that it should not matter how one regroups the events of a given set. The entropy of the set should stay the same. To see the meaning of Equation 10.3, consider an experiment whose outcome is known with absolute certainty. Then, the corresponding probability density is a Kronecker
1

146

Bayesian versus Frequentist Methods of Inference

delta pi = δiq , where q is the certain event, and H(p) = 0 (no uncertainty). If there are two equally likely outcomes, then H(p) = −1/2 log(1/2) − 1/2 log(1/2) = log 2. Whereas if one of the events has a probability 1/10 and one has probability 9/10 then the entropy is H(p) = −1/10 log(1/10) − 9/10 log(9/10) = log 10 − .9 log 9, which is about half that in the equally likely case. For M equally likely events H(p) = log M . And in the limit that M goes to inﬁnity, then the uncertainty must too. The usefulness of the deﬁnition 10.3 depends on the deﬁnition of the probability p. If the p is a 1D distribution associated with the frequency of outcomes of the possible events, then p is not aﬀected by the correlation of the points. E.g., if we sample 10 points pseudo-randomly from a probability distribution with two equally likely outcomes 0 and 1, we might see something like the following 0010110110. As luck would have it there are 5 0’s and 5 1’s. Now sort these outcomes in increasing order 0000011111. There are still 5 0’s and 5 1’s but we certainly would not regard the latter experiment as representing the same degree of uncertainty as the former. Similarly, if we sample 1000 points independently from a Gaussian we’ll see a nice bell-shaped curve. But the frequencies of the binned events are independent of their order. So, once again, sorting them into monotonic order will not change the entropy. of course, the dependence (correlation) among events is not being considered by an unidimensional probability distribution. Once multidimensional probabilities are used in deﬁnition (10.3), such correlation can be accounted for in entropy calculations. Now, in inverse theory we are always comparing one state of information to another— what we know before relative to what we know after. So it is more appropriate to measure the relative information of one probability compared to another. Let us therefore revise the original deﬁnition of discrete entropy (Equation (10.3)), and introduce the concept of relative entropy. Thus, H[pi ; qi ] = − pi log
i

10.3 Priors in High Dimensional Spaces: The Curse of Dimensionality 147 which is ﬁnite. We will adopt Equation (10.4) as the deﬁnition of relative entropy in the discrete case, and, as commonly done, the last expression of Equation (10.6) as the deﬁnition of relative entropy in the continuous case. q(x), or qi , represents a state of information against which we make comparisons. Finally it is worth mentioning that the negative of the quantity H[p(x); q(x)], known as cross-entropy, was ﬁrst deﬁned by Kullback [Kul59] as the directed divergence. This quantity deﬁnes the amount of information of the probability density p(x) with respect to q(x). See also [SJ81]. Convinced that entropy is a suitable measure for the uncertainty of a probability distribution, Jaynes [Jay57] showed that a useful tool for conservatively assigning probabilities was to maximize the entropy of the unknown distribution subject to constraints on its moments. Mathematically this variational problem can be expressed by maximizing Equation (10.4) subjected to the normalization of the distribution
N

p(xi ) = 1,
i=1

(10.6)

and to other constraints given in the form of expectations
N

wk (x) =
i=1

wk (xi )p(xi ),

k = 1, ..., K.

(10.7)

This is equivalent to the unconstrained problem, given by
N

S(p; λ, q) = −

p (xi ) ln
i=1 N

p (xi ) q (xi ) p (xi ) − 1 (10.8)

−(λ0 − 1)
K

i=1 N

−

λk
k=1 i=1

wk(xi ) p (xi ) − µk ,

where µk are sample estimates of wk (x) and the λk are the Lagrange multipliers associated with the constraints. Note that the term (λ0 − 1) is just a redeﬁnition of the zero-order Lagrange multiplier introduced for convenience. If we take the ﬁrst variation of the functional S(p; λ, q) with respect to the probabilities, we get that δS(p; λ, q) equals N K ∂H − (λ0 − 1) − λk wk(xi ) δp(xi ), (10.9) i=1 ∂p(xi ) k=1 with p (xi ) ∂H = − ln +1 . ∂p(xi ) q (xi )
1

(10.10)

148

BIBLIOGRAPHY

The solution to the problem can be found in the usual way by letting δS(p) = 0, which yields
K

However, to determine the maximum-entropy probability function we still need to ﬁnd the values for the other Lagrange multipliers and to specify the prior probability q(xi ). Techniques and examples of maximum entropy calculations of this sort are described in [GMS01].

Chapter 11 Iterative Linear Solvers
We have seen throughout the course that least squares problems are ubiquitous in inverse theory for two main reasons. First, least squares gives rise to linear problems, which are relatively easy to deal with. And secondly, ﬁnding the maximum of a Gaussian distribution is a least squares problem. That means that if the ﬁnal a posteriori probability on the models is Gaussian, then ﬁnding the maximum a posteriori (MAP) model amounts to solving a weighted least squares problem. For both reasons, least squares is very important and a number of specialized numerical techniques have been developed. In this chapter we digress and discuss a very useful class of iterative algorithms for solving linear systems. These methods are at the core of most large-scale inverse calculations.

11.1

Classical Iterative Methods for Large Systems of Linear Equations

A direct method for solving linear systems involves a ﬁnite sequence of steps, the number of which is known in advance and does not depend on the matrix involved. Usually nothing can be gained by halting a direct method early; it’s all or nothing. If the matrix is sparse, direct methods will almost always result in intermediate ﬁll, the creation of new nonzero matrix elements during the solution. Fill can usually be mitigated by carefully ordering operations and/or the matrix elements. Even so, direct methods for sparse linear systems require a certain amount of sophistication and careful programming. On the other hand, iterative methods start with some approximation, perhaps random, to the solution of the linear system and reﬁne it successively until some “stopping criterion” is satisﬁed. Most iterative methods do not require the matrix to be explicitly deﬁned; it often suﬃces to know the action of the matrix (and perhaps its transpose) on arbitrary vectors. As a result, ﬁll does not occur and the data structures necessary to store and manipulate the matrix elements can be quite simple. The general
1

152

Iterative Linear Solvers

subject of iterative methods is vast and in no sense will a survey be attempted. The aim of the ﬁrst section is simply to get the ball rolling and introduce a few classical methods before diving into conjugate gradient. In addition, the classical iterative methods are mostly based on matrix “splitting” which plays a key role in preconditioned conjugate gradient. This brief discussion is patterned on Chapter 8 of [SB80] and Chapter 4 of [You71]. Young is the pioneer in computational linear algebra and his book is the standard reference in the ﬁeld. Let A be a nonsingular n × n matrix and x = A−1 h be the exact solution of the system Ax = h. A general class of iterative methods is of the form xi+1 = Φ(xi ), i = 0, 1, 2, . . . (11.2) (11.1)

In order for this to work one must be able to solve (11.5). Further, the closer B is to A, the smaller the moduli of the eigenvalues of I − B −1 A will be, and the more rapidly will (11.6) converge. Many of the common iterative methods can be illustrated with the following splitting. A=D−E−F (11.7) where D = diag(A), −E is the lower triangular part of A and −F is the upper triangular part of A. Now, using the abbreviations L ≡ D −1 E, U ≡ D −1 F, J ≡ L + U, H ≡ (I − L)−1 U (11.8)

and assuming ai,i = 0 ∀i, one has
a

The spectral radius of an operator is the least upper bound of its spectrum σ: ρ(Φ) ≡ sup | λ |
λ∈σ(Φ)

where the subscript in parentheses refers to the iteration number. Jacobi’s method is also called the “total step method.” To get the “single step” or Gauss-Seidel method choose B to be the lower triangular part of A including the diagonal: Algorithm 2 Gauss-Seidel Method B = D − E, aj,k xk(i+1) + aj,j xj(i+1) +
k<j k>j

More generally still, one may consider using a class of splitting matrices B(ω) depending on a parameter ω, and choosing ω in such a way as to make the spectral radius of I − B −1 (ω)A as small as possible. The “relaxation” methods are based on the following choice for B: Algorithm 3 Relaxation Methods 1 D(I − ωL) ω = (B(ω) − A)xi + h i = 0, 1, . . . B(ω) = (11.13) (11.14)

B(ω)xi+1

For ω > 1 this is called overrelaxation, while for ω < 1 it is called underrelaxation. For ω = 1 (11.14) reduces to Gauss-Seidel. The rate of convergence of this method is determined by the spectral radius of I − B −1 (ω)A = (I − ωL)−1 [(1 − ω)I + ωU ] (11.15)

In particular, the Gauss-Seidel method (ω = 1) converges for positive deﬁnite matrices. For a proof of this result, see [SB80], pages 547–548. This result can be considerably sharpened for what Young calls type-A matrices or the “consistently ordered” matrices (see, for example, [You71], chapter 5).
b

A matrix A is positive if (x, Ax) ≥ 0 for all x. It is positive deﬁnite if the inequality is strict.
1

154

Iterative Linear Solvers

11.2

Conjugate Gradient

Conjugate gradient is by far the most widely used iterative method for solving large linear systems. In its simplest forms it is easy to program and use, yet retains the ﬂexibility to tackle some very demanding problems. Theoretically, CG is a descendant of the method of steepest descent, which is where the discussion begins. But ﬁrst, a few deﬁnitions.

11.2.1

Inner Products

We will assume that vectors lie in ﬁnite dimensional Cartesian spaces such as Rn . An inner product is a scalar-valued function on Rn × Rn , whose values are denoted by (x, y), which has the following properties:

where A ∈ Rn×n ; h, x ∈ Rn ; and c is a constant. The quadratic form is said to be symmetric, positive, or positive deﬁnite, according to whether the matrix A has these properties. The gradient of a symmetric quadratic form f is f (x) = Ax − h. (11.22)

The fact that solutions of Ax = h can be maxima or saddle points complicates things slightly. We will use the concept of positivity for a matrix. A matrix A is said to be positive if (x, Ax) ≥ 0 for all x. So one must make a few assumptions which are clariﬁed by the following lemma. Lemma 2 Suppose that z is a solution of the system Az = h, A is positive and symmetric, and f (x) is the quadratic form associated with A, then 1 f (x) = f (z) + ((x − z), A(x − z)). 2 (11.23)

This means that z must be a minimum of the quadratic form since the second term on the right is positive. Thus the value of f at an arbitrary point x must be greater than its value at z. To prove this, let x = z + p where Az = h. Then 1 f (x) = f (z + p) = ((z + p), A(z + p)) − (h, (z + p)) + c. 2 1 = f (z) + {(z, Ap) + (p, Az) + (p, Ap)} − (h, p). 2 If A is symmetric, the ﬁrst two terms in brackets are equal, hence: 1 f (x) = f (z) + (p, Ap) + (Az, p) − (h, p). 2 But by assumption Az = h, so that 1 1 f (x) = f (z) + (p, Ap) = f (z) + ((x − z), A(x − z)) 2 2 which completes the proof. As a corollary one observes that if A is positive deﬁnite as well as symmetric, then z is the unique minimum of f (z) since in that case the term ((x − z), A(x − z)) is equal to zero if and only if x = z. It will be assumed, unless otherwise stated, that the matrices are symmetric and positive deﬁnite. The level surfaces of a positive deﬁnite quadratic form (i.e, the locus of points for which f (x) is constant) is an ellipsoid centered about the global minimum. And the semiaxes of this ellipsoid are related to the eigenvalues of the deﬁning matrix. The negative gradient of any function points in the direction of steepest descent of the function. Calling this direction r one has r = −f (x) = h − Ax = A(z − x) (11.24)

since Az = h. The idea behind the method of steepest descents is to repeatedly minimize f along lines deﬁned by the residual vector. A prescription for this is given by the following lemma.
1

In other words, there exists a constant α such that by moving by an amount 2α along the residual, one ends up on the other side of the ellipsoid f (x) = constant. And further, if one moves to the midpoint of this line, one is assured of being closer to (or at the very least, not farther away from) the global minimum. The proof of this assertion is by construction. From the deﬁnition of f one has for arbitrary x, α f (x + 2αr) = 1 ((x + 2αr), A(x + 2αr)) − (h, (x + 2αr)) + c 2 1 = f (x) + {(2αr, Ax) + (2αr, A2αr) + (x, A2αr)} − (h, 2αr) 2 = f (x) + 2α(r, Ax) + 2α2 (r, Ar) − 2α(h, r) = f (x) − 2α(r, r) + 2α2 (r, Ar)

is a monotone sequence which is bounded below by the unique minimum f (z). That such a sequence must converge is intuitively clear and indeed follows from the Monotone Convergence Theorem. The proof of this theorem relies on a surprisingly deep property of real numbers: any nonempty set of real numbers which has a lower bound, has a greatest lower bound (called the inﬁnum). Having thus established the convergence of f (xk ) to f (z), the convergence of xk to z follows from Lemma 2 and the properties of inner products: 1 f (z) − f (xk ) = − (xk − z, A(xk − z)) → 0 ⇒ xk − z → 0 2 since A is positive deﬁnite. There is a drawback to steepest descent, which occurs when the ratio of the largest to the smallest eigenvalue (the condition number κ) is very large; the following result quantiﬁes ill-conditioning for quadratic minimization problems. Theorem 13 Let λmax and λmin be the largest and smallest eigenvalues of the symmetric positive deﬁnite matrix A. Let z be the minimum of f (x) and r the residual associated with an arbitrary x. Then r r ≤ f (x) − f (z) ≤ 2λmax 2λmin where x
2

(11.29)

(11.30)

≡ (x, x) is the Euclidean norm.

If all the eigenvalues of A were the same, then the level surfaces of f would be spheres, and the steepest descent direction would point towards the center of the sphere for any initial vector x. Similarly, if there are clusters of nearly equal eigenvalues, then steepest descent will project out the spherical portion of the level surfaces associated with those eigenvalues nearly simultaneously. But if there are eigenvalues of very diﬀerent magnitude then the portions of the level surfaces associated with them will be long thin ellipsoids. As a result, the steepest descent direction will not point towards the quadratic minimum. Depending upon the distribution of eigenvalues, steepest descent has a tendency to wander back and forth across the valleys, with the residual changing very little from iteration to iteration. The proof of this result is as follows. Let x = z+p where x is arbitrary. From Lemma 2, 1 f (x) − f (z) = (p, Ap) 2 Now, p = −A−1 r, so that 1 1 −1 (A r, r) = (p, Ap) 2 2
1

where A = diag(10, 1) and h = (1, −1). If we start the steepest descent iterations with x0 = (0, 0) then the ﬁrst few residuals vectors are: (1, −1) , (−9/11, −9/11), (81/121, 81/121) and so on. In general the even residuals are proportional to (1, −1) and the odd ones are proportional to (−1, −1). The √ coeﬃcients are (9/11)n , so the matrix were norm of the residual vector at the i-th step is ri = 2(9/11)i . If the √ A = diag(100, 1) instead, the norm of the i-th residual would be ri = 2(99/101)i : steepest descent would be very slow to converge. This can be seen graphically from a plot of the solution vector as a function of iteration superposed onto a contour plot of the quadratic form associated with the matrix A, shown in Figure (11.1). It is not a coincidence that the residuals at each step of steepest descent are orthogonal to the residuals before and after. We can prove this generally: rk = h − Axk = h − A(xk−1 + αk rk−1 ) = rk−1 − αk Ark−1 .
1

So the residuals are pairwise orthogonal. The question naturally arises, is convergence always asymptotic? Is there ever a situation in which SD terminates in exact arithmetic? Using the above expression rk = rk−1 − αk Ark−1 (11.36) we see that rk = 0 if and only if rk−1 = αk Ark−1 . But this just means that the residual at the previous step must be an eigenvector of the matrix A. We know that the eigenvectors of any symmetric matrix are mutually orthogonal, so this means that unless we start the steepest descent iteration so that the ﬁrst residual lies along one of the principal axes of the quadratic form, convergence is not exact.

Symbolic arithmetic packages, including Mathematica, can solve the problem in exact arithmetic. This is very useful for analysis of the eﬀects of rounding errors. Figure out geometrically what steepest descent is doing. Does SD ever converge in ﬁnitely many steeps on this problem in exact arithmetic? In this case you should be able to derive an analytic expression for the residual vector. Make plots showing the level curves of the quadratic form associated with A. Then plot the solution vector as a function of iteration. The changes should always be normal to the contours. Under what circumstances can the residual vector be exactly zero? What is the geometrical interpretation of this? Do your conclusions generalize to symmetric non-diagonal matrices? What happens if you change the matrix from diag(10,1) to diag(100,1)?

11.2.5

The Method of Conjugate Directions

The problem with steepest descent (SD) is that for ill-conditioned matrices the residual vector doesn’t change much from iteration to iteration. A simple scheme for improving its performance goes back to Fox, Husky, and Wilkinson [FHW49] and is called the conjugate direction (CD) method. Instead of minimizing along the residual vector, as in SD, minimize along “search vectors” pk which are assumed (for now) to be orthogonal with respect to the underlying matrix. This orthogonality will guarantee convergence to the solution in at most n steps, where n is the order of the matrix. So replace the step xk = xk−1 + αk rk−1 with xk = xk−1 + αk pk−1 where p is to be deﬁned. As in SD the idea is to minimize f along these lines. The scale factors α, as in SD, are determined by the minimization. Using the proof of Lemma 2, 1 f (xk + αpk ) = f (xk ) + (αpk , Aαpk ) − (αpk , rk ) 2 1 2 = f (xk ) + α (pk , Apk ) − α(pk , rk ) 2 Setting
∂f (xk +αpk ) ∂α

So, provided the scale factors αk satisfy the last equation, one is guaranteed to minimize the residual along the search vector pk . The conditions necessary for the search vectors are given by the following theorem. Theorem 14 Conjugate Direction Theorem Suppose that the search vectors are chosen such that (pi , Apj ) = 0 if i = j (A-orthogonality), then the CD method converges to the exact solution in at most n steps. Proof. Using the CD iteration xk = xk−1 + αk pk−1 , one has by induction xk − x0 = α1 p0 + α2 p1 + · · · + αk pk−1 for any x0 chosen. Since the p vectors are A-orthogonal, it follows that (pk , A(xk − x0 )) = 0. The A-orthogonality also implies that the p vectors must be linearly independent. Thus any vector in Rn can be represented as an expansion in the {pk }n−1 . In particular, the k=0 unknown solution z of the linear system can be written z = γ0 p0 + · · · + γn−1 pn−1 . Taking the inner product of this equation with ﬁrst A and then pi , and using the A-orthogonality gives (pi , Az) = γi (pi , Api ) ⇒ γi = (pi , Az) . (pi , Api ) k = 1, 2, . . .

The idea of the proof is to show that these numbers, namely the γi , are precisely the coeﬃcients of the CD algorithm; that would automatically yield convergence since by proceeding with CD we would construct this expansion of the solution. Just as an arbitrary vector x can be expanded in terms of the linearly independent search vectors, so can z − x0 where x0 is still the initial approximation. Thus,
n−1

The A-orthogonality can be seen to arise geometrically from the fact that the vector which points from the current location x to the global minimum of the quadratic form z must be A-orthogonal to the tangent plane of the quadratic form. To see this observe that since since the residual r must be normal to the surface, a tangent t must satisfy (t, r) = 0. Therefore 0 = (t, Ax − h) = (t, Ax − Az) = (t, Ap), where p = x − z. So far, all this shows is that if n vectors, orthogonal with respect to the matrix A can be found, then the conjugate direction algorithm will give solutions to the linear systems of the matrix. One can imagine applying a generalized form of Gram-Schmidt orthogonalization to an arbitrary set of linearly independent vectors. In fact Hestenes and Stiefel [HS52] show that A-orthogonalizing the n unit vectors in Rn and using them in CD leads essentially to Gaussian elimination. But this is no real solution since Gram-Schmidt requires O(n3 ) operations, and the search vectors, which will generally be dense even when the matrix is sparse, must be stored. The real advance to CD was made by Hestenes and Stiefel, who showed that A-orthogonal search vectors could be computed on the ﬂy. This is the conjugate gradient method.

11.2.6

The Method of Conjugate Gradients

Using the machinery that has been developed, it is a relatively easy task to describe the conjugate gradient (CG) algorithm as originally proposed by Hestenes and Stiefel
1

11.2 Conjugate Gradient

163

[HS52].c In going from steepest descent to conjugate directions, minimization along the residuals was replaced by minimization along the search vectors. So it makes sense to consider computing the search vectors iteratively from residuals. Suppose we make the ansatz p0 = r0 and pk+1 = rk+1 + βk+1 pk . (11.41) Can the coeﬃcients β be chosen so as to guarantee the A-orthogonality of the p vectors? Using (11.41), one has (pk , Apk+1 ) = (pk , Ark+1 ) + βk+1 (pk , Apk ). If one chooses βk+1 = − (pk , Ark+1 ) (pk , Apk ) (11.42)

As a ﬁnal touch, notice that the residuals can be calculated recursively; by induction rk+1 ≡ = = = The result of all this work is: Algorithm 6 Method of Conjugate Gradients Choose x0 . Put p0 = r0 = h−Ax0 . Then for k = 0, 1, 2, . . . (r ) αk+1 = (pkk ,rk k ) ,Ap xk+1 = xk + αk+1 pk rk+1 = rk − αk+1 Apk (11.43) (rk+1 ,rk+1 ) βk+1 = (rk ,rk ) pk+1 = rk+1 + βk+1 pk The α coeﬃcients are the same as in the CD algorithm, whereas the β coeﬃcients arise from the CG ansatz: p0 = r0 , pk+1 = rk+1 + βk+1 pk . From a computational point of view, note the simplicity of the algorithm. It involves nothing more than: • The inner product of a matrix and a vector; and only one per iteration since Apk can be calculated once and stored.
c The method was invented independently by M. Hestenes [Hes51] and E. Stiefel [Sti52], who later collaborated on the famous paper of 1952.

h − Axk+1 h − A(xk + αk+1 pk ) (h − Axk ) − αk+1 Apk rk − αk+1 Apk .

1

164 • The inner product of two vectors. • The sum of a vector and a scalar times a vector.

Iterative Linear Solvers

Since most of the calculation in CG will be taken up by the matrix-vector products, it is ideally suited for use on sparse matrices. Whereas a dense matrix-vector inner product takes O(n2 ) ﬂoating point operations, if the matrix is sparse, this can be reduced to O(nzero), where nzero is the number of nonzero matrix elements. To close this section a number of related details for the CD and CG algorithms will be shown. Lemma 4 (ri , pj ) = 0 for 0 ≤ j < i ≤ n (ri , pi ) = (ri , ri ) for i ≤ n (ri , rj ) = 0 for 0 ≤ i < j ≤ n − (rk+1 , rk+1 ) (rk+1 , Apk ) = (pk , Apk ) (rk , rk ) (pk , rk ) (rk , rk ) = (pk , Apk ) (pk , Apk ) (11.44) (11.45) (11.46)

(11.47) (11.48)

Proof. (11.44), (11.45), and (11.46) are by induction on n. (11.47) and (11.48) then follow immediately from this. Details are left as an exercise. Equation (11.45) arises interestingly if we ask under what circumstances the conjugate gradient residual is exactly zero. It can be shown that ri+1 = 0 if and only if (ri , pi ) = (ri , ri ). As a ﬁnal consideration, notice that although the gradient algorithms guarantee that the error z − xk is reduced at each iteration, it is not the case that the residual h−Axk is also reduced. Of course, the overall trend is for the residual to be reduced, but from step to step, relatively large ﬂuctuations may be observed. There are several generalizations of the basic Hestenes-Stiefel CG algorithm, known as residual reducing methods, which are guaranteed to reduce the residual at each step. For more details see Paige and Saunders [PS82] and Chandra [Cha78].

11.2.7

Finite Precision Arithmetic

The exact convergence implied by the Conjugate Direction Theorem is never achieved in practice with CG since the search vectors are computed recursively and tend to loose their A-orthogonality with time. CD methods were originally conceived as being “direct” in the sense of yielding the “exact” solution after a ﬁnite sequence of steps,
1

11.2 Conjugate Gradient

165

the number of which was known in advance. It soon became apparent that CG could be used as an iterative method. One can show that [Cha78]: x − xk
A≤

x − x0

A

√ 1− κ √ 1+ κ

2k

(11.49)

where κ ≡ λmax /λmin is the condition number of the matrix and x A ≡ (x, Ax). If √ √ the condition number is very nearly one, then (1 − κ)/(1 + κ) is very small and the iteration converges rapidly. On the other hand if κ = 105 it may take several hundred iterations to get a single digit’s improvement in the solution. But (11.49) is only an upper bound and probably rather pessimistic unless the eigenvalues of the matrix are well separated. For some problems, a comparatively small number of iterations will yield acceptable accuracy. And in any event, the convergence can be accelerated by a technique known as preconditioning. The idea behind preconditioning is to solve a related problem having a much smaller condition number, and then transform the solution of the related problem into the one you want. If one is solving Ax = h, then write this instead as Ax = h AC Cx = h Ax = h
−1

(11.50)

where A ≡ AC −1 and Cx ≡ x . To be useful, it is necessary that • κ(A ) κ(A)

• Cx = h should be easily solvable. In this case, CG will converge much more rapidly to a solution of A x = h than of Ax = h and one will be able to recover x by inverting Cx = x . Alternatively, one could write the preconditioned equations as Ax = h DAx = Dh Ax = h where DA ≡ A and Dh ≡ h . The most eﬀective preconditioner would be the inverse of the original matrix, since then CG would converge in a single step. At the other extreme, the simplest preconditioner from a computational standpoint would be a diagonal matrix; whether any useful preconditioning can be obtained from so simple a matrix is another matter. Between these two extremes lies a vast array of possible methods many of which are based upon an approximate factorization of the matrix. For example one could imagine doing a
1

(11.51)

166

Iterative Linear Solvers

Cholesky decomposition of the matrix and simply throwing away any nonzero elements which appear where the original matrix had a zero. In other words, one could enforce the sparsity pattern of A on its approximate factorization. For details on these “incomplete factorization” methods see [Man80],[Ker78], and [GvL83], for example.

11.2.8

CG Methods for Least-Squares

Conjugate gradient can be extended to the least squares solution of arbitrary linear systems. Solutions of the normal equations AT Ax = AT h are critical points of the function Ax − h
2

(11.52)

≡ ((Ax − h), (Ax − h)).

(11.53)

Note that AT A is always symmetric and nonnegative. The basic facts for least-squares solutions are these: if the system Ax = h is overdetermined, i.e., if there are more rows than columns, and if the columns are linearly independent, then there is a unique leastsquares solution. On the other hand, if the system is underdetermined or if some of the columns are linearly dependent then the least-squares solutions are not unique. (For a complete discussion see the book by Campbell and Meyer [CM79].) In the latter case, the solution to which CG converges will depend on the initial approximation. Hestenes [Hes75] shows that if x0 = 0, the usual case, then CG converges to the least-squares solution of smallest Euclidean norm. In applying CG to the normal equations avoid explicitly forming the products AT A. This is because the matrix AT A is usually dense even when A is sparse. But CG does not actually require the matrix, only the action of the matrix on arbitrary vectors. So one could imagine doing the matrix-vector vector multiplies AT Ax by ﬁrst doing Ax and then dotting AT into the resulting vector. Unfortunately, since the condition number of AT A is the square of the condition number of A, this results in slowly convergent iteration if κ(A) is reasonably large. The solution to this problem is contained, once again, in Hestenes’ and Stiefel’s original paper [HS52]. The idea is to apply CG to the normal equations, but to factor terms of the form AT h − AT Ax into AT (h − Ax), doing the subtraction before the ﬁnal matrix multiplication. The result is

[Cha78] shows that factoring the matrix multiplications in this way results in improved rounding behavior. For a more detailed discussion of the applications of CGLS see [HS52], [L¨u59], [Hes75], a [Law73], and [Bj¨75]. Paige and Saunders [PS82] present a variant of CGLS called o LSQR which is very popular since it is freely available through the Transactions on Mathematical Software. [PS82] also has a very useful discussion of stopping criteria for least squares problems. Our experience is that CGLS performs just as well as LSQR and since the CG code is so easy to write, it makes sense to do this in order to easily take advantage of the kinds of weighting and regularization schemes that will be discussed later.

11.2.9

Computer Exercise: Conjugate Gradient

Write a program implementing CG for symmetric, positive deﬁnite matrices. Consider the following matrix, right-hand side, and initial approximation:

n=6; A = Table[1/(i+j-1),{i,n},{j,n}]; h = Table[1,{i,nx}]; x = Table[0,{i,nx}]; To switch to ﬂoating point arithmetic, use i+j −1. instead of i+j −1 in the speciﬁcation of the matrix. The ﬁrst step is to familiarize yourself with CG and make sure your code is working. First try n = 4 or n = 5. On a PC, ﬂoating point arithmetic should work nearly perfectly in the sense that you get the right answer in n iterations. Now go to n = 6. You should begin to see signiﬁcant discrepancies between the exact and ﬂoating point answers if you use only n iterations. On other machines, these particular values of n may be diﬀerent, but the trend will always be the same.
1

168

Iterative Linear Solvers

Try to assess what’s going on here in terms of the numerical loss of A-orthogonality of the search vectors. You’ll need to do more than look at adjacent search vectors. You might try comparing p0 with all subsequent search vectors. Now see if you can ﬁx this problem simply by doing more iterations. If you get the right answer ultimately, why? What are the search vectors doing during these extra iterations. This is a subtle problem. Don’t be discouraged if you have trouble coming up with a compelling answer. Are the residuals monotonically decreasing? Should they be? What’s the condition number of the matrix for n = 6?

11.3
11.3.1

Practical Implementation
Sparse Matrix Data Structures

Clearly one needs to store all of the nonzero elements of the sparse matrix and enough additional information to be able to unambiguously reconstruct the matrix. But these two principles leave wide latitude for data structures.d It would seem that the more sophisticated a data structure, and hence the more compact its representation of the matrix, the more diﬃcult are the matrix operations. Probably the simplest scheme is to store the row and column indices in separate integer arrays. Calling the three arrays elem (a real array containing the nonzero elements of A), irow and icol, one has elem(i) = A(irow(i), icol(i)) i = 1, 2, . . . , N Z (11.55)

where N Z is the number of nonzero elements in the matrix. Thus if the matrix is


1 0 0 4    3 −2 0 0  0 0 −1 0 then elem = (1, 4, 3, −2, −1), irow = (1, 1, 2, 2, 3), and icol = (1, 4, 1, 2, 3). The storage requirement for this scheme is nzero real words plus 2 × nzero integer words. But clearly there is redundant information in this scheme. For example, instead of storing all of the row indices one could simply store a pointer to the beginning of each new row within elem. Then irow would be (1, 3, 5, 6). The 6 is necessary so that one knows how many nonzero elements there are in the last row of A. The storage requirement for this scheme (probably the most common in use) is nzero real words plus nzero + nrow integer words, where nrow is the number of rows in A. In the ﬁrst scheme, call it the
d A well-written and thorough introduction to sparse matrix methods is contained in Serge Pissanetsky’s book Sparse Matrix Technology [Pis84].

It is left as an exercise to construct similar operation within the row-pointer scheme. The matrix-vector inner product in the row-pointer scheme amounts to taking the inner product of each sparse row of the matrix with the vector and adding them up. If the rows are long enough, this way of doing things amounts to a substantial savings on a vectorizing computer since each row-vector inner product vectorizes with gatherscatter operations. At the same time, the long vector length would imply a substantial memory economy in this scheme. On the other hand, if the calculation is done on a scalar computer, and if memory limitations are not an issue, the full-index scheme is very eﬃcient in execution since partial sums of the individual row-vector inner products are accumulated simultaneously. For the same reason, a loop involving the result vector will be recursive and hence not easily vectorized.

11.3.2

Data and Parameter Weighting

For inverse problems one is usually interested in weighted calculations: weights on data both to penalize(reward) bad(good) data and to eﬀect a dimensionless stopping criterion such as χ2 , and weights on parameters to take into account prior information on model space. If the weights are diagonal, they can be incorporated into the matrix–vector multiply routines via: ∀k y(icol(k)) = y(icol(k)) + elem(k) ∗ x(irow(k)) ∗ W 1(irow(k)) (11.58)

for column or parameter weighting. Here, W 1 and W 2 are assumed to contain the diagonal elements of the weighting matrices.

11.3.3

Regularization

Just as most real inverse calculations involve weights, most real inverse calculations must be regularized somehow. This is because in practice linear least squares calculations usually involve singular matrices or matrices that are numerically singular (have very small eigenvalues). Regularization is the process by which these singularities are
1

170

Iterative Linear Solvers

tamed. We saw two diﬀerent examples of regularization in Chapter 5. The ﬁrst was truncating the SVD. We can throw away zero or small singular values and that regularizes the problem. But we also found that in the presence of a model null space it was useful to be able to penalize the size of the solution as well as the data misﬁt. In other words we replaced the minimization problem min with min Ax − h
2

Ax − h +

2

(11.60)
2

x

.

(11.61)

The ﬁrst term is the data misﬁt, and the second is the “regularization” term. As shown here, the two aspects of the minimization (make the data misﬁt small, make the model norm small) get equal weight. Now let us go two steps beyond this. First, let us introduce a fudge factor λ to control the tradeoﬀ between the two terms. Next, let us consider the possibility of minimizing not the norm of the model itself, but the norm of some linear function of the model Rh. min Ax − h 2 +λ Rx 2 (11.62) If R = I, then we’re back to our familiar regularization. But now suppose that R ≡ ∂ n , n = 0, 1, 2, . . . and ∂ n is an n − th order discrete diﬀerence operator. In this case the term Rx 2 penalizes the slope, roughness, or higher order derivative of the model. Penalizing roughness would be useful if we want a smooth solution. The “normal equations” associated with this generalized objective function, obtained by setting the derivative of (11.62) equal to zero, are (AT A + λRT R)x = AT h. (11.63)

This sort of regularization is straightforward to implement in a sparse matrix framework by augmenting the matrix with the regularization term: ˜ A≡ √A λR .

From this you can tell right away that R must have the same number of columns as A. But in principle it can have any number of rows. For example, we might use 1 −1 0 · · ·    0 1 −1 · · ·      . . R=  .
     

which is singular but has the same null space as the continuous derivative operator; i.e., it maps constant vectors into 0. We must also augment the right hand side, with a number of zeros equal to the number of rows in the regularization matrix. We write the augmented right hand side ˜ y≡ y 0 .

˜ ˜ ˜ ˜ ˜ ˜ Since AT y = AT y and AT A = AT A + λRT R, the least squares solutions of Ax = y satisfy (AT A + λRT R)x = AT h. So to incorporate any regularization of the form of (11.62) all one has to do is augment the sparse matrix. Most commonly this means either damping, in which case R is diagonal, or second-diﬀerence smoothing, in which case R is tridiagonal.

11.3.4

Jumping Versus Creepinge

The pseudo-inverse A† itself has something of a smoothness condition built in. If the matrix A has full column rank and the number of rows is greater than or equal to the number of columns (in which case the system is said to be overdetermined) then the least squares solution is unique. But if the system is underdetermined, the least squares solution is not unique since A has a nontrivial null space. All of the least squares solutions diﬀer only by elements of the null space of A. Of all of these, the pseudo-inverse solution is the one of smallest norm. That is, x† ≤ x for every x such that AT Ax = AT y, as we saw in Chapter 5. This means, for example, that in a nonlinear least squares problem, where we perturb about a reference model and compute this perturbation at each step by solving a linear least squares problem, then the size of the steps will be minimized if the pseudo-inverse is used. This has led to the term “creeping” being used for this sort of inversion. On the other hand, if at each nonlinear step we solve for the unknown model directly, then using the pseudo-inverse will smallest norm will enforce the smalleste norm property on the model itself, not the perturbation of this model about the background. This is called “jumping” since the size of the change in the solution between nonlinear iterations is not constrained to be small. The terms creeping and jumping are due to Parker [Par94].
e

This section and the next are taken from [SDG90].
1

172

Iterative Linear Solvers

This point merits a brief digression since the eﬀects of damping or smoothing will be diﬀerent according as one is doing jumping or creeping. Suppose the nonlinear inverse problem is: given y, ﬁnd x such that y − F (x) is minimized in some sense. Expanding the forward problem F to ﬁrst order in a Taylor series about some model x0 gives y = y0 + F (x0 )(x − x0 ) (11.64)

Expressing the initial model in terms of its components in the row space and null space of A, x0 = xrow + xnull (11.68) 0 0 and noting that xrow = A† Ax0 0 then xj = xrow + A† (y − y0 ) 0 and (11.67) becomes xj − xc = −xnull . 0 (11.71) Thus, the creeping and jumping solutions diﬀer by the component of the initial model that lies in the null space of A: some remnants of the initial model that appear in xc are not present in xj . Only if A is of full column rank (giving A† A = I ) will the two solutions be the same for any initial guess. In the next sections it will be seen that this analysis must be modiﬁed when regularization is employed. (11.70) (11.69)

11.3.5

How Smoothing Aﬀects Jumping and Creeping

In the absence of regularization, the jumping and creeping solutions diﬀer only by the component of the initial model in the null space of the Jacobian matrix. Regularization changes things somewhat since the matrix associated with the regularized forward
1

As in (11.71), the diﬀerence between the two solutions depends on the initial model. But when smoothing is applied, the creeping solution possesses components related to the slope of x0 (ﬁrst diﬀerence smoothing) or to the roughness of x0 (second diﬀerence smoothing) which are not present in the jumping solution. An important corollary of this result is that for smooth initial models, jumping and creeping will give the same results when roughness penalties are employed to regularize the calculation. Examples illustrating the comparative advantages of jumping and creeping are contained in [SDG90].
f If the columns of the regularization operator are linearly independent, then the columns of the augmented matrix are too.

1

174

Iterative Linear Solvers

11.4

Sparse Singular Value Calculationsg

The singular value decomposition is one of the most useful items in the inverter’s toolkit. With the SVD one can compute the pseudo-inverse solution of rectangular linear systems, analyze resolution (within the linear and Gaussian assumptions), study the approximate null space of the forward problem, and more. The now classical GolubReinsch approach to SVD [GR70] begins by reducing the matrix to block bidiagonal form via a sequence of transformations known as Householder transformations. The Householder transformations annihilate matrix elements below the diagonal, one column at a time. Unfortunately, after each transformation has been applied, the sparsity pattern in the remaining lower triangular part of the matrix is the union of the sparsity pattern of the annihilated column and the rest of the matrix. After a very few steps, one is working with nearly full intermediate matrices. This makes conventional SVD unsuitable for large, sparse calculations. On the other hand, for some problems, such as studying the approximate null space of the forward problem, one doesn’t really need the entire SVD; it suﬃces to compute the singular vectors associated with the small singular values (“small” here is deﬁned relative the level of noise in the data). Or perhaps from experience one knows that one must iterate until all those eigenvectors down to a certain eigenvalue level have been included in the solution. Conventional SVD gives no choice in this matter, it’s all or nothing. In this section we shall consider the use of iterative methods such as conjugate gradient for computing some or all singular value/singular vector pairs.

In other words, just by doing CG one gets a symmetric tridiagonalization of the matrix for free. Needless to say, computing the eigenvalues of a symmetric tridiagonal matrix is vastly simpler and less costly than extracting them from the original matrix. For rectangular matrices, simply apply the least squares form of CG and use the α and β scale factors in (11.83) and (11.84), to get a symmetric tridiagonalization of the normal equations. Then, just take their positive square roots to get the singular values. The calculation of the eigenvalues of symmetric tridiagonal matrices is the subject of a rather large literature. See [Sca89] for details. The following example illustrates the idea of iterative eigenvalue computation. We will 1 consider the Hilber matrix, whose i − j element is i+j+1 . This matrix arises in the theory of approximation and is known to be highly ill-conditioned.h The matrix in question is an eighth-order Hilbert matrix:
h A simple explanation for this was contributed to the Usenet news group sci.math by Zdislav V. Kovarik. The idea is you can interpret the i − j element as the inner product of xi and xj on the 1 interval [0, 1]. Now, the cosine of the angle between xk and x( k + 1) is just 2∗k+2 . So you can see that as k increases, this matrix. which consists of the scalar products of these almost linearly dependent vectors, is bound to be nearly singular.

1

176 nx =8; A = Table[1/(i+j-1.),{i,nx},{j,nx}];

Iterative Linear Solvers

The condition number of this matrix is 1010 . The exact solution to the system Ax = y, where y consists of all ones is: (−8, 504, −7560, 46200, −138600, 216216, −168168, 51480). After just 5 iterations, using 16 digits of precision, CG produces the following solution: (0.68320, −4.01647, −127.890, 413.0889, −19.3687, −498.515, −360.440, 630.5659) which doesn’t look very promising. However, even after only 5 iterations we have excellent approximations to the ﬁrst 4 eigenvalues. The progression towards these eigenvalues is illustrated in the following table, which shows the fractional error in each eigenvalue as a function of CG iterations. Even after only one iteration, we’ve already got the ﬁrst eigenvalue to within 12%. After 3 iterations, we have the ﬁrst eigenvalue to within 1 part in a million and the second eigenvalue to within less than 1%. Eigenvalue Fractional error in CG-computed eigenvalues 1.6959389 0.2981252 0.0262128 0.0014676 0.0000543 Iteration 0.122 1 0.015 0.52720 2 1.0 10−5 0.006 1.284 3 −12 −7 9.0 10 1.9 10 0.002 1.184 4 0.0 7.3 10−15 1.13 10−8 8.0 10−4 1.157 5

11.4.2

Finite Precision Arithmetic

Using CG or Lanczos methods to compute the spectrum of a matrix, rather than simply solving linear systems, gives a close look at the very peculiar eﬀects of rounding error on these algorithms. Intuitively one might think that the main eﬀects of ﬁnite precision arithmetic would be a general loss of accuracy of the computed eigenvalues. This does not seem to be the case. Instead, “spurious” eigenvalues are calculated. These spurious eigenvalues fall into two categories. First, there are numerically multiple eigenvalues; in other words duplicates appear in the list of computed eigenvalues. Secondly, and to a much lesser extent, there are extra eigenvalues. The appearance of spurious eigenvalues is associated with the loss of orthogonality of the CG search vectors. A detailed explanation of this phenomenon, which was ﬁrst explained by Paige [Pai71] is beyond the scope of this discussion. For an excellent review see ([CW85], Chapter 4). In practice, the duplicate eigenvalues are not diﬃcult to detect and remove. Various strategies have been developed for identifying the extra eigenvalues. These rely either on changes in the T matrix from iteration to iteration (in other words, on changes in Tm as m increases), or diﬀerences in the spectra between T (at a given iteration) and the principle submatrix of T formed by deleting its ﬁrst row and column. An extensive
1

11.4 Sparse SVD

177

discussion of the tests used for detecting spurious eigenvalues is given by [CW85]. It is also not obvious how many iterations of CG are necessary to generate a given number of eigenvalues. At best it appears that for “large enough [number of iterations] m, every distinct eigenvalue of A is an eigenvalue of Tm ”—the Lanczos phenomenon[CW80]. On the other hand, the spurious eigenvalues crop up because one has been content to let the search vectors lose orthogonality: computing a lot of iterations, throwing away a lot of duplicate eigenvalues, and relying on the Lanczos phenomenon to assure that eventually one will calculate all of the relevant eigenvalues. The examples in [CW85] and the example that will be shown presently would seem to indicate that that is not an unreasonable goal. On the other hand, many (perhaps most) practitioners of the Lanczos method advocate some sort of partial or selective reorthogonalization. In other words, orthogonalize by hand the current search vector with respect to the last, say, N vectors, which then must be stored. Some examples of reorthogonalized Lanczos are given by [Par80]. It is diﬃcult to do justice to the controversy which surrounds this point; suﬃce it to say, whether one uses reorthogonalized methods or not, care must be taken to insure, on the one hand, that spurious eigenvalues are not mistakenly included, and on the other, that reorthogonalization is suﬃciently selective that the speed of the method is not completely lost. Here is a simple example of the use of CG-tridiagonalization from [Sca89]. The problem is a small, 1500 or so rays, travel time inversion of reﬂection seismic data. The model has about 400 unknown elastic parameters. In the table below are listed the ﬁrst 40 singular values of the Jacobian matrix computed with an SVD (on a Cray X-MP) and using Conjugate Gradient. Duplicate singular values have been removed. The results are extremely close except for the three spurious singular values 7, 24, and 38. In all I was able to compute about half of the nonzero singular values without diﬃculty. Most of these were accurate to at least 6 or 7 decimal places.

Finally, we point out a clever result of Hestenes which seems to have been largely ignored. In the paper [Hes75] he proves the following. Let r be the rank of A an arbitrary matrix, and let p and q be the CGLS search vectors, and let x0 = 0. Then A† = p0 p0 p1 p1 pr−1 pr−1 + +··· AT (q0 , q0 ) (q1 , q1 ) (qr−1 , qr−1 ) (11.85)

is the generalized pseudo-inverse of A. A generalized pseudo-inverse satisﬁes only two of the four Penrose conditions, to wit: A† AA† = A† AA† A = A To illustrate this result, consider the following least squares problem:
              

(11.86) (11.87)

1 2 −4 5 −1 3 2 −7

x y

=

5 6 5 −12

  . 

The column rank of the matrix is 2. It is straightforward to show that AT A Therefore the pseudo-inverse is A† = AT A
−1 −1

3. Prove Lemma 4. 4. With steepest descent, we saw that in order for the residual vector to be exactly zero, it was necessary for the initial approximation to the solution to lie on one of the principle axes of the quadratic form. Show that with CG, in order for the residual vector to be exactly zero we require that (ri , pi ) = (ri , ri ) which is always true by virtue of Lemma 3.

where A is the derivative of the forward problem (an n by m matrix), m is a vector of unknown model parameters, and d contains the observed data. m is a vector in R m . However, it represents a discretization of the model slowness s(r), which is a scalar function deﬁned on a closed subset Ω of RD , D ∈ (1, 2, 3 . . .). It will be assumed that the set of all possible models lies in some linear function space M. It is useful to introduce an orthonormal basis of functions (we will use our old friends the pixel functions) which span the model space M. Suppose that Ω is completely covered by m closed, convex, mutually disjoint sets (cells) σ ∈ R D : Ω = ∪ σi such that σi ∩ σj = if i = j. The basis functions are then deﬁned to be hi (r) = if r ∈ σi νi 0 otherwise
−1/2 −1/2

is made where νi is the volume of the ith cell. The choice of the normalization νi −1/2 can be to remove bias introduced by cell size. If a constant cell size is adopted, νi replaced with 1. Given the deﬁnition of hi , it is clear that
Ω

hi (r)hj (r)dD r = δij .
∞ i=1

Thus an arbitrary function can be written as an expansion in hi m(r) = mi hi (r) ≡ m · h(r).
1

(12.2)

184

More on the Resolution-Variance Tradeoﬀ

In practice this sum will usually be truncated at a ﬁnite number of terms. For inﬁnite dimensional vectors it is more traditional to write the inner product as (m, h(r)) but we will continue to use the dot notation. A basic took of BG theory is the point spread function A. The point spread function (PSF) is deﬁned in a formal manner by noting that a local average of the model m(r) can be obtained at any point r0 by integrating the model and a locally deﬁned unimodular function m(r0 ) = A(r, r0 )m(r)dD r. (12.3)
Ω

Unimodular means that the function integrates to 1.
Ω

A(r, r0 )dD r = 1.

Also, it is assumed that A(r, r0 ) ∈ M for each r0 ∈ Ω and that the support of A is concentrated at the point r0 . Naturally, the more accurately the model is determined at each point, the more closely the PSF resembles a delta function – at that point. So estimating the PSF is equivalent to estimating a local average of the model. The more delta function-like is the PSF, the more precise our estimate of the model. Like any other function in M, the PSF can be expanded in terms of hi . A(r, r0 ) = Thus (12.2) and (12.3) imply that m(r0 ) =
∞ i=1 ∞ i=1

ai (r0 ) hi (r) ≡ a(r0 ) · h(r).

ai (r0 ) mi = a (r0 ) · m.

(12.4)

It is clear that one can construct a PSF which will yield a local average of the model – any approximation to a delta function will do. Unfortunately, there is a tradeoﬀ between the sharpness of the PSF and the variance, or RMS error, of the solution. To show this, BG assume that the local average of the model is a linear function of the data n m(r0 ) =
i=1

bi (r0 ) di = b (r0 ) · d = (b (r0 ) , Am + e)

(12.5)

where b is to be determined. For the moment let’s neglect the noise–for zero mean noise we can just take expectations. Comparing (12.4) and (12.5) one sees that the expansion coeﬃcients of the PSF are simply a (r0 ) = AT b (r0 ) (12.6)

Now, how one measures the “width” of the PSF is largely a matter of taste. Nolet [Nol85] makes the following natural choice W (r0 ) = cD
Ω

A (r, r0 )2 |r − r0 |D+1 dD r
1

12.2 Using the SVD

185

where cD is a scale factor chosen to make W have as simple a form as possible. For example, Nolet chooses c3 = 56π/9. Plugging the deﬁnition of the pixel functions and of the PSF into (12.6) it follows that W (r0 ) =
i,j,k

fi Aji Aki bj (r) bk (r0 )

(12.7)

where fi = c D
Ω

|r − r0 |D+1 hi (r)2 dD r ≈ cD |ˆi − r0 |D+1 r

and where ri is the centroid of the ith cell. Equation (12.7) shows how the width of the ˆ PSF depends on the b coeﬃcients. Now all one needs is a similar expression for the error of the average model value m (r0 )
n

σ 2 = Var m (r0 ) =
i,j=1

bi (r0 ) bj (r0 ) Cov(di , dj ) = b (r0 ) · b (r0 ) .

The last equality follows since if one assumes that the data are uncorrelated, then weights can always be chosen such that Cov(di , dj ) = δij . Thus, it has been shown that both the width of the PSF and the variance of the solution depend on b; Thus one cannot tighten up the PSF without aﬀecting the variance. The solution, at least formally, is to introduce a tradeoﬀ parameter, say w, and jointly minimize J(w, r0 ) ≡ W (r0 ) + w 2 σ 2 (r0 ). This last problem is straightforwardly solved but note that to compute the coeﬃcients b which jointly minimize the variance and the width of the PSF requires the solution of a (large) least squares problem at each point in the model where the resolving power is desired. For large, sparse operators A, a far more eﬃcient approach would be using the conjugate gradient methods outlined in Chapter 11.

12.2

Using the SVD

Now let us look at this tradeoﬀ for a ﬁnite-dimensional problem using the SVD. Let A be the forward modeling operator, now assumed to map Rm to Rn : d = Am + e ˆ where e is an n-vector of random errors. The least squares estimated model m is given † † by A d, where A is the pseudo-inverse of A. ˆ ˆ ˆ The covariance of m is E[mmT ].a We can get a simple result for this matrix using the singular value decomposition. The singular value decomposition of A is A = U ΛV T
a

where U is an orthogonal matrix of “data” eigenvectors (i.e., they span R n ) and V is an orthogonal matrix of “model” eigenvectors (they span R m ). Λ is the n × m diagonal matrix of singular values λi . The pseudo-inverse of A is A† = V Λ−1 U T where Λ−1 where denotes the m × n diagonal matrix obtained by inverting the nonzero singular values. To keep things simple, let’s assume that the covariance of the data errors is just the identity matrix. This will let us look at the structure of the covariance ˆ of m as a function of the forward operator alone. It is easy to see that in this case ˆ ˆ ˆ Cov(m) = E[mmT ] = A† Cov(d)A† = V Λ−2 V T =
T m i=1 T λ−2 vi vi . i

The last term on the right is the sum of the outer products of the columns of V (these are the model space eigenvectors). So the covariance can be seen as a weighted projection operator onto the row space of A, with weights given by the inverse-square of the singular values. ˆ With this it is not diﬃcult to see that the j-th diagonal element of Cov(m), which is the variance of the j-th model parameter is
m

ˆ Var(mj ) =
i=1

λ−2 (vi )2 i j

where (vi )j is the j-th component of the i-th eigenvector. If the rank of A is less than m, say r, then all of the sums involving the pseudo-inverse are really only over the r eigenvectors/eigenvalues. In particular
r

ˆ Var(mj ) =
i=1

λ−2 (vi )2 . i j

This is because A = U ΛV T = Ur Λr VrT where the subscript r means that we have eliminated the terms associated with zero singular values. Now suppose we decide not to use all the r model eigenvectors spanning the row space of A?b For example we might need only p eigenvectors to actually ﬁt the data. Let ˆ ˆ us denote by mp the resulting estimate of m (which is obviously conﬁned to the pm dimensional subspace of R spanned by the ﬁrst p model singular vectors):
p

ˆ m ≡

p

vi
i=1

uT d i λi

where ui is the i-th column of U (i.e., the i-th data eigenvector). Using the result above ˆ for the variance of the j-th component of m we can see that
p

ˆj Var(mp )

=
i=1

λ−2 (vi )2 . i j

b Remember that if a vector is in the null space of a matrix, then it is orthogonal to all the rows of the matrix. Hence the row space and the null space are orthogonal complements of one another.

1

BIBLIOGRAPHY

187

ˆ It follows that the variance of the j-th component of mp is monotonically nondecreasing with p. So, while we can formally decrease the variance by using fewer eigenvectors, we end up with a less resolution because we won’t have enough structure in the remaining eigenvectors to characterize the model. In general we cannot compute the bias for an estimate of the true model without taking into account the discretization, but let’s neglect this for the moment and assume that A represents the exact forward problem and that the true model lies within R m . The bias ˆ of mr is the component of the true model in the row space of A, assuming zero-mean errors.c So, apart from the component of the true model in the row space of A, the bias ˆ of mp is r uT d ˆ ˆ ˆ bias(mp ) = E[mp − mr ] = vi i i . λ i=p+1