MIMS EPrints
Not a member yet
2151 research outputs found
Sort by
Information distance estimation between mixtures of multivariate Gaussians
There are efficient software programs for
extracting from image sequences certain mixtures of distributions, such as
multivariate Gaussians, to
represent the important features needed for accurate document retrieval from databases.
This note describes a method to use information geometric methods to measure distances
between distributions in mixtures of multivariate Gaussians.
There is no general analytic solution for the information
geodesic distance between two k-variate Gaussians,
but for many purposes the absolute information distance is not essential and comparative
values suffice for proximity testing.
For two mixtures of multivariate Gaussians
we must resort to approximations to incorporate the weightings.
In practice, the relation between
a reasonable approximation and a true geodesic distance is likely to be monotonic, which
is adequate for many applications. Here we compare several choices for the incorporation of
weightings in distance estimation and provide illustrative results from simulations of
differently weighted mixtures of multivariate Gaussians
Bounds for the Distance to the Nearest Correlation Matrix
In a wide range of practical problems correlation matrices are formed in such a way that, while symmetry and a unit diagonal are assured, they may lack semidefiniteness. We derive a variety of new upper bounds for the distance from an arbitrary symmetric matrix to the nearest correlation matrix. The bounds are of two main classes: those based on the eigensystem and those based on a modified Cholesky factorization. Bounds from both classes have a computational cost of flops for a matrix of order but are much less expensive to evaluate than the nearest correlation matrix itself. For unit diagonal with for all the eigensystem bounds are shown to overestimate the distance by a factor at most . We show that for a collection of matrices from the literature and from practical applications the eigensystem-based bounds are often good order of magnitude estimates of the actual distance; indeed the best upper bound is never more than a factor larger than a related lower bound. The modified Cholesky bounds are less sharp but also less expensive, and they provide an efficient way to test for definiteness of the putative correlation matrix. Both classes of bounds enable a user to identify an invalid correlation matrix relatively cheaply and to decide whether to revisit its construction or to compute a replacement, such as the nearest correlation matrix
The RKFIT algorithm for nonlinear rational approximation
The RKFIT algorithm outlined in [M. Berljafa and S. Guettel, Generalized rational Krylov decompositions with an application to rational approximation, SIAM J. Matrix Anal. Appl., 2015] is a Krylov-based approach for solving nonlinear rational least squares problems. This paper puts RKFIT into a general framework, allowing for its extension to nondiagonal rational approximants and a family of approximants sharing a common denominator. Furthermore, we derive a strategy for the degree reduction of the approximants, as well as methods for their conversion to partial fraction form, for the efficient evaluation, and root-finding. We also discuss commons and differences of RKFIT and the popular vector fitting algorithm. A MATLAB implementation of RKFIT is provided and numerical experiments, including the fitting of a MIMO dynamical system and an optimization problem related to exponential integration, demonstrate its applicability
High Speed Imaging of Dynamic Processes with a Switched Source X-ray CT System
Conventional x-ray computed tomography (CT) scanners are limited in their scanning speed by the mechanical constraints of their rotating gantries, and as such do not provide the necessary temporal resolution for imaging of fast-moving dynamic processes, such as moving fluid flows. The Real Time Tomography (RTT) system is a family of fast cone beam CT scanners which instead use multiple fixed discrete sources, and complete rings of detectors in an offset geometry. We demonstrate the potential of this system for use in the imaging of such high speed dynamic processes, and give results using simulated and real experimental data. The unusual scanning geometry results in some challenges in image reconstruction, which are overcome using algebraic iterative reconstruction techniques and explicit regularisation. Through the use of a simple temporal regularisation term, and by optimising the source firing pattern, we show that temporal resolution of the system may be increased at the expense of spatial resolution, which may be advantageous in some situations. Results are given showing temporal resolution of approximately 500μs with simulated data, and 3ms with real experimental data
An Efficient Bound for the Condition Number of the Matrix Exponential
A new bound for the condition number of the matrix exponential
is presented. Using the bound, we propose an efficient
approximation to the condition number, denoted by
, that \emph{avoids} the computation of the
Fr\'echet derivative of the matrix exponential that underlies
condition number estimation in the existing algorithms. We
exploit the identity for a nonnegative
integer with the properties of the Fr\'echet
derivative operator to obtain the bound.
Our cost analysis reveals that considerable computational
savings are possible since estimating the condition number by
the existing algorithms requires several invocation of the Fr\'echet
derivative of the matrix exponential whose single invocation
costs as twice as the cost of the matrix exponential itself.
The bound and hence only involve Fr\'echet
derivative of a monomial of degree , which can be computed
exactly in matrix multiplications. We propose two versions
of the scaling and squaring algorithm that implement
.
%that represents the squaring phase of
%the scaling and squaring method.
Our numerical experiments show that captures
the behavior of the condition number and moreover outperforms
the condition number in the estimation of relative forward
errors for a wide range of problems
New Algorithms for Computing the Matrix Sine and Cosine Separately or Simultaneously
Several existing algorithms for computing the matrix cosine employ polynomial or rational approximations combined with scaling and use of a double angle formula. Their derivations are based on forward error bounds. We derive new algorithms for computing the matrix cosine, the matrix sine, and both simultaneously, that are backward stable in exact arithmetic and behave in a forward stable manner in floating point arithmetic. Our new algorithms employ both Pad\'e approximants of and new rational approximants to and obtained from Pad\'e approximants to . The amount of scaling and the degree of the approximants are chosen to minimize the computational cost subject to backward stability in exact arithmetic. Numerical experiments show that the new algorithms have backward and forward errors that rival or surpass those of existing algorithms and are particularly favorable for triangular matrices
Bifurcation from stable fixed point to -dimensional attractor in the border collision normal form
The -dimensional border collision normal form describes bifurcations of piecewise smooth systems. It is shown that there is an open set of parameters such that on one side of the bifurcation the map has a stable fixed point and on the other an attractor with Hausdorff dimension . For generic
parameters this attractor contains open sets and hence has
topological dimension equal to
An Algorithm for the Matrix Lambert W Function
An algorithm is proposed for computing primary matrix Lambert functions of a square matrix , which are solutions of the matrix equation . The algorithm employs the Schur decomposition and blocks the triangular form in such a way that Newton's method can be used on each diagonal block, with a starting matrix depending on the block. A natural simplification of Newton's method for the Lambert function is shown to be numerically unstable. By reorganizing the iteration a new Newton variant is constructed that is proved to be numerically stable. Numerical experiments demonstrate that the algorithm is able to compute the branches of the matrix Lambert function in a numerically reliable way
Parallel Adaptive Importance Sampling
Markov chain Monte Carlo methods are a powerful and commonly used family of numerical methods for sampling from complex probability distributions. As applications of these methods increase in size and complexity, the need for efficient methods which can exploit the parallel architectures which are prevalent in high performance computing increases. In this paper, we aim to develop a framework for scalable parallel MCMC algorithms. At each iteration, an importance sampling proposal distribution is formed using the current states of all of the chains within an ensemble. Once weighted samples have been produced from this, a state-of-the-art resampling method is then used to create an evenly weighted sample ready for the next iteration. We demonstrate that this parallel adaptive importance sampling (PAIS) method outperforms naive parallelisation of serial MCMC methods using the same number of processors, for low dimensional problems, and in fact shows better than linear improvements in convergence rates with respect to the number of processors used