Acta Numerica 1998 (Volume 7)

381 Pages·1998·18.099 MB·English

by  Arieh Iserles (Editor)

A c t a 
Numerica 
O  Q ) 
Volume  7 
0  o 
D 
01  •  
ooo o  ool  •  oio O  oi I O 
oooo •  oool  •
Acta Numeric a 199 8
Managing editor
A. Iserles
DAMPT, University of Cambridge, Silver Street
Cambridge CBS 9EW, England
Editorial Board
C. de Boor, University of Wisconsin, Madison, USA
F. Brezzi, Instituto di Analisi Numerica del CNR, Italy
J. C. Butcher, University of Auckland, New Zealand
P.G  . Ciarlet, Universite Paris VI, France
G.H. Golub, Stanford University, USA
H.B  . Keller, California Institute of Technology, USA
H.O­ . Kreiss, University of California, Los Angeles, USA
K.W  . Morton, University of Oxford, England
M.J  .D . Powell, University of Cambridge, England
R. Temam, Universite Paris-Sud, France
eta
umerica
1998
Volume 7
CAMBRIDGE
UNIVERSITY PRES S
Published by the Press Syndicate of the University of Cambridge
The Pitt Building, Trumpington Street, Cambridge CB2 1RP
40 West 20th Street, New York, NY 100114­ 211 ,U  SA
10 Stamford Road, Oakleigh, Melbourne 3166, Australia
© Cambridge University Press 1998
First published 1998
Printed in Great Britain at the University Press, Cambridge
Library of Congress cataloguing in publication data available
A catalogue record for this book is available from the British Library
ISBN 05­ 216­ 4316­ 3 hardback
ISSN 09624­ 92 9
Contents
Monte Carlo an d quasi­Mont e Carl o method s  1
Russel E. Caflisch
Nonlinear approximatio n  51
Ronald A. DeVore
Relative perturbatio n result s for matri x  eigenvalue s
and singula r value s  151
Use C. F. Ipsen
Stability for time­dependen t  differentia l  equation s  203
Heinz-Otto Kreiss and Jens Lorenz
Direct searc h algorithm s for optimizatio n calculation s  287
M. J. D. Powell
Choice of norm s for dat a  fitting  an d
function approximatio n  337
G. A. Watson
Ada Numerica (1998),  pp .1  4­ 9  ©  Cambridge University Press, 1998
Monte Carl o an d quasiM­ ont e  Carl o
methods
Russel E . Caflisch*
Mathematics Department, UCLA,
Los Angeles, CA 90095-1555, USA
E-mail: [email protected]
Monte Carlo iso  neo  ft he most versatile and widely used numerical methods.
Its convergence rate, O(N~1^2),  isi ndependent ofd  imension, which shows
Monte Carlo t ob  ev  ery robust bu ta  lso slow. This article presents a ni ntro­
duction t o Monte Carlo methods fori ntegration problems, including conver­
gence theory, sampling methods and variance reduction techniques. Acceler­
ated convergence forM  onte Carlo quadrature isa  ttained using quasir­ ando m
(also called lowd­ iscrepancy ) sequences, which area  d  eterministic alternative
to random orp  seudor­ ando m sequences.  Th ep  oints i na   quasir­ ando ms  e­
quence are correlated to provide greater uniformity. The resulting quadrature
method, called quasiM­ ont e Carlo, hasa   convergence rate ofa  pproximately
O((log N^N'1).  For quasiM­ ont e Carlo, both theoretical error estimates and
practical limitations arep  resented.  Although th ee  mphasis i nt his articlei  s
on integration, Monte Carlo simulation ofr arefied gasd  ynamics isa  lso dis­
cussed. In the limit of small mean free path (that is, the fluid dynamic limit),
Monte Carlol osesi ts effectiveness because the collisional distance ism  uch less
than the fluid dynamic length scale. Computational examples ar ep  resented
throughout the text t o illustrate the theory. A number ofo  pen problemsa  re
described.
Research supported in part by the Army Research Office under grant number DAAH04­
951­ 0­ 155 .
R. E.  CAFLISCH
CONTENTS
1  Introductio n  2
2  Monte Carlo integratio n  4
3  Generation and sampling method s  8
4  Variance reduction  13
5  Quasi­rando m numbers  23
6  QuasiM­ ont e Carlo techniques  33
7  Monte Carlo method s forr  arefied gas dynamics  42
References  46
1.  Introductio n
Monte Carlo provides a d  irect method forp  erforming simulation and integ­
ration.  Because i t is simple an dd  irect , Monte Carlo is easy t o use.  I ti  s
also robust , since it sa  ccuracy depends o no  nly th ec  rudest measure oft  h e
complexity of th e problem. For example, Monte Carlo integration converges
at a  rat e O{N~1/2)  tha t is independent oft  h e dimension oft  h e integral.
For thi s reason, Monte Carlo ist  h e only viable method fora   wide range of
highd­ imensiona l problems, ranging from atomic physics t o finance.
The price fori  t sr  obustness ist  ha t Monte Carlo ca nb  ee  xtremely slow.
The order O(N~1^2)  convergence rat e is decelerating, since a n additional
factor of4  i ncrease in computationa l effort only provides an additional factor
of 2i  mprovement i na  ccuracy. Th e result oft hi s combination ofe  ase of use,
wider  ange ofa  pplicability and slowc  onvergence ist  ha t ane  normous amount
of computer time iss  pent o nM  onte Carlo computations .
This represents a  great opportunit y forr  esearchers i nc  omputational sci­
ence.  Even modest  improvements i n th e Monte Carlo method ca nh  ave
substantial impact o nt  h e efficiency  an d range of applicability for Monte
Carlo methods .  Indeed, much oft  h e effort  i n th e development of Monte
Carlo ha s been i n construction ofv  ariance reduction methods which speed
up th ec  omputation .  A  description ofs  ome oft  h e most common variance
reduction method s isg  iven i nS  ection4  .
Variance reduction method s accelerate th ec  onvergence rat e byr  educing
the constan t i n front oft h e O(N~1/2)  forM  onte Carlo methods using ran­
dom o rp  seudo­rando m sequences. A na  lternativ e approach t o acceleration
ist  oc  hange th e choice ofs  equence. QuasiM­ ont e Carlo methods use quasi­
random  (also known a s lowd­ iscrepancy )  sequences instead of random o r
pseudo­random .  Unlike pseudo­rando m sequences, quasir­ ando m sequences
do no ta  ttemp t t oi  mitat e th e behaviour ofr  andom sequences. Instead ,t  h e
elements ofa  q  uasi­rando m sequence are correlated t om  ake them more uni­
form tha n rando m sequences. For thi s reason, Monte Carloi  ntegration using
MONTE CARLO AND QUASI­MONT E CARLO  3
quasir­ ando m points converges more rapidly, at a rat e O(N~1 (log N)k),  for
somec  onstant k. Quasir­ ando m sequences are described in Sections 5a  nd 6.
In spite of their  importance in applications,  Monte Carlo methods re­
ceive relatively littl e attentio n from numerical analysts and applied math ­
ematicians. Instead, it is number theorist s and statistician s who design th e
pseudor­ andom ,  quasir­ ando m  and other types of sequence tha t  are used
in Monte Carlo, while th e innovations in Monte Carlo techniques are de­
veloped mainly by practitioners , including physicists, systems engineers and
statisticians.
The reasons for th e near neglect of Monte Carlo in numerical analysis and
applied mathematics are related t o it s robustness. First , Monte Carlo meth­
ods require less sophisticated mathematics tha n other numerical methods .
Finite difference  and finite element methods, for example, require  careful
mathematical analysis because of possible stability problems, bu t stability
is not an issue for Monte Carlo.  Instead , Monte Carlo nearly always gives
an answer tha t  is qualitatively correct, bu t  acceleration  (error reduction)
is always needed.  Second, Monte Carlo methods are often phrased in non­
mathematical terms .  In rarefied  gas dynamics, for example, Monte Carlo
allows for direct simulation of th e dynamics of th e gas of particles, as de­
scribed in Section 7. Finally, it iso  ften difficult t o obtain definitive resultso  n
Monte Carlo, because of th e random noise.  Thu s computationa l improve­
ments often  come more from experience tha n from a particula r  insightful
calculation.
Thisa  rticle isi ntended t op  rovide ani  ntroduction t o MonteC  arlo methods
for numerical analysts and applied mathematicians .  In spite of th e reasons
cited above, there are ample opportunitie s for thi s community t o make sig­
nificant contributions t o Monte Carlo.  Firs t of all, any improvements can
have a big impact, because of th e prevalence of Monte Carlo computations .
Second, the methodology of numerical analysis and applied mathematics ,
including well controlled computational experiments on canonical problems,
is needed for Monte Carlo.  Finally, ther e are some outstandin g problems
on which a numerical analysis or applied mathematic s viewpoint is clearly
needed; for example:
  design ofM  onteC  arlos  imulation for transpor t problems int  h e diffusion
limit (Section 7)
  formulation  of  a  quasiM­ ont e  Carlo  method  for  th e  Metropolis  al­
gorithm (Section 6)
  explanation of why quasiM­ ont e  Carlo behaves like standar d  Monte
Carlo when th e dimension is large and th e number of simulation is of
moderate size (Section 6).
Someo  lder, bu t still very good, general references on Monte Carlo are Kalos
and Whitlock (1986) and Hammersley and Handscomb (1965).
4  R. E. CAFLISCH
The focus oft  hi s article iso  n Monte Carlof  or integration problems. Integ­
ration problems ar es  imply stated , bu t they ca nb  ee  xtremely challenging.
In addition, integration problems contain most of th e difficulties  tha t ar e
found i n more general Monte Carlo computations, such a s simulationa  n d
optimization.
The next section formulates th eM  onte Carlo method fori ntegration and
describes it s convergence.  Section 3  describes random number generators
and sampling methods .  Variance reduction methods ar ed  iscussed i nS  ec­
tion 4a  n dq  uasiM­ ont e Carlo methods i n Section 5. Effective use ofq  uasi­
Monte Carlor  equiress  omem  odification ofs  tandar d Monte Carlot echniques,
as described i n Section 6. Monte Carlo methods for rarefied ga sd  ynamics
are described i n Section 7, with emphasis on th e loss of effectiveness for
Monte Carlo i nt h e fluid dynamic limit.
2.  Mont e Carl o integratio n
The integral ofa  L  ebesgue integrable function f(x) can b ee  xpressed as th e
average or expectation of th e function  /   evaluated a t a  random location.
Consider a ni ntegral ont h eo  ned­ imensiona l unit interval
nn= f1 f(x)dx = f.  (2.i )
Jo
Let a; bea  r  andom variable tha t isu  niformly distribute d ont  h e unit interval.
Then
I[f] = E[f(x)].  (2.2 )
For a ni ntegral ont h eu  nit cube Id = [0, l]di  n dd  imensions,
)dx,  (2.3)
in which a;i  sa   uniformly distribute d vector i nt  h eu  nit cube.
The Monte Carlo quadratur e formula isb  ased ont h ep  robabilistic inter­
pretation ofa  n integral. Consider as  equence {x } sampled from th e uniform
n
distribution. Then a ne  mpirical approximation t ot  h ee  xpectationi  s
M/] = ^  £/(*n) .  (2.4)
n=l
According t ot  h eS  trong Law ofL  arge Numbers (Feller 1971), thi s approx­
imation i sc  onvergent with probability one; tha ti  s ,
(2.5)
In addition, i ti  su  nbiased, which means tha t th e average of Iff[  /] ise  xactly
/[/] for any N; tha ti  s ,
E[I [f]] = I[f],  (2.6)
N