Search

We consider the convolution equation $F*X=B$, where $F\in\mathbb{R}^{3\times 3}$ and $B\in\mathbb{R}^{m\times n}$ are given, and $X\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 analyze the unique solvability of the convolution equation.

This paper considers the computation of the matrix exponential $\mathrm{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 paper, to utilize the DE formula, we analyze the truncation error and propose two algorithms. The first one approximates $\mathrm{e}^A$ with the fixed mesh size which is a parameter in the DE formula affecting the accuracy. Second one computes $\mathrm{e}^A$ based on the first one with automatic selection of the mesh size depending on the given error tolerance.

The choice of relaxation parameter in the projected successive overrelaxation (PSOR) method for nonnegative quadratic programming problems is problem-dependent. We present novel adaptive PSOR algorithms that adaptively control the relaxation parameter using the Wolfe conditions. The method and its variants can be applied to various problems without requiring additional assumptions, barring the positive semidefiniteness concerning the matrix that defines 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 matrix $A$, vector $\boldsymbol{b}$ and function $f$, the computation of vector $f(A)\boldsymbol{b}$ arises in many scientific computing applications. We consider the problem of obtaining quantum state $\lvert f \rangle$ corresponding to vector $f(A)\boldsymbol{b}$. There is a quantum algorithm to compute state $\lvert f \rangle$ using eigenvalue estimation that uses phase estimation and Hamiltonian simulation $\mathrm{e}^{\mathrm{{\bf i}} A t}$. However, the algorithm based on eigenvalue estimation needs $\textrm{poly}(1/\epsilon)$ runtime, where $\epsilon$ is the desired accuracy of the output state. Moreover, if matrix $A$ is not Hermitian, $\mathrm{e}^{\mathrm{{\bf i}} A t}$ is not unitary and we cannot run eigenvalue estimation. In this paper, we propose a quantum algorithm that uses Cauchy's integral formula and the trapezoidal rule as an approach that avoids eigenvalue estimation. We show that the runtime of the algorithm is $\mathrm{poly}(\log(1/\epsilon))$ and the algorithm outputs state $\lvert f \rangle$ even if $A$ is not Hermitian., Comment: 23 pages, 1 figure

The matrix functions can be defined by Cauchy's integral formula and can be approximated by the linear combination of inverses of shifted matrices using a quadrature formula. In this paper, we show a concrete construction of a framework to implement the linear combination of the inverses on quantum computers and propose a quantum algorithm for matrix functions based on the framework. Compared with the previous study [S. Takahira, A. Ohashi, T. Sogabe, and T.S. Usuda, Quant. Inf. Comput., 20, 1&2, 14--36, (Feb. 2020)] that proposed a quantum algorithm to compute a quantum state for the matrix function based on the circular contour centered at the origin, the quantum algorithm in the present paper can be applied to a more general contour. Moreover, the algorithm is described by the block-encoding framework. Similarly to the previous study, the algorithm can be applied even if the input matrix is not a Hermitian or normal matrix., Comment: 23 pages, 5 figures

Two quadrature-based algorithms for computing the matrix fractional power $A^\alpha$ are presented in this paper. These algorithms are based on the double exponential (DE) formula, which is well-known for its effectiveness in computing improper integrals as well as in treating nearly arbitrary endpoint singularities. The DE formula transforms a given integral into another integral that is suited for the trapezoidal rule; in this process, the integral interval is transformed to the infinite interval. Therefore, it is necessary to truncate the infinite interval into an appropriate finite interval. In this paper, a truncation method, which is based on a truncation error analysis specialized to the computation of $A^\alpha$, is proposed. Then, two algorithms are presented -- one computes $A^\alpha$ with a fixed number of abscissas, and the other computes $A^\alpha$ adaptively. Subsequently, the convergence rate of the DE formula for Hermitian positive definite matrices is analyzed. The convergence rate analysis shows that the DE formula converges faster than the Gaussian quadrature when $A$ is ill-conditioned and $\alpha$ is a non-unit fraction. Numerical results show that our algorithms achieved the required accuracy and were faster than other algorithms in several situations., Comment: The title of the manuscript was changed. The former title is "Computing the matrix fractional power based on the double exponential formula"

We develop K$\omega$, an open-source linear algebra library for the shifted Krylov subspace methods. The methods solve a set of shifted linear equations $(z_k I-H)x^{(k)}=b\, (k=0,1,2,...)$ for a given matrix $H$ and a vector $b$, simultaneously. The leading order of the operational cost is the same as that for a single equation. The shift invariance of the Krylov subspace is the mathematical foundation of the shifted Krylov subspace methods. Applications in materials science are presented to demonstrate the advantages of the algorithm over the standard Krylov subspace methods such as the Lanczos method. We introduce benchmark calculations of (i) an excited (optical) spectrum and (ii) intermediate eigenvalues by the contour integral on the complex plane. In combination with the quantum lattice solver $\mathcal{H} \Phi$, K$\omega$ can realize parallel computation of excitation spectra and intermediate eigenvalues for various quantum lattice models., Comment: 12 pages, 5 figures

A thick-restart Lanczos type algorithm is proposed for Hermitian $J$-symmetric matrices. Since Hermitian $J$-symmetric matrices possess doubly degenerate spectra or doubly multiple eigenvalues with a simple relation between the degenerate eigenvectors, we can improve the convergence of the Lanczos algorithm by restricting the search space of the Krylov subspace to that spanned by one of each pair of the degenerate eigenvector pairs. We show that the Lanczos iteration is compatible with the $J$-symmetry, so that the subspace can be split into two subspaces that are orthogonal to each other. The proposed algorithm searches for eigenvectors in one of the two subspaces without the multiplicity. The other eigenvectors paired to them can be easily reconstructed with the simple relation from the $J$-symmetry. We test our algorithm on randomly generated small dense matrices and a sparse large matrix originating from a quantum field theory., Comment: 22 pages, 13 figures, svjour3.cls is used. The ambiguity of the term "J-symmetric" was removed. Another test case was added. Japan J. Indust. Appl. Math. (2020)

We consider applying the Strang splitting to semilinear parabolic problems. The key ingredients of the Strang splitting are the decomposition of the equation into several parts and the computation of approximate solutions by combining the time evolution of each split equation. However, when the Dirichlet boundary condition is imposed, order reduction could occur due to the incompatibility of the split equations with the boundary condition. In this paper, to overcome the order reduction, a modified Strang splitting procedure is presented for the one-dimensional semilinear parabolic equation with first-order spatial derivatives, like the Burgers equation.

The T-congruence Sylvester equation is the matrix equation $AX+X^{\mathrm{T}}B=C$, where $A\in\mathbb{R}^{m\times n}$, $B\in\mathbb{R}^{n\times m}$, and $C\in\mathbb{R}^{m\times m}$ are given, and $X\in\mathbb{R}^{n\times m}$ is to be determined. Recently, Oozawa et al. discovered a transformation that the matrix equation is equivalent to one of the well-studied matrix equations (the Lyapunov equation); however, the condition of the transformation seems to be too limited because matrices $A$ and $B$ are assumed to be square matrices ($m=n$). In this paper, two transformations are provided for rectangular matrices $A$ and $B$. One of them is an extension of the result of Oozawa et al. for the case $m\ge n$, and the other is a novel transformation for the case $m\le n$.

We propose a Fourier pseudo-spectral scheme for the space-fractional nonlinear Schr\"odinger equation. The proposed scheme has the following features: it is linearly implicit, it preserves two invariants of the equation, its unique solvability is guaranteed without any restrictions on space and time step sizes. The scheme requires solving a complex symmetric linear system per time step. To solve the system efficiently, we also present a certain variable transformation and preconditioner.

We consider the computation of the matrix logarithm by using numerical quadrature. The efficiency of numerical quadrature depends on the integrand and the choice of quadrature formula. The Gauss--Legendre quadrature has been conventionally employed; however, the convergence could be slow for ill-conditioned matrices. This effect may stem from the rapid change of the integrand values. To avoid such situations, we focus on the double exponential formula, which has been developed to address integrands with endpoint singularity. In order to utilize the double exponential formula, we must determine a suitable finite integration interval, which provides the required accuracy and efficiency. In this paper, we present a method for selecting a suitable finite interval based on an error analysis as well as two algorithms, and one of these algorithms addresses error control.

Because the expense of estimating the optimal value of the relaxation parameter in the successive over-relaxation (SOR) method is usually prohibitive, the parameter is often adaptively controlled. In this paper, new adaptive SOR methods are presented that are applicable to a variety of symmetric positive definite linear systems and do not require additional matrix-vector products when updating the parameter. To this end, we regard the SOR method as an algorithm for minimising a certain objective function, which yields an interpretation of the relaxation parameter as the step size following a certain change of variables. This interpretation enables us to adaptively control the step size based on some line search techniques, such as the Wolfe conditions. Numerical examples demonstrate the favourable behaviour of the proposed methods.

The iterative nature of many discretisation methods for continuous dynamical systems has led to the study of the connections between iterative numerical methods in numerical linear algebra and continuous dynamical systems. Certain researchers have used the explicit Euler method to understand this connection, but this method has its limitation. In this study, we present a new connection between successive over-relaxation (SOR)-type methods and gradient systems; this connection is based on discrete gradient methods. The focus of the discussion is the equivalence between SOR-type methods and discrete gradient methods applied to gradient systems. The discussion leads to new interpretations for SOR-type methods. For example, we found a new way to derive these methods; these methods monotonically decrease a certain quadratic function and obtain a new interpretation of the relaxation parameter. We also obtained a new discrete gradient while studying the new connection.

We consider computing the $k$-th eigenvalue and its corresponding eigenvector of a generalized Hermitian eigenvalue problem of $n\times n$ large sparse matrices. In electronic structure calculations, several properties of materials, such as those of optoelectronic device materials, are governed by the eigenpair with a material-specific index $k.$ We present a three-stage algorithm for computing the $k$-th eigenpair with validation of its index. In the first stage of the algorithm, we propose an efficient way of finding an interval containing the $k$-th eigenvalue $(1 \ll k \ll n)$ with a non-standard application of the Lanczos method. In the second stage, spectral bisection for large-scale problems is realized using a sparse direct linear solver to narrow down the interval of the $k$-th eigenvalue. In the third stage, we switch to a modified shift-and-invert Lanczos method to reduce bisection iterations and compute the $k$-th eigenpair with validation. Numerical results with problem sizes up to 1.5 million are reported, and the results demonstrate the accuracy and efficiency of the three-stage algorithm.

Incremental Newton (IN) iteration, proposed by Iannazzo, is stable for computing the matrix $p$th root, and its computational cost is $\mathcal{O}(n^3p)$ flops per iteration. In this paper, a cost-efficient variant of IN iteration is presented. The computational cost of the variant well agrees with $\mathcal{O} (n^3 \log p)$ flops per iteration, if $p$ is up to at least 100.

We consider the numerical integration of the Hunter--Saxton equation, which models the propagation of weakly nonlinear orientation waves. For the equation, we present two weak forms and their Galerkin discretizations. The Galerkin schemes preserve the Hamiltonian of the equation and can be implemented with cheap $H^1$ elements. Numerical experiments confirm the effectiveness of the schemes.

This paper is dedicated to compute Pfaffian and determinant of one type of skew centrosymmetric matrices in terms of general number sequence of second order.

We consider the T-congruence Sylvester equation $AX+X^{\rm T}B=C$, where $A\in \mathbb R^{m\times n}$, $B\in \mathbb R^{n\times m}$ and $C\in \mathbb R^{m\times m}$ are given, and matrix $X \in \mathbb R^{n\times m}$ is to be determined. The T-congruence Sylvester equation has recently attracted attention because of a relationship with palindromic eigenvalue problems. For example, necessary and sufficient conditions for the existence and uniqueness of solutions, and numerical solvers have been intensively studied. In this note, we will show that, under a certain condition, the T-congruence Sylvester equation can be transformed into the Lyapunov equation. This may lead to further properties and efficient numerical solvers by utilizing a great deal of studies on the Lyapunov equation., Comment: 7 pages

The present paper gives a review of our recent progress and latest results for novel linear-algebraic algorithms and its application to large-scale quantum material simulations or electronic structure calculations. The algorithms are Krylov-subspace (iterative) solvers for generalized shifted linear equations, in the form of (zS-H)x=b,in stead of conventional generalized eigen-value equation. The method was implemented in our order-$N$ calculation code ELSES (http://www.elses.jp/) with modelled systems based on ab initio calculations. The code realized one-hundred-million-atom, or 100-nm-scale, quantum material simulations on the K computer in a high parallel efficiency with up to all the built-in processor cores. The present paper also explains several methodological aspects, such as use of XML files and 'novice' mode for general users. A sparse matrix data library in our real problems (http://www.elses.jp/matrix/) was prepared. Internal eigen-value problem is discussed as a general need from the quantum material simulation. The present study is a interdisciplinary one and is sometimes called 'Application-Algorithm-Architecture co-design'. The co-design will play a crucial role in exa-scale scientific computations., Comment: 13 pages, 6 figures

We explains the shifted COCG method which can solve a series of the linear equations generated by numbers of scaler shifts, without time consuming matrix-vector operations, except at the only one reference energy. This is a family of the CG method and sharing the robustness and the capability of the accuracy estimation. Then shifted COCG is quite useful to calculate the Green's function of the many-electron Hamiltonian which have very large dimension. We applied it to the double orbital extended Hubbard model with twelve electrons on the periodic sqrt(8) x sqrt(8) site system, the dimension of the Hamiltonian equals to 64,128,064, and found the ground state is insulator. We also explained the crucial points of the shifted COCG algorithm for reducing the amount of required memory., Comment: 9 pages, 4 figures