Abstract: To harness the full power of computer technology, economists need to use a broad range of mathematical techniques. In this book, Kenneth Judd presents techniques from the numerical analysis and applied mathematics literatures and shows how to use them in economic analyses. The book is divided into five parts. Part I provides a general introduction. Part II presents basics from numerical analysis on R^n, including linear equations, iterative methods, optimization, nonlinear equations, approximation methods, numerical integration and differentiation, and Monte Carlo methods. Part III covers methods for dynamic problems, including finite difference methods, projection methods, and numerical dynamic programming. Part IV covers perturbation and asymptotic solution methods. Finally, Part V covers applications to dynamic equilibrium analysis, including solution methods for perfect foresight models and rational expectation models. A web site contains supplementary material including programs and answers to exercises.

Abstract: We present a blind deconvolution algorithm based on the total variational (TV) minimization method proposed by Acar and Vogel (1994). The motivation for regularizing with the TV norm is that it is extremely effective for recovering edges of images as well as some blurring functions, e.g., motion blur and out-of-focus blur. An alternating minimization (AM) implicit iterative scheme is devised to recover the image and simultaneously identify the point spread function (PSF). Numerical results indicate that the iterative scheme is quite robust, converges very fast (especially for discontinuous blur), and both the image and the PSF can be recovered under the presence of high noise level. Finally, we remark that PSFs without sharp edges, e.g., Gaussian blur, can also be identified through the TV approach.

Abstract: We have examined an optimization approach to iterative control design. The important ingredient is that the gradient of the design criterion is computed from measured closed loop data. The approach is thus not model-based. The scheme converges to a stationary point of the design criterion under the assumption of boundedness of the signals in the loop. From a practical viewpoint, the scheme offers several advantages. It is straightforward to apply. It is possible to control the rate of change of the controller in each iteration. The objective can be manipulated between iterations in order to tighten or loosen performance requirements. Certain frequency regions can be emphasized if desired. This direct optimal tuning algorithm is particularly well suited for the tuning of the basic control loops in the process industry, which are typically PID loops. These primary loops are often very badly tuned, making the application of more advanced (for example, multivariable) techniques rather useless. A first requirement in the successful application of advanced control techniques is that the primary loops be tuned properly. This new technique appears to be a very practical way of doing this, with an almost automatic procedure.

Abstract: The need for solving multivariate optimization problems is pervasive in engineering and the physical and social sciences. The simultaneous perturbation stochastic approximation (SPSA) algorithm has recently attracted considerable attention for challenging optimization problems where it is difficult or impossible to directly obtain a gradient of the objective function with respect to the parameters being optimized. SPSA is based on an easily implemented and highly efficient gradient approximation that relies on measurements of the objective function, not on measurements of the gradient of the objective function. The gradient approximation is based on only two function measurements (regardless of the dimension of the gradient vector). This contrasts with standard finite-difference approaches, which require a number of function measurements proportional to the dimension of the gradient vector. This paper presents a simple step-by-step guide to implementation of SPSA in generic optimization problems and offers some practical suggestions for choosing certain algorithm coefficients.

Abstract: This paper is devoted to difference of convex functions (d.c.) optimization: d. c. duality, local and global optimality conditions in d. c. programming, the d. c. algorithm (DCA), and its application to solving the trust-region problem. The DCA is an iterative method that is quite different from well-known related algorithms. Thanks to the particular structure of the trust-region problem, the DCA is very simple (requiring only matrix-vector products) and, in practice, converges to the global solution. The inexpensive implicitly restarted Lanczos method of Sorensen is used to check the optimality of solutions provided by the DCA. When a nonglobal solution is found, a simple numerical procedure is introduced both to find a feasible point having a smaller objective value and to restart the DCA at this point. It is shown that in the nonconvex case, the DCA converges to the global solution of the trust-region problem, using only matrix-vector products and requiring at most 2m+2 restarts, where m is the number of distinct negative eigenvalues of the coefficient matrix that defines the problem. Numerical simulations establish the robustness and efficiency of the DCA compared to standard related methods, especially for large-scale problems.

Abstract: A family of new iteration methods is presented for designing quantum optimal controls to manipulate the transition probability Theoretical analysis shows that these new methods exhibit quadratic and monotonic convergence Numerical calculations verify that for these new methods, within very few steps, the optimized objective functional comes close to its convergent limit

Abstract: A new iteration method is presented for achieving quantum optimal control over the expectation value of a positive definite operator. Theoretical analysis shows that this new algorithm exhibits quadratic and monotonic convergence. Numerical calculations verify that for this new algorithm, within a few steps, the optimized objective functional comes close to its converged limit.

Abstract: We introduce a well-developed Newton iterative (truncated Newton) algorithm for solving large-scale nonlinear systems. The framework is an inexact Newton method globalized by backtracking. Trial steps are obtained using one of several Krylov subspace methods. The algorithm is implemented in a Fortran solver called NITSOL that is robust yet easy to use and provides a number of useful options and features. The structure offers the user great flexibility in addressing problem specificity through preconditioning and other means and allows easy adaptation to parallel environments. Features and capabilities are illustrated in numerical experiments.

Abstract: A fast, second-order accurate iterative method is proposed for the elliptic equation \[ \grad\cdot(\beta(x,y) \grad u) =f(x,y) \] in a rectangular region $\Omega$ in two-space dimensions. We assume that there is an irregular interface across which the coefficient $\beta$, the solution u and its derivatives, and/or the source term f may have jumps. We are especially interested in the cases where the coefficients $\beta$ are piecewise constant and the jump in $\beta$ is large. The interface may or may not align with an underlying Cartesian grid. The idea in our approach is to precondition the differential equation before applying the immersed interface method proposed by LeVeque and Li [ SIAM J. Numer. Anal., 4 (1994), pp. 1019--1044]. In order to take advantage of fast Poisson solvers on a rectangular region, an intermediate unknown function, the jump in the normal derivative across the interface, is introduced. Our discretization is equivalent to using a second-order difference scheme for a corresponding Poisson equation in the region, and a second-order discretization for a Neumann-like interface condition. Thus second-order accuracy is guaranteed. A GMRES iteration is employed to solve the Schur complement system derived from the discretization. A new weighted least squares method is also proposed to approximate interface quantities from a grid function. Numerical experiments are provided and analyzed. The number of iterations in solving the Schur complement system appears to be independent of both the jump in the coefficient and the mesh size.

Abstract: The standard incomplete LU (ILU) preconditioners often fail for general sparse indefinite matrices because they give rise to "unstable" factors L and U. In such cases, it may be attractive to approximate the inverse of the matrix directly. This paper focuses on approximate inverse preconditioners based on minimizing ||I-AM||F, where AM is the preconditioned matrix. An iterative descent-type method is used to approximate each column of the inverse. For this approach to be efficient, the iteration must be done in sparse mode, i.e., with "sparse-matrix by sparse-vector" operations. Numerical dropping is applied to maintain sparsity; compared to previous methods, this is a natural way to determine the sparsity pattern of the approximate inverse. This paper describes Newton, "global," and column-oriented algorithms, and discusses options for initial guesses, self-preconditioning, and dropping strategies. Some limited theoretical results on the properties and convergence of approximate inverses are derived. Numerical tests on problems from the Harwell--Boeing collection and the FIDAP fluid dynamics analysis package show the strengths and limitations of approximate inverses. Finally, some ideas and experiments with practical variations and applications are presented.

Abstract: Restoration of images that have been blurred by the effects of a Gaussian blurring function is an ill-posed but well-studied problem. Any blur that is spatially invariant can be expressed as a convolution kernel in an integral equation. Fast and effective algorithms then exist for determining the original image by preconditioned iterative methods. If the blurring function is spatially variant, however, then the problem is more difficult. In this work we develop fast algorithms for forming the convolution and for recovering the original image when the convolution functions are spatially variant but have a small domain of support. This assumption leads to a discrete problem involving a banded matrix. We devise an effective preconditioner and prove that the preconditioned matrix differs from the identity by a matrix of small rank plus a matrix of small norm. Numerical examples are given, related to the Hubble Space Telescope (HST) Wide-Field/Planetary Camera. The algorithms that we develop are applicable to other ill-posed integral equations as well.

Abstract: Dissecting Euclidean d -space with the power diagram of n weighted point sites partitions a given m -point set into clusters, one cluster for each region of the diagram. In this manner, an assignment of points to sites is induced. We show the equivalence of such assignments to constrained Euclidean least-squares assignments. As a corollary, there always exists a power diagram whose regions partition a given d -dimensional m -point set into clusters of prescribed sizes, no matter where the sites are placed. Another consequence is that constrained least-squares assignments can be computed by finding suitable weights for the sites. In the plane, this takes roughly O(n 2 m) time and optimal space O(m) , which improves on previous methods. We further show that a constrained least-squares assignment can be computed by solving a specially structured linear program in n+1 dimensions. This leads to an algorithm for iteratively improving the weights, based on the gradient-descent method. Besides having the obvious optimization property, least-squares assignments are shown to be useful in solving a certain transportation problem, and in finding a least-squares fitting of two point sets where translation and scaling are allowed. Finally, we extend the concept of a constrained least-squares assignment to continuous distributions of points, thereby obtaining existence results for power diagrams with prescribed region volumes. These results are related to Minkowski's theorem for convex polytopes. The aforementioned iterative method for approximating the desired power diagram applies to continuous distributions as well.

Abstract: We introduce an economical Gram--Schmidt orthogonalization on the extended Krylov subspace originated by actions of a symmetric matrix and its inverse. An error bound for a family of problems arising from the elliptic method of lines is derived. The bound shows that, for the same approximation quality, the diagonal variant of the extended subspaces requires about the square root of the dimension of the standard Krylov subspaces using only positive or negative matrix powers. An example of an application to the solution of a 2.5-D elliptic problem attests to the computational efficiency of the method for large-scale problems.

Abstract: This paper considers simultaneous gate and wire sizing for general very large scale integrated (VLSI) circuits under the Elmore delay model. We present a fast and exact algorithm which can minimize total area subject to maximum delay bound. The algorithm can be easily modified to give exact algorithms for optimizing several other objectives (e.g., minimizing maximum delay or minimizing total area subject to arrival time specifications at all inputs and outputs). No previous algorithm for simultaneous gate and wire sizing can guarantee exact solutions for general circuits. Our algorithm is an iterative one with a guarantee on convergence to global optimal solutions. It is based on Lagrangian relaxation and "one-gate/wire-at-a-time" greedy optimizations, and is extremely economical and fast. For example, we can optimize a circuit with 27648 gates and wires in 11.53 min using under 23 Mbytes memory on a PC with a 333-MHz Pentium II processor.

Abstract: A new optimization-based iterative learning control algorithm is proposed and its properties derived. An important characteristic of this algorithm is that it uses present and future predicted errors to compute the current control, in a similar manner to model-based predictive control using a receding horizon. In particular, it enables the algorithm designer to achieve good control over convergence rate. The actual implementation has a multimodel structure but uses standard linear quadratic regulator methods for a causal formulation (in the iterative learning sense) of what is originally a non-causal algorithm. The results are illustrated by simulations.

Abstract: Gradient adaptation is a useful technique for adjusting a set of parameters to minimize a cost function. While often easy to implement, the convergence speed of gradient adaptation can be slow when the slope of the cost function varies widely for small changes in the parameters. In this paper, we outline an alternative technique, termed natural gradient adaptation, that overcomes the poor convergence properties of gradient adaptation in many cases. The natural gradient is based on differential geometry and employs knowledge of the Riemannian structure of the parameter space to adjust the gradient search direction. Unlike Newton's method, natural gradient adaptation does not assume a locally-quadratic cost function. Moreover, for maximum likelihood estimation tasks, natural gradient adaptation is asymptotically Fisher-efficient. A simple example illustrates the desirable properties of natural gradient adaptation.

Abstract: An iterative algorithm to determine the dynamic user equilibrium with respect to link costs defined by a traffic simulation model is presented. Each driver's route choice is modeled by a discrete probability distribution which is used to select a route in the simulation. After each simulation run, the probability distribution is adapted to minimize the travel costs. Although the algorithm does not depend on the simulation model, a queuing model is used for performance reasons. The stability of the algorithm is analyzed for a simple example network. As an application example, a dynamic version of Braess's paradox is studied.

Abstract: In bearings-only tracking, observer maneuver is critical to ensure observability and to obtain an accurate target localization. Here, optimal control theory is applied to the determination of the course of a constant speed observer that minimizes an accuracy criterion deduced from the Fisher information matrix (FIM). Necessary conditions for optimal maneuver (Euler equations) are established and resolved, partly by analytical means and partly by an iterative numerical procedure. Examples of optimal observer maneuvers are presented and discussed.

Abstract: A fast accurate iterative reconstruction (FAIR) method suitable for low-statistics positron volume imaging has been developed. The method, based on the expectation maximization-maximum likelihood (EM-ML) technique, operates on list-mode data rather than histogrammed projection data and can, in just one pass through the data, generate images with the same characteristics as several ML iterations. Use of list-mode data preserves maximum sampling accuracy and implicitly ignores lines of response (LORs) in which no counts were recorded. The method is particularly suited to systems where sampling accuracy can be lost by histogramming events into coarse LOR bins, and also to sparse data situations such as fast whole-body and dynamic imaging where sampling accuracy may be compromised by storage requirements and where reconstruction time can be wasted by including LORs with no counts. The technique can be accelerated by operating on subsets of list-mode data which also allows scope for simultaneous data acquisition and iterative reconstruction. The method is compared with a standard implementation of the EM-ML technique and is shown to offer improved resolution, contrast and noise properties as a direct result of using improved spatial sampling, limited only by hardware specifications.

Abstract: . The Picard and modified Picard iteration schemes are often used to numerically solve the nonlinear Richards equation governing water flow in variably saturated porous media. While these methods are easy to implement, they are only linearly convergent. Another approach to solve the Richards equation is to use Newton's iterative method. This method, also known as Newton–Raphson iteration, is quadratically convergent and requires the computation of first derivatives. We implemented Newton's scheme into the mixed form of the Richards equation. As compared to the modified Picard scheme, Newton's scheme requires two additional matrices when the mixed form of the Richards equation is used and requires three additional matrices, when the pressure head-based form is used. The modified Picard scheme may actually be viewed as a simplified Newton scheme. Two examples are used to investigate the numerical performance of different forms of the 1D vertical Richards equation and the different iterative solution schemes. In the first example, we simulate infiltration in a homogeneous dry porous medium by solving both, the h based and mixed forms of Richards equation using the modified Picard and Newton schemes. Results shows that, very small time steps are required to obtain an accurate mass balance. These small times steps make the Newton method less attractive. In a second test problem, we simulate variable inflows and outflows in a heterogeneous dry porous medium by solving the mixed form of the Richards equation, using the modified Picard and Newton schemes. Analytical computation of the Jacobian required less CPU time than its computation by perturbation. A combination of the modified Picard and Newton scheme was found to be more efficient than the modified Picard or Newton scheme.

Abstract: We present a new approach to the least-squares design of stable infinite impulse response (IIR) digital filters. The design is accomplished by using an iterative scheme in which the denominator polynomial obtained from the preceding iteration is treated as a part of the weighting function, and each iteration is carried out by solving a standard quadratic programming problem that yields a stable rational function. When the iteration converges, a stable and truly least-squares solution is obtained. The method is then extended to address the least-squares design of stable IIR two-dimensional (2-D) filters. Examples are included to illustrate the proposed design techniques.

Abstract: A convergence rate is established for nonstationary iterated Tikhonov regularization, applied to ill-posed problems involving closed, densely defined linear operators, under general conditions on the iteration parameters. It is also shown that an order-optimal accuracy is attained when a certain a posteriori stopping rule is used to determine the iteration number.

Abstract: We present a finite-element-based algorithm for reconstruction of fluorescence lifetime and yield in turbid media, using frequency-domain data. The algorithm is based on a set of coupled diffusion equations that describe the propagation of both excitation and fluorescent emission light in multiply scattering media. Centered on Newton's iterative method, we implemented our algorithm by using a synthesized scheme of Marquardt and Tikhonov regularizations. A low-pass spatial filter is also incorporated into the algorithm for enhancing image reconstruction. Simulation studies using both noise-free and noisy data have been performed with the nonzero photon density boundary conditions. Our results suggest that quantitative images can be produced in terms of fluorescent lifetime and yield values and location, size, and shape of heterogeneities within a circular background region.

Abstract: SUMMARY The direct transcription method for solving optimal control problems involves the use of a discrete approximation to the original problem. This paper describes a technique for changing the discretization in order to improve the accuracy of the approximation. An integer programming technique is used to minimize the maximum error during the refinement iterations. The eƒciency of the method is illustrated for an application with path inequality constraints. ( 1998 John Wiley & Sons, Ltd.