NATIONAL BUREAU OF STANDARDS REPORT

2365

"MONTE CARLO" METHODS EOR THE ITERATION OE LINEAR OPERATORS

by

J. H. Curtiss

U. S. DEPARTMENT OF COMMERCE

NATIONAL BUREAU OF STANDARDS
U. S. DEPARTMENT OF COMMERCE

Sinclair Weeks, Secretary

NATIONAL BUREAU OF STANDARDS

A. V. Astin, Director THE NATIONAL BUREAU OF STANDARDS

The scope of activities of the National Bureau of Standards is suggested in the following listing of

the divisions and sections engaged in technical work. In general, each section is engaged in special*

ized research, development, and engineering in the field indicated by its title. A brief description

of the activities, and of the resultant reports and publications, appears on the inside of the nack

cover of this report.

Electricity. Resistance Measurements. Inductance and Capacitance. Electrical Instruments. Magnetic Measurements. Applied Electricity. Electrochemistry.

Optics and Metrology. Photometry and Colorimetry. Optical Instruments. Photographic

Technology. Length. Gage.

Heat and Power. Temperature Measurements. Thermodynamics. Cryogenics. Engines and

Lubrication. Engine Fuels. Cryogenic Engineering.

Atomic and Radiation Physics. Spectroscopy. Radiometry. Mass Spectrometry. Solid State

Physics. Electron Physics. Atomic Physics. Neutron Measurements. Infrared Spectroscopy.

Nuclear Physics. Radioactivity. X-Rays. Betatron. Nucleonic Instrumentation. Radiological Equipment. Atomic Energy Commission Instruments Branch.

Chemistry. Organic Coatings. Surface Chequstry. Organic Chemistry. Analytical Chemistry.

Inorganic Chemistry. Electrodeposition. Gas Chemistry. Physical Chemistry. Thermochemistry. Spectrochemistry. Pure Substances.

Mechanics. Sound. Mechanical Instruments. Aerodynamics. Engineering Mechanics. Hydraulics. Mass. Capacity, Density, and Fluid Meters. Organic and Fibrous Materials. Rubber. Textiles. Paper. Leather. Testing and Specifications. Polymer Structure. Organic Plastics. Dental Research.

Metallurgy. Thermal Metallurgy. Chemical Metallurgy. Mechanical Metallurgy. Corrosion.

Mineral Products. Porcelain and Pottery. Glass. Refractories. Enameled Metals. Concreting Materials. Constitution and Microstructure. Chemistry of Mineral Products.

Building Technology. Structural Engineering. Fire Protection. Heating and Air Conditioning. Floor, Roof, and Wall Coverings. Codes and Specifications.

Applied Mathematics. Numerical Analysis. Computation. Statistical Engineering. Machine

Development.

Electronics. Engineering Electronics. Electron Tubes. Electronic Computers. Electronic

Instrumentation.

Radio Propagation. Upper Atmosphere Research. Ionospheric Research. Regular Propagation Services. Frequency Utilization Research. Tropospheric Propagation Research. High

Frequency Standards. Microwave Standards.

Ordnance Development. These three divisions are engaged in a broad program of research

Electromechanical Ordnance, and development in advanced ordnance. Activities include

Ordnance Electronics. basic and applied research, engineering, pilot production, field testing, and evaluation of a wide variety of ordnance materiel. Special skills and facilities of other

NBS divisions also contribute to this program. The activity is sponsored by the Department of

Defense.

Missile Development. Missile research and development: engineering, dynamics, intelligence,

instrumentation, evaluation. Combustion in jet engines. These activities are sponsored by the

Department of Defense. • Office of Basic Instrumentation • Office of Weights and Measures.
NATIONAL BUREAU OF STANDARDS REPORT

NBS PROJECT

i 101-10-5100 March 19? 1953

NBS REPORT

2365

"MONTE CARLO'8 METHODS FOR THE ITERATION OF LINEAR OPERATORS

by

Jo Curtiss

Harvard University

The National Bureau of Standards

The publication, repr Approved for public release by the

unless permission is oi Director of the National Institute of

25 ,d,c. Suchpermi Standards and Technology (NIST)

cal ly prepared If tha on October 9, 2015

part, Is prohibited

dards, Washington

rt has been speclflirt for Its own use,

"MONTE CARLO1

' METHODS FOR THE ITERATION OF LINEAR OPERATORS

by Jo H® Curtiss

Harvard University aod

The National Bureau of Standards

1 o Int rodue tion 0 A purely formal description of the type of

problem to be dealt with in this paper is as follows 0 Let

c = c(x) and Uq = Uq(x) be real-valued functions defined on a coordinate space R, which may b© multidimensional » Let

L = L(f) be a linear transf ©raa tion defined ©o the space of all

real-valued functions f whose arguments belong to R® Required,

to calculate the sequence ©f functions u^ ? Ug? 000 ? defined by

the recursion formula

( 1 o 1 ) U

1 ” L ^ “ Ojl 5 ^J o o o

This problem arises in many contexts in both pure aod

applied mathematics a Much attention has recently been focused

on the numerical aspects of it by the nuclear physicists, because the recursion formula (id) is obtained when time-dependent diffusion and transport problems are formulated in a discrete form*, The methods to be discussed here for estimating

the solution are an outgrowth ©f a novel stochastic attack

suggested during the late war by Von Nuemann and Ul®m in connection with diffusion problems®* Their idea was t© bypass

th© mathematical formulation and set up a computing

procedure with various random decisions in it which more or

less closely imitated the physical phenomenon under study®

^Various practical aspects ©f the stochastic estimation are

presented in [l4].
2

This type of approach to distribution problems has long

been known to statisticians under the name of '“model sampling."

The physicists have thought up a new name for it that seems

likely to sticks the "Monte Carlo Method."

The formal solution of (1.1) is the truncated Neumann

series

(1.2) ~ c + L(c) + L^(e) + . . .+ ^ (c) + L^(uq), N > 0,

K where L means the K-th iterate of L. Especial interest, of

course, lies in the case in which the corresponding infinite

series converges. If it does, it represents a function u which

satisfies the equation

(1.3) u = L(u) + c o

The error estimate is then provided by

(1.4-) - u = L^(uq - u).

Equation (1 »2) shows that it makes no difference in the long

run how Uq was chosen, but Eg. (1 .4) suggests that the nearer

Uq is to u, the faster the convergence will be.

The stochastic approach to estimating the solution

of (1.1) will now be described in correspondingly general terms.

Consider the space R over which the functions c, Uq, u^ , . . .

,

are defined. Let x be a point in it at which it is desired

to calculate u^. A probability distribution is now set up on

each of the spaces R, R X R , R X R X R,o.., where R X R denotes the Cartesian product of R and R. Random variables
3

Z*qs*\ Q Z 4 6)> ®^2t ?

® ° ® ;e re t spaces®

The chain of probability distributions and of random variables

is such that the sequence of conditional mean values*

vQ

(x) = E(Z

Q

z

n

= Uq(x)) w

1: (x) = E(Z^ z,0

u“0n(x) j

,

v

2

ix,

'o - Uq(x) satisfies (1.1). Since Vq =

it follows that v^(x) = n^Cx),

The computational problem then becomes one of calculating repeated realizations of and combining them into an

appropriate statistical estimator ©f v^„

The stochastic method is particularly well adapted to

the case in which the value of u^(x) or u(x) is to be estimated

at only ©n© point x®

In the case in which the interest lies in estimating

the solution of (1®3) ? it is quite possible to carry the sequence 9 Z 9 , 0S8 ©n to infinity and define a random variable Z

on the infinite product-space EX EX ®oo But this has no significance for actual practice^ and curiously ®n@ugh ? it turns

out to be theoretically disadvantageous® That arises from the

through absolutely convergent integrals and sums® This in turn

*We shali^lsse the symbol "FrTaj b) to denote the conditional

probability of the event a 9 given that b has occurred , the

symbol E(Y) to denote the mean value of the random variable

Y 9 and’ the symbol E(Y

j

b) to denote the mean value of the

conditional distribution of Y 9 given that b has occurred®
4

imposes some irrelevant restrictions on L if is to be identified with u 9 at least in the cases t© be considered in this

paper** Therefore the attitude here in the ^steady state”

situation will be that a suitably large but finite N will be

chosen once and for all and held fixed during the sampling*

The mathematical material preceding Section 7 is in

the main a rearrangement and mathematical formulation of known

procedures^ presented so as to show up their relationships*

The method of error analysis proposed in Section 4 and most of

the material in Sections 7? 8 § 9? and 10 are believed to be new*

However^ a good deal ©f work on the Monte Carlo Method has been

^published” in privately circulated , sometimes classified 9 reports ? and one can never be quite sure of a priority under such

eireums taoces

A word of caution to the reader may be in order* The

Monte Carlo method as a computational procedure has had its

chief successes in problems which had natural stochastic bases

and which were at the same time so complicated that they were

inaccessible to ordinary analytic or numerical methods* This

paper makes no pretense of putting the method into competition

with the standard numerical practices 9 especially for the

*As applied to the solution of simultaneous linear equations,

the methods of the present paper were anticipated by Forsythe

and Leibler [10] and Was aw [16]* Both of these papers deal

directly with the infinite product space* See also Curtiss

[ 4] o
5

simpler type of problem for which good methods already exist*,

The idea here is merely t© present some theory which may be of

some interest for itself ©lone ? and which unifies and clarifies

certain of the Monte Carlo devices which have been proposed ? and which lays the groundwork for further numerical experimentation aimed at investigating the limits of usefulness of the

method for non=»§ tochas tie problems®

We propose to stop short of describing the practical

computational details® Except for the process of generating

the necessary random elements 9 they can all be supplied easily

by the experts in computing® The process of correctly generating

the random elements and making sure that they are appropriately

random could be the subject of another article at least as long

as this ©n@o Iff it were well d@ne § it might contain substantial contributions to the philosophy of probability®
6

2. Specialization of the problem* The foregoing introduction

has been very vague as to the nature of the operator L 0 Actually* the stochastic method to be presented here seems to be

automatically confined to operators L of a Lebdsgue-S tie! t jes

type *

( 2 . 1 ) L(u) u(y) dyk(x *y)

,

in which k is of bounded variation and h has integrability

properties which will permit the iteration* We shall not pursue

the question of generability any further* however* Instead we

shall present the theory for two special eases? that in which

L is an ordinary integral transform**

(2*2) L(u) = K

J

h(x*y) u(y) dy* ' R

and that in which L is a matrix and u is a vector*

The second case is the special case of (2*1) in which

k is a step-function* Most of the exposition except for that

in Section 4- will be directed toward this case* It is the more

fundamental one in numerical analysis and the majority of the

results can be carried over so readily to the continuous operator (2*2) that no comment on the matter will be necessary***

It is convenient to introduce new notation for the

matrix case* The space R will be thought of in this case as

*The integral can be taken in either the Riemann or Lebesque

sense in the sequel*

**See Cutkosky [5 J*
consisting of a finite discrete set of points s x ?5 ,,. ,x n

, The function u = u(x) will be represented by the vector

u = (o^ , 02 ?o..jU

jj ) = (u(x^},,.,,u{x

o

)), Similarly «, we write

U®, — ( U«,(Xj 3 } UmyiX/^} J O e O 5 U M(X„ j ) “ ( UjUi 5 ?

® ® ® ? N— OjJjSj

‘N N 1 ' ? N 2

? ? N n and c = (c(x.j), c (xg) ? • « • »c (x Q# , .-

1

- - - /

n the function h(x ? y) is represented by the matrix

N1 ? N2 ?

° • * ? Nn‘

C -a <J>Co9©©©9^.

H = [h(x

i? y^ ij

uation (1 ®1} becomes in matrix notation

(2o3)

N+'

Hu^ + c o

The problem connected with this eq

more components of the vector e^ ?

nation is to calculate one or

given H

^ uQ

=(

u

Q1 ,Uq 2

, . . • ?uQn )

The series solution (1 ®2) becomes

(2.^) u N

= c + He + fi 2

® +ooo+ H

N=1

c + E

N

nQ

= (I=fi

N

) (I-H)”1 c

+ H^Uq ? M = 1 j2 ? o o ® ®

Here the symbol 1 stands for the unit matrix® Of course the

third member of the equation can be written down only if I - H

if non-singular

®

It is well known that the necessary and sufficient

m

condition for lira H = 0 is that all the eigenvalues of fi N

’~J> w

(that is s the roots ©f the deterirtinaetal equation |Xl — H J = 0

)

must lie within the unit circle of the complex plane®* Iff this
8

ii the case—and we shall always assume that it is whenever we

are discussing the situation as N~“£- oo «=» then 9 (I=H) ^ exists^

=*1 and u = lira uN

= (I - H) c, which satisfies the linear equaN-* 00 N

tions

(2o5) u = Hu + ®o

3o The solution of Au = bo

We digress for a moment, her© to note the relationship

between (2o3)> (2o**), and (2o5)j and the important problem of

solving the system off linear equations

( 3 o 1 } Au = b

,

where b = (b.0 » e o 9 b ) is an arbitrary vector,, Choose the matrix

and also a new one M s so that H + MA = I 9 and choose c = Mb 0 Then (2 0 5) reduces to MAu = Mbo Ii M is non-singular 9 then

this system is precisely equivalent to (3°1) in the sense that

each solution u ©ff (2„5) ii a solution off (3»1) 9 and vice versa,.

If there is more than ©n@ solution t® (3°1) 9 then A is singular^

and in this case it Is easily checked that one of the eigenvalues of H = I = MA ii unity,, This means that in (

2

0 3)

and (2o 4) cannot converge to the solution,. We shall therefore

exclude this ease and assume that A is non- singular,,

With M and A both n@n=siogular 9 and with H having its

eigenvalues all inside the unit circle^ (2o3) becomes what is

known in the theory of linear algebraic systems as a stationary
9

linear iterative process,, ©r Wittmeyer process., for solving

Au - bo There are obviously an infinite number ©f ways of

ehooiiog M and H so that the ©©editions are fulfilled® One

standard method is t© split up A into the difference of two

matrices V and W 9 where V is easily inverted® That is 9 let

A = V - W 0 Then take M = 9 H = V~^W® It is easily checked

that H + MA = I, and it is not hard to arrange things so that

the eigenvalues ©f H are all sufficiently small in meduluso*

The calculation ©f A

=1

9 if that is the problem,, can of

course be accomplished by specializing b appropriately® Alternatively) the Wittraeyer process can be directly adapted to

this problem by replacing € by M, and Uq 9 9 Up ? ®o®«, by

square matrices Uq 9 9 Ug?*®® « The importance of these remarks in the present context

is merely that they show that with a little preliminary preparation! any system An = b with a non-singular matrix can be solved

by the Monte Carl© methods ^ t© be described hereinafter®

*A method

J d¥e

=tT@^l©rHT

=Ti^e [7? PP® 132-133]) takes V as a triangular matrix obtained by replacing all elements of A above

the principal diagonal by zeros® The iterative method of Jacobi takes ¥ as a principal diagonal matrix whose diagonal is

that of A® The well-known method of Seidel is another Wittmeyer process® The so-called relaxation method is not a linear

process® So far as the author is awa re ? the extension of

Monte Carl© methods t© non-linear processes has not yet been

accomplished 9 and may be impossible® For an interesting and

scholarly classification of o@o-s tochastic methods of solving

linear equations ? see Forsythe [9]°
10

The stochastic methods for N = 1 » The basic idea in the

Monte Carlo attack on problems (Id) and (1 ®3) will now be described for the case of the zero-th iteration of the recursion

formula., To avoid dealing with a trivial problem, we present

the ideas here with L taken as the continuous operator (2 0 2)

rather than as a discrete operator..

The problem then is to estimate the numerical value of

u

1 (x) = J

h(x,y) Uq

(

y) dy + c(x).

* R

The function c(x) will play no significant role in the present

discussion, nor will K ? so we confine ourselves to the estimation of

I(x) h(x,y) u(y) dy»

Now let functions z(x,y) and p(x,y) be chosen so that*

(1 ) zp = h

,

(2) p £ 0,

(3) ) p(x,y) dy = 1

„

s R

Then for each x, p may be regarded as a probability density on

R» Let X be a vector random variable with the probability distribution defined by p. Consider the random variable

Z = z(x, X) u(X)» Clearly

E(Z) = ft z(x,y) u(y) p(x,y) dy

= h(x,y) u (y) dy = I(x)»

*The stochastic mechanism which is being arranged here has been

proposed by various writers.. See Kahn [11], and also [l4-] 0
11

A$ ao example, let R be one dimensional , and let

. ^v.:4w .w h(x,y) = e“Xy

? 0

y > °

y < 0

Then I(x) — J

e u(y) dy, the Laplace transform of u(y) c A

natural choice for p(x,y) would be the Pearson Type III density

function, e“x

^

x

0 Then z(x 9 y) = x, and Z = x u(X)«

One of the standard statistical procedures for estimating the mean of the distribution of a random variable Z is to

make J independent observations on the random variable and then

take their average Z as the estimator,. If the standard deviation*

^

of the distribution of the random variable is finite,

then the probability distribution of Z is asymptotically Gaussian, or "normal®11 , with (of course) the same average, and a standard deviation equal to

^ //>/2, (This is a special case

of the famous Central Limit Theorem; see Cramer [3, pp 0 21 5=217]

The approximation is usually very close for J 2 30„)

More than 99 per cent of" a normal distribution lies

within the intervals [mean + 3 X standard deviation],. From

this we can easily calculate the sample size J theoretically

necessary to achieve with this level of certainty a given statistical accuracy in using Z as an estimator of I(x)»

The standard deviation of Y is defined t© be ^E[ (Y-E(Y) )^]J

^ ^ (e(y2

) - [e(y)] 2

J

1/2

. The square of the standard deviation is

called the variance; we shall write it as Var(Y)„
12

Let us say that we wish to be almost sure that Z will

lie in the interval I + AI» That is, we want to arrange things

so that Pr ( | Z - I | ^ AI) > „99° This means that AI >

from which we get

(4ol) 9

(AI ) 2

9 Var(Z)

(AI) 2

The formula shows that if AI is to be small, say 0o005,

then unless Var(Z) , is also small, f will be well up in the

hundreds of thousands., The two redeeming traits of the stochastic method are that (1) the sampling error and necessary

sample size are independent of dimensionality and (2) they are

independent of how locally smooth the integrand of I(x) is

«

(It will be recalled that the error terras in the standard quadrature formulas involve high-order deriva tives »

)

But it is clearly worth while to consider methods of

reducing the sample size« The statisticians have a number of

devices for increasing the accuracy of sampling surveys, and

almost all of them are applicable here. In the present paper

we shall study only the procedures known as "'sampling with

probabilities in proportion to size," or "importance sampling"’ .*

They take advantage of the fact that we have an infinite number

of ways to choose z and p for any given problem, and an astute

choice may decrease the variance by a surprising amount.

*See Deraing [6 pp 0 92-93]°
13

The variance of Z (and we henceforth assume that it is

finite) is given by the formula

= E (Z2

) - LE(i:)]

2 = f 2 hu2dy - I2

0 y R

If the answer I were known in advance , and z were

chosen to be I (p would then have to be hij/l ? which would

mean that the integrand hu was non-negative) , then clearly

2 = Oo This leads us to propose the following arrangement as

a guide in choosing z and p„ First choose an integrable, nonnegative function p

8 such that the function 6=6 (x,y) defined by

(4»2) hu - p

8 = p

8 £

is as small as possible in absolute value,, The function p

8 is,

moreover, to be chosen so that the integral p” dy can be

obtained numerically without too much trouble,. Let p = p*/$

and let z = h/p = J(1 + i )/u„ Then it is easily seen that z and p satisfy the conditions (1 )

( 3) listed above,, The estimator Z is given by Z = J(x) (1 + 6(x 9 X))»

The variance of ©< + (3> Y, where c< and (3 are constants

2 and Y is a random variable, is simply (3' Var(Y)o Thus

( 4 . 3 ) 6\ 2 = Var(Z) = J

2 Var [ € (x,X)] .

If e(x) is the least upper bound of

J £ (x,y) f for Y on then

it might be safe to presume that Var(c) i e2/2„ (If the distribution of 6 were rectangular with range 2© — an unfavorable
case — then Var (6) = e /3« This appraisal gives us the formula

(4.4) / £

4. 5 J

2

e 2

(AI)

2

for the sample size theoretically required to achieve an accuracy of + AI with at least 99 °/o certainty .*

It might be noted that (*+<>3) and (4.4) are quite independent of whether or not j€| remains small for all y. If it does,

then there is an implication that hu cannot go very far in the

negative direction, since p

9 cannot be negative. These restrictions on the usefulness of the arrangement (4.2) can be circumvented to some extent, but we shall not go into the matter here.

We have been carrying the somewhat superfluous parameter

x throughout the above discussion mainly to emphasize the link

between quadrature and other problems. The usual quadrature

problem would of course be presented with h(x,y) - h(y),

u(y) = 1 . But it is perhaps worthwhile to observe that if there

is a parameter x in the problem, then by choosing p so .that it

is dependent only on x - y, a set of ^determinations of X can

be gotten once and for all, from which the statistics for Z

can be computed over and over again for as many different values

of x as may be desired.

One final remark which applies to all the remaining

sections of this paper as well as to the present section is this:

*In any case, certainly Var(6) <> e2
15

Th® statistical error was the only kind of error under consideration here, but of course there would be many other possible

sources of error in the actual numerical applications a For

example, there would be possible round-off errors, mistakes,

and systematic errors of one kind or another® In accordance

with th© resolution expressed at the end of the Introduction,

non© of these non-§ tatis tical errors will be discussed®

5° The Case N > 1 . Suppose that the method ©£ the preceding

section had been applied to the trivial problem of estimating

n n

u li

= ^ h(x i? x,)o(xJ

j=1 J J JET h. .u. r— i ] j J=1 J J Hu] i

?

for some fixed i, where the notation is that of Section 3*

®

The procedure would have been to select numbers • = z(Xj ,Xj)

and p 4 j = p(x^,Xj) such that

(1) Zij P ij

= h

i j

’ i j j » = 1 , o o o , n ,

(2) Pij i °> 3 ? j j " I,©©®, n,

(3) 2?ij = 1> x — 1 , 000 , n|

j

then to let X b® a vector random variable with the probability

distribution given by Fr(X = x.) = p. and to form the random

^ ^ J variable Z = z(x^,X)u(X) as before,, Its mean value would be

E(Z) = 21 z(x. ,x .)u(x .) p(x. 9 x .) = u 1

. o j x J J x J

1

A

*By bj

j

, where b is Tvector, we mean the i-th component of b®
16

We shall now show how to extend this stochastic scheme

to the iterations of H*

We continue to define z. • and p. . as above* Consider a

random walk on R defined as follows s the starting point Xq

is a random variable with a probability distribution given by

Pi = p(x.), where p^ > 0, P^ = 1 > but p i is otherwise ari

bitrary** Thereafter the successive positions or states visited by the random walk are random variables X, , whose

distributions are given by the formula Pr(X^

+

-= x . Xk=x.}=

i , J — 1 ) O » O )Ho

These directions have the effect of unambiguously specif yin g a probability distribution on the product space

RXRXoooXR for each N„ The typical point is assigned the

probability p.p. . p .

. p . . 1 1J 1 j

1 J 2

J 2

J

o o p The chain of random

n-1

variables Xq ,Xp ,

<> o o so defined is a simple Markov chain***

Now define a new chain of random variables os followss

Zq = Zq[ U ] = U (Xq)

Z-j = Z

1 [u,z] = z(Xq ? X-j ) u(X^ )

^2 ~ ^2 L u , z ] — z Xq ?X^ } z (X

1 5X2 ) u ( X2 5

©o —• o Zj^ — Z

N [ U , Z ] Z (Xq ,X.j ) « 0 a Z CXN__ 1 ?Xj^) U a

We shall call Zq, Z., ? Z

2

, = * c , an m -chain (m for mu 1 tipi

i

cative) 0

*This distribution does not play an intrinsic role in our discussion and is introduced only for logical completeness

„

**See Feller [8, Chap 0 1 5]° Actually our specifications uniquely

assigns a probability distribution to the infinite product

space RXRX°°° (see [12, pp 28-33]), but we shall not make use
17

Let the vectors

v o = v 01 o

VN±VN

Vq 5 , o o

o

9 be defined as follows?

(Xjju) = u(x

i

) ? (x^u^z) = E(Zn

[o ? 2] X

q

= X

A

), N=1 ,2, o

Consider now a typical path of the random walk represented by the sequence of random variables ,X^ 9 • ®

.

? X^o Let this

path be x.,x. 9 x. ?000? x. ? x.® The conditional probability

1 J 2

' J N-1 J that the random walk actually takes this path is p. . 000

1J i p j,j 2

p. -o Therefore by the definition of the concept of mean

J N-1 J value in the theory of probability s n n n n

¥»,<> = Ni X

jl” 1 z: oo

j 2“"^

J'ssm" 1 j=1

Zi

Ji

2

j

1 j 2

° ° °

2 • .u .p . . p . . © (, .p . • J N-1 J j 1 3

-J J1 J 2

J N~1 J

n n n N -i zl o o o 2 h 4

. h . . ,> ® ® h . .u. =HuJ..

j^1 j 2”"^ j=1 1J 1 J

1 J 2

J N-1 J J 1

Thus

( 5 . la) ¥

N+1 Hv

n

? N = 0 ? 1 y 2 j 9 » o 9

( 5 . 1b) VN

= HN

u s N = o 0 0 o o CM

The question of estimating u^ in the recursion relation

( 2 . 3 ) can now be resolved ratS^er easily® It is clear that

u Ni “ ( c +

^

H

N=1

c + H

N

u0

)] i = V0

(xi? €) + V^X^CjZ) + ...

+ vn-1 + vN

(x

i ;u0? 2g)
18

- E (Zq[ c J + Z^ [c ) z]+ o o L c J+Zj^[ Uq ,zj

| Xq — x.^} ,

~ e ^£n-i^ c?z ^ + 2

r/

u

o ?z ^

I x

o

” x

c^ *

where JT N

[u ? z] = Z

Q

[u] + Z^ [u,z]+o 0 o+Z

n

[u,z]o This shows that

if we define vectors w

Q

,w^ ,w2

, « • « , by

w

0i

= W’u) = u(xP’

w

Ni

= w

N

( Xi ;c ? u ? z) ~ E( 12Tj^ii= -| [ C , Z J + Z

N

[u ? z] X

Q “ x

i) 9 N” 1

, 2 , o o o

,

then with u = Uq, these vectors satisfy the relations

(5.2a) WNM = HwN

+ c ’ Z

1! O CM

(5.2b) w

N

= (I-H N

)(I-H)

-1

c+HN

u0

, N = 1 ,2,0

Furthermore, if we assume that the eigenvalues of H are

in modulus all less than unity,? then the vector w = Lim w.T OO M v N N —

>

00

exists and is the solution of the system u = Hu + c,

Thus the conditional mean values of the chain of random

variables 2 ^

, + Z,%T , N = 1 ,2,o«« provide a solution to our basic

problem 0 The actual computation would consist in making a large

number of realizations of the vector random variable (Xq ,X^ , »

»

o ,XN

),

calculating ,/t N

, + Z^ for each realization, and averaging the

results, or otherwise combining them into a statistical estimator of the mean value Q Several remarks are now in order 0

(1) This method of solving linear equations was proposed by Wasowkl 6] 0 He considered only the infinite product

space RXRXooo, and was thereby forced to use the restriction

that the matrix of the absolute values of the elements of H
19

had its eigenvalues in the unit circle® The pedigree of Wasow's

suggestion goes back to a well-known paper of Gourant, Friedrichs , and Dewy [l], in which the Green 9 ® function of an elliptic difference equation is identified with the mean number

of visits to the points of the lattice made by a particle performing a certain random walk on the lattice® (See the remark

(3) below®)

(2) Wasow and other writers describe the process of

realising Z

N

as a weighted random walk in which the particle

starts at x^ with a "mass" unity, which ii then multiplied by

the factor z(x. ,x. ), then multiplied by z(x . ,x. ), and so

1

^ J 2

forth® He deals in [16] with matrix inversion and chooses

u0

= 0, so his c = (6 k j, j = 1 9 2,...,n), where 6^ = 1, j = k,

6rj = 0, j

j* k® With this specialization,^^ then becomes in

Wasow 8 s words the total amount of mass carried through the

point X|. during the N steps of the random walk®

(3) If in the f unction IE ^[c 9 z] we take . = 1,

c = (6^j), then ^

is the total number of visits t© x^ made

by the N-step random walk starting at x^, counting in one

visit for the starting point if i = k; and w|s^q 9 ^jl)

is the mean number of visits to X|„® (In making this observation we are of course no longer tying ourselves down to the

requirement that z. .p. . = h. . where h . . is given in advance,

1 J X j X J X J but instead we are assuming that the p^. are given a priori®)

<*P
20

(*0 The chain Zq, Z, , represents a type of branching processo It is a Markov chain which may be of any order

from one to infinity, depending on the choice of the z®s 0 Such

chains have been studied by Montroll [153* together with the

companion type obtained by replacing multiplications by additions o They have numerous applications to theoretical physics

and physical chemistry,

(5) In the solution of u = Hu + c , or in the problem of

inverting of I H , it should be noted that once a large set of

realizations of (Xq ,X i X^T ) has been obtained with all

the random walks starting from some fixed x. , then by proper

JL bookkeeping procedure they can be used a number of times* Not

only can they be used to get the statistics for a number of

different vectors c-~= say, all the columns of the identity matrix

I, which would be the procedure for inverting I-H— but also,

by considering a visit to x., j ^ i, as starting a new random

walk, they can be used to find components of u^, N* < N, or of

u, other than the i-th componento Nevertheless a peculiarity

of the stochastic method here presented is that it seems to appear to the best advantage in comparison with the standard deterministic methods when the problem is to find only one component of Uj^j or of u, or one row of „ This fact has already been commented upon once before, in the Introduction*

(6) A final remark of minor importance is that the

estimator used in the quadrature problem of Section 4- was
21

analogous to the Z . in the present section, not to Zq Q Thus

to achieve strict parallelism between the two sections, we

should have assigned an appropriate distribution to RXR in

Section 4 rather than just to R, and then considered a condi=

tional distribution in R/R under a hypothesis relating to R 0

6 0 The m-chain method Q The development of the preceding section can easily be modified so that an m-chain does the work

of the random sequence^

q

+ Z ; 2-j + Z ? , 0 .. o To accomplish this,

we adjoin n additional points , x n+2 ? ° °

?x2n t0 and set

up the 2n x 2n partitioned matrix

h*(x

i ,x.)

i"* 1 $ © o o ^ 2n

j'3* 1 ^ © o o

c i = ®< x

i> = c ^ xi+n^ » 1 = 1 *

(6o1

)

HN

H

1

1 I

1 0

f

‘ I

1 _

o o o ?o0n ?

Do It

H

2

o o o h

n

) (I -

[N=1

-1

0

Now let zj . and pj. be chosen so that .p| . = h| . ,pt .£0

,

sj J ar s) O 2n

5- p*

j=1 ij

tice x^

,

ties p* . ij

1

o o o

o

o Set up the random walk X*q, X*., , 00 o, on the lat-

,x

n

j x

n+1 ? ° °

?x -n ? using the transition probabiliThen construct the new m-chairio
22

Z*0

=

Z* . =

z* = ^ N

Z*0 [(u:c)] = u(X*0

), x

*o

= 9

0 0 O ,X R

,

c(X*Q

),

II 0 X* x

n+1 5OO05X2n°-

= z*(X*0

,X*

1 )x|

u(X* 1 ) , X*

1 =x

c (X*

1 ) , X*^ =x 9000 9 X 2n

Z* N [ (use) ,z*

z* (x* 0

?x*

1 ) ooe z*(X

N==1

,X* N

) X

u i X* j^) ? X* ^

joo

»

c (X* N

) 5 X*^ ~ 9000 ?^ 2 n'

Let the vectors v*Q,v*p 00 o, be defined as follows

v* oi = v* 0 (Xi5 (use) ) = U^ , 1 ““ 1 O O O J

G« , 1 n^- 1 5 n+ 2 1 o o o «2 do 1—

n

7 IT?

v*Ni = v*n C x

i> (usc),z*) = E(Z*[ (use) ,z*]

j Xg = x

±

), N= 1 , 2 ,oo

Then according to ( 5 ° 1 ) these vectors satisfy the relations

(60 2a ) v

j!sB-1 ~ VN J ^ = }2)ooo|

(60 2b) v* = H* N

(usc), N= 0 , 1 , 2 , o,o o

The components of these vectors for i = n+ 1 ,ooo, 2n deserve some attention,, From the definition of vg, vg

n+ ^~ c^,

i = 1 ,ooo ? n„ As for v| , it will be seen by looking at H* that

*lj Pjj = i > n, j > n Q This means that vf ±

= v

| ?n+i

vN,n+i " c i’ 1 ~ 1 * o o o ^ n o

Using (6„1) to (6,2) and these facts, we obtain, letting u = Uq
23

(6.3a) vg

i = fu0i

c .

(6.3b) v*^. = jHvft]. + c.,

» i=o+1 ,.,o,2n

(6.3c) vj^ = j

hN“0 + (I-H N)(I-H)

_1

c] i i=1 ,n

c i-n’ X” n+ 1 j

«

o » j ki ii

.

N“ 1 , d.

, ® O 9

(In forming the product of the non-eonf ©rmafole factors II and v* ? only the first n components of v* are to be used,)

These relations show that the first n components of v*

satisfy the recursion relation ( 2 ® 3 J with the initial vector Uq 0 M

If j-j then v* = Lira vf, exists and its first n components

CO SB T w I >1 N ~->©o

satisfy the linear equations (2,5)°

Although the vector w^ of the preceding section and the

vector v^ of this section must be identical, they are based on

quite different stochastic processes,. For one thing, in the

n previous section it was understood that p„ . always equaled

j=i 1J

unity, which implied that the random walk remained forever

on x,,, 9 o,x

n? whereas now the corresponding sum of probabilities must be less than unity, and the random walk will not nee-

©ssarily remain forever ©n the first n points of the lattice.

In fact, if we choose p$ . = 6^j for i > n ? then the

new points ’ x

ni? ,eS9,x2n become trap states or absorbing

states for the random walk X| , , ®

®

• It is now a random walk

with absorbing barriers 0
To increase the comparability of the processes of this

section and of the preceding one for solving (2„3) and (2,5)

we can set up the solution of the preceding section with H* replacing H, z* replacing z, p* replacing p, (uQ

!o) = (u

Q1

,uQ2 , . „

,

u0n’ °>**o>0) replacing uQ

, and(csO) = (c

1 ,e 2

> • . . ,c n?

0 ? . . , ,0)

replacing c, The vectors (u^sO) and (csO) have the effect of

N

annulling the right-hand blocks of H* • Then the process

Z^CCc'COjZ*] + z*[(u

Q sO),z*], where^> * = ZJ + <,,<>+ Zjjj, will ^

N

have the same conditional mean values for XJ = i <> n, as

*r Zfc[ (

u

n % c )z*

]

* The random variable) + _f, is identically equal

to zero if the random walk ever visits the new part of the lattice (that is, the point-set x

n+ i

’ x

n+2 ’ ° ° 0 ,x 2n^ -*- n t*ie course

of the first N steps.

The m-chain method of inverting matrices was first proposed by Forsythe and Leibler [10], following a suggestion of

von Neumann, They considered only the infinite product space

RXR)( and thus had to make a restriction on H similar to

that which Wasow used.

7« Variances , The variance of the random variable Z^[u,z] of

section 5 is easily derived. Let K = [z. . h. .] = [h. ./p.

-tj

Then by the argument which led to our fundamental equations

(5»1 )

,

we find that

(7.1) 2

VNi 2v

N^ xi’ ~ ^

f ^N^ u,z ^ X

0=x± = K

N

u

2

].,

where u means the vector (u^ , U

2 ,o,.,u

n ), (We shall frequently
25

use this exponential notation to represent the operation of

squaring the components of a vector,)

The conditional variance of Z^? given that Xq = x^? is

2 of course the vector 0 vaT— vM . 2 N N

Turning to the random variable) [c ? z] + ZN

[u?z]? we

N-1

first find by using the method that led to (5°1) that for § £ t

,

E(Z

s [c,z] Zt u,z

|

X

Q

= x

A

) = K

sCH t_s u] i

,

where € is a principal diagonal matrix whose diagonal elements

are ? C 2 }o.»)C^ a We apply this relation and (7»1) to the

expansion of the square in the last member of

2

w

Mi = 2

w

N

(x i ;c ? u,z) = e|( [c-,z] + Z

N

[u?z]) 2

j XQ

= x

i J

= E {

(ZQ

[ c ] + Z

1 [c?z]+» « „+ Z N_ 1 [c ,z] + Z^Lu?^]) 2

| X

Q

=x

By going through some matrix algebra? the following rather formidable (but equivalent) formulas can be reached?

(7<>2a)

2

w

n

= (I - SO”1 ®

2

+ K

N

[u 2

- (1“K)”'

! J

+ 2(I-K)"1 [C(I-HPM )(I-H)"1 He - S^c]

+ 23 N-1 Ho *

(7* 2b)

2

w

N

= K

N==1 e 2

+ K

N

u 2 + (I-K)^1 (I-K^1 )[2C(I=fi)“ 1 c-c

2

]

- (CH N“ 1 +S

n=2 )(I-H)“1 He + 2S n=1

Hu ?

where = KN

CH + K

N==1 CH 2+ OO0 + K2CHN~ 1 + KCHN

.
26

The variance is of course to be obtained by subtracting

2 the vector from the vector 2

w

N°

The limiting forms of these formulas will now be derived

,

under the hypothesis that the eigenvalues of both K and H lie

inside the unit circle in the complex plane® For this we need

the following simple result®

Lemma » If A is any n x n matrix with complex elements , and with

the property that all of its eigenvalues lie in the unit circle , then there exist constants m = m( A) > 0 and r = r(A) > 1 which

are independent of N, and are such that a

i ? j ” i ;2j o o o )R j N = 1 , 2

,

(N)

ij < m/rN

, (N) where a . . is the element in the

N

i-th row and the j«=th column off A

For the proof, we first observe that the eigenvalues of

rA, where r is now any scalar, r f 0, are those of A multiplied

by r, because if y is an eigenvector for the eigenvalue A ? then the equations Ay = y and rAy = r A y are equivalent® This

means that if A satisfies the hypothesis of the Lemma, then a real number r exists, r > 1 ,such that the eigenvalues of rA

N

are still less than unity in modulus® Therefore Lira (rA) = 0®

N oo Therefore there is a positive number ra, independent of N, such

that for all N, 1 rNa.^|< tj m< The result follows at once from

this

We now apply the Lemma to S

N‘

Let c = max Then*

^We are using the notation EL . for the i,j-th element of the

matrix B«
27

J(K

sC)

i jj

< m(K)c/[ r(K)

]

S

, and

| (KS CH N“ S+1

) i j

| < nra(K)in(ii)c/[r(K)] S [r(H )]

N" s+1 ’ 0 *

S - N ’

With m ~ max[m(K), m(H)] ? F = min[r(K), r(H)], this becomes

j (KSCH^ I ^ nra ^c/r « From this we get the fact that

if the eigenvalue of K and 11 are all less than unity in modulus , then there exist numbers m = m(K,H,C) > 0 and r = r(K,H) > 1 independent of N such that

(7*3) j (S N) ^

. < n N m/ r M, i , j = 1 , . . . ,ra; N =

(This estimate could be refined so that it is

of n, using for example the methods described in [2,

note)], but it is sufficiently precise as it stands

purposes •

)

1 ^ 2 ^

e 9 j «

independent

p.16 (footfor present

The inequality (7°3) implies, of course, that Lira S w

=0®

N-> oo Therefore the limit of the vector pW^, which we shall

denote by 2W

00 ? hy the formulas

(7Aa)

(7Ab)

2w = (I-K) ’c

2 + 2(I-K) (I—H) 'He,

2

w = (I-K)

-1

L 2C (I-H)

-1

c - c 2

].

It is of interest to specialize the m-ehain variance

given through (7.1) to the augmented 2n x 2n matrix H* discussed

in Section 6. Let z* . and pi.. be chosen as previously, but

3 <3 with the proviso now that pi . = 0 whenever hi .= 0 , i=1 ,

2

, e » «,

,

n , J J j=n+1 ,n+2, . . „ ,2n, and also that pi. = 6^., i > n, j > n, Then

the matrix K* = [z* . hi.] 5:138 t5ie appearance
28

K* =

i=1

(6»2) that if we define

2

u

i

2

zl? . h . . 1 ij ij i/pf 0

1 -L“ 1 5 © © & 0 | l/p|

j=1 9. . .n 1

i

0

©

© •1/p

0 I

n - X efr It follows as

- v5i

-

C .i-n j

i = 1 , • .

9 n,

i — o+ 1 j

* » » j 2n j

vfo = e{ (Z^L (use) ,z*])j

2

I Xg = xJ , N = 1 ,2,...

then letting u = Uq ?

(7o5a) v* = i ^2v

Np ^i

2V

N+1 , i ) ^

c i-n ?

Qc i — 1 i ~ n+ 1 jooo^£i.n

N=0 , 1 , * * . ,

(7»5b) ovti 2w N,i = i kN°0

2

^1 + (I“0 (I-K) , Qc^] i

U i=1

2-

2

C o « 1“IS ? xTM n+

1

1 ^ 2 ? ® ® ® « Here Q denotes the principal diagonal matrix in the upper righthand corner of K* , and K means the n x n matrix* [z|. h. .]•

or The limiting form of (7«5) is

(7.6a) 2vSoi = K( 2v

oo + Qc

i’

1 = 1 >•••>"

(7 • fob) 2vw i = ^ ' Qc ]j_s

*To use K in this way is slightly inconsistent with the previous

definition of K; perhaps K** should be used; but the author felt

that in deference to good taste and the British readers

,

the

page should not be made to look any more star-spangled than it

already does®
29

It was noted in section 6 that the process

* y [ (e;0) ,z*] + Zi, t(an

.

: 0),z*] ^

N-1 TM u

has the same conditional mean value as Z^[(uie) ? z*] for = x^,

i ^n Q With K redefined as above to mean z| . h . . , (7°2a)

(7 o2b) 9 (7Aa) , and (7® L*-b) give the components for i < n @f the

'\— aje second moment of the random variable) + Z* exactly as they

N-1

stand now, The components with i > n are all zero.

Several remarks will now be made,

(1)

The formula (7®6b) was derived by Forsythe and Leibler [16] as the second moment of the conditional distribution of

the random variable defined on the infinite product-space

RXRX . © o ©

(2)

The formulas for the second moments suggest that to

keep down the variances® the numbers z» . and p. . or zf . and pf . ij ij ij ij

should be chosen so that on the average (speaking intuitively)

the elements of the matrix K should be as small in absolute value

as possible. One way not to achieve this end is by letting

p o . op pf . be positive when h . . = 0 or h* . = 0 , because the

unnecessary positive values of the p

8 s in a given row could be

portioned out to the other elements of K or K* in that row so

as to make the elements h. ./p. . or h$ ./p* . a little smaller,

(We have anticipated this remark to some extent by choosing

p±j = 0 when h^j = 0 for j > nj

(3)

A corollary of our results of some slight interest for matrix theory is this. If Lh. .] is a real n x m

J
30

matrix with all of its eigenvalues inside the unit circle, and

if there exists a set of positive numbers p^ . , i ? j = 1

? ,,, 9 n

such that I>_ p. . = 1

5 i = 1 ?

, ,

,

9 n and the eigenvalues of

3

J 2 [hf are all inside the omit circle 3 then, the diagonal

elements of (I-=H)”"® are all greater than 1/2, This follows

from (7°4b) with c specialised to c = ( 6 ^, j= 1

? oo, ? n) ?

k c=" *!<aooo^Oo

(4-) The natural statistical estimation of the mean values

VN

? W

N

9 VN 9 are course the arithmetic averages of many determinations of the respective random variables of which these

are the theoretical mean values. Clearly the variances of

these estimators always exist for finite matrices and finite

values of N ? no matter where the eigenvalues of H or H* lie.

Therefore in particular^ the formula (4,1 ) for the number of

random walks which must be performed to attain a given theoretical accuracy are here valid ? with AI replaced by Aw^? etCo

,

as the case may be.

80 Comparison of the two methods of inverting matrices . We

shall now use the formulas of the preceding section to effect

a comparison of the statistical error of the ra-chain method of

inverting matrices (Forsythe-Leibler) to that of the method

based on/ + Z M (Wasow) given in Section 5, u N-1 w

It will be impossible to compare these methods unless

the z

8 s and p

8 s are chosen comparably. Therefore for both
31

methods we shall suppose that the arrangement described in Section 6 has been set upj that is , the one which uses H* , 7\y

i ? 'N' nd Z*„ The comparison will be made by compari ng

only the limits of the variances as oo , seems to be rather int rae tible <=

finite case

for

With K = [

o h. i ? j = 1 and Q defined as

1 J 1 J

j the formulas to compare are

2wo©= (I-K)” 1 L2CCI-H)" 1 c - e 2

]

and

2 oo = ( 1= =1

The specialization to th

accomplished

random walks

by letting c =

start from x =

( 0

k j

9 x . , we

i

case of matrix inversion is

= 1 ?e ® OJ n )5 then if the

shall be estimating (X-H)^ k

®

( 8 . 1 ) 2W o©i = (X-K) ik [2(MI) k

1k - 1],

2voo i = (l~K)7k/p| .

The ratio of these second moments is

(8.2) 2 coi _ _* /t -si “7^— - Pf 12 d“B) kk - 1J

2 ooi

Since (X-4I) kk is fixed by the terms of the problem

,

we have at our disposal only the numbers p| in this formula®

It is clear that for a given H these can always be chosen so

that the ra-ehain method (Forsy the-Leibler ) is poorer than the

other method (Wasow). And for an H such that the diagonal ele

i i meets of (1=41) are small, the Wasow method might have a
32

smaller variance than the other method ? no matter how the random

walk was set up®

These negative ©bs

in an attempt to make the

(3 0 2), then the numbers p

the factors z* . will have ij

will make the elements of

of both 0w_ and 0 v_ will 2 o© 2 o©

able® On the other hand ?

ervations come far from telling the

are adjusted so as to be quite large ? m=ehaim method show up favorably in

$ . ? i < n 9 j < n 9 will be small and

1 j

• S

t© be proportionately larger® This

and (I-K)”' large ? and the values

be large 9 which of course is undesirif the pf Is are chosen so that the xj

values of p| are small ? then the mean duration of the random walk

becomes very long (see Section 9 below) ? and this is perhaps

more undesirable for the Wasow method then for the m-chaio

method because the former requires a little more computing per

The question of an optimum method from the point of view

of minimizing computing has not as yet been settled ® Probably

there is no such things because as in the case N = 1 treated

in Section it will turn out that the more one knows about

the solution., the better one can do® This phenomenon will

again be in evidence in Section 10®

9° The duration of the random walk ® We shall now apply the

formulas of the preceding section to obtain theoretical mean"

values and dispersions in two classical problems connected with

random walks , which h®v© a bearing ©n the usefulness of the
J J

Monte Carlo methods discussed in this paper®

We shall be considering the random walk X^ , X| , * , . , on

the lattice x^ ?

® oe ,x

n

, x

n+1 ,.„, 5 Xp

n

, Section 6® However., in

considering the associated matrices H* = [zfj an<^ 2 K* = [z#

u PjjJ> (the latter was introduced in Section 7)?

the point of view will be that the transition probabilities

pjj and the factors zjj are chosen first, and the values of

h| . are determined thereby, instead of the other way around®

We suppose that p$ . = b . . , i > n« j > n* ij ij

In each of the problems to be considered, z. . = 1.

t J This means that K* = H* = pj . , K = H = Pj

j

? i, j = 1 ,

. •

.

,n.

The first of these matrices is a stochastic matrix® We denote

the second one by P®

The first problem is that of determining the mean

value and variance of the number of visits to x^ if the random

walk starts at Xj ® In accordance with remark (3) in Section 5?

the total number of visits to in an N-step random walk starting at x^ is

Z N

= 6

k (x8

) + 6

k (x

T

) + 5 k (x2

} + --+ 6

k (x^>

where 6^(Xj) = 6^j ? j = 1,2,...,2n, The answers to the problem

are thus given by (5.2) and (7.2a) or (7®2b)®

Let the mean number of visits to x^ if the random

(N) walk started at x^ be denoted by . Then

’kj

[(X_pN+ 1 )(I_p)-1 ] ik , i < n, k < n,

G

ik

} = ^ki’

^ X

6

= *i>

h. k > n
3^

Assuming that lim F = 0 ? we have

N —$> oo

i < n ? k < n ?

(9.1b) G

( ~} = 6 ki , k > n«

Then by substituting into (7®^b), we obtain*

lim Var(J [ (o. .) , 1 ] X* = x.) = 2G

N ~H» ©0 N

J '

The first two terms of the formula are of course just

a special case of (8.1).

The other problem concerns the duration of the random

walk Xg 9 X| ? o • • We define the duration to mean the total number

of visits made to each of the transient states ^2 $ • • ® before absorption takes place® In counting the visits , we count

in the starting point, of the walk as one visit®

If the walk is limited t© N + 1 visits (that is 9 to N

steps) ? then the duration is clearly

(IsO)

*By Var(Yj b) ' we mean the variance of the conditional distribution of the random variable Y ? given that the event b has

occurred

.
The vector d N satisfies the recursion relation

i *'•- n 9 N — 0 , 1

, 2 ,oo®

with d0 i “ 1 ° Its limiting value is

n (9-2) <*00 1

= ru-r>:l, i<n.

The variance of the duration can be obtained by substituting appropriately into (7® 2a) or (7®2b)„, We shall again write

out only the limiting case. This is obtained easily from

( 7 °h-b)| the i-th component is

where d = max d

0 oo i

1

The duration of non-truncated random walks in an infinite product space has recently been investigated in very generstudied at some length by statistical theorists in connection

with the problem on the mean length of a sequential test® (See

[*+] for further results and references,,)

o

al cases by Wasow [17] Previously, special cases had been

The conditional probability that Xj^ falls on one of
36

n the states x.g ? ® 0 « jX^ ? given that = x^? is clearly CP

1*)

j=1

1J

This is the probability that the walk lasts for N steps without

falling into a trap state,, To assert that lira = 0 is

N—» 00

therefore equivalent t© saying that the probability that the

walk lasts more than N steps approaches zero with N„ In fact?

if we consider the infinite produc t-spaee EXKX o» e for a moment 9 lira pN = 0 means that the walk is ^almost certainly”

N -> oo of finite duration? and conversely^ It is known (see for

example Curtiss [*+? Sfectioe 11])? that if from each one of

the states ?x 0 ? . ® » ?x n it is possible to reach a trap state

over a path in RXRX°°° with moo-zero probability? then the

walk is almost certainly of finite duration,.

This provides a sufficient condition for the validity

of (9.1) @nd (9o2)® If it is satisfied? then all the eigenTM

values of F be less than unity in absolute value,,

We mentioned in Section 8 that if the number p| = 1 n

K - Pf . were small? then the mean duration would tend t©

j=i J be large® This follows intuitively from the fact that the

smaller the absorption probabilities are? the longer the walk

will go on „ A somewhat more rigorous demonstration can be

given by considering the dependence of the eigenvalues of P

on the row-sums of P 0 If pj.. > 0? i = 1?<,® s ?n? j=1?® 0 ®?n?

then as the minimum row-sum approaches unity? the eigenvalue

of maximum absolute value* approaches unity? and thus the eon-

*Xt happens to be real and positive? since F is a matrix

of positive elements

.
37

N

vergenee t© zero of P becomes slower and slower,

10® Importance sampling ,, In this section we shall discuss

the problem of the control of the statistical error® The

treatment will be analogous to the one in the latter part of

Section 4-® It will pertain only to the estimators Z^ and Z*,

and not to the estimator/ + We shall assume throughout

N-1 W

that the problem is to estimate the solution of (2,3) or (2*5)

»

Therefore the factors and P^j? or and P|j? will always

be related in the usual way to the elements of the matrix H

or H* ®

It is worth while first to inquire into the conditions

under which an m-ehain can be a zero-variance estimator® It

will suffice here to examine the situation only for the function

Zj^[u ? z] of Section 5°

If did indeed have a conditional variance of zero,

given that the random walk starts at x^ , then the value assumed

on every path with positive probability would depend only on x^®

One way to insure this would be to choose = Xn^/Uj, where

is a constant; then

'N

N u(Xn ) u(X1 ) u(XM 1 ) N

- ^ . __L_ ... N 1 u(X„) = \

Nu(X0

).

u(X, ) u(X2

) u(XN

)

Conversely, it can be shown that this is roughly the

most general choice of the factors . which will insure the

t J desired result® We shall not try to give an accurate formulation of the theorem here®
38

Now if z. . = \u./u., then the requirement that z. .p. . = h. . t J t j

implies that the choice of p. . must be h.-u./u. • This in turn

requires not only that all of these quotients must be real and

non-negative, but also that ^ h

^

.u ./ x u ^ - 1 0 The latter relaj J J tion states that u must be an eigenvector of H corresponding to

the eigenvalue X °

If we stop insisting on stationary transition probabilities for our random walk, and also permit the factors z^j to vary

from step to step, then a zero-variance in-chain estimator of

U N

“ H u

o can vefy ©asily be constructed for any Uq, provided

tha t h. . > 0 and u. > 0, i,j = 1 We simply choose

J J ,0) «

h HUN-1 ,.i

ij »N>i

,(2) .

h ii “N-2,1

(i) _

u Na

ij u N“1 ,j

ij uN-1 , i

p : .

,(2 ) _

u n-i a

lj UN-2,j

Oo

o (N)_ u li

ij uOj

(K) It is easily checked that formed with these factors z ^j

and

with u = Uq

,

has the constant value u^

^ for any random walk

starting at x^, For each K, = H Vj,, where v^ E(ZjJXq=x^),

but the chain is not a zero-variance chain unless K = N»

The argument can be extended so as to allow zero elements

in H and Uq by giving a little attention to the undefined quotients

It is not difficult to proceed from here to a practical

arrangement whereby given only approximate values of Uq,u^,

u. , o N

! a chain can be set up in which the variance can be
39

exhibited as a function of the errorj of the approximations, as was

done in the case N = 1 in Section Because of limitations of

space we shall not pursue the matter further* Instead we shall

study a special class of m chains which use stationary transition

probabilities and for which the statistical error analysis is

usually easy to make* These m=-ehain§, however, have the slight

disadvantage that they are connected with a special type of

matrix H*

The main problem is, as usual, to estimate the vector u^

in the recursion relation

(10.1) “ ^ UN

+ e e

We shall also be interested in the problems of estimating u in

the equation

(10® 2) u = Hu + c ,

but as elsewhere in this paper this problem will be considered as

the limiting case of the first problem as N becomes infinite*

Whenever (10*2) is in view, we shall as usual assume that all

the eigenvalues of H lie inside the unit circle*

We now impose certain conditions on They are that

> 0, and that a vector Uq with positive components shall

exist such that the components of (I - H)uq are all positive*

The existence of such a vector Oq follows automatically from

known results on matrices if the eigenvalues of H all lie inside

the unit circle and H is non-singular • We write Uq = HuQ + c

+ £ I thenc. + 6.>0 ? j = 1 9 . G * ? n

.

J: <3 If the problem in view is to solve Au = h, then
theoretically speaking, matrices fi and M can always be selected

so that the equivalent system u = Hu + Mb, with H + MA = I,

satisfies the conditions imposed here. For example, choose H

to be a principal diagonal matrix with elements lying between 0

and 1 o Then its eigenvalues lie in the unit circle, and the

existence of a vector uQ with the required property is assured,

(Of coarse, any particular a priori choice of H such as this us-

—1

ually means that M must be determined from M = A (I-H) , which

-I requires a knowledge of A , In practice therefore it is M

that will probably be chosen first, not H„)

We note that 6 (or rather -£) is the K residual” in

the classical theory of solutions of the system (I - H)u = c,

It is easy to show that if the successive approximations u^ , u^jooo to the solution u of (10,2) are defined by the recursion relation (10,1), then

u N

- u = H

N

(I - H)

-1

£ o

The second member of this equation is the truncation error (as

opposed to the statistical error) in the statistical solution

about to be proposed.

At this point we reformulate the problem in terms of

he partitioned matrix and vectors introduced in Section 6, The

equation (10,1) can be rewritten as

(10o3) (uN+1 :c) =H*X(uN

sc)

and the equation (10,2) becomes
(10.4) (uic) = H* X (use)

41

The problem ©f solving (10*2) then becomes one of finding an

eigenvec torof H* for the eigenvalue unity, given preassigned

values c. ? C 2

? oo ® 9 c

ri for the last n of the 2n components of

the eigenvec tor <, This formulation of the ^steady state’8 problem

as an eigenvalue problem permits us to make an approach to the

zero-variance sampling situation discussed earlier in the section*

We define correspondences between the points x

n+ ., » x n+2 ? 0 ° ° ? x 2n ? an<^ coraP on€nts °f 6 and c respectively by

c (x« ) = c. , € (x Un ) = ^ i = and we introduce

the vector e* whose components are given by

I Up ( X o ) , i “ 1 ? k J 8 e 3 c* = c*(x.) = ] j c ( x » ) + £ (Xj ) , i = n+1 ,

• o • ,2n.

That is f c* = (uq§c+£ ). The n equations represented by

«q = Hup + c + £ which define the components of £ are equivalent t© the 2n equations represented by

(10.5) c* = H*c*.

We are now ready to set up the basic random walk,

and the corresponding m-chain whose mean value is the solution

of (10.1)*

The transition probabilities of the random walk will

be given by

c

3$*; p* . = h* . -*

F xj ij c*

It is to be noted th § t . > 0 ? and because of (10.5)? i j st

i = 1 , • • * ,2n

2n

& pli = 1 ’
Letting the random variables XJ, X|

,

e » « denote as usual the

successive states visited by a random walk having these transition probabilities , and substituting into Z*[ (uQ

1 c) ,z*] , we

obtain

C *( X

N-1 ^

’ «*(X*) X

c(X*) ,if X* = x

i ,i=n+1

,

n+ 2 ,

<> o o , 2n

o

Cancelling terms and using the definition of c*, we find

for X^ on ,ooo,x

n € (Xfi)

(10.6) Z* = uQ (XJ) [1 - &(Xjjp

c (X^) + € (Xft)^

’

that

where

1 — 1 ^

o 9 > ) n

x n+ 1 ,o, g ,2n.

It is known from the results of Section 6 that

E(Z*

| Xq = x^) = , i 1 n, where u^ is defined by ( 1

0

0 1 ) above

The practical procedure implied here for estimating u^

is simply thiss Start a random walk at x^, using the transition

probabilities p* If absorption has not taken place after N

•I steps (that is, if the random walk has not reached any of the

points x

n+ ^

, x

n+2’ 3 3 3 ?x 2n^ ? then record Uq^ 0 If absorption

does take place during the N steps, stop the walk then and
there, note the index i of the last point x. touched before

X

absorption, and record u

Qj_ c j/( cj_ + Do this for many

random walks, and then average up the recorded values*

We now consider the statistical error of the procedure*

F rom (10*6)*

i 2 6 (Xfc) 1 Var(ZN I X

5

- x

i

) - u0i Var(6(X

fc) T(XfiTT £ (X*3'I

X

6

= x

i>’

1 - 1 , O <3 *> 4 1 1 ©

(We note in passing that if £ = 0, then the variance vanishes,

as it should). Proceeding as we did in Section h-, and letting

e = max

i

*i

C . + . 1 X

we get

Var(Z* X

6

= x.) < i = 1

, n

,

which for practical purposes is about as satisfactory as appraisal

of the variance as it seems possible to obtain* It contains

none of the unknown quantities in the problem.

By using the formulas and methods of Section 7j it is

easy to get an explicit formula for the variance. The result

is

*By the symbol Var (Yj b) , we mean the variance of the conditional

distribution of Y, given that the event b has occurred.
1+lf

Var(zfi| x

8 = xi>

= u oi (X-H N

)( I-H ) _1 f ']

± -

{(I-H) -1

?].]

2

,

i = 1 , ,n,

where

£

8 =

i c . + £• . ? i e i

i = 1, O }

One somewhat interesting conclusion that can be drawn

from this formula is that if Uq is chosen so as to have constant

components, then the variance will have one standard limiting

vector as N —> oo , no matter what the magnitude of the components

of Uq may be „ It turns out that this limiting vector is

(I - H)

1

u2

,

where u is the solution of ( 1 0 o 2) and

2

C! =

j=1 1J

i = 1 , °

»

n,
45

REFERENCES

1 o Courant, R 0 , Friedrichs , K0 , and Lewy, Ho, Uber die Partiellen

Dif f erenzengleichungen der Mathema tischen Physik,

Math 0 Annalen, vol„ 100(1928), pp 0 32-74

o

2 0 Courant, R 0

, and Hilbert, Do, Methoden der Mathema tischen

Physik, Berlin, 1931°

3o Cramer, Ho, Mathematical Methods of Statistics, Princeton,

1 946 o 4 0 Curtiss, J 0 H 0 , Sampling Methods Applied to Differential and

Difference Equations 0 Proceedings of a Seminar on

Scientific Computation, held by the International

Business Machines Corporation, Endicott, New York,

November, 1949° PP» 87-109®

5° Cutkosky, R 0 E», A Monte Carlo method for solving a class of

integral equations. Journal of Research of the National Bureau of StandardSo

6 9 Demin, W 9 Edwards, Some Theory of Sampling, New York, 1950°

7 0 Frazer, R 0 A*, Duncan, W„ J e , and Collar, A„ R», Elementary

Matrices, Cambridge, 1 950 0

8o Feller, W 0 , An Introduction to Probability Theory and its

Applications, New York, 1950 o

9o Forsythe, G„ E„, Tentative Classification of methods and

bibliography on solving systems of linear equations 0 To appear in Simultaneous Equations and the Determination of Eigenvalues— Proceedings of a National

Bureau of Standards Symposium held in Los Angeles,

Augus t , 1 951

o

10o Forsythe, Go E« , and Leibler, R„ A» } Matrix Inversion by a Monte Carlo Method, Mathematical Tables and Other

Aids to Computation, vol 0 *+(1950), pp» 127-129®

1 1 o Kahn, Ho, Stochastic (Monte Carlo) Attenuation Analysis,

RAM) Corpo , Santa Monica, Calif ornia 0 Report No,

R-163, June 14, 1949°
REFERENCES (Continued)

12. Kolmogorov, A. N» , Foundations of the Theory of Probability,

New York, 1950.

1 3o MacDuffee, C. Co, The Theory of Matrices, New York, 19*+6.

14 0 Proceedings of a Symposium on the Monte Carlo Method,

National Bureau of Standards Applied Mathematics

Series No. 12, Washington, 1 951

o

1 5o Montroll, E. W Q , On the theory of Markov Chains. Annals of

MathoStatis tics , vol. 1 8(1 9*+7)> pp« 1 o—3^> <*

16 0 Wasow, W., On the inversion of matrices by random walks,

Mathematical Tables and Other Aids to Computation,

volo 6(1951)? PP« 78—81 .

1 7» Wasow, W 0 , On the mean duration of random walks, Journal of

Research of the National Bureau of Standards, vol.

**6(1951), PPo 462-472o
THE NATIONAL BUREAU OF STANDARDS

Functions and Activities

The functions of the National Bureau of Standards are set forth in the Act

of Congress, March 3, 1901, as amended hy Congress in Public Law 619, 1950.

These include the development and maintenance of the national standards of

measurement and the provision of means and methods for making measurements

consistent with these standards; the determination of physical constants and

properties of materials; the development of methods and instruments for testing

materials, devices, and structures; advisory services to Government Agencies on

scientific and technical problems; invention and development of devices to serve

special needs of the Government; and the development of standard practices,

codes, and specifications. The Avork includes basic and applied research, development, engineering, instrumentation, testing, evaluation, calibration services

and various consultation and information services. A major portion of the

Bureau’s work is performed for other Government Agencies, particularly the

Department of Defense and the Atomic Energy Commission. The scope of activities is suggested by the fisting of divisions and sections on the inside of the

front cover. Reports and Publications

The results of the Bureau’s work take the form of either actual equipment and

devices or published papers and reports. Reports are issued to the sponsoring

agency of a particular project or program. Published papers appear either in the

Bureau’s own series of publications or in the journals of professional and scientific

societies. The Bureau itself publishes three monthly periodicals, available from

the Government Printing Office: The Journal of Research, which presents complete papers reporting technical investigations; the Technical News Bulletin, which

presents summary and preliminary reports on work in progress; and Basic Radio

Propagation Predictions, which provides data for determining the best frequencies

to use for radio communications throughout the world. There are also five series

of nonperiodical publications: The Applied Mathematics Series, Circulars, Handbooks, Building Materials and Structures Reports, and Miscellaneous Publications.

Information on the Bureau’s publications can be found in NBS Circular 460,

Publications of the National Bureau of Standards ($1.00). Information on cali- bration services and fees can be found in NBS Circular 483, Testing by the National

Bureau of Standards (25 cents). Both are available from the Government Printing Office. Inquiries regarding the Bureau’s reports and publications should he

addressed to the Office of Scientific Publications, National Bureau of Standards,

Washington 25, D. C.