Stanford NA: 1963-1984 From Pine (Polya?) Hall to Serra House

1. Stanford NA: 1963-1984From Pine (Polya?) Hall to Serra House Presented by: Victor Pereyra Weidlinger Associates Inc. Mountain View, CA Presented at: Numerical Analysis History Series Stanford University February 2007

2. September 7, 1963 J. B. Rosen

3. Calculating the Pseudoinverse CS 13, Computer Science Dept., Stanford University (1964). ftp://reports.stanford.edu/pub/cstr/reports/cs/tr/64/13/CS-TR-64-13.pdf

4. Difference Corrections • "The Difference Correction Method for Nonlinear Two-Point Boundary Value Problems", CS18, Comp. Science Dept., Stanford University (1965). • Fox (1947) Henrici (1960) M. Lees (1964) Gene VP "The Difference Correction Method for Nonlinear Two Point Boundary Value Problems of Class M", Rev. U.M.A. 22:184-201 (1965)

5. Transition • From Stanford to Madison, Wisconsin (1965)

6. PhD Thesis (Univ. of Wisconsin, 1967) • "Highly Accurate Discrete Methods for Nonlinear Problems", Univ. Microfilms, Ann Arbor, Michigan, Tech. Rep. 728, Math. Res. Center, University of Wisconsin, Madison, 117 pp. (1967). • "Accelerating the Convergence of Discretization Algorithms", SIAM J. Numer. Anal. 4:508-533 (1967). • "Iterated Deferred Corrections for Nonlinear Operator Equations", Numer. Math. 10:316-323 (1967). • "Iterated Deferred Corrections for Initial Value Problems", with J. W. Daniel and L. L. Schumaker, Acta Cient. Venezolana 19:128-135 (1968). • "Iterated Deferred Corrections for Nonlinear Boundary Value Problems", Numer. Math. 11:111-125 (1968).

7. Deferred Corrections in a Nutshell • General nonlinear problem: F(x) = 0 • Low Order Discretization: Fh2(xh2) = 0 • Residual (defect, local truncation error), calculated with higher order approximation: Rh2 = Dh(F(x)) - Fh4(xh2)) • Correction (it uses the same low order operator): Fh2(xh4) = -Rh2 • Also observe that ||Rh2||is a local error estimator and || xh2 - xh4|| is a global error estimator, so that the door was open for an adaptive h-p method for BV and other problems (invented around 1972), see: • "Variable Order, Variable Step Finite Difference Method for Nonlinear Boundary Value Problems", Lecture Notes in Mathematics, 363:118-133, Springer (1973). • A whole series of public domain programs for 2PBVP (PASVAn, n=1,…,4) were written between 1970 and 1984. • A good view of the state of the art circa 1984 can be found in theSpecial Issue of Computing:“Defect Correction Methods. Theory and Applications”. Editors: K. Bohmer and H. Stetter. Springer (1984)

8. Modern Work • Bertil Gustafson and his students at Uppsalahave made good use of DC for Elliptic problems (Domain Decomposition) and for hyperbolic systems (I.,e., initial/boundary value problems) in recent years: “Deferred Correction Methods for Initial Boundary Value Problems”. W. Kress and B. Gustafsson, Journal of Scientific Computing17 : 241 - 251 (2002). “Implicit high-order difference methods and domain decomposition for hyperbolic problems”. Bertil Gustafsson, and Lina Hemmingsson-Franden, Applied Numerical Math. 33:493-500 (2000). • Some contributions from England to Computational Fluid Mechanics: “Adjoint and defect error bounding and correction for functional estimates”.M. B. Giles and N. L. Pierce.Journal of Computational Physics200:769-794 (2004). • Spectral deferred corrections: A. Dutt, L. Greengard, V. Rokhlin, Spectral deferred correction methods for ordinary differential equations, BIT 40:241-266 (2000) . `` High-Order Multi-Implicit Spectral Deferred Correction Methods for Problems of Reactive Flows.'' Anne Bourlioux, Anita T. Layton, and Michael L. Minion, Journal of Computational Physics, 189: 651-675 (2003).

9. Another Transition • From Madison, Wisconsin to Caracas, Venezuela (1967).

10. Visiting Stanford • Visiting Professor, Computer Science Department. Stanford University, California. - January - March 1973; - August 1974; - July/August 1976,(Caltech 1974-1978), 1979, 1980, 1981, 1982; - January 1981; - January 1983. • In 1984 I moved to the area for good and became a full time employee of WAI.

11. The Caracas-Stanford Axis (1967-1984) • Fast Vandermonde solvers (1D): "On the Construction of Discrete Approximations to Linear Differential Expressions", with C. Ballester, Math. Comp. 21: 297-302 (1967). "Solution of Vandermonde Systems of Equations", with Ake Bjorck, Math. Comp. 24:893-904 (1970). "Solving Confluent Vandermonde Systems of Hermite Type", with G. Galimberti, Numer. Math. 18:44-60 (1971). Bjorck (Sweden) Golub (Stanford) VP (Caracas) • Carefully analyzed by N. Higham, “Error analysis of the Bjorck-Pereyra algorithms for solving Vandermonde systems”. Numer. Math. 50:613-632 (1987). • Many extensions from Stanford IEEE byTom Kailah and V. Olshevsky. See:“Fast Algorithms for Structured Matrices: Theory and Applications”. V. Olshevsky Editor. Contemporary Mathematics 3, Siam Pub. (2003).

12. Extensions to Multiple Dimensions "Numerical Differentiation and the Solution of Multidimensional Vandermonde Systems", with G. Galimberti, Math. Comp. 24:357-362 (1970). "Efficient Computer Manipulation of Tensor Products with Applications in Multidimensional Approximation", with G. Scherer, Math. Comp. 27:595-605 (1973). • These were all tools to automatically produce high order numerical differentiation formulas on unevenly spaced meshes for their use in iterated deferred correction algorithms. • Extension to least squares problems by Serra House alumniEric Grosse: “Tensor spline approximation.” Linear Algebra and Its Applications, 34:29-41 (1980). • See also: “Efficient Computer Manipulation of Tensor Products.” C. de Boor, ACM Trans. Math. Software 5: 173-182 (1979).

13. Fast Tensor (Kronecker) Product Operations

14. Variable Projections (VP) • "Least Squares Estimation for a Class of Nonlinear Models", with I. Guttman and H. D. Scolnik, Technometrics 15: 209-218 (1973). • "The Differentiation of Pseudoinverses and Nonlinear Least Squares Problems Whose Variables Separate", with G. H. Golub, SIAM J. Numer. Anal. 10:413-432 (1973). • "Differentiation of Pseudoinverses, Separable Nonlinear Least Squares Problems, and Other Tales", with G. H. Golub, Proc. Math. Res. Center Seminar on Generalized Inverses and Their Applications, Univ. Wisconsin, Madison, (Ed. Z. Nashed), 303-324, Academic Press, New York (1976). • "A Method for Separable Nonlinear Least Squares with Separable Nonlinear Equality Constraints", with L. Kaufman, SIAM J. Numer. Anal. 15:12-20 (1978).

15. Separable nonlinear least squares • A separable model F(t;a, w)=ai I(t;w) • To fit in the LSQ sense the data {tj,yj}: Minw ||F(t;a, w) - y||22 • Or Minw || (t;w) a - y|| 22, Where (t;w)={I(tj,w)}.

16. Variable Projections

17. VARPRO • A public domain program carefully implementing these ideas was written (in Venezuela and Uruguay). • Serra House students John Bolstad, Linda Kaufman and Randy LeVeque contributed to its generalization and beautification. • An important extension: G.H. Golub and R.J. LeVeque, “Extensions and uses of the variable projection algorithm for solving nonlinear least squares problems.” In Proceedings of the 1979 Army Numerical Analysis and Computers Conference 1-12, (1979). • Amodern version is available from nanet (implemented at Bell Labs. by L. Kaufman and D. Gay). • The National Bureau of Standards was an early enthusiastic adopter and they produced a version with a friendly user interface and extensive documentation. • At the Universidad de Barcelona, Spain, there is a center with hundreds of users for NMRS.

18. 30 Years Down the Line … • Many applications of VP. See "Separable Nonlinear Least Squares: the Variable Projection Method and its Applications". With G. Golub. Inverse Problems 19:R1-R26 (2003). • Most significant: It is a basic tool for in vivo Nuclear Magnetic Resonance Spectroscopy (NMRS). • Most amusing: What can Lattice Quantum Chromo-Dynamics theorists learn from NMR spectroscopysts?, G. T. Fleming (2005). • The answer: VARPRO

19. 2PBVP Course (1973) • "High Order Finite Difference Solution of Differential Equations", STAN-CS-73-348, Computer Science Dept., Stanford University, California, 86 pp. (1973). • It has some unpublished nuggets, such as U2DCG: an Universal 2PBVP Deferred Correction Generator (FORTRAN code included). • At this time we had codes using uniform meshes and it became clear that they would not work well for problems with sharp gradients, boundary layers and such. That lead naturally to the next stage.

20. Equidistributed Meshes (1974) "Mesh Selection for Discrete Solution of Boundary Problems in Ordinary Differential Equations", with E. G. Sewell, Numer. Math. 23:261-268 (1975). • This concept is now central to adaptive meshes. • At Stanford it was picked up by Oliger and Berger (1984) for finite difference methods. • For Finite Element methods Babuska and Rheinboldt (h-method) in the early 80’s used the same principle. • The actual origin of these ideas comes from splines technology: “Splines with Optimal Knots Are Better.” J
H Burchard, Applicable Math, 1974. “Good approximation by splines with variable knots.” C. de Boor, Lect. Notes Math. 363:12-20 (1973).

21. The unwritten Concus, Golub, Pereyra paper • During my extended visit to Serra House in 1973 I did some additional work that is mentioned in: “Variable Order Variable Step Finite Difference Methods for Nonlinear Boundary Value Problems”. Lecture Notes in Math. 363:118-133 (1973). • Section 4.3 of this paper deals with Poisson, Helmholtz and mildly nonlinear elliptic equations in rectangular regions. There I combine Buneman fast Poisson solver with deferred corrections. For the Poisson equation -u=50sin(5(x+y)), discretized in the unit square with a 32x32 uniform grid I obtained typical results as shown next.

22. Deferred Corrections Results By the way, using Moore’s Law and without any improvements in the algorithm these results could be obtained today in a high end PC in 3 miliseconds!

23. Computational Geophysics (1976-present) • In the summer of 1976, while I was a Professor at Caltech I visited Stanford and met Willie H. K. Lee (Seismologist, USGS) through Gene. • We started working on seismic ray tracing and earth tomography, where I could combine my two-point BVP and nonlinear least squares skills. • That lasted for the next 30 years. After learning that Herb Keller and his brother Joe had a lot to say about Ray Theory I started some early contributions with Herb: "Computational Methods for Inverse Problems in Geophysics: Inversion of Travel Time Observations", with H. B. Keller and W. H. K. Lee, Physics of the Earth and Planetary Interiors 21:120-125 (1980). "Two-Point Ray Tracing in Heterogeneous Media and the Inversion of Travel Time Data", Comp. Methods Applied Science and Engineering (Eds. R. Glowinski and J. L. Lions), North-Holland Pub. Co., Amsterdam, 553-570 (1980).

24. And Here is Some of the Developed Software for Modeling and 3D Ray Tracing on the Gulf of Mexico (INTEGRA)

25. Conclusion • My research life has been closely intertwined with Stanford’s NA group and with its main engine and facilitator Gene H. Golub. I deeply thank him as one of the many that have benefited from his generosity and foresight. • If I can leave a last thought as moral for the students: There is a social obligation to publish your work and for those that have the ability and the inclination, there is a significant payoff to implementing your results in usable software. • For more details, please visit: • http://homepage.mac.com/vpereyra/pereyra-bio.html