67 research outputs found
Rational Krylov and ADI iteration for infinite size quasi-Toeplitz matrix equations
No full textWe consider a class of linear matrix equations involving semi-infinite matrices which have a quasi-Toeplitz structure. These equations arise in different settings, mostly connected with PDEs or the study of Markov chains such as random walks on bidimensional lattices. We present the theory justifying the existence of the solution in an appropriate Banach algebra which is computationally treatable, and we propose several methods for computing them. We show how to adapt the ADI iteration to this particular infinite dimensional setting, and how to construct rational Krylov methods. Convergence theory is discussed, and numerical experiments validate the proposed approaches
Rank-structured qr for chebyshev rootfinding
No full textWe consider the computation of roots of polynomials expressed in the Chebyshev basis. We extend the QR iteration presented in [Y. Eidelman, L. Gemignani, and I. Gohberg, Numer. Algorithms, 47 (2008), pp. 253-273] introducing an aggressive early deflation strategy, and showing that the rank-structure allows one to parallelize the algorithm, avoiding data dependencies which would be present in the unstructured QR. We exploit the particular structure of the colleague linearization to achieve quadratic complexity and linear storage requirements. The (unbalanced) QR iteration used for Chebyshev rootfinding does not guarantee backward stability on the polynomial coefficients, unless the vector of coefficients satisfy ||p|| ≈ 1, a hypothesis which is almost never verified for polynomials approximating smooth functions. Even though the presented method is mathematically equivalent to the unbalanced QR algorithm, we show that exploiting the rank structure allows one to guarantee a small backward error on the polynomial, up to an explicitly computable amplification factor γ1(p), which depends on the polynomial under consideration. We show that this parameter is almost always of moderate size, making the method accurate on several numerical tests, in contrast with what happens in the unstructured unbalanced QR
Computing performability measures in Markov chains by means of matrix functions
No full textWe discuss the efficient computation of performance, reliability, and availability measures for Markov chains; these metrics — and the ones obtained by combining them, are often called performability measures.
We show that this computational problem can be recasted as the evaluation of a bilinear form induced by appropriate matrix functions, and thus solved by leveraging the fast methods available for this task.
We provide a comprehensive analysis of the theory required to translate the problem from the language of Markov chains to the one of matrix functions. The advantages of this new formulation are discussed, and it is shown that this setting allows to easily study the sensitivities of the measures with respect to the model parameters.
Numerical experiments confirm the effectiveness of our approach; the tests we have run show that we can outperform the solvers available in state of the art commercial packages on a representative set of large scale examples
A Multilinear Nyström Algorithm for Low-Rank Approximation of Tensors in Tucker Format
No full textThe Nyström method offers an effective way to obtain low-rank approximation of SPD matrices and has been recently extended and analyzed to nonsymmetric matrices (leading to the generalized Nyström method). It is a randomized, single-pass, streamable, cost-effective, and accurate alternative to the randomized SVD, and it facilitates the computation of several matrix low-rank factorizations. In this paper, we take these advancements a step further by introducing a higher-order variant of Nyström's methodology tailored to approximating low-rank tensors in the Tucker format: the multilinear Nyström technique. We show that, by introducing appropriate small modifications in the formulation of the higher-order method, strong stability properties can be obtained. This algorithm retains the key attributes of the generalized Nyström method, positioning it as a viable substitute for the randomized higher-order SVD algorithm
Mixed Precision Recursive Block Diagonalization for Bivariate Functions of Matrices
No full textVarious numerical linear algebra problems can be formulated as evaluating bivariate function of matrices. The most notable examples are the Fr'echet derivative along a direction, the evaluation of (univariate) functions of Kronecker-sum-structured matrices, and the solution of, Sylvester matrix equations. In this work, we propose a recursive block diagonalization algorithm for computing bivariate functions of matrices of small to medium size, for which dense linear algebra is appropriate. The algorithm combines a blocking strategy, as in the Schur-Parlett scheme, and an evaluation procedure for the diagonal blocks. We discuss two implementations of the latter. The first is a natural choice based on Taylor expansions, whereas the second is derivative-free and relies on a multiprecision perturb-and-diagonalize approach. In particular, the appropriate use of multiprecision guarantees backward stability without affecting the efficiency in the generic case. This makes the second approach more robust. The whole method has cubic complexity, and it is closely related to the well-known Bartels-Stewart algorithm for Sylvester matrix equations when applied to f(x, y) = x+1y . We validate the performances of the proposed numerical method on several problems with different conditioning properties
Rational Krylov for Stieltjes matrix functions: convergence and pole selection
No full textEvaluating the action of a matrix function on a vector, that is x= f(M) v, is an ubiquitous task in applications. When M is large, one usually relies on Krylov projection methods. In this paper, we provide effective choices for the poles of the rational Krylov method for approximating x when f(z) is either Cauchy–Stieltjes or Laplace–Stieltjes (or, which is equivalent, completely monotonic) and M is a positive definite matrix. Relying on the same tools used to analyze the generic situation, we then focus on the case M= I⊗ A- BT⊗ I, and v obtained vectorizing a low-rank matrix; this finds application, for instance, in solving fractional diffusion equation on two-dimensional tensor grids. We see how to leverage tensorized Krylov subspaces to exploit the Kronecker structure and we introduce an error analysis for the numerical approximation of x. Pole selection strategies with explicit convergence bounds are given also in this case
Sampling the eigenvalues of random orthogonal and unitary matrices
No full textWe develop an efficient algorithm for sampling the eigenvalues of random matrices distributed according to the Haar measure over the orthogonal or unitary group. Our technique samples directly a factorization of the Hessenberg form of such matrices, and then computes their eigenvalues with a tailored core-chasing algorithm. This approach requires a number of floating-point operations that is quadratic in the order of the matrix being sampled, and can be adapted to other matrix groups. In particular, we explain how it can be used to sample the Haar measure over the special orthogonal and unitary groups and the conditional probability distribution obtained by requiring the determinant of the sampled matrix be a given complex number on the complex unit circle
Low-rank tensor structure preservation in fractional operators by means of exponential sums
No full textThe use of fractional differential equations is a key tool in modeling non-local phenomena. Often, an efficient scheme for solving a linear system involving the discretization of a fractional operator is computing inverse fractional powers of the standard discretized Laplace operator. In this work, an exponential sum approximation for such fractional powers is derived. It is accurate over all positive real numbers larger than one, and allows to efficiently approximate the action of such operators on tensors stored in a variety of low-rank formats (CP, TT, Tucker). The results are relevant from a practical and theoretical perspective, as they predict the low-rank approximability of the solutions of these linear systems in low-rank tensor formats
Fast Solvers for Two-Dimensional Fractional Diffusion Equations Using Rank Structured Matrices
No full textWe consider the discretization of time-space diffusion equations with fractional derivatives in space and either one-dimensional (1D) or 2D spatial domains. The use of an implicit Euler scheme in time and finite differences or finite elements in space leads to a sequence of dense large scale linear systems describing the behavior of the solution over a time interval. We prove that the coefficient matrices arising in the 1D context are rank structured and can be efficiently represented using hierarchical formats (scrH -matrices, HODLR). Quantitative estimates for the rank of the off-diagonal blocks of these matrices are presented. We analyze the use of HODLR arithmetic for solving the 1D case and we compare this strategy with existing methods that exploit the Toeplitz-like structure to precondition the GMRES iteration. The numerical tests demonstrate the convenience of the HODLR format when at least a reasonably low number of time steps is needed. Finally, we explain how these properties can be leveraged to design fast solvers for problems with 2D spatial domains that can be reformulated as matrix equations. The experiments show that the approach based on the use of rank-structured arithmetic is particularly effective and outperforms current state of the art techniques
RANDOMIZED SKETCHED TT-GMRES FOR LINEAR SYSTEMS WITH TENSOR STRUCTURE*
No full textIn the past decade, tensors have shown their potential as valuable tools for various
tasks in numerical linear algebra. While most of the research has been focusing on how to compress a
given tensor in order to maintain information as well as reducing the storage demand for its allocation,
the solution of linear tensor equations is a less explored venue. Even if many of the routines available
in the literature are based on alternating minimization schemes (ALS), we pursue a different path
and utilize Krylov methods instead. The use of Krylov methods in the tensor realm is not new.
However, these routines often turn out to be rather expensive in terms of computational cost, and ALS
procedures are preferred in practice. We enhance Krylov methods for linear tensor equations with
a panel of diverse randomization-based strategies which remarkably increase the efficiency of these
solvers, making them competitive with state-of-the-art ALS schemes. The up-to-date randomized
approaches we employ range from sketched Krylov methods with incomplete orthogonalization and
structured sketching transformations to streaming algorithms for tensor rounding. The promising
performance of our new solver for linear tensor equations is demonstrated by many numerical results
- …
