Search

In the last 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.

The randomized singular value decomposition proposed in [12] has certainly become one of the most well-established randomization-based algorithms in numerical linear algebra. The key ingredient of the entire procedure is the computation of a subspace which is close to the column space of the target matrix $\mathbf{A}$ up to a certain probabilistic confidence. In this paper we propose a modification to the standard randomized SVD procedure which leads, in general, to better approximations to $\text{Range}(\mathbf{A})$ at the same computational cost. To this end, we explicitly construct information from the row space of $\mathbf{A}$ enhancing the quality of our approximation. We also observe that very few pieces of information from $\text{Range}(\mathbf{A}^T)$ are indeed necessary. We thus design a variant of our algorithm equipped with a subsampling step which largely increases the efficiency of our procedure while attaining competitive accuracy records. Our findings are supported by both theoretical analysis and numerical results., Comment: 24 pages, 6 figures

Sequences of parametrized Lyapunov equations can be encountered in many application settings. Moreover, solutions of such equations are often intermediate steps of an overall procedure whose main goal is the computation of $\text{trace}(EX)$ where $X$ denotes the solution of a Lyapunov equation and $E$ is a given matrix. We are interested in addressing problems where the parameter dependency of the coefficient matrix is encoded as a low-rank modification to a \emph{seed}, fixed matrix. We propose two novel numerical procedures that fully exploit such a common structure. The first one builds upon the Sherman-Morrison-Woodbury (SMW) formula and recycling Krylov techniques, and it is well-suited for small dimensional problems as it makes use of dense numerical linear algebra tools. The second algorithm can instead address large-scale problems by relying on state-of-the-art projection techniques based on the extended Krylov subspace. We test the new algorithms on several problems arising in the study of damped vibrational systems and the analyses of output synchronization problems for multi-agent systems. Our results show that the algorithms we propose are superior to state-of-the-art techniques as they are able to remarkably speed up the computation of accurate solutions.

Of all the possible projection methods for solving large-scale Lyapunov matrix equations, Galerkin approaches remain much more popular than minimal-residual ones. This is mainly due to the different nature of the projected problems stemming from these two families of methods. While a Galerkin approach leads to the solution of a low-dimensional matrix equation per iteration, a matrix least-squares problem needs to be solved per iteration in a minimal-residual setting. The significant computational cost of these least-squares problems has steered researchers towards Galerkin methods in spite of the appealing properties of minimal-residual schemes. In this paper we introduce a framework that allows for modifying the Galerkin approach by low-rank, additive corrections to the projected matrix equation problem with the two-fold goal of attaining monotonic convergence rates similar to those of minimal-residual schemes while maintaining essentially the same computational cost of the original Galerkin method. We analyze the well-posedness of our framework and determine possible scenarios where we expect the residual norm attained by two low-rank-modified variants to behave similarly to the one computed by a minimal-residual technique. A panel of diverse numerical examples shows the behavior and potential of our new approach.

Thanks to its great potential in reducing both computational cost and memory requirements, combining sketching and Krylov subspace techniques has attracted a lot of attention in the recent literature on projection methods for linear systems, matrix function approximations, and eigenvalue problems. Applying this appealing strategy in the context of linear matrix equations turns out to be far more involved than a straightforward generalization. These difficulties include analyzing well-posedness of the projected problem and deriving possible error estimates depending on the sketching properties. Further computational complications include the lack of a natural residual norm estimate and of an explicit basis for the generated subspace. In this paper we propose a new sketched-and-truncated polynomial Krylov subspace method for Sylvester equations that aims to address all these issues. The potential of our novel approach, in terms of both computational time and storage demand, is illustrated with numerical experiments. Comparisons with a state-of-the-art projection scheme based on rational Krylov subspaces are also included.

Adaptive cubic regularization methods for solving nonconvex problems need the efficient computation of the trial step, involving the minimization of a cubic model. We propose a new approach in which this model is minimized in a low dimensional subspace that, in contrast to classic approaches, is reused for a number of iterations. Whenever the trial step produced by the low-dimensional minimization process is unsatisfactory, we employ a regularized Newton step whose regularization parameter is a by-product of the model minimization over the low-dimensional subspace. We show that the worst-case complexity of classic cubic regularized methods is preserved, despite the possible regularized Newton steps. We focus on the large class of problems for which (sparse) direct linear system solvers are available and provide several experimental results showing the very large gains of our new approach when compared to standard implementations of adaptive cubic regularization methods based on direct linear solvers. Our first choice as projection space for the low-dimensional model minimization is the polynomial Krylov subspace; nonetheless, we also explore the use of rational Krylov subspaces in case where the polynomial ones lead to less competitive numerical results.

Among randomized numerical linear algebra strategies, so-called sketching procedures are emerging as effective reduction means to accelerate the computation of Krylov subspace methods for, e.g., the solution of linear systems, eigenvalue computations, and the approximation of matrix functions. While there is plenty of experimental evidence showing that sketched Krylov solvers may dramatically improve performance over standard Krylov methods, many features of these schemes are still unexplored. We derive a new sketched Arnoldi-type relation that allows us to obtain several different new theoretical results. These lead to an improvement of our understanding of sketched Krylov methods, in particular by explaining why the frequently occurring sketched Ritz values far outside the spectral region of A do not negatively influence the convergence of sketched Krylov methods for f (A)b. Our findings also help to identify, among several possible equivalent formulations, the most suitable sketched approximations according to their numerical stability properties. These results are also employed to analyze the error of sketched Krylov methods in the approximation of the action of matrix functions, significantly contributing to the theory available in the current literature.

Time-parallel time integration has received a lot of attention in the high performance computing community over the past two decades. Indeed, it has been shown that parallel-in-time techniques have the potential to remedy one of the main computational drawbacks of parallel-in-space solvers. In particular, it is well-known that for large-scale evolution problems space parallelization saturates long before all processing cores are effectively used on today's large scale parallel computers. Among the many approaches for time-parallel time integration, ParaDiag schemes have proved themselves to be a very effective approach. In this framework, the time stepping matrix or an approximation thereof is diagonalized by Fourier techniques, so that computations taking place at different time steps can be indeed carried out in parallel. We propose here a new ParaDiag algorithm combining the Sherman-Morrison-Woodbury formula and Krylov techniques. A panel of diverse numerical examples illustrates the potential of our new solver. In particular, we show that it performs very well compared to different ParaDiag algorithms recently proposed in the literature.

Algorithms for data assimilation try to predict the most likely state of a dynamical system by combining information from observations and prior models. Variational approaches, such as the weak-constraint four-dimensional variational data assimilation formulation considered in this problem, can ultimately be interpreted as a minimization problem. One of the main challenges of such a formulation is the solution of large linear systems of equations which arise within the inner linear step of the adopted nonlinear solver. Depending on the adopted approach, these linear algebraic problems amount to either a saddle point linear system or a symmetric positive definite (SPD) one. Both formulations can be solved by means of a Krylov method, like GMRES or CG, that needs to be preconditioned to ensure fast convergence in terms of the number of iterations. In this paper we illustrate novel, efficient preconditioning operators which involve the solution of certain Stein matrix equations. In addition to achieving better computational performance, the latter machinery allows us to derive tighter bounds for the eigenvalue distribution of the preconditioned linear system for certain problem settings. A panel of diverse numerical results displays the effectiveness of the proposed methodology compared to current state-of-the-art approaches.

One of the most computationally expensive steps of the low-rank ADI method for large-scale Lyapunov equations is the solution of a shifted linear system at each iteration. We propose the use of the extended Krylov subspace method for this task. In particular, we illustrate how a single approximation space can be constructed to solve all the shifted linear systems needed to achieve a prescribed accuracy in terms of Lyapunov residual norm. Moreover, we show how to fully merge the two iterative procedures in order to obtain a novel, efficient implementation of the low-rank ADI method, for an important class of equations. Many state-of-the-art algorithms for the shift computation can be easily incorporated into our new scheme, as well. Several numerical results illustrate the potential of our novel procedure when compared to an implementation of the low-rank ADI method based on sparse direct solvers for the shifted linear systems.

We identify a relationship between the solutions of a nonsymmetric algebraic T-Riccati equation (T-NARE) and the deflating subspaces of a palindromic matrix pencil, obtained by arranging the coefficients of the T-NARE. The interplay between T-NARE and palindromic pencils allows one to derive both theoretical properties of the solutions of the equation, and new methods for its numerical solution. In particular, we propose methods based on the (palindromic) QZ algorithm and the doubling algorithm, whose effectiveness is demonstrated by several numerical tests

We introduce a very weak space-time variational formulation for the wave equation, prove its well-posedness (even in the case of minimal regularity) and optimal inf-sup stability. Then, we introduce a tensor product-style space-time Petrov-Galerkin discretization with optimal discrete inf-sup stability, obtained by a non-standard definition of the trial space. As a consequence, the numerical approximation error is equal to the residual, which is particularly useful for a posteriori error estimation. For the arising {discrete linear systems} in space and time, we introduce efficient numerical solvers that appropriately exploit the equation structure, either at the preconditioning level or in the approximation phase by using a tailored Galerkin projection. This Galerkin method shows competitive behavior concerning {wall-clock} time, accuracy and memory as compared with a standard time-stepping method in particular in low regularity cases. Numerical experiments with a 3D (in space) wave equation illustrate our findings., Comment: 25 pages, 4 Tabes, 4 figures

The Loewner framework is one of the most successful data-driven model order reduction techniques. If $N$ is the cardinality of a given data set, the so-called Loewner and shifted Loewner matrices $\mathbb{L}\in\mathbb{C}^{N\times N}$ and $\mathbb{S}\in\mathbb{C}^{N\times N}$ can be defined by solely relying on information encoded in the considered data set and they play a crucial role in the computation of the sought rational model approximation. In particular, the singular value decomposition of a linear combination of $\mathbb{S}$ and $\mathbb{L}$ provides the tools needed to construct accurate models which fulfill important approximation properties with respect to the original data set. However, for highly-sampled data sets, the dense nature of $\mathbb{L}$ and $\mathbb{S}$ leads to numerical difficulties, namely the failure to allocate these matrices in certain memory-limited environments or excessive computational costs. Even though they do not possess any sparsity pattern, the Loewner and shifted Loewner matrices are extremely structured and, in this paper, we show how to fully exploit their Cauchy-like structure to reduce the cost of computing accurate rational models while avoiding the explicit allocation of $\mathbb{L}$ and $\mathbb{S}$. In particular, the use of the \emph{hierarchically semiseparable} format allows us to remarkably lower both the computational cost and the memory requirements of the Loewner framework obtaining a novel scheme whose costs scale with $N \log N$.

Rational Krylov subspaces have become a reference tool in dimension reduction procedures for several application problems. When data matrices are symmetric, a short-term recurrence can be used to generate an associated orthonormal basis. In the past this procedure was abandoned because it requires twice the number of linear system solves per iteration compared with the classical long-term method. We propose an implementation that allows one to obtain the rational subspace reduced matrices at lower overall computational costs than proposed in the literature by also conveniently combining the two system solves. Several applications are discussed where the short-term recurrence feature can be exploited to avoid storing the whole orthonormal basis. We illustrate the advantages of the proposed procedure with several examples.

The nonsymmetric T-Riccati equation is a quadratic matrix equation where the linear part corresponds to the so-called T-Sylvester or T-Lyapunov operator that has previously been studied in the literature. It has applications in macroeconomics and policy dynamics. So far, it presents an unexplored problem in numerical analysis, and both, theoretical results and computational methods, are lacking in the literature. In this paper we provide some sufficient conditions for the existence and uniqueness of a nonnegative minimal solution and its efficient computation is deeply analyzed. Both the small-scale and the large-scale setting are addressed and Newton-Kleinman-like methods are derived. The convergence of these procedures to the minimal solution is proved and several numerical results illustrate the computational efficiency of the proposed methods.

Block Krylov subspace methods (KSMs) comprise building blocks in many state-of-the-art solvers for large-scale matrix equations as they arise, e.g., from the discretization of partial differential equations. While extended and rational block Krylov subspace methods provide a major reduction in iteration counts over polynomial block KSMs, they also require reliable solvers for the coefficient matrices, and these solvers are often iterative methods themselves. It is not hard to devise scenarios in which the available memory, and consequently the dimension of the Krylov subspace, is limited. In such scenarios for linear systems and eigenvalue problems, restarting is a well explored technique for mitigating memory constraints. In this work, such restarting techniques are applied to polynomial KSMs for matrix equations with a compression step to control the growing rank of the residual. An error analysis is also performed, leading to heuristics for dynamically adjusting the basis size in each restart cycle. A panel of numerical experiments demonstrates the effectiveness of the new method with respect to extended block KSMs.

Variational formulations of time-dependent PDEs in space and time yield $(d+1)$-dimensional problems to be solved numerically. This increases the number of unknowns as well as the storage amount. On the other hand, this approach enables adaptivity in space and time as well as model reduction w.r.t. both type of variables. In this paper, we show that matrix oriented techniques can significantly reduce the computational timings for solving the arising linear systems outperforming both time-stepping schemes and other solvers.

Low-rank Krylov methods are one of the few options available in the literature to address the numerical solution of large-scale general linear matrix equations. These routines amount to well-known Krylov schemes that have been equipped with a couple of low-rank truncations to maintain a feasible storage demand in the overall solution procedure. However, such truncations may affect the convergence properties of the adopted Krylov method. In this paper we show how the truncation steps have to be performed in order to maintain the convergence of the Krylov routine. Several numerical experiments validate our theoretical findings.

We show that the discrete operator stemming from the time and space discretization of evolutionary partial differential equations can be represented in terms of a single Sylvester matrix equation. A novel solution strategy that combines projection techniques with the full exploitation of the entry-wise structure of the involved coefficient matrices is proposed. The resulting scheme is able to efficiently solve problems with a tremendous number of degrees of freedom while maintaining a low storage demand as illustrated in several numerical examples.

Galerkin and Petrov-Galerkin methods are some of the most successful solution procedures in numerical analysis. Their popularity is mainly due to the optimality properties of their approximate solution. We show that these features carry over to the (Petrov-)Galerkin methods applied for the solution of linear matrix equations. Some novel considerations about the use of Galerkin and Petrov-Galerkin schemes in the numerical treatment of general linear matrix equations are expounded and the use of constrained minimization techniques in the Petrov-Galerkin framework is proposed.

The numerical solution of the algebraic Riccati equation is a challenging task especially for very large problem dimensions. In this paper we present a new algorithm that combines the very appealing computational features of projection methods with the convergence properties of the inexact Newton-Kleinman procedure equipped with a line search. In particular, the Newton scheme is completely merged in a projection framework with a single approximation space so that the Newton-Kleinman iteration is only implicitly performed. Moreover, the line search turns out to be exact in our setting, i.e., the existence of a local minimum of the Riccati residual norm along the current search direction is guaranteed and the corresponding minimizer is chosen as step-size. This property determines a monotone decrease of the Riccati residual norm under some mild assumptions. Several numerical results are reported to illustrate the potential of our novel approach.

We consider the problem of efficiently solving Sylvester and Lyapunov equations of medium and large scale, in case of rank-structured data, i.e., when the coefficient matrices and the right-hand side have low-rank off-diagonal blocks. This comprises problems with banded data, recently studied by Haber and Verhaegen in "Sparse solution of the Lyapunov equation for large-scale interconnected systems", Automatica, 2016, and by Palitta and Simoncini in "Numerical methods for large-scale Lyapunov equations with symmetric banded data", SISC, 2018, which often arise in the discretization of elliptic PDEs. We show that, under suitable assumptions, the quasiseparable structure is guaranteed to be numerically present in the solution, and explicit novel estimates of the numerical rank of the off-diagonal blocks are provided. Efficient solution schemes that rely on the technology of hierarchical matrices are described, and several numerical experiments confirm the applicability and efficiency of the approaches. We develop a MATLAB toolbox that allows easy replication of the experiments and a ready-to-use interface for the solvers. The performances of the different approaches are compared, and we show that the new methods described are efficient on several classes of relevant problems.

The numerical solution of large-scale Lyapunov matrix equations with symmetric banded data has so far received little attention in the rich literature on Lyapunov equations. We aim to contribute to this open problem by introducing two efficient solution methods, which respectively address the cases of well conditioned and ill conditioned coefficient matrices. The proposed approaches conveniently exploit the possibly hidden structure of the solution matrix so as to deliver memory and computation saving approximate solutions. Numerical experiments are reported to illustrate the potential of the described methods.

We consider generalizations of the Sylvester matrix equation, consisting of the sum of a Sylvester operator and a linear operator $\Pi$ with a particular structure. More precisely, the commutator of the matrix coefficients of the operator $\Pi$ and the Sylvester operator coefficients are assumed to be matrices with low rank. We show (under certain additional conditions) low-rank approximability of this problem, i.e., the solution to this matrix equation can be approximated with a low-rank matrix. Projection methods have successfully been used to solve other matrix equations with low-rank approximability. We propose a new projection method for this class of matrix equations. The choice of subspace is a crucial ingredient for any projection method for matrix equations. Our method is based on an adaption and extension of the extended Krylov subspace method for Sylvester equations. A constructive choice of the starting vector/block is derived from the low-rank commutators. We illustrate the effectiveness of our method by solving large-scale matrix equations arising from applications in control theory and the discretization of PDEs. The advantages of our approach in comparison to other methods are also illustrated.

In the numerical treatment of large-scale Sylvester and Lyapunov equations, projection methods require solving a reduced problem to check convergence. As the approximation space expands, this solution takes an increasing portion of the overall computational effort. When data are symmetric, we show that the Frobenius norm of the residual matrix can be computed at significantly lower cost than with available methods, without explicitly solving the reduced problem. For certain classes of problems, the new residual norm expression combined with a memory-reducing device make classical Krylov strategies competitive with respect to more recent projection methods. Numerical experiments illustrate the effectiveness of the new implementation for standard and extended Krylov subspace methods.

We are interested in the numerical solution of nonsymmetric linear systems arising from the discretization of convection-diffusion partial differential equations with separable coefficients and dominant convection. Preconditioners based on the matrix equation formulation of the problem are proposed, which naturally approximate the original discretized problem. For certain types of convection coefficients, we show that the explicit solution of the matrix equation can effectively replace the linear system solution. Numerical experiments with data stemming from two and three dimensional problems are reported, illustrating the potential of the proposed methodology.