Computation Of Matrix-vector Products Research Articles

PDMM: A Novel Primal-Dual Majorization-Minimization Algorithm for Poisson Phase-Retrieval Problem

In this paper, we introduce a novel iterative algorithm for the problem of phase-retrieval where the measurements consist of only the magnitude of linear function of the unknown signal, and the noise in the measurements follow Poisson distribution. The proposed algorithm is based on the principle of majorization-minimization (MM); however, the application of MM here is very novel and distinct from the way MM has been usually used to solve optimization problems in the literature. More precisely, we reformulate the original minimization problem into a saddle point problem by invoking Fenchel dual representation of the <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"><tex-math notation="LaTeX">$\log (\cdot)$</tex-math></inline-formula> term in the Poisson likelihood function. We then propose tighter surrogate functions over both primal and dual variables resulting in a double-loop MM algorithm, which we have named as Primal-Dual Majorization-Minimization (PDMM) algorithm. The iterative steps of the resulting algorithm are simple to implement and involve only computing matrix vector products. We also extend our algorithm to handle various <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"><tex-math notation="LaTeX">$\ell _{1}$</tex-math></inline-formula> regularized Poisson phase-retrieval problems (which exploit sparsity). The proposed algorithm is compared with previously proposed algorithms such as wirtinger flow (WF), MM (conventional), and alternating direction methods of multipliers (ADMM) for the Poisson data model. The simulation results under different experimental settings show that PDMM is faster than the competing methods, and its performance in recovering the original signal is at par with the state-of-the-art algorithms.

“Short-Dot”: Computing Large Linear Transforms Distributedly Using Coded Short Dot Products

We consider the problem of computing a matrix-vector product $Ax$ using a set of $P$ parallel or distributed processing nodes prone to “straggling,” i.e. , unpredictable delays. Every processing node can access only a fraction $({s}/{N})$ of the $N$ -length vector $x$ , and all processing nodes compute an equal number of dot products. We propose a novel error correcting code-that we call “Short-Dot”–that introduces redundant, shorter dot products such that only a subset of the nodes’ outputs are sufficient to compute $Ax$ . To address the problem of straggling in computing matrix-vector products, prior work uses replication or erasure coding to encode parts of the matrix $A$ , but the length of the dot products computed at each processing node is still $N$ . The key novelty in our work is that instead of computing the long dot products as required in the original matrix-vector product, we construct a larger number of redundant and short dot products that only require a fraction of $x$ to be accessed during the computation. Short-Dot is thus useful in a communication-constrained scenario as it allows for only a fraction of $x$ to be accessed by each processing node. Further, we show that in the particular regime where the number of available processing nodes is greater than the total number of dot products, Short-Dot has lower expected computation time under straggling under an exponential model compared to existing strategies, e.g. replication, in a scaling sense. We also derive fundamental limits on the trade-off between the length of the dot products and the recovery threshold, i.e., the required number of processing nodes, showing that Short-Dot is near-optimal.

PFOSM: An efficient algorithm for aerodynamic robust design based on continuous adjoint and matrix-vector products

This article proposes a novel robust design tool based on the Method of Moments (MoM), demonstrated in problems governed by laminar flows of incompressible fluids. The objective function of the robust design problem consists of the weighted sum of the mean value and standard deviation of the Quantity of Interest (QoI), computed based on the value and gradient of the latter with respect to the uncertain variables. First-order derivatives with respect to the uncertain and design variables are computed using the continuous adjoint method. Using a first-order MoM to compute the objective function and by minimizing it with a gradient-based approach leads to the need of computing up to second-order derivatives of the QoI with respect to the design and uncertain variables. This can lead to a high computational cost in problems with a great number of uncertain variables. The proposed robust design algorithm avoids computing second-order derivatives by computing matrix-vector products instead, using a combination of adjoint and direct differentiation. This makes the cost per gradient computation independent of the numbers of both the design and uncertain variables. The proposed method is demonstrated through the constrained shape optimization of an isolated 2D airfoil under uncertain farfield flow conditions.

Efficient computation of matrix power-vector products: Application for space-fractional diffusion problems

A novel algorithm is proposed for computing matrix–vector products Aαv, where A is a symmetric positive semidefinite sparse matrix and α>0. The method can be applied for the efficient implementation of the matrix transformation method to solve space-fractional diffusion problems. The performance of the new algorithm is studied in a comparison with the conventional MATLAB subroutines to compute matrix powers.

Using monomer vibrational wavefunctions to compute numerically exact (12D) rovibrational levels of water dimer.

We compute numerically exact rovibrational levels of water dimer, with 12 vibrational coordinates, on the accurate CCpol-8sf ab initio flexible monomer potential energy surface [C. Leforestier et al., J. Chem. Phys. 137, 014305 (2012)]. It does not have a sum-of-products or multimode form and therefore quadrature in some form must be used. To do the calculation, it is necessary to use an efficient basis set and to develop computational tools, for evaluating the matrix-vector products required to calculate the spectrum, that obviate the need to store the potential on a 12D quadrature grid. The basis functions we use are products of monomer vibrational wavefunctions and standard rigid-monomer basis functions (which involve products of three Wigner functions). Potential matrix-vector products are evaluated using the F matrix idea previously used to compute rovibrational levels of 5-atom and 6-atom molecules. When the coupling between inter- and intra-monomer coordinates is weak, this crude adiabatic type basis is efficient (only a few monomer vibrational wavefunctions are necessary), although the calculation of matrix elements is straightforward. It is much easier to use than an adiabatic basis. The product structure of the basis is compatible with the product structure of the kinetic energy operator and this facilitates computation of matrix-vector products. Compared with the results obtained using a [6 + 6]D adiabatic approach, we find good agreement for the inter-molecular levels and larger differences for the intra-molecular water bend levels.

A tree-parenchyma coupled model for lung ventilation simulation.

In this article, we develop a lung ventilation model. The parenchyma is described as an elastic homogenized media. It is irrigated by a space-filling dyadic resistive pipe network, which represents the tracheobronchial tree. In this model, the tree and the parenchyma are strongly coupled. The tree induces an extra viscous term in the system constitutive relation, which leads, in the finite element framework, to a full matrix. We consider an efficient algorithm that takes advantage of the tree structure to enable a fast matrix-vector product computation. This framework can be used to model both free and mechanically induced respiration, in health and disease. Patient-specific lung geometries acquired from computed tomography scans are considered. Realistic Dirichlet boundary conditions can be deduced from surface registration on computed tomography images. The model is compared to a more classical exit compartment approach. Results illustrate the coupling between the tree and the parenchyma, at global and regional levels, and how conditions for the purely 0D model can be inferred. Different types of boundary conditions are tested, including a nonlinear Robin model of the surrounding lung structures.

A fast algorithm for solving the space–time fractional diffusion equation

In this paper, we propose a fast algorithm for efficient and accurate solution of the space–time fractional diffusion equations defined in a rectangular domain. The spatial discretization is done by using the central finite difference scheme and matrix transfer technique. Due to its nonlocality, numerical discretization of the spectral fractional Laplacian (−Δ)sα/2 results in a large dense matrix. This causes considerable challenges not only for storing the matrix but also for computing matrix–vector products in practice. By utilizing the compact structure of the discrete system and the discrete sine transform, our algorithm avoids to store the large matrix from discretizing the nonlocal operator and also significantly reduces the computational costs. We then use the Laplace transform method for time integration of the semi-discretized system and a weighted trapezoidal method to numerically compute the convolutions needed in the resulting scheme. Various experiments are presented to demonstrate the efficiency and accuracy of our method.

On the structured backward error of inexact Arnoldi methods for (skew)-Hermitian and (skew)-symmetric eigenvalue problems

In this paper, we present the inexact structure preserving Arnoldi methods for structured eigenvalue problems. They are called structure preserving because the computed eigenvalues and eigenvectors can preserve the desirable properties of the structures of the original matrices, even with large errors involved in the computation of matrix-vector products. A backward error matrix is called structured if it has the same structure as the original matrix. We derive a common form for the structured backward errors that can be used for different structure preserving processes, and prove the derived form has the minimum Frobenius norm among all possible backward errors. Furthermore, we employ the derived backward errors for some specific structure preserving processes to estimate the accuracy of the solutions obtained by inexact Arnoldi methods for eigenvalue problems. We aim to give, wherever possible, formulae that are inexpensive to compute so that they can be used routinely in practice. Numerical experiments are provided to support the theoretical results.

Automatic generation of fast algorithms for matrix–vector multiplication

ABSTRACTThis paper describes the methods for finding fast algorithms for computing matrix–vector products including the procedures based on the block-structured matrices. The proposed methods involve an analysis of the structural properties of matrices. The presented approaches are based on the well-known optimization techniques: the simulated annealing and the hill-climbing algorithm along with its several extensions. The main idea of the proposed methods consists in finding a decomposition of the original matrix into a sparse matrix and a matrix corresponding to an appropriate block-structured pattern. The main criterion for optimizing is a reduction of the computational cost. The methods presented in this paper can be successfully implemented in many digital signal processing tasks.

A Distributed-Memory Package for Dense Hierarchically Semi-Separable Matrix Computations Using Randomization

We present a distributed-memory library for computations with dense structured matrices. A matrix is considered structured if its off-diagonal blocks can be approximated by a rank-deficient matrix with low numerical rank. Here, we use Hierarchically Semi-Separable (HSS) representations. Such matrices appear in many applications, for example, finite-element methods, boundary element methods, and so on. Exploiting this structure allows for fast solution of linear systems and/or fast computation of matrix-vector products, which are the two main building blocks of matrix computations. The compression algorithm that we use, that computes the HSS form of an input dense matrix, relies on randomized sampling with a novel adaptive sampling mechanism. We discuss the parallelization of this algorithm and also present the parallelization of structured matrix-vector product, structured factorization, and solution routines. The efficiency of the approach is demonstrated on large problems from different academic and industrial applications, on up to 8,000 cores. This work is part of a more global effort, the STRUctured Matrices PACKage (STRUMPACK) software package for computations with sparse and dense structured matrices. Hence, although useful on their own right, the routines also represent a step in the direction of a distributed-memory sparse solver.

Performance of Preconditioned Nonstationary Methods for Electromagnetic Scattering From One Dimensional Dielectric Rough Surfaces

For the analysis of electromagnetic scattering from one dimensional dielectric randomly rough surfaces, we apply the right-preconditioning technique to three nonstationary iterative algorithms, namely, the generalized minimal residual (GMRES), conjugate gradient squared (CGS), and biconjugate gradient stabilized (BiCGSTAB). The spectral acceleration technique is also utilized to expedite computation of matrix-vector product. Among the nonstationary methods considered, with the right-preconditioning (RP) applied, GMRES-RP is the most robust, followed by BiCGSTAB-RP, with CGS-RP falling behind. In terms of number of iterations to achieve convergence, for cases when they all converge, GMRES-RP requires the most whereas BiCGSTAB-RP the least, with CGS-RP slightly more demanding than BiCGSTAB-RP. However, in terms of run time, GMRES-RP is the most efficient, followed by BiCGSTAB-RP, and CGS-RP is the least efficient. This relative performance is in contrast to the unpreconditioned case as reported in the literature, which speaks of the impact of the adopted right-preconditioning on the spectral distribution of these nonstationary algorithms.

Solving Scattering by Multilayer Dielectric Objects Using JMCFIE-DDM-MLFMA

In this letter, a hybrid approach combining nonconformal, nonoverlapping domain decomposition method (DDM) and electric/magnetic current combined-field integral equation (JMCFIE) is presented to solve scattering problems of three-dimensional multilayer dielectric objects. The original structure can be decomposed into some closed subdomains according to material distribution, and JMCFIE is formulated in each subdomain. The major ingredients in the proposed approach are the use of the Robin-type transmission conditions (TCs) and nonconformal, nonoverlapping domain partitioning, which remarkably increases the flexibility of mesh discretization and enhances the efficiency of solution. The multilevel fast multipole algorithm (MLFMA) is also applied to realize fast computation of matrix-vector product. Finally, some numerical results are given to demonstrate the validity of the JMCFIE-DDM-MLFMA.

Performance assessment of exponential Rosenbrock methods for large systems of ODEs

This article studies time integration methods for stiff systems of ordinary differential equations of large dimension. For such problems, implicit methods generally outperform explicit methods because the step size is usually less restricted by stability constraints. Recently, however, a family of explicit methods, called exponential integrators, have become popular for large stiff problems due to their favourable stability properties and the rapid convergence of non-preconditioned Krylov subspace methods for computing matrix-vector products involving exponential-like functions of the Jacobian matrix. In this article, we implement the so-called exponential Rosenbrock methods using Krylov subspaces. Numerical experiments on a challenging real-world test problem reveal that these methods are a promising preconditioner-free alternative to well-established approaches based on preconditioned Newton--Krylov implementations of the backward differentiation formulas. References E. J. Carr, T. J. Moroney and I. W. Turner. Efficient simulation of unsaturated flow using exponential time integration. Appl. Math. Comput., 217(14): 6587--6596, 2011. doi:10.1016/j.amc.2011.01.041 E. J. Carr, I. W. Turner and M. Ilic. Krylov subspace approximations for the exponential Euler method: error estimates and the harmonic Ritz approximant. ANZIAM Journal, 52: C612--C627, 2011. http://journal.austms.org.au/ojs/index.php/ANZIAMJ/article/view/3938 E. J. Carr, I. W. Turner and P. Perre. A variable step size Jacobian-free exponential integrator for simulating transport in heterogeneous porous media: application to wood drying. J. Comput. Phys., 233: 66-82, 2013. doi:10.1016/j.jcp.2012.07.024 Sundials user documentation. https://computation.llnl.gov/casc/sundials/documentation/documentation.html M. Hochbruck, C. Lubich and H. Selhofer. Exponential integrators for large systems of differential equations. SIAM J. Sci. Comput., 19(5): 1552--1574, 1998. doi:10.1137/S1064827595295337 M. Hochbruck, A. Ostermann and J. Schweitzer. Exponential Rosenbrock-type methods. SIAM J. Numer. Anal., 47(1): 786--803, 2009. doi:10.1137/080717717 C. T. Kelley. Solving nonlinear equations with Newton's method. SIAM, USA, 2003. D. A. Knoll and D. E. Keyes. Jacobian-free Newton--Krylov methods: a survey of approaches and applications. J. Comput. Phys., 193(2): 357--397, 2004. doi:10.1016/j.jcp.2003.08.010 M. Tokman. Efficient integration of large stiff systems of ODEs with exponential propagation iteration (EPI) methods. J. Comput. Phys., 213: 748--776, 2006. doi:10.1016/j.jcp.2005.08.032

Parameter estimation in high dimensional Gaussian distributions

In order to compute the log-likelihood for high dimensional Gaussian models, it is necessary to compute the determinant of the large, sparse, symmetric positive definite precision matrix. Traditional methods for evaluating the log-likelihood, which are typically based on Cholesky factorisations, are not feasible for very large models due to the massive memory requirements. We present a novel approach for evaluating such likelihoods that only requires the computation of matrix-vector products. In this approach we utilise matrix functions, Krylov subspaces, and probing vectors to construct an iterative numerical method for computing the log-likelihood.

A Robust Preconditioned GMRES Method for Electromagnetic Scattering From Dielectric Rough Surfaces

In this paper, we propose a generalized minimal residual procedure which is right preconditioned with the forward-backward method (FBM) for the analysis of electromagnetic scattering from 1-D dielectric randomly rough surfaces. The spectral acceleration (SA) technique is also combined to expedite the computation of matrix-vector product. For lossy dielectric, the Green's function is attenuative and suitably truncated, without the need to apply the SA technique for the lower medium. The applied preconditioning transforms the original linear system from near singular to stable with a good condition number. The spectrum of the preconditioned matrix is found condensed in the vicinity of the point 1 in the complex plane, which speaks of the good approximation quality of the preconditioner to the original matrix. Moreover, the construction of the preconditioner does not require the knowledge of the distribution of the impedance matrix spectrum, which means that the proposed method can be used as a general purpose iterative solver. The proposed method has demonstrated the desirable properties of a numerical algorithm: robust and efficient. Regarding robustness, the proposed method significantly improves convergence over the FBM-SA, in particular, for exponentially correlated rough surfaces. Regarding efficiency, the proposed method runs two to four times faster than the FBM-SA for horizontal (HH) polarization, and is about twice faster for vertical (VV) polarization. These features indicate that the proposed method can be effectively used to analyze electromagnetic scattering from 1-D dielectric Gaussian surfaces with both Gaussian and exponential spectra.

Improved Model-Based Parameter Estimation Approach for Accelerated Periodic Method of Moments Solutions With Application to the Analysis of Convoluted Frequency Selected Surfaces and Metamaterials

An improved ¿smart¿ interpolation approach known as model-based parameter estimation (MBPE) is applied to the wide-band interpolation of periodic method of moments (PMM) impedance matrices for normal and oblique incidence cases. Prior to interpolation, easy to calculate but hard to interpolate, phase terms are removed from the impedance matrices. An efficient spectral-domain PMM formulation is introduced for the accelerated analysis of frequency selective surface (FSS) problems with a large number of unknowns, employing a one dimensional O(N log N) FFT-based method to speed up the computation of matrix-vector products within the bi-conjugate gradient (BCG) iterative solver, which is made possible by the asymmetric multilevel block-Toeplitz structure of the impedance-matrix. The MBPE interpolation algorithm provides a faster matrix fill time than the brute force method and is comparable or even faster than the 2-D FFT-based method for a large number of unknowns. It also has the advantage that it can be applied to non-uniform gridding cases. The accuracy and efficiency of the proposed techniques for large FSS problems are demonstrated by several design examples for both the normal and oblique incidence cases. We also apply this efficient analysis tool to the design of multiband single-layer FSS filters and artificial magnetic conductors (AMC) comprised of a 2-D periodic arrangement of convoluted metallic strips in the shape of a Hilbert curve. The multiband properties of the Hilbert curve FSS filters are studied for different iteration orders (i.e., different degrees of space-filling).

Some properties of generalized K-centrosymmetric H-matrices

Every n × n generalized K-centrosymmetric matrix A can be reduced into a 2 × 2 block diagonal matrix (see [Z. Liu, H. Cao, H. Chen, A note on computing matrix–vector products with generalized centrosymmetric (centrohermitian) matrices, Appl. Math. Comput. 169 (2) (2005) 1332–1345]). This block diagonal matrix is called the reduced form of the matrix A. In this paper we further investigate some properties of the reduced form of these matrices and discuss the square roots of these matrices. Finally exploiting these properties, the development of structure-preserving algorithms for certain computations for generalized K-centrosymmetric H-matrices is discussed.

A note on computing matrix–vector products with generalized centrosymmetric (centrohermitian) matrices

To reduce the costs of computing matrix–vector product Ax related to a centrosymmetric matrix A as compared to the case of an arbitrary matrix A, two algorithms were proposed recently, one was designed by Melman [A. Melman, Symmetric centrosymmetric matrix–vector multiplication, Linear Algebra Appl. 320 (2000) 193–198] for symmetric centrosymmetric matrices, another was presented by Fassbender and Ikramov [H. Fassbender, K.D. Ikramov, Computing matrix–vector products with centrosymmetric and centrohermitian matrices, Linear Algebra Appl. 364 (2003) 235–241] for general centrosymmetric matrices. In this note we further discuss this topic of computing Ax, where A is a generalized centrosymmetric matrix. We firstly investigate the reducibility of a generalized centrosymmetric matrix and then provide an algorithm which can be viewed as a generalization of one [H. Fassbender, K.D. Ikramov, Computing matrix–vector products with centrosymmetric and centrohermitian matrices, Linear Algebra Appl. 364 (2003) 235–241]. We show that our algorithm is suitable for computing many matrix–vector products with the same matrix. Furthermore, we show that similar results can be obtained for certain subclasses of generalized centrosymmetric matrices, as compared to the case of a centrosymmetric matrix, and analogue gains are available for generalized skew-centrosymmetric or generalized centrohermitian matrices.

Low-Frequency Model-Order Reduction of Electromagnetic Fields Without Matrix Factorization

In this paper, we develop a reduced-order modeling technique, which is based on a low-frequency expansion of the electromagnetic field. The expansion can be written in terms of the pseudoinverse of a so-called system matrix. This pseudoinverse is given explicitly, and it is shown that it satisfies a reciprocity relation. Moreover, we show that computing matrix-vector products with this pseudoinverse essentially amounts to repeatedly solving Poisson's equation. The latter two properties allow us to efficiently compute reduced-order models via a Lanczos-type algorithm. The proposed method is illustrated by a number of numerical examples.

Multipole-based preconditioners for large sparse linear systems

Dense operators for preconditioning sparse linear systems have traditionally been considered infeasible due to their excessive computational and memory requirements. With the emergence of techniques such as block low-rank approximations and hierarchical multipole approximations, the cost of computing and storing these preconditioners has reduced dramatically. This paper describes the use of multipole operators as parallel preconditioners for sparse linear systems. Hierarchical multipole approximations of explicit Green’s functions are effective preconditioners due to their bounded-error properties. By enumerating nodes in proximity preserving order, one can achieve high parallel efficiency in computing matrix–vector products with these dense preconditioners. The benefits of the approach are illustrated on the Poisson problem and the generalized Stokes problem arising in incompressible fluid flow simulations. Numerical experiments show that the multipole-based techniques are effective preconditioners that can be parallelized efficiently on multiprocessing platforms.