Search

In this paper, we apply the Anderson acceleration technique to the existing relaxation fixed‐point iteration for solving the multilinear PageRank. In order to reduce computational cost, we further consider the periodical version of the Anderson acceleration. The convergence of the proposed algorithms is discussed. Numerical experiments on synthetic and real‐world datasets are performed to demonstrate the advantages of the proposed algorithms over the relaxation fixed‐point iteration and the extrapolated shifted fixed‐point method. In particular, we give a strategy for choosing the quasi‐optimal parameters of the associated algorithms when they are applied to solve the test problems with different sizes but the same structure. [ABSTRACT FROM AUTHOR]

The goal of the present paper is the design of embeddings of a general sparse graph into a set of points in ℝd for appropriate d ≥ 2. The embeddings that we are looking at aim to keep vertices that are grouped in communities together and keep the rest apart. To achieve this property, we utilize coarsening that respects possible community structures of the given graph. We employ a hierarchical multilevel coarsening approach that identifies communities (strongly connected groups of vertices) at every level. The multilevel strategy allows any given (presumably expensive) graph embedding algorithm to be made into a more scalable (and faster) algorithm. We demonstrate the presented approach on a number of given embedding algorithms and large‐scale graphs and achieve speed‐up over the methods in a recent paper. [ABSTRACT FROM AUTHOR]

In this paper, enlightened by the effectiveness of Anderson acceleration (AA), a well‐established technique for accelerating fixed‐point solvers, we first present the Anderson accelerating modulus‐based matrix splitting (AAMS) algorithms for a class of horizontal nonlinear complementarity problems. Then, by introducing the strong semi‐smoothness of the absolute value function, we establish the local convergence theory of the AAMS algorithms. More importantly, we provide the optimal parameter for the AAMS algorithms, which is independent of iteration and applicable to other reduced MS algorithms. Finally, the numerical experiments are conducted to clarify the efficiency and practicability of the AAMS algorithms and its optimal parameter, and the effects of two parameters (both derived from AA) of the AAMS algorithms are analyzed. [ABSTRACT FROM AUTHOR]

This paper presents the two‐parameter scaling memoryless Broyden–Fletcher–Goldfarb–Shanno (BFGS) method for solving a system of monotone nonlinear equations. The optimal values of the scaling parameters are obtained by minimizing the measure function involving all the eigenvalues of the memoryless BFGS matrix. The optimal values can be used in the analysis of the quasi‐Newton method for ill‐conditioned matrices. This algorithm can also be described as a combination of the projection technique and memoryless BGFS method. Global convergence of the method is provided. For validation and efficiency of the scheme, some test problems are computed and compared with existing results. [ABSTRACT FROM AUTHOR]

This paper is concerned with improving the empirical convergence speed of block‐coordinate descent algorithms for approximate nonnegative tensor factorization (NTF). We propose an extrapolation strategy in‐between block updates, referred to as heuristic extrapolation with restarts (HER). HER significantly accelerates the empirical convergence speed of most existing block‐coordinate algorithms for NTF, in particular for challenging computational scenarios, while requiring a negligible additional computational budget. [ABSTRACT FROM AUTHOR]

Summary: This paper presents some recent advances for parallel‐in‐time methods applied to linear elasticity. With recent computer architecture changes leading to stagnant clock speeds, but ever increasing numbers of cores, future speedups will be available through increased concurrency. Thus, sequential algorithms, such as time stepping, will suffer a bottleneck. This paper explores multigrid reduction in time (MGRIT) for an important application area, linear elasticity. Previously, efforts at parallel‐in‐time for elasticity have experienced difficulties, for example, the beating phenomenon. As a result, practical parallel‐in‐time algorithms for this application area currently do not exist. This paper proposes some solutions made possible by MGRIT (e.g., slow temporal coarsening and FCF‐relaxation) and, more importantly, a different formulation of the problem that is more amenable to parallel‐in‐time methods. Using a recently developed convergence theory for MGRIT and Parareal, we show that the changed formulation of the problem avoids the instability issues and allows the reduction of the error using two temporal grids. We then extend our approach to the multilevel case, where we demonstrate how slow temporal coarsening improves convergence. The paper ends with supporting numerical results showing a practical algorithm enjoying speedup benefits over the sequential algorithm. [ABSTRACT FROM AUTHOR]

Summary: This paper studies a low‐communication algorithm for solving elliptic PDEs on high‐performance machines, the nested iteration with range decomposition (NIRD) algorithm. Previous work has shown that NIRD converges to a high level of accuracy within a small, fixed number of iterations (usually one or two) when applied to simple elliptic problems. This paper makes some improvements to the NIRD algorithm (including the addition of adaptivity during preprocessing, wider choice of partitioning functions, and modified error measurement) that enhance the method's accuracy and scalability, especially on more difficult problems. In addition, an updated convergence proof is presented based on heuristic assumptions that are supported by numerical evidence. Furthermore, a new performance model is developed that shows increased performance benefits for NIRD when problems are more expensive to solve using traditional methods. Finally, extensive testing on a variety of elliptic problems provides additional insight into the behavior of NIRD and additional evidence that NIRD achieves excellent convergence on a wide class of elliptic PDEs and, as such, should be a very competitive method for solving PDEs on large parallel computers. [ABSTRACT FROM AUTHOR]

Tensor decompositions such as the canonical format and the tensor train format have been widely utilized to reduce storage costs and operational complexities for high‐dimensional data, achieving linear scaling with the input dimension instead of exponential scaling. In this paper, we investigate even lower storage‐cost representations in the tensor ring format, which is an extension of the tensor train format with variable end‐ranks. Firstly, we introduce two algorithms for converting a tensor in full format to tensor ring format with low storage cost. Secondly, we detail a rounding operation for tensor rings and show how this requires new definitions of common linear algebra operations in the format to obtain storage‐cost savings. Lastly, we introduce algorithms for transforming the graph structure of graph‐based tensor formats, with orders of magnitude lower complexity than existing literature. The efficiency of all algorithms is demonstrated on a number of numerical examples, and in certain cases, we demonstrate significantly higher compression ratios when compared to previous approaches to using the tensor ring format. [ABSTRACT FROM AUTHOR]

In the present paper, we describe an adaptive modified rational global Lanczos algorithm for model-order reduction problems using multipoint moment matching-based methods. The major problem of these methods is the selection of some interpolation points. We first propose a modified rational global Lanczos process and then we derive Lanczos-like equations for the global case. Next, we propose adaptive techniques for choosing the interpolation points. Second-order dynamical systems are also considered in this paper, and the adaptive modified rational global Lanczos algorithm is applied to an equivalent state space model. Finally, some numerical examples will be given. [ABSTRACT FROM AUTHOR]

Summary: The tensor SVD (t‐SVD) for third‐order tensors, previously proposed in the literature, has been applied successfully in many fields, such as computed tomography, facial recognition, and video completion. In this paper, we propose a method that extends a well‐known randomized matrix method to the t‐SVD. This method can produce a factorization with similar properties to the t‐SVD, but it is more computationally efficient on very large data sets. We present details of the algorithms and theoretical results and provide numerical results that show the promise of our approach for compressing and analyzing image‐based data sets. We also present an improved analysis of the randomized and simultaneous iteration for matrices, which may be of independent interest to the scientific community. We also use these new results to address the convergence properties of the new and randomized tensor method as well. [ABSTRACT FROM AUTHOR]

Summary: In this paper, we study the nearest stable matrix pair problem: given a square matrix pair (E,A), minimize the Frobenius norm of (ΔE,ΔA) such that (E+ΔE,A+ΔA) is a stable matrix pair. We propose a reformulation of the problem with a simpler feasible set by introducing dissipative Hamiltonian matrix pairs: A matrix pair (E,A) is dissipative Hamiltonian if A=(J−R)Q with skew‐symmetric J, positive semidefinite R, and an invertible Q such that QTE is positive semidefinite. This reformulation has a convex feasible domain onto which it is easy to project. This allows us to employ a fast gradient method to obtain a nearby stable approximation of a given matrix pair. [ABSTRACT FROM AUTHOR]

Summary: In this paper, we propose a fast algorithm for computing the spectral radii of symmetric nonnegative tensors. In particular, by this proposed algorithm, we are able to obtain the spectral radii of weakly reducible symmetric nonnegative tensors without requiring the partition of the tensors. As we know, it is very costly to determine the partition for large‐sized weakly reducible tensors. Numerical results are reported to show that the proposed algorithm is efficient and also able to compute the spectral radii of large‐sized tensors. As an application, we present an algorithm for testing the positive definiteness of Z‐tensors. By this algorithm, it is guaranteed to determine the positive definiteness for any Z‐tensor. [ABSTRACT FROM AUTHOR]

The contour integral-based eigensolvers are the recent efforts for computing the eigenvalues inside a given region in the complex plane. The best-known members are the Sakurai-Sugiura method, its stable version CIRR, and the FEAST algorithm. An attractive computational advantage of these methods is that they are easily parallelizable. The FEAST algorithm was developed for the generalized Hermitian eigenvalue problems. It is stable and accurate. However, it may fail when applied to non-Hermitian problems. Recently, a dual subspace FEAST algorithm was proposed to extend the FEAST algorithm to non-Hermitian problems. In this paper, we instead use the oblique projection technique to extend FEAST to the non-Hermitian problems. Our approach can be summarized as follows: (a) construct a particular contour integral to form a search subspace containing the desired eigenspace and (b) use the oblique projection technique to extract desired eigenpairs with appropriately chosen test subspace. The related mathematical framework is established. Comparing to the dual subspace FEAST algorithm, we can save the computational cost roughly by a half if only the eigenvalues or the eigenvalues together with their right eigenvectors are needed. We also address some implementation issues such as how to choose a suitable starting matrix and design-efficient stopping criteria. Numerical experiments are provided to illustrate that our method is stable and efficient. [ABSTRACT FROM AUTHOR]

SUMMARY In a recent paper an algorithm for large-scale Tikhonov regularization in standard form called GKB-FP was proposed and numerically illustrated. In this paper, further insight into the convergence properties of this method is provided, and extensions to general-form Tikhonov regularization are introduced. In addition, as alternative to Tikhonov regularization, a preconditioned LSQR method coupled with an automatic stopping rule is proposed. Preconditioning seeks to incorporate smoothing properties of the regularization matrix into the computed solution. Numerical results are reported to illustrate the methods on large-scale problems. Copyright © 2013 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

SUMMARY In 1907, Erhard Schmidt published a paper in which he introduced an orthogonalization algorithm that has since become known as the classical Gram-Schmidt process. Schmidt claimed that his procedure was essentially the same as an earlier one published by J. P. Gram in 1883. The Schmidt version was the first to become popular and widely used. An algorithm related to a modified version of the process appeared in an 1820 treatise by P. S. Laplace. Although related algorithms have been around for almost 200 years, it is the Schmidt paper that led to the popularization of orthogonalization techniques. The year 2007 marked the 100th anniversary of that paper. In celebration of that anniversary, we present a comprehensive survey of the research on Gram-Schmidt orthogonalization and its related QR factorization. Its application for solving least squares problems and in Krylov subspace methods are also reviewed. Software and implementation aspects are also discussed. Copyright © 2012 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

The ChebFilterCG algorithm, proposed by Golub, Ruiz, and Touhami [SIAM J. Matrix Anal. Appl. 29 (2007), pp. 774-795] is an iterative method that combines Chebyshev filter and conjugate gradient for solving symmetric positive definite linear systems with multiple right-hand sides. The Chebyshev filter is used to produce initial residuals rich in eigenvectors corresponding to the smallest eigenvalues, which are then used in the initial phase of the conjugate gradient. This paper presents a convergence analysis of ChebFilterCG. In particular, it is shown theoretically and numerically that the algorithm yields an approximation of the invariant subspace associated with the smallest eigenvalues that can be recycled for solving several linear systems with the same matrix and different right-hand sides. A refined error bound when solving these systems is also given. The choice and influence of the Chebyshev filtering steps is discussed. Numerical experiments are described to illustrate that the Chebyshev filter does not degrade the distribution of the smallest eigenvalues and highlight the effect of rounding errors when large outlying eigenvalues are present. Finally, it is shown that the method may become more effective when an additional Chebyshev filtering step is used in the initialization phase of ChebFilterCG. [ABSTRACT FROM AUTHOR]

SUMMARY This paper describes the development of explicit time stepping methods for linear PDEs that are specifically designed to cope with the stiffness of the system of ODEs that results from spatial discretization. As stiffness is caused by the contrasting behavior of coupled components of the solution, it is proposed to adopt a componentwise approach in which each coefficient of the solution in an appropriate basis is computed using an individualized approximation of the solution operator. This has been accomplished by Krylov subspace spectral (KSS) methods, which use techniques from 'matrices, moments and quadrature' to approximate bilinear forms involving functions of matrices via block Gaussian quadrature rules. These forms correspond to coefficients with respect to the chosen basis of the application of the solution operator of the PDE to the solution at an earlier time. In this paper, it is proposed to substantially enhance the efficiency of KSS methods through the prescription of quadrature nodes on the basis of asymptotic analysis of the recursion coefficients produced by block Lanczos iteration for each Fourier coefficient as a function of frequency. The potential of this idea is illustrated through numerical results obtained from the application of the modified KSS methods to diffusion equations and wave equations. Copyright © 2012 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

In this paper, two accelerated divide-and-conquer (ADC) algorithms are proposed for the symmetric tridiagonal eigenvalue problem, which cost O(N2r) flops in the worst case, where N is the dimension of the matrix and r is a modest number depending on the distribution of eigenvalues. Both of these algorithms use hierarchically semiseparable (HSS) matrices to approximate some intermediate eigenvector matrices, which are Cauchy-like matrices and are off-diagonally low-rank. The difference of these two versions lies in using different HSS construction algorithms, one (denoted by ADC1) uses a structured low-rank approximation method and the other (ADC2) uses a randomized HSS construction algorithm. For the ADC2 algorithm, a method is proposed to estimate the off-diagonal rank. Numerous experiments have been carried out to show their stability and efficiency. These algorithms are implemented in parallel in a shared memory environment, and some parallel implementation details are included. Comparing the ADCs with highly optimized multithreaded libraries such as Intel MKL, we find that ADCs could be more than six times faster for some large matrices with few deflations. [ABSTRACT FROM AUTHOR]

We consider the nonlinear eigenvalue problem M(λ)x= 0, where M(λ) is a large parameter-dependent matrix. In several applications, M(λ) has a structure where the higher-order terms of its Taylor expansion have a particular low-rank structure. We propose a new Arnoldi-based algorithm that can exploit this structure. More precisely, the proposed algorithm is equivalent to Arnoldi's method applied to an operator whose reciprocal eigenvalues are solutions to the nonlinear eigenvalue problem. The iterates in the algorithm are functions represented in a particular structured vector-valued polynomial basis similar to the construction in the infinite Arnoldi method [Jarlebring, Michiels, and Meerbergen, Numer. Math., 122 (2012), pp. 169- 195]. In this paper, the low-rank structure is exploited by applying an additional operator and by using a more compact representation of the functions. This reduces the computational cost associated with orthogonalization, as well as the required memory resources. The structure exploitation also provides a natural way in carrying out implicit restarting and locking without the need to impose structure in every restart. The efficiency and properties of the algorithm are illustrated with two large-scale problems. [ABSTRACT FROM AUTHOR]

SUMMARY Bootstrap algebraic multigrid (BAMG) is a multigrid-based solver for matrix equations of the form Ax = b. Its aim is to automatically determine the interpolation weights used in algebraic multigrid by locally fitting a set of test vectors that have been relaxed as solutions to the corresponding homogeneous equation, Ax = 0. This paper studies an improved form of BAMG, called relaxation-corrected bootstrap algebraic multigrid ( rBAMG), that involves adding scaled residuals of the test vectors to the least-squares equations. The basic rBAMG scheme was introduced in an earlier paper [1] and analyzed on a simple model problem. The purpose of the current paper is to further develop this algorithm by incorporating several new critical components and to systematically study its performance on an interesting model problem from quantum chromodynamics. Whereas the earlier paper introduced a new least-squares principle involving the residuals of the test vectors, a simple extrapolation scheme is developed here to accurately estimate the convergence factors of the evolving algebraic multigrid solver. Such a capability is essential to the effective development of a fast solver, and the approach introduced here is shown numerically to be much more effective than the conventional approach of just observing successive error reduction factors. Another component of the setup process developed here is an adaptive cycling process. This component assesses the effectiveness of the V-cycle constructed in the initial rBAMG phase by applying it to the homogeneous equation. When poor convergence is observed, the set of test vectors is enhanced with the resulting error, enabling the subsequent least-squares fit of interpolation to produce an improved V-cycle. A related component is the scaling and recombination Ritz process that targets the so-called weak approximation property in an attempt to reveal the important elements of these evolving error and test vector spaces. The aim of the numerical study documented here is to provide insight into the various design choices that arise in the development of an rBAMG algorithm. With this in mind, the results for quantum chromodynamics focus on the behavior of rBAMG in terms of the number of initial test vectors used, the number of relaxation sweeps applied to them, and the size of the target matrices. Copyright © 2012 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

SUMMARY This special issue contains a selection of papers from the Sixth International Workshop on the Numerical Solution of Markov Chains, held in Williamsburg, Virginia on September 16-17, 2010. The papers cover a broad range of topics including perturbation theory for absorbing chains, bounding techniques, steady-state and transient solution methods, multilevel algorithms, preconditioning, and applications. Copyright © 2011 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

SUMMARY In this paper, we revisit the steady-state solution method for Markov Regenerative Processes (MRP) proposed in the work by German. This method solves the embedded Markov chain P of the MRP without storing the matrix P explicitly. We address three issues left open in German's Work: 1) the solution method is restricted to Power method; 2) it has been defined only for ergodic MRPs; and 3) no preconditioning is available to speed-up the computation. This paper discusses how to lift these limitations by extending the algorithm to preconditioned Krylov-subspace methods and by generalizing it to the non-ergodic case. An MRP-specific preconditioner is also proposed, which is built from a sparse approximation of the MRP matrix, computed via simulation. An experimental assessment of the proposed preconditioner is then provided. Copyright © 2011 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

SUMMARY For a Markov reward process, where upper and lower bounds for the transition rates and rewards are known, a new approach to bound the expected reward is presented. Based on a previous paper where sharp bounds have been defined for the problem, but only an inefficient and unstable algorithm is proposed, this paper presents algorithms to compute the bounds by interpreting the problem as a Markov Decision Process. In this way, the well known value and policy iteration algorithms can be adopted to compute reward bounds in a stable and fairly efficient way. Different numerical techniques are presented for computing the reward bounds. Copyright © 2011 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

This paper is concerned with the tight enclosure of matrix multiplication AB for two floating-point matrices A and B. The aim of this paper is to compute component-wise upper and lower bounds of the exact result C of the matrix multiplication AB by floating-point arithmetic. Namely, an interval matrix enclosing C is obtained. In this paper, new algorithms for enclosing C are proposed. The proposed algorithms are designed to mainly exploit the level 3 operations in BLAS. Although the proposed algorithms take around twice as much costs as a standard algorithm promoted by Oishi and Rump, the accuracy of the result by the proposed algorithms is better than that of the standard algorithm. At the end of this paper, we present numerical examples showing the efficiency of the proposed algorithms. Copyright © 2010 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

This paper concerns the LBM factorization of unsymmetric tridiagonal matrices, where L and M are unit lower triangular matrices and B is block diagonal with 1×1 and 2×2 blocks. In some applications, it is necessary to form this factorization without row or column interchanges while the tridiagonal matrix is formed. Bunch and Kaufman proposed a pivoting strategy without interchanges specifically for symmetric tridiagonal matrices, and more recently, Bunch and Marcia proposed pivoting strategies that are normwise backward stable for linear systems involving such matrices. In this paper, we extend these strategies to the unsymmetric tridiagonal case and demonstrate that the proposed methods both exhibit bounded growth factors and are normwise backward stable. Copyright © 2009 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

Preconditioners based on various multilevel extensions of two-level finite element methods (FEM) lead to iterative methods which have an optimal order computational complexity with respect to the size of the system. Such methods were first presented in Axelsson and Padiy (SIAM. J. Sci. Stat. Comp. 1990; 20:1807) and Axelsson and Vassilevski (Numer. Math. 1989; 56:157), and are based on (recursive) two-level splittings of the finite element space. The key role in the derivation of optimal convergence rate estimates is played by the constant γ in the so-called Cauchy–Bunyakowski–Schwarz (CBS) inequality, associated with the angle between the two subspaces of the splitting. It turns out that only existence of uniform estimates for this constant is not enough but accurate quantitative bounds for γ have to be found as well. More precisely, the value of the upper bound for γ∈(0,1) is part of the construction of various multilevel extensions of the related two-level methods. In this paper, an algebraic two-level preconditioning algorithm for second-order elliptic boundary value problems is constructed, where the discretization is done using Crouzeix–Raviart non-conforming linear finite elements on triangles. An important point to make is that in this case the finite element spaces corresponding to two successive levels of mesh refinements are not nested. To handle this, a proper two-level basis is considered, which enables us to fit the general framework for the construction of two-level preconditioners for conforming finite elements and to generalize the method to the multilevel case. The major contribution of this paper is the derived estimates of the related constant γ in the strengthened CBS inequality. These estimates are uniform with respect to both coefficient and mesh anisotropy. To our knowledge, the results presented in the paper are the first such estimates for non-conforming FEM systems. Copyright © 2004 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

Helmholtz equations with variable coefficients are known to be hard to solve both analytically and numerically. In this paper, we introduce a numerical multigrid solver for one-dimensional Helmholtz eigenvalue problems with periodic potentials and solutions. The solvers employ wave–ray methodology suggested by Brandt, Livshits for Helmholtz equations with constant coefficients. The paper concludes with numerical experiments and the discussion of future efforts for solving two-dimensional problems. Copyright © 2004 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

Bramble et al. (Math. Comput. (2000); 69 (230) : 667–689) discussed the Uzawa type algorithms for iteratively solving nonsymmetric stable saddle point problems. In this paper, we consider the Uzawa type algorithms for iteratively solving generalized saddle point problems. Our main concern is how to accelerate the convergence of the inexact Uzawa type algorithms. The contributions of the paper are that a new non-linear Uzawa type algorithm is presented and its convergence is analysed. Applications to Navier–Stokes equations by mixed finite element discretization with an unstable pair of approximation spaces, i.e. Q1-P0 pair of approximation spaces, are discussed and the results of the numerical experiments of applying our new algorithm are presented. Copyright © 2003 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

The Jacobi-Davidson (JD) algorithm is considered one of the most efficient eigensolvers currently available for non-Hermitian problems. It can be viewed as a coupled inner-outer iteration, where the inner one expands the search subspace and the outer one reduces the eigenpair residual. One of the difficulties in the JD efficient use stems from the definition of the most appropriate inner tolerance, so as to avoid useless extra work and keep the number of outer iterations under control. To this aim, the use of an efficient preconditioner for the inner iterative solver is of paramount importance. The present paper describes a fresh implementation of the JD algorithm with controlled inner iterations and block factorized sparse approximate inverse preconditioning for non-Hermitian eigenproblems in a parallel computational environment. The algorithm performance is investigated by comparison with a freely available software package such as SLEPc. The results show that combining the inner tolerance control with an efficient preconditioning technique can allow for a significant improvement of the JD performance, preserving a good scalability. Copyright © 2016 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

This paper deals with studying some of well-known iterative methods in their tensor forms to solve a Sylvester tensor equation. More precisely, the tensor form of the Arnoldi process and full orthogonalization method are derived by using a product between two tensors. Then tensor forms of the conjugate gradient and nested conjugate gradient algorithms are also presented. Rough estimation of the required number of operations for the tensor form of the Arnoldi process is obtained, which reveals the advantage of handling the algorithms based on tensor format over their classical forms in general. Some numerical experiments are examined, which confirm the feasibility and applicability of the proposed algorithms in practice. Copyright © 2016 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

In this paper, we study a class of weakly nonlinear complementarity problems arising from the discretization of free boundary problems. By reformulating the complementarity problems as implicit fixed-point equations based on splitting of the system matrices, we propose a class of modulus-based matrix splitting algorithms. We show their convergence by assuming that the system matrix is positive definite. Moreover, we give several kinds of typical practical choices of the modulus-based matrix splitting iteration methods based on the different splitting of the system matrix. Numerical experiments on two model problems are presented to illustrate the theoretical results and examine the numerical effectiveness of our modulus-based matrix splitting algorithms. Copyright © 2016 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

Among numerous iterative methods for solving the minimal nonnegative solution of an M-matrix algebraic Riccati equation, the structure-preserving doubling algorithm (SDA) stands out owing to its overall efficiency as well as accuracy. SDA is globally convergent and its convergence is quadratic, except for the critical case for which it converges linearly with the linear rate 1/2. In this paper, we first undertake a delineatory convergence analysis that reveals that the approximations by SDA can be decomposed into two components: the stable component that converges quadratically and the rank-one component that converges linearly with the linear rate 1/2. Our analysis also shows that as soon as the stable component is fully converged, the rank-one component can be accurately recovered. We then propose an efficient hybrid method, called the two-phase SDA, for which the SDA iteration is stopped as soon as it is determined that the stable component is fully converged. Therefore, this two-phase SDA saves those SDA iterative steps that previously have to have for the rank-one component to be computed accurately, and thus essentially, it can be regarded as a quadratically convergent method. Numerical results confirm our analysis and demonstrate the efficiency of the new two-phase SDA. Copyright © 2015 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

We revisit the shift-and-invert Arnoldi method proposed in [S. Lee, H. Pang, and H. Sun. Shift-invert Arnoldi approximation to the Toeplitz matrix exponential, SIAM J. Sci. Comput., 32: 774-792, 2010] for numerical approximation to the product of Toeplitz matrix exponential with a vector. In this approach, one has to solve two large-scale Toeplitz linear systems in advance. However, if the desired accuracy is high, the cost will be prohibitive. Therefore, it is interesting to investigate how to solve the Toeplitz systems inexactly in this method. The contribution of this paper is in three regards. First, we give a new stability analysis on the Gohberg-Semencul formula (GSF) and define the GSF condition number of a Toeplitz matrix. It is shown that when the size of the Toeplitz matrix is large, our result is sharper than the one given in [M. Gutknecht and M. Hochbruck. The stability of inversion formulas for Toeplitz matrices, Linear Algebra Appl., 223/224: 307-324, 1995]. Second, we establish a relation between the error of Toeplitz systems and the residual of Toeplitz matrix exponential. We show that if the GSF condition number of the Toeplitz matrix is medium-sized, then the Toeplitz systems can be solved in a low accuracy. Third, based on this relationship, we present a practical stopping criterion for relaxing the accuracy of the Toeplitz systems and propose an inexact shift-and-invert Arnoldi algorithm for the Toeplitz matrix exponential problem. Numerical experiments illustrate the numerical behavior of the new algorithm and show the effectiveness of our theoretical results. Copyright © 2015 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

The use of matchings is a powerful technique for scaling and ordering sparse matrices prior to the solution of a linear system A x = b. Traditional methods such as implemented by the HSL software package MC64 use the Hungarian algorithm to solve the maximum weight maximum cardinality matching problem. However, with advances in the algorithms and hardware used by direct methods for the parallelization of the factorization and solve phases, the serial Hungarian algorithm can represent an unacceptably large proportion of the total solution time for such solvers. Recently, auction algorithms and approximation algorithms have been suggested as alternatives for achieving near-optimal solutions for the maximum weight maximum cardinality matching problem. In this paper, the efficacy of auction and approximation algorithms as replacements for the Hungarian algorithm is assessed in the context of sparse symmetric direct solvers when used in problems arising from a range of practical applications. High-cardinality suboptimal matchings are shown to be as effective as optimal matchings for the purposes of scaling. However, matching-based ordering techniques require that matchings are much closer to optimality before they become effective. The auction algorithm is demonstrated to be capable of finding such matchings significantly faster than the Hungarian algorithm, but our [ABSTRACT FROM AUTHOR]

In this paper, we consider the solution of linear systems of saddle point type by correcting the Uzawa algorithm, which has been proposed in [K. Arrow, L. Hurwicz, H. Uzawa, Studies in nonlinear programming, Stanford University Press, Stanford, CA, 1958]. We call this method as corrected Uzawa (CU) method. The convergence of the CU method is analyzed for solving nonsingular saddle point problem as well as the semi-convergence for the singular case. First, the corrected model for the Uzawa algorithm is established, and the CU algorithm is presented. Then we study the geometric meaning of the CU model. Moreover, we introduce the overall reduction coefficient α to measure the effect of the CU process. It is shown that the CU method converges faster than the Uzawa method and several other methods if the overall reduction coefficient α satisfies certain conditions. Numerical experiments are presented to illustrate the theoretical results and examine the numerical effectiveness of the CU method. Copyright © 2015 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

Alternating least squares (ALS) is often considered the workhorse algorithm for computing the rank- R canonical tensor approximation, but for certain problems, its convergence can be very slow. The nonlinear conjugate gradient (NCG) method was recently proposed as an alternative to ALS, but the results indicated that NCG is usually not faster than ALS. To improve the convergence speed of NCG, we consider a nonlinearly preconditioned NCG (PNCG) algorithm for computing the rank- R canonical tensor decomposition. Our approach uses ALS as a nonlinear preconditioner in the NCG algorithm. Alternatively, NCG can be viewed as an acceleration process for ALS. We demonstrate numerically that the convergence acceleration mechanism in PNCG often leads to important pay-offs for difficult tensor decomposition problems, with convergence that is significantly faster and more robust than for the stand-alone NCG or ALS algorithms. We consider several approaches for incorporating the nonlinear preconditioner into the NCG algorithm that have been described in the literature previously and have met with success in certain application areas. However, it appears that the nonlinearly PNCG approach has received relatively little attention in the broader community and remains underexplored both theoretically and experimentally. Thus, this paper serves several additional functions, by providing in one place a concise overview of several PNCG variants and their properties that have only been described in a few places scattered throughout the literature, by systematically comparing the performance of these PNCG variants for the tensor decomposition problem, and by drawing further attention to the usefulness of nonlinearly PNCG as a general tool. In addition, we briefly discuss the convergence of the PNCG algorithm. In particular, we obtain a new convergence result for one of the PNCG variants under suitable conditions, building on known convergence results for non-preconditioned NCG. Copyright © 2014 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

In this paper, we present a new algorithm to accelerate the Chambolle gradient projection method for total variation image restoration. The new proposed method considers an approximation of the Hessian based on the secant equation. Combined with the quasi-Cauchy equations and diagonal updating, we can obtain a positive definite diagonal matrix. In the proposed minimization method model, we use the positive definite diagonal matrix instead of the constant time stepsize in Chambolle's method. The global convergence of the proposed scheme is proved. Some numerical results illustrate the efficiency of this method. Moreover, we also extend the quasi-Newton diagonal updating method to solve nonlinear systems of monotone equations. Performance comparisons show that the proposed method is efficient. A practical application of the monotone equations is shown and tested on sparse signal reconstruction in compressed sensing. Copyright © 2015 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

SUMMARY This paper introduces a new strategy for setting the regularization parameter when solving large-scale discrete ill-posed linear problems by means of the Arnoldi-Tikhonov method. This new rule is essentially based on the discrepancy principle, although no initial knowledge of the norm of the error that affects the right-hand side is assumed; an increasingly more accurate approximation of this quantity is recovered during the Arnoldi algorithm. Some theoretical estimates are derived in order to motivate our approach. Many numerical experiments performed on classical test problems as well as image deblurring problems are presented. Copyright © 2014 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

SUMMARY Novel memory-efficient Arnoldi algorithms for solving matrix polynomial eigenvalue problems are presented. More specifically, we consider the case of matrix polynomials expressed in the Chebyshev basis, which is often numerically more appropriate than the standard monomial basis for a larger degree d. The standard way of solving polynomial eigenvalue problems proceeds by linearization, which increases the problem size by a factor d. Consequently, the memory requirements of Krylov subspace methods applied to the linearization grow by this factor. In this paper, we develop two variants of the Arnoldi method that build the Krylov subspace basis implicitly, in a way that only vectors of length equal to the size of the original problem need to be stored. The proposed variants are generalizations of the so-called quadratic Arnoldi method and two-level orthogonal Arnoldi procedure methods, which have been developed for the monomial case. We also show how the typical ingredients of a full implementation of the Arnoldi method, including shift-and-invert and restarting, can be incorporated. Numerical experiments are presented for matrix polynomials up to degree 30 arising from the interpolation of nonlinear eigenvalue problems, which stem from boundary element discretizations of PDE eigenvalue problems. Copyright © 2013 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

SUMMARY In this paper, a rapid iterative algorithm is proposed to find robust approximations for the inverse of nonsingular matrices. The analysis of convergence reveals that this high-order method possesses eighth-order convergence. The interesting point is that, this rate is attained using less number of matrix-by-matrix multiplications in contrast to the existing methods of the same type in the literature. The extension of the method for finding Moore-Penrose inverse of singular or rectangular matrices is also presented. Numerical comparisons will be given to show the applicability, stability and consistency of the new scheme by paying special attention on the computational time. Copyright © 2013 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

SUMMARY In this paper, some important spectral characterizations of symmetric nonnegative tensors are analyzed. In particular, it is shown that a symmetric nonnegative tensor has the following properties: (i) its spectral radius is zero if and only if it is a zero tensor; (ii) it is weakly irreducible (respectively, irreducible) if and only if it has a unique positive (respectively, nonnegative) eigenvalue-eigenvector; (iii) the minimax theorem is satisfied without requiring the weak irreducibility condition; and (iv) if it is weakly reducible, then it can be decomposed into some weakly irreducible tensors. In addition, the problem of finding the largest eigenvalue of a symmetric nonnegative tensor is shown to be equivalent to finding the global solution of a convex optimization problem. Subsequently, algorithmic aspects for computing the largest eigenvalue of symmetric nonnegative tensors are discussed. Copyright © 2013 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

SUMMARY In this paper, we address the solution of three-dimensional heterogeneous Helmholtz problems discretized with second-order finite difference methods with application to acoustic waveform inversion in geophysics. In this setting, the numerical simulation of wave propagation phenomena requires the approximate solution of possibly very large indefinite linear systems of equations. For that purpose, we propose and analyse an iterative two-grid method acting on the Helmholtz operator where the coarse grid problem is solved inaccurately. A cycle of a multigrid method applied to a complex shifted Laplacian operator is used as a preconditioner for the approximate solution of this coarse problem. A single cycle of the new method is then used as a variable preconditioner of a flexible Krylov subspace method. We analyse the properties of the resulting preconditioned operator by Fourier analysis. Numerical results demonstrate the effectiveness of the algorithm on three-dimensional applications. The proposed numerical method allows us to solve three-dimensional wave propagation problems even at high frequencies on a reasonable number of cores of a distributed memory computer. Copyright © 2012 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

SUMMARY The main purpose of this paper is the extension of the classical spectral direct and inverse analysis of Jacobi matrices for the non-self-adjoint setting. Matrices of this class appear in the context of non-Hermitian quantum mechanics. The reconstruction of a pseudo-Jacobi matrix from its spectrum and the spectra of two complementary principal matrices is investigated in the context of indefinite inner product spaces. An existence and uniqueness theorem is given, and a strikingly simple algorithm, based on the Euclidean division algorithm, to reconstruct the matrix from the spectral data is presented. A result of Friedland and Melkman stating a necessary and sufficient condition for a real sequence to be the spectrum of a non-negative Jacobi matrix is revisited and generalized. Namely, it is shown that a suitable set of prescribed eigenvalues defines a unique non-negative pseudo-Jacobi matrix, which is J-Hermitian for a fixed J. Copyright © 2012 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

SUMMARY In this paper, an inner-outer preconditioner for symmetric positive definite matrices based on hierarchically semi-separable (HSS) matrix representation is presented. A sequence of new HSS algorithms are developed, including some ULV-type HSS methods and triangular HSS methods. During the construction of this preconditioner, off-diagonal blocks are compressed in parallel, and an approximate HSS form is used as the outer HSS representation. In the meantime, the internal dense diagonal blocks are approximately factorized into HSS Cholesky factors. Such inner factorizations are guaranteed to exist for any approximation accuracy. The inner-outer preconditioner combines the advantages of both direct HSS and triangular HSS methods, and are both scalable and robust. Systematic complexity analysis and discussions of the errors are presented. Various tests on some practical numerical problems are used to demonstrate the efficiency and robustness. In particular, the effectiveness of the preconditioner in iterative solutions has been shown for some ill-conditioned problems. This work also gives a practical way of developing inner-outer preconditioners using direct rank structured factorizations (other than iterative methods). Copyright © 2012 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

SUMMARY The choice of the preconditioner is a key factor to accelerate the convergence of eigensolvers for large-size sparse eigenproblems. Although incomplete factorizations with partial fill-in prove generally effective in sequential computations, the efficient preconditioning of parallel eigensolvers is still an open issue. The present paper describes the use of block factorized sparse approximate inverse (BFSAI) preconditioning for the parallel solution of large-size symmetric positive definite eigenproblems with both a simultaneous Rayleigh quotient minimization and the Jacobi-Davidson algorithm. BFSAI coupled with a block diagonal incomplete decomposition proves a robust and efficient parallel preconditioner in a number of test cases arising from the finite element discretization of 3D fluid-dynamical and mechanical engineering applications, outperforming FSAI even by a factor of 8 and exhibiting a satisfactory scalability. Copyright © 2011 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

SUMMARY An iterative method for finding the largest eigenvalue of a nonnegative tensor was proposed by Ng, Qi, and Zhou in 2009. In this paper, we establish an explicit linear convergence rate of the Ng-Qi-Zhou method for essentially positive tensors. Numerical results are given to demonstrate linear convergence of the Ng-Qi-Zhou algorithm for essentially positive tensors. Copyright © 2011 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

SUMMARY In this paper, we consider the stabilization of the nonstationary incompressible Navier-Stokes equations around a stationary solution by a boundary linear feedback control. The feedback operator is obtained from the solution of the algebraic Bernoulli equation associated with the penalized linearized Navier-Stokes equations around an unstable stationary solution and is used to locally stabilize the original nonlinear equations. We give the explicit factorized form of the stabilizing solution of the algebraic Bernoulli equation. The numerical effectiveness of this approach is demonstrated by stabilizing the vortex shedding behind a circular obstacle. Copyright © 2011 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

SUMMARY In this paper, a practical two-term acceleration algorithm is proposed, the interval of the parameter which guarantees the convergence of the acceleration algorithm is analyzed in detail. Further, the acceleration ratio of the new acceleration algorithm is obtained in advance. The new acceleration algorithm is less sensitive to the parameter than the Chebyshev semi-iterative method. Finally, some numerical examples show that the accelerated algorithm is effective. Copyright © 2011 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

SUMMARY The aim of the nuclear norm minimization problem is to find a matrix that minimizes the sum of its singular values and satisfies some constraints simultaneously. Such a problem has received more attention largely because it is closely related to the affine rank minimization problem, which appears in many control applications including controller design, realization theory, and model reduction. In this paper, we first propose an exact version alternating direction method for solving the nuclear norm minimization problem with linear equality constraints. At each iteration, the method involves a singular value thresholding and linear matrix equations which are solved exactly. Convergence of the proposed algorithm is followed directly. To broaden the capacity of solving larger problems, we solve approximately the subproblem by an iterative method with the Barzilai-Borwein steplength. Some extensions to the noisy problems and nuclear norm regularized least-square problems are also discussed. Numerical experiments and comparisons with the state-of-the-art method FPCA show that the proposed method is effective and promising. Copyright © 2011 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]

SUMMARY Although there have been significant advances in robust algebraic multigrid (AMG) methods in recent years, numerical studies and emerging hardware architectures continue to favor structured-grid approaches. Specifically, implementations of logically structured robust variational multigrid algorithms, such as the black box multigrid (BoxMG) solver, have been shown to be 10 times faster than AMG for three-dimensional heterogeneous diffusion problems on structured grids. BoxMG offers important features such as operator-induced interpolation for robustness, while taking advantage of direct data access and bounded complexity in the Galerkin coarse-grid operator. Moreover, because BoxMG uses a variational framework, it can be used to explore advances of modern adaptive AMG approaches in a structured setting. In this paper, we show how to extend the adaptive multigrid methodology to the BoxMG setting. This extension not only retains the favorable properties of the adaptive framework but also sheds light on the relationship between BoxMG and AMG. In particular, we show how classical BoxMG can be viewed as a special case of classical AMG and how this viewpoint leads to a richer family of adaptive BoxMG approaches. We present numerical results that explore this family of adaptive methods and compare its robustness and efficiency to the classical BoxMG solver.Copyright © 2012 John Wiley & Sons, Ltd. [ABSTRACT FROM AUTHOR]