A numerical comparison of solvers for large-scale, continuous-time algebraic Riccati equations and LQR problems

In this paper, we discuss numerical methods for solving large-scale continuous-time algebraic Riccati equations. These methods have been the focus of intensive research in recent years, and significant progress has been made in both the theoretical understanding and efficient implementation of various competing algorithms. There are several goals of this manuscript: first, to gather in one place an overview of different approaches for solving large-scale Riccati equations, and to point to the recent advances in each of them. Second, to analyze and compare the main computational ingredients of these algorithms, to detect their strong points and their potential bottlenecks. And finally, to compare the effective implementations of all methods on a set of relevant benchmark examples, giving an indication of their relative performance

Evaluating parametric holonomic sequences using rectangular splitting

We adapt the rectangular splitting technique of Paterson and Stockmeyer to the problem of evaluating terms in holonomic sequences that depend on a parameter. This approach allows computing the $n$-th term in a recurrent sequence of suitable type using $O(n^{1/2})$ "expensive" operations at the cost of an increased number of "cheap" operations. Rectangular splitting has little overhead and can perform better than either naive evaluation or asymptotically faster algorithms for ranges of $n$ encountered in applications. As an example, fast numerical evaluation of the gamma function is investigated. Our work generalizes two previous algorithms of Smith.Comment: 8 pages, 2 figure

Computing hypergeometric functions rigorously

We present an efficient implementation of hypergeometric functions in arbitrary-precision interval arithmetic. The functions ${}_0F_1$, ${}_1F_1$, ${}_2F_1$ and ${}_2F_0$ (or the Kummer $U$-function) are supported for unrestricted complex parameters and argument, and by extension, we cover exponential and trigonometric integrals, error functions, Fresnel integrals, incomplete gamma and beta functions, Bessel functions, Airy functions, Legendre functions, Jacobi polynomials, complete elliptic integrals, and other special functions. The output can be used directly for interval computations or to generate provably correct floating-point approximations in any format. Performance is competitive with earlier arbitrary-precision software, and sometimes orders of magnitude faster. We also partially cover the generalized hypergeometric function ${}_pF_q$ and computation of high-order parameter derivatives.Comment: v2: corrected example in section 3.1; corrected timing data for case E-G in section 8.5 (table 6, figure 2); adjusted paper siz

Stochastic approximation of score functions for Gaussian processes

We discuss the statistical properties of a recently introduced unbiased stochastic approximation to the score equations for maximum likelihood calculation for Gaussian processes. Under certain conditions, including bounded condition number of the covariance matrix, the approach achieves $O(n)$ storage and nearly $O(n)$ computational effort per optimization step, where $n$ is the number of data sites. Here, we prove that if the condition number of the covariance matrix is bounded, then the approximate score equations are nearly optimal in a well-defined sense. Therefore, not only is the approximation efficient to compute, but it also has comparable statistical properties to the exact maximum likelihood estimates. We discuss a modification of the stochastic approximation in which design elements of the stochastic terms mimic patterns from a $2^n$ factorial design. We prove these designs are always at least as good as the unstructured design, and we demonstrate through simulation that they can produce a substantial improvement over random designs. Our findings are validated by numerical experiments on simulated data sets of up to 1 million observations. We apply the approach to fit a space-time model to over 80,000 observations of total column ozone contained in the latitude band $40^{\circ}$-$50^{\circ}$N during April 2012.Comment: Published in at http://dx.doi.org/10.1214/13-AOAS627 the Annals of Applied Statistics (http://www.imstat.org/aoas/) by the Institute of Mathematical Statistics (http://www.imstat.org

Information Theoretic Study of Gaussian Graphical Models and Their Applications

In many problems we are dealing with characterizing a behavior of a complex stochastic system or its response to a set of particular inputs. Such problems span over several topics such as machine learning, complex networks, e.g., social or communication networks; biology, etc. Probabilistic graphical models (PGMs) are powerful tools that offer a compact modeling of complex systems. They are designed to capture the random behavior, i.e., the joint distribution of the system to the best possible accuracy. Our goal is to study certain algebraic and topological properties of a special class of graphical models, known as Gaussian graphs. First, we show that how Gaussian trees can be used to determine a particular complex system\u27s random behavior, i.e., determining a security robustness of a public communication channel characterized by a Gaussian tree. We show that in such public channels the secrecy capacity of the legitimate users Alice and Bob, in the presence of a passive adversary Eve, is strongly dependent on the underlying structure of the channel. This is done by defining a relevant privacy metric to capture the secrecy capacity of a communication and studying topological and algebraic features of a given Gaussian tree to quantify its security robustness. Next, we examine on how one can effectively produce random samples from such Gaussian tree. The primary concern in synthesis problems is about efficiency in terms of the amount of random bits required for synthesis, as well as the modeling complexity of the given stochastic system through which the Gaussian vector is synthesized. This is done through an optimization problem to propose an efficient algorithm by which we can effectively generate such random vectors. We further generalize the optimization formulation from Gaussian trees to Gaussian vectors with arbitrary structures. This is done by introducing a new latent factor model obtained by solving a constrained minimum determinant factor analysis (CMDFA) problem. We discuss the benefits of factor models in machine learning applications and in particular 3D image reconstruction problems, where our newly proposed CMDFA problem may be beneficial

Two Iterative algorithms for the matrix sign function based on the adaptive filtering technology

In this paper, two new efficient algorithms for calculating the sign function of the large-scale sparse matrix are proposed by combining filtering algorithm with Newton method and Newton Schultz method respectively. Through the theoretical analysis of the error diffusion in the iterative process, we designed an adaptive filtering threshold, which can ensure that the filtering has little impact on the iterative process and the calculation result. Numerical experiments are consistent with our theoretical analysis, which shows that the computational efficiency of our method is much better than that of Newton method and Newton Schultz method, and the computational error is of the same order of magnitude as that of the two methods.Comment: 18 pages,12 figure