A reduced conjugate gradient basis method for fractional diffusion

This work is on a fast and accurate reduced basis method for solving discretized fractional elliptic partial differential equations (PDEs) of the form $\mathcal{A}^su=f$ by rational approximation. A direct computation of the action of such an approximation would require solving multiple (20$\sim$30) large-scale sparse linear systems. Our method constructs the reduced basis using the first few directions obtained from the preconditioned conjugate gradient method applied to one of the linear systems. As shown in the theory and experiments, only a small number of directions (5$\sim$10) are needed to approximately solve all large-scale systems on the reduced basis subspace. This reduces the computational cost dramatically because: (1) We only use one of the large-scale problems to construct the basis; and (2) all large-scale problems restricted to the subspace have much smaller sizes. We test our algorithms for fractional PDEs on a 3d Euclidean domain, a 2d surface, and random combinatorial graphs. We also use a novel approach to construct the rational approximation for the fractional power function by the orthogonal greedy algorithm (OGA)

Solving large linear algebraic systems in the context of integrable non-abelian Laurent ODEs

Abstract The paper reports on a computer algebra program LSSS (Linear Selective Systems Solver) for solving linear algebraic systems with rational coefficients. The program is especially efficient for very large sparse systems that have a solution in which many variables take the value zero. The program is applied to the symmetry investigation of a non-abelian Laurent ODE introduced recently by M. Kontsevich. The computed symmetries confirmed that a Lax pair found for this system earlier generates all first integrals of degree at least up to 14

Krylov methods for large-scale modern problems in numerical linear algebra

Large-scale problems have attracted much attention in the last decades since they arise from different applications in several fields. Moreover, the matrices that are involved in those problems are often sparse, this is, the majority of their entries are zero. Around 40 years ago, the most common problems related to large-scale and sparse matrices consisted in solving linear systems, finding eigenvalues and/or eigenvectors, solving least square problems or computing singular value decompositions. However, in the last years, large-scale and sparse problems of different natures have appeared, motivating and challenging numerical linear algebra to develop effective and efficient algorithms to solve them. Common difficulties that appear during the development of algorithms for solving modern large-scale problems are related to computational costs, storage issues and CPU time, given the large size of the matrices, which indicate that direct methods can not be used. This suggests that projection methods based on Krylov subspaces are a good option to develop procedures for solving large-scale and sparse modern problems. In this PhD Thesis we develop novel and original algorithms for solving two large-scale modern problems in numerical linear algebra: first, we introduce the R-CORK method for solving rational eigenvalue problems and, second, we present projection methods to compute the solution of T-Sylvester matrix equations, both based on Krylov subspaces. The R-CORK method is an extension of the compact rational Krylov method (CORK) [104] introduced to solve a family of nonlinear eigenvalue problems that can be expressed and linearized in certain particular ways and which include arbitrary polynomial eigenvalue problems, but not arbitrary rational eigenvalue problems. The R-CORK method exploits the structure of the linearized problem by representing the Krylov vectors in a compact form in order to reduce the cost of storage, resulting in a method with two levels of orthogonalization. The first level of orthogonalization works with vectors of the same size as the original problem, and the second level works with vectors of size much smaller than the original problem. Since vectors of the size of the linearization are never stored or orthogonalized, R-CORK is more efficient from the point of view of memory and orthogonalization costs than the classical rational Krylov method applied to the linearization. Moreover, since the R-CORK method is based on a classical rational Krylov method, the implementation of implicit restarting is possible and we present an efficient way to do it, that preserves the compact representation of the Krylov vectors. We also introduce in this dissertation projection methods for solving the TSylvester equation, which has recently attracted considerable attention as a consequence of its close relation to palindromic eigenvalue problems and other applications. The theory concerning T-Sylvester equations is rather well understood, and before the work in this thesis, there were stable and efficient numerical algorithms to solve these matrix equations for small- to medium- sized matrices. However, developing numerical algorithms for solving large-scale T-Sylvester equations was a completely open problem. In this thesis, we introduce several projection methods based on block Krylov subspaces and extended block Krylov subspaces for solving the T-Sylvester equation when the right-hand side is a low-rank matrix. We also offer an intuition on the expected convergence of the algorithm based on block Krylov subspaces and a clear guidance on which algorithm is the most convenient to use in each situation. All the algorithms presented in this thesis have been extensively tested, and the reported numerical results show that they perform satisfactorily in practice.Adicionalmente se recibió ayuda parcial de los proyectos de investigación: “Structured Numerical Linear Algebra: Matrix Polynomials, Special Matrices, and Conditioning” (Ministerio de Economía y Competitividad de España, Número de proyecto: MTM2012-32542) y “Structured Numerical Linear Algebra for Constant, Polynomial and Rational Matrices” (Ministerio de Economía y Competitividad de España, Número de proyecto: MTM2015-65798-P), donde el investigador principal de ambos proyectos fue Froilán Martínez Dopico.Programa Oficial de Doctorado en Ingeniería MatemáticaPresidente: José Mas Marí.- Secretario: Fernando de Terán Vergara.- Vocal: José Enrique Román Molt

Sparse Gr\"obner Bases: the Unmixed Case

Toric (or sparse) elimination theory is a framework developped during the last decades to exploit monomial structures in systems of Laurent polynomials. Roughly speaking, this amounts to computing in a \emph{semigroup algebra}, \emph{i.e.} an algebra generated by a subset of Laurent monomials. In order to solve symbolically sparse systems, we introduce \emph{sparse Gr\"obner bases}, an analog of classical Gr\"obner bases for semigroup algebras, and we propose sparse variants of the $F_5$ and FGLM algorithms to compute them. Our prototype "proof-of-concept" implementation shows large speed-ups (more than 100 for some examples) compared to optimized (classical) Gr\"obner bases software. Moreover, in the case where the generating subset of monomials corresponds to the points with integer coordinates in a normal lattice polytope $\mathcal P\subset\mathbb R^n$ and under regularity assumptions, we prove complexity bounds which depend on the combinatorial properties of $\mathcal P$. These bounds yield new estimates on the complexity of solving $0$-dim systems where all polynomials share the same Newton polytope (\emph{unmixed case}). For instance, we generalize the bound $\min(n_1,n_2)+1$ on the maximal degree in a Gr\"obner basis of a $0$-dim. bilinear system with blocks of variables of sizes $(n_1,n_2)$ to the multilinear case: $\sum n_i - \max(n_i)+1$. We also propose a variant of Fr\"oberg's conjecture which allows us to estimate the complexity of solving overdetermined sparse systems.Comment: 20 pages, Corollary 6.1 has been corrected, ISSAC 2014, Kobe : Japan (2014

Solving Sparse Integer Linear Systems

We propose a new algorithm to solve sparse linear systems of equations over the integers. This algorithm is based on a $p$-adic lifting technique combined with the use of block matrices with structured blocks. It achieves a sub-cubic complexity in terms of machine operations subject to a conjecture on the effectiveness of certain sparse projections. A LinBox-based implementation of this algorithm is demonstrated, and emphasizes the practical benefits of this new method over the previous state of the art