This study considers quadrature-based algorithms to compute A^αb, the action of a real power of a Hermitian positive-definite matrix A on a vector b. In these algorithms, the computation of an integral representation of A^α b is reduced to solving several tens or hundreds of shifted linear systems. Current approaches usually analyze the quadrature discretization error, but rarely take into account the additional error introduced by solving these shifted linear systems with iterative solvers. Here, we bound this error with the residual of the approximated solution of these linear systems. This allows the derivation of a stopping criterion for iterative solvers to keep the error of A^αb below a prescribed error tolerance. Numerical results demonstrate that the proposed criterion enables the computation of A^αb within prescribed tolerance limits.
Krylov subspace methods, such as the Conjugate Gradient (CG) and BiCGSTAB methods, are widely used in scientific computing for solving linear systems. In this study, we propose a new framework for solving large Sylvester equations in a low-rank format by reconstructing matrix-oriented Krylov subspace methods. The framework realizes efficient algorithms that are mathematically equivalent to the matrix-oriented Krylov subspace methods by exploiting the mathematical properties of the Sylvester operator and the low-rank structure of the right-hand side. Specifically, by leveraging these properties, approximate solutions can be expressed in a low-rank factorized form, enabling efficient computation and reduced memory requirements. The effectiveness of our algorithms is demonstrated through numerical experiments.
The generalized product-type method based on the biconjugate gradient (hereafter referred to as the GPBiCG method) has been recognized as an efficient Krylov subspace method for solving nonsymmetric linear systems. In this paper, tensor form of the GPBiCG method and its preconditioned variant are presented for solving the Stein tensor equations. The effectiveness of the proposed methods is verified by numerical experiments.
This paper focuses on integration-type methods that employ the double-exponential (DE) integral formula, DE-type methods, for computing matrix functions. DE-type methods exhibit high convergence rate thanks to the highly efficient DE integral formula; however, it is known that existence of eigenvalues near a specific region on the complex plane adversely affects the convergence. This paper proposes a deflation approach that improves the convergence of the DE-type methods by extracting and separating the eigenspaces corresponding to the eigenvalues within such a region. To extract these interior eigenspaces, the proposed approach uses a complex moment-based eigensolver. The components of the matrix function corresponding to the interior eigenspace are directly obtained from the eigenpairs, whereas those corresponding to the exterior eigenspace are computed through the DE-type method with deflation. Error analysis and numerical experiments demonstrate the effectiveness of the proposed approach.
Several methods for computing the action of the matrix exponential e^Ab are expressed by substituting A into a rational approximation of the scalar exponential function. The error of such methods can be estimated using the numerical range of A, which enables the computation of e^Ab with a prescribed accuracy. However, when the input matrix has the structure A = τM^-1K, this approach is challenging because computing the bounding box of numerical range is difficult and the numerical range may be too large to construct rational approximations on it. In this paper, focusing on the case where M is a well-conditioned symmetric positive definite matrix, we propose considering the numerical range of a similarity transformed matrix of A. The numerical range of transformed matrix is not only numerically computable but can also be theoretically bounded depending on properties of K. Numerical experiments confirm that the computations can be performed within the prescribed error tolerance.
The Stein tensor equation with the Einstein product arises in multidimensional signal processing, control theory, and data reconstruction. However, its inherent nonlinearity and high dimensionality present computational challenges. In this work, we design a predefined-time zeroing neural network model equipped with a novel activation function to efficiently solve the time-varying Stein tensor equation. Theoretical analysis establishes the model’s convergence and robustness, demonstrating that it not only converges to the exact solution within a predefined time but also exhibits strong resilience to two types of noise. Numerical experiments further confirm the superior performance of the proposed model in terms of accuracy, convergence speed, and noise robustness.
Variational quantum algorithms (VQAs) have attracted attention as quantum algorithms employing on noisy intermediate-scale quantum devices. VQAs for solving the Poisson equation were recently proposed, as this method transforms the linear equation with a symmetric positive definite matrix obtained by discretizing the Poisson equation into a minimization problem. In this study, we propose a VQA for second-order linear differential equations including the Poisson equation on the basis of the referenced study, where the coefficient matrix obtained through discretization is a non-symmetric matrix. Furthermore, in our VQAs, we succeeded in decomposing the coefficient matrices arising from the periodic and the Dirichlet boundary condition into linear combinations of measurement operators and unitary matrices.
This letter considers the double exponential (DE) formula for computing the matrix function -A log(A), where A is a Hermitian positive semidefinite matrix and tr(A) = 1. The motivation of this work is to utilize the DE formula without selecting parameters. In order to accomplish this, we present a method for truncating the infinite interval transformed by the DE transformation based on an error analysis to achieve the required accuracy. We also discuss techniques to select the number of abscissas. In addition, we report that an appropriate choice of a parameter in the DE transformation will improve the convergence.
Toward the end of the 20th century, S.-L. Zhang constructed the so-called Zhang’s framework that successfully incorporated the CGS and Bi-CGSTAB methods and subsequently introduced the GPBi-CG method derived from the framework. While the GPBi-CG method often converges faster than the CGS and Bi-CGSTAB methods, there are still cases where the CGS method performs the best. This observation has motivated us to revisit the framework to find a hybrid algorithm, combining the GPBi-CG and CGS methods. In this paper, we introduce the hybrid algorithm, show its efficiency in some numerical experiments, and discuss a possible reason behind its effectiveness.
It is well-known that a multilinear system with a nonsingular M-tensor and a positive right-hand side has a unique positive solution. Tensor splitting methods generalizing the classical iterative methods for linear systems have been proposed for finding the unique positive solution. The Alternating Anderson-Richardson (AAR) method is an effective method to accelerate the classical iterative methods. In this study, we apply the idea of AAR for finding the unique positive solution quickly. We first present a tensor Richardson method based on tensor regular splittings, then apply Anderson acceleration to the tensor Richardson method and derive a tensor Anderson-Richardson method, finally, we periodically employ the tensor Anderson-Richardson method within the tensor Richardson method and propose a tensor AAR method. Numerical experiments show that the proposed method is effective in accelerating tensor splitting methods.
The shifted LOPBiCG method can circumvent a problem of the optimal choice for the initial seed system through a seed-switching technique when solving shifted linear systems, thanks to its natural collinearity of residual vectors. However, LOPBiCG, employing the same one-degree accelerating polynomial as that in the Bi-CGSTAB method, may not perform well for linear systems whose real coefficient matrix has complex eigenvalues with relatively large imaginary parts. Here, we utilize an & lscr;-degree accel-erating polynomial to extend LOPBiCG to LOPBiCG(& lscr;) and apply LOPBiCG(& lscr;) to solve shifted linear systems, which is referredto asshifted LOPBiCG(& lscr;). Since the natural collinearity for LOPBiCG(& lscr;) is also satisfied, we can incorporate the seed-switchingtechnique into shifted LOPBiCG(& lscr;). Numerical experiments demonstrate that our method effectively solves linear systems witha complex spectrum and is free from the issue associated with choosing the initial seed when solving shifted linear systems
This note considers the computation of the logarithm of symmetric positive definite matrices using the Gauss-Legendre (GL) quadrature. The GL quadrature becomes slow when the condition number of the given matrix is large. In this note, we propose a technique dividing the matrix logarithm into two matrix logarithms, where the condition numbers of the divided logarithm arguments are smaller than that of the original matrix. Although the matrix logarithm needs to be computed twice, each computation can be performed more efficiently, and it potentially reduces the overall computational cost. It is shown that the proposed technique is effective when the condition number of the given matrix is approximately between 130 and 3.0 x 105.
This article considers the computation of the matrix exponential e A {{\rm{e}}}<^>{A} with numerical quadrature. Although several quadrature-based algorithms have been proposed, they focus on (near) Hermitian matrices. In order to deal with non-Hermitian matrices, we use another integral representation including an oscillatory term and consider applying the double exponential (DE) formula specialized to Fourier integrals. The DE formula transforms the given integral into another integral whose interval is infinite, and therefore, it is necessary to truncate the infinite interval. In this article, to utilize the DE formula, we analyze the truncation error and propose two algorithms. The first one approximates e A {{\rm{e}}}<^>{A} with the fixed mesh size, which is a parameter in the DE formula affecting the accuracy. The second one computes e A {{\rm{e}}}<^>{A} based on the first one with automatic selection of the mesh size depending on the given error tolerance.
When solving shifted linear systems using shifted Krylov subspace methods, selecting a seed system is necessary, and an unsuitable seed may result in many shifted systems being unsolved. To avoid this problem, a seed-switching technique has been proposed to help switch the seed system to another linear system as a new seed system without losing the dimension of the constructed Krylov subspace. Nevertheless, this technique requires collinear residual vectors when applying Krylov subspace methods to the seed and shifted systems. Since the product-type shifted Krylov subspace methods cannot provide such collinearity, these methods cannot use this technique. In this article, we propose a variant of the shifted BiCGstab method, which possesses the collinearity of residuals, and apply the seed-switching technique to it. Some numerical experiments show that the problem of choosing the initial seed system is circumvented.
We consider the convolution equation F*X=BF* X=B, where F∈R3×3F\in {{\mathbb{R}}}^{3\times 3} and B∈Rm×nB\in {{\mathbb{R}}}^{m\times n} are given and X∈Rm×nX\in {{\mathbb{R}}}^{m\times n} is to be determined. The convolution equation can be regarded as a linear system with a coefficient matrix of special structure. This fact has led to many studies including efficient numerical algorithms for solving the convolution equation. In this study, we show that the convolution equation can be represented as a generalized Sylvester equation. Furthermore, for some realistic examples arising from image processing, we show that the generalized Sylvester equation can be reduced to a simpler form, and we analyze the unique solvability of the convolution equation.
The optimal value of the projected successive overrelaxation (PSOR) method for nonnegative quadratic programming problems is problem-dependent. We present a novel adaptive PSOR algorithm that adaptively controls the relaxation parameter using the Wolfe conditions. The method and its variants can be applied to various problems without requiring a specific assumption regarding the matrix defining the objective function, and the cost for updating the parameter is negligible in the whole iteration. Numerical experiments show that the proposed methods often perform comparably to (or sometimes superior to) the PSOR method with a nearly optimal relaxation parameter.
For fourth-order geometric evolution equations for planar curves with the dissipation of the bending energy,including the Willmore and the Helfrich flows,we consider a numerical approach.In this study,we construct a structure-preserving method based on a discrete variational derivative me-thod.Furthermore,to prevent the vertex concentration that may lead to numer-ical instability,we discretely introduce Deckelnick's tangential velocity.Here,a modification term is introduced in the process of adding tangential velocity.This modified term enables the method to reproduce the equations'properties while preventing vertex concentration.Numerical experiments demonstrate that the proposed approach captures the equations'properties with high accu-racy and avoids the concentration of vertices.
The tensor biconjugate gradient (TBiCG) method has recently been proposed for solving Sylvester tensor equations. The TBiCG method is based on the BiCG method that may exhibit irregular convergence behavior. To overcome the limitations, product-type methods, such as BiCGSTAB and GPBiCG, have been proposed. In this study, we apply the idea of product-type methods to solve Sylvester tensor equations and propose tensor GPBiCG and BiCGSTAB methods. Furthermore, we consider preconditioned algorithms of the tensor GPBiCG and BiCGSTAB methods using the nearest Kronecker product preconditioner. Numerical experiments illustrate that the proposed methods are competitive with some existing methods. & COPY; 2023 Elsevier Inc. All rights reserved.