Abstract: Large linear systems of saddle point type arise in a wide variety of applications throughout computational science and engineering. Due to their indefiniteness and often poor spectral properties, such linear systems represent a significant challenge for solver developers. In recent years there has been a surge of interest in saddle point problems, and numerous solution techniques have been proposed for this type of system. The aim of this paper is to present and discuss a large selection of solution methods for linear systems in saddle point form, with an emphasis on iterative methods for large and sparse problems.

Abstract: We consider the stabilization problem for a class of networked control systems in the discrete-time domain with random delays. The sensor-to-controller and controller-to-actuator delays are modeled as two Markov chains, and the resulting closed-loop systems are jump linear systems with two modes. The necessary and sufficient conditions on the existence of stabilizing controllers are established. It is shown that state-feedback gains are mode-dependent. An iterative linear matrix inequality (LMI) approach is employed to calculate the state-feedback gains.

Abstract: We present an iterative linear-quadratic-Gaussian method for locally-optimal feedback control of nonlinear stochastic systems subject to control constraints. Previously, similar methods have been restricted to deterministic unconstrained problems with quadratic costs. The new method constructs an affine feedback control law, obtained by minimizing a novel quadratic approximation to the optimal cost-to-go function. Global convergence is guaranteed through a Levenberg-Marquardt method; convergence in the vicinity of a local minimum is quadratic. Performance is illustrated on a limited-torque inverted pendulum problem, as well as a complex biomechanical control problem involving a stochastic model of the human arm, with 10 state dimensions and 6 muscle actuators. A Matlab implementation of the new algorithm is availabe at www.cogsci.ucsd.edu//spl sim/todorov.

Abstract: A maximum likelihood (ML) acoustic source location estimation method is presented for the application in a wireless ad hoc sensor network. This method uses acoustic signal energy measurements taken at individual sensors of an ad hoc wireless sensor network to estimate the locations of multiple acoustic sources. Compared to the existing acoustic energy based source localization methods, this proposed ML method delivers more accurate results and offers the enhanced capability of multiple source localization. A multiresolution search algorithm and an expectation-maximization (EM) like iterative algorithm are proposed to expedite the computation of source locations. The Crame/spl acute/r-Rao Bound (CRB) of the ML source location estimate has been derived. The CRB is used to analyze the impacts of sensor placement to the accuracy of location estimates for single target scenario. Extensive simulations have been conducted. It is observed that the proposed ML method consistently outperforms existing acoustic energy based source localization methods. An example applying this method to track military vehicles using real world experiment data also demonstrates the performance advantage of this proposed method over a previously proposed acoustic energy source localization method.

Abstract: In this correspondence, we consider the problem of maximizing sum rate of a multiple-antenna Gaussian broadcast channel (BC). It was recently found that dirty-paper coding is capacity achieving for this channel. In order to achieve capacity, the optimal transmission policy (i.e., the optimal transmit covariance structure) given the channel conditions and power constraint must be found. However, obtaining the optimal transmission policy when employing dirty-paper coding is a computationally complex nonconvex problem. We use duality to transform this problem into a well-structured convex multiple-access channel (MAC) problem. We exploit the structure of this problem and derive simple and fast iterative algorithms that provide the optimum transmission policies for the MAC, which can easily be mapped to the optimal BC policies.

Abstract: The problem of computing low rank approximations of matrices is considered. The novel aspect of our approach is that the low rank approximations are on a collection of matrices. We formulate this as an optimization problem, which aims to minimize the reconstruction (approximation) error. To the best of our knowledge, the optimization problem proposed in this paper does not admit a closed form solution. We thus derive an iterative algorithm, namely GLRAM, which stands for the Generalized Low Rank Approximations of Matrices. GLRAM reduces the reconstruction error sequentially, and the resulting approximation is thus improved during successive iterations. Experimental results show that the algorithm converges rapidly. We have conducted extensive experiments on image data to evaluate the effectiveness of the proposed algorithm and compare the computed low rank approximations with those obtained from traditional Singular Value Decomposition (SVD) based methods. The comparison is based on the reconstruction error, misclassification error rate, and computation time. Results show that GLRAM is competitive with SVD for classification, while it has a much lower computation cost. However, GLRAM results in a larger reconstruction error than SVD. To further reduce the reconstruction error, we study the combination of GLRAM and SVD, namely GLRAM + SVD, where SVD is preceded by GLRAM. Results show that when using the same number of reduced dimensions, GLRAM + SVD achieves significant reduction of the reconstruction error as compared to GLRAM, while keeping the computation cost low.

Abstract: Preface. 1. MATLAB Usage and Computational Errors. 1.1 Basic Operations of MATLAB. 1.1.1 Input/Output of Data from MATLAB Command Window. 1.1.2 Input/Output of Data Through Files. 1.1.3 Input/Output of Data Using Keyboard. 1.1.4 2-D Graphic Input/Output. 1.1.5 3-D Graphic Output. 1.1.6 Mathematical Functions. 1.1.7 Operations on Vectors and Matrices. 1.1.8 Random Number Generators. 1.1.9 Flow Control. 1.2 Computer Errors Versus Human Mistakes. 1.2.1 IEEE 64-bit Floating-Point Number Representation. 1.2.2 Various Kinds of Computing Errors. 1.2.3 Absolute/Relative Computing Errors. 1.2.4 Error Propagation. 1.2.5 Tips for Avoiding Large Errors. 1.3 Toward Good Program. 1.3.1 Nested Computing for Computational Efficiency. 1.3.2 Vector Operation Versus Loop Iteration. 1.3.3 Iterative Routine Versus Nested Routine. 1.3.4 To Avoid Runtime Error. 1.3.5 Parameter Sharing via Global Variables. 1.3.6 Parameter Passing Through Varargin. 1.3.7 Adaptive Input Argument List. Problems. 2. System of Linear Equations. 2.1 Solution for a System of Linear Equations. 2.1.1 The Nonsingular Case (M = N). 2.1.2 The Underdetermined Case (M N): Least-Squares Error Solution. 2.1.4 RLSE (Recursive Least-Squares Estimation). 2.2 Solving a System of Linear Equations. 2.2.1 Gauss Elimination. 2.2.2 Partial Pivoting. 2.2.3 Gauss-Jordan Elimination. 2.3 Inverse Matrix. 2.4 Decomposition (Factorization). 2.4.1 LU Decomposition (Factorization): Triangularization. 2.4.2 Other Decomposition (Factorization): Cholesky, QR, and SVD. 2.5 Iterative Methods to Solve Equations. 2.5.1 Jacobi Iteration. 2.5.2 Gauss-Seidel Iteration. 2.5.3 The Convergence of Jacobi and Gauss-Seidel Iterations. Problems. 3. Interpolation and Curve Fitting. 3.1 Interpolation by Lagrange Polynomial. 3.2 Interpolation by Newton Polynomial. 3.3 Approximation by Chebyshev Polynomial. 3.4 Pade Approximation by Rational Function. 3.5 Interpolation by Cubic Spline. 3.6 Hermite Interpolating Polynomial. 3.7 Two-dimensional Interpolation. 3.8 Curve Fitting. 3.8.1 Straight Line Fit: A Polynomial Function of First Degree. 3.8.2 Polynomial Curve Fit: A Polynomial Function of Higher Degree. 3.8.3 Exponential Curve Fit and Other Functions. 3.9 Fourier Transform. 3.9.1 FFT Versus DFT. 3.9.2 Physical Meaning of DFT. 3.9.3 Interpolation by Using DFS. Problems. 4. Nonlinear Equations. 4.1 Iterative Method Toward Fixed Point. 4.2 Bisection Method. 4.3 False Position or Regula Falsi Method. 4.4 Newton(-Raphson) Method. 4.5 Secant Method. 4.6 Newton Method for a System of Nonlinear Equations. 4.7 Symbolic Solution for Equations. 4.8 A Real-World Problem. Problems. 5. Numerical Differentiation/Integration. 5.1 Difference Approximation for First Derivative. 5.2 Approximation Error of First Derivative. 5.3 Difference Approximation for Second and Higher Derivative. 5.4 Interpolating Polynomial and Numerical Differential. 5.5 Numerical Integration and Quadrature. 5.6 Trapezoidal Method and Simpson Method. 5.7 Recursive Rule and Romberg Integration. 5.8 Adaptive Quadrature. 5.9 Gauss Quadrature. 5.9.1 Gauss-Legendre Integration. 5.9.2 Gauss-Hermite Integration. 5.9.3 Gauss-Laguerre Integration. 5.9.4 Gauss-Chebyshev Integration. 5.10 Double Integral. Problems. 6. Ordinary Differential Equations. 6.1 Euler's Method. 6.2 Heun's Method: Trapezoidal Method. 6.3 Runge-Kutta Method. 6.4 Predictor-Corrector Method. 6.4.1 Adams-Bashforth-Moulton Method. 6.4.2 Hamming Method. 6.4.3 Comparison of Methods. 6.5 Vector Differential Equations. 6.5.1 State Equation. 6.5.2 Discretization of LTI State Equation. 6.5.3 High-Order Differential Equation to State Equation. 6.5.4 Stiff Equation. 6.6 Boundary Value Problem (BVP). 6.6.1 Shooting Method. 6.6.2 Finite Difference Method. Problems. 7. Optimization. 7.1 Unconstrained Optimization [L-2, Chapter 7]. 7.1.1 Golden Search Method. 7.1.2 Quadratic Approximation Method. 7.1.3 Nelder-Mead Method [W-8]. 7.1.4 Steepest Descent Method. 7.1.5 Newton Method. 7.1.6 Conjugate Gradient Method. 7.1.7 Simulated Annealing Method [W-7]. 7.1.8 Genetic Algorithm [W-7]. 7.2 Constrained Optimization [L-2, Chapter 10]. 7.2.1 Lagrange Multiplier Method. 7.2.2 Penalty Function Method. 7.3 MATLAB Built-In Routines for Optimization. 7.3.1 Unconstrained Optimization. 7.3.2 Constrained Optimization. 7.3.3 Linear Programming (LP). Problems. 8. Matrices and Eigenvalues. 8.1 Eigenvalues and Eigenvectors. 8.2 Similarity Transformation and Diagonalization. 8.3 Power Method. 8.3.1 Scaled Power Method. 8.3.2 Inverse Power Method. 8.3.3 Shifted Inverse Power Method. 8.4 Jacobi Method. 8.5 Physical Meaning of Eigenvalues/Eigenvectors. 8.6 Eigenvalue Equations. Problems. 9. Partial Differential Equations. 9.1 Elliptic PDE. 9.2 Parabolic PDE. 9.2.1 The Explicit Forward Euler Method. 9.2.2 The Implicit Backward Euler Method. 9.2.3 The Crank-Nicholson Method. 9.2.4 Two-Dimensional Parabolic PDE. 9.3 Hyperbolic PDE. 9.3.1 The Explicit Central Difference Method. 9.3.2 Two-Dimensional Hyperbolic PDE. 9.4 Finite Element Method (FEM) for solving PDE. 9.5 GUI of MATLAB for Solving PDEs: PDETOOL. 9.5.1 Basic PDEs Solvable by PDETOOL. 9.5.2 The Usage of PDETOOL. 9.5.3 Examples of Using PDETOOL to Solve PDEs. Problems. Appendix A: Mean Value Theorem. Appendix B: Matrix Operations/Properties. Appendix C: Differentiation with Respect to a Vector. Appendix D: Laplace Transform. Appendix E: Fourier Transform. Appendix F: Useful Formulas. Appendix G: Symbolic Computation. Appendix H: Sparse Matrices. Appendix I: MATLAB. References. Subject Index. Index for MATLAB Routines. Index for Tables.

Abstract: In this note, we apply a hierarchical identification principle to study solving the Sylvester and Lyapunov matrix equations. In our approach, we regard the unknown matrix to be solved as system parameters to be identified, and present a gradient iterative algorithm for solving the equations by minimizing certain criterion functions. We prove that the iterative solution consistently converges to the true solution for any initial value, and illustrate that the rate of convergence of the iterative solution can be enhanced by choosing the convergence factor (or step-size) appropriately. Furthermore, the iterative method is extended to solve general linear matrix equations. The algorithms proposed require less storage capacity than the existing numerical ones. Finally, the algorithms are tested on computer and the results verify the theoretical findings.

Abstract: Numerical crack propagation schemes were augmented in an elegant manner by the X-FEM method. The use of special tip enrichment functions, as well as a discontinuous function along the sides of the crack allows one to do a complete crack analysis virtually without modifying the underlying mesh, which is of industrial interest, especially when a numerical model for crack propagation is desired. This paper improves the implementation of the X-FEM method for stress analysis around cracks in three ways. First, the enrichment strategy is revisited. The conventional approach uses a 'topological' enrichment (only the elements touching the front are enriched). We suggest a 'geometrical' enrichment in which a given domain size is enriched. The improvements obtained with this enrichment are discussed. Second, the conditioning of the X-FEM both for topological and geometrical enrichments is studied. A preconditioner is introduced so that 'off the shelf' iterative solver packages can be used and perform as well on X-FEM matrices as on standard FEM matrices. The preconditioner uses a local (nodal) Cholesky based decomposition. Third, the numerical integration scheme to build the X-FEM stiffness matrix is dramatically improved for tip enrichment functions by the use of an ad hoc integration scheme. A 2D benchmark problem is designed to show the improvements and the robustness.

Abstract: Shuffled versions of iterative decoding of low-density parity-check codes and turbo codes are presented. The proposed schemes have about the same computational complexity as the standard versions, and converge faster. Simulations show that the new schedules offer better performance/complexity tradeoffs, especially when the maximum number of iterations has to remain small.

Abstract: For the augmented system of linear equations, Golub, Wu and Yuan recently studied an SOR-like method (BIT 41(2001)71---85). By further accelerating it with another parameter, in this paper we present a generalized SOR (GSOR) method for the augmented linear system. We prove its convergence under suitable restrictions on the iteration parameters, and determine its optimal iteration parameters and the corresponding optimal convergence factor. Theoretical analyses show that the GSOR method has faster asymptotic convergence rate than the SOR-like method. Also numerical results show that the GSOR method is more effective than the SOR-like method when they are applied to solve the augmented linear system. This GSOR method is further generalized to obtain a framework of the relaxed splitting iterative methods for solving both symmetric and nonsymmetric augmented linear systems by using the techniques of vector extrapolation, matrix relaxation and inexact iteration. Besides, we also demonstrate a complete version about the convergence theory of the SOR-like method.

Abstract: In recent years, kernel principal component analysis (KPCA) has been suggested for various image processing tasks requiring an image model such as, e.g., denoising or compression. The original form of KPCA, however, can be only applied to strongly restricted image classes due to the limited number of training examples that can be processed. We therefore propose a new iterative method for performing KPCA, the kernel Hebbian algorithm, which iteratively estimates the kernel principal components with only linear order memory complexity. In our experiments, we compute models for complex image classes such as faces and natural images which require a large number of training examples. The resulting image models are tested in single-frame super-resolution and denoising applications. The KPCA model is not specifically tailored to these tasks; in fact, the same model can be used in super-resolution with variable input resolution, or denoising with unknown noise characteristics, in spite of this, both super-resolution and denoising performance are comparable to existing methods.

Abstract: Recent work has looked at extending the k-Means algorithm to incorporate background information in the form of instance level must-link and cannot-link constraints We introduce two ways of specifying additional background information in the form of δ and constraints that operate on all instances but which can be interpreted as conjunctions or disjunctions of instance level constraints and hence are easy to implement We present complexity results for the feasibility of clustering under each type of constraint individually and several types together A key finding is that determining whether there is a feasible solution satisfying all constraints is, in general, NP-complete Thus, an iterative algorithm such as k-Means should not try to find a feasible partitioning at each iteration This motivates our derivation of a new version of the k-Means algorithm that minimizes the constrained vector quantization error but at each iteration does not attempt to satisfy all constraints Using standard UCI datasets, we find that using constraints improves accuracy as others have reported, but we also show that our algorithm reduces the number of iterations until convergence Finally, we illustrate these benefits and our new constraint types on a complex real world object identification problem using the infra-red detector on an Aibo robot

Abstract: We consider the stabilization problem for a kind of networked control systems in discrete-time domain with random delays. The sensor-to-controller and controller-to-sensor delays are modeled as two Markov chains, and the resulting closed-loop systems are jump linear systems with two modes. The necessary and sufficient conditions on the existence of the stabilizing controllers are established. It is shown that the state-feedback gains are different with different modes. An iterative linear matrix inequality (LMI) approach is employed to calculate the state-feedback gains.

Abstract: Error-propagation phenomena and computational complexity of the filters' design are important drawbacks of existing decision-feedback equalizers (DFE) for dispersive channels. In this paper, we propose a new iterative block DFE (IBDFE) which operates iteratively on blocks of the received signal. Indeed, a suitable data-transmission format must be used to allow an efficient implementation of the equalizer in the frequency domain, by means of the discrete Fourier transform. Two design methods are considered. In the first method, hard detected data are used as input of the feedback, and filters are designed according to the correlation between detected and transmitted data. In the second method, the feedback signal is directly designed from soft detection of the equalized signal at the previous iteration. Estimators of the parameters involved in the IBDFE design are also derived. From performance simulations on a wireless dispersive fading channel, we observed that the IBDFE outperforms existing DFEs. Moreover, the IBDFE exhibits a reduction of the computational complexity when compared against existing schemes, both in signal processing and in filter design.

Abstract: This paper addresses the compensation of the dynamics-coupling effect in piezoscanners used for positioning in atomic force microscopes (AFMs). Piezoscanners are used to position the AFM probe, relative to the sample, both parallel to the sample surface (x and y axes) and perpendicular to the sample surface (z axis). In this paper, we show that dynamics-coupling from the scan axes (x and y axes) to the perpendicular z axis can generate significant positioning errors during high-speed AFM operation, i.e., when the sample is scanned at high speed. We use an inversion-based iterative control approach to compensate for this dynamics-coupling effect. Convergence of the iterative approach is investigated and experimental results show that the dynamics-coupling-caused error can be reduced, close to the noise level, using the proposed approach. Thus, the main contribution of this paper is the development of an approach to substantially reduce the dynamics-coupling-caused error and thereby, to enable high-speed AFM operation.

Abstract: Initial conditions, or initial resetting conditions, play a fundamental role in all kinds of iterative learning control methods. In this note, we study five different initial conditions, disclose the inherent relationship between each initial condition and corresponding learning convergence (or boundedness) property. The iterative learning control method under consideration is based on Lyapunov theory, which is suitable for plants with time-varying parametric uncertainties and local Lipschitz nonlinearities.

Abstract: The Lanczos method is an iterative procedure to compute an orthogonal basis for the Krylov subspace generated by a symmetric matrix $A$ and a starting vector $v$. An interesting application of this method is the computation of the matrix exponential $\exp(-\tau A)v$. This vector plays an important role in the solution of parabolic equations where $A$ results from some form of discretization of an elliptic operator. In the present paper we will argue that for these applications the convergence behavior of this method can be unsatisfactory. We will propose a modified method that resolves this by a simple preconditioned transformation at the cost of an inner-outer iteration. A priori error bounds are presented that are independent of the norm of $A$. This shows that the worst case convergence speed is independent of the mesh width in the spatial discretization of the elliptic operator. We discuss, furthermore, a posteriori error estimation and the tuning of the coupling between the inner and outer iteration. We conclude with several numerical experiments with the proposed method.

Abstract: In this paper, we present a sequence of iterative methods improving Newton's method for solving nonlinear equations. The Adomian decomposition method is applied to an equivalent coupled system to construct the sequence of the methods whose order of convergence increases as it progresses. The orders of convergence are derived analytically, and then rederived by applying symbolic computation of Maple. Some numerical illustrations are given.

Abstract: Algebraic multigrid (AMG) is a very efficient iterative solver and preconditioner for large unstructured sparse linear systems. Traditional coarsening schemes for AMG can, however, lead to computational complexity growth as problem size increases, resulting in increased memory use and execution time, and diminished scalability. Two new parallel AMG coarsening schemes are proposed that are based solely on enforcing a maximum independent set property, resulting in sparser coarse grids. The new coarsening techniques remedy memory and execution time complexity growth for various large three-dimensional (3D) problems. If used within AMG as a preconditioner for Krylov subspace methods, the resulting iterative methods tend to converge fast. This paper discusses complexity issues that can arise in AMG, describes the new coarsening schemes, and examines the performance of the new preconditioners for various large 3D problems.

Abstract: In this correspondence, we consider the problem of maximizing sum rate of a multiple-antenna Gaussian broadcast channel (BC). It was recently found that dirty-paper coding is capacity achieving for this channel. In order to achieve capacity, the optimal transmission policy (i.e., the optimal transmit covariance structure) given the channel conditions and power constraint must be found. However, obtaining the optimal transmission policy when employing dirty-paper coding is a computationally complex nonconvex problem. We use duality to transform this problem into a well-structured convex multiple-access channel (MAC) problem. We exploit the structure of this problem and derive simple and fast iterative algorithms that provide the optimum transmission policies for the MAC, which can easily be mapped to the optimal BC policies.

Abstract: This paper presents a reduced-complexity soft-input soft-output detection scheme, called iterative tree search detection, for multiple-input multiple-output wireless communication systems employing turbo processing at the receiver. In this scheme, a reduced search space is selected with the aid of the M-algorithm, and QAM signal constellations with block partitionable labels are used in order to make the detection complexity per bit almost independent of the modulation order, as well as asymptotically linear in the number of transmit antennas. Results from computer simulations are presented which demonstrate the capability of the scheme to approach optimal performance at considerably reduced complexity.

Abstract: An effective method for optimizing the performance of a fixed current distribution, uniformly spaced antenna array has been to adjust its element positions to provide performance improvement. In comparison with the default uniform structure, this approach yields performance improvements such as smaller sidelobe levels or beamwidth values. Additionally, it provides practical advantages such as reductions in size, weight and number of antenna elements. The objective of this paper is to describe a unified mathematical approach to nonlinear optimization of multidimensional array geometries. The approach utilizes a class of limiting properties of sinusoidal, Bessel or Legendre functions that are dictated by the array geometry addressed. The efficacy of the method is demonstrated by its generalized application to synthesis of rectangular, cylindrical and spherical arrays. The unified mathematical approach presented below is a synthesis technique founded on the mathematical transformation of the desired field pattern, followed by the application of limiting forms of the transformation, and resulting in the development of a closed form expression for the element positions. The method offers the following advantages over previous techniques such as direct nonlinear optimization or genetic algorithms. First, it is not an iterative, searching algorithm, and provides element spacing values directly in a single run of the algorithm, thereby saving valuable CPU time and memory storage. Second, It permits the array designer to place practical constraints on the array geometry, (e.g., the minimum/maximum spacing between adjacent elements)

Abstract: The boundary element method has become a popular tool for the solution of Maxwell's equations in electromagnetism. From a linear algebra point of view, this leads to the solution of large dense complex linear systems, where the unknowns are associated with the edges of the mesh defined on the surface of the illuminated object. In this paper, we address the iterative solution of these linear systems via preconditioned Krylov solvers. Our primary focus is on the design of an efficient parallelizable preconditioner. In that respect, we consider an approximate inverse method based on the Frobenius-norm minimization. The preconditioner is constructed from a sparse approximation of the dense coefficient matrix, and the patterns both for the preconditioner and for the coefficient matrix are computed a priori using geometric information from the mesh. We describe how such a preconditioner can be naturally implemented in a parallel code that implements the multipole technique for the matrix-vector product calculation. We investigate the numerical scalability of our preconditioner on realistic industrial test problems and show that it exhibits some limitations on very large problems of size close to one million unknowns. To improve its robustness on those large problems we propose an embedded iterative scheme that combines nested GMRES solvers with different fast multipole computations. We show through extensive numerical experiments that this new scheme is extremely robust at affordable memory and CPU costs for the solution of very large and challenging problems.