This workshop is the eighth in a series of workshops on matrix equations and tensor techniques. As in the previous meetings, the focus will be on the latest developments in the theory, computation and applications of linear and nonlinear matrix equations and tensor equations.
Sparse tensor decompositions have become increasingly popular in the literature due to their capability to naturally model high dimensional sparse data with many features, and glean hidden relations owing to underlying low-rank structures within the data. They have been successfully employed in many application settings including recommender systems, graph analytics, healthcare data analytics, web search, cyber security, and many others. The aptitude of tensor methods for such big data analysis applications solicited the development of efficient parallel tensor factorization algorithms capable of handling datasets of billions of entries, which has been among the most trending areas of research in the HPC community in the recent past.
In this talk, we will discuss parallelization techniques for sparse tensor factorization together with various partitioning strategies for balancing computation/memory costs and reducing communication. We will compare advantages and limitations of various approaches and touch upon outstanding challenges for better parallel scalability. We will conclude the talk with an overview of the capabilities of the PACOS library (PArtitioning and COmmunication framework for Sparse irregular applications) that enable devising efficient and scalable parallel sparse tensor factorization kernels as well as partitioning routines, and facilitate reproducibility of scalability results.
We propose guaranteed and fully computable upper bound on the energy norm of the error in low rank Tensor Train (TT) approximate solutions of (possibly) high dimensional reaction-diffusion problems. The error bound is obtained from Euler-Lagrange equations for a complementary flux reconstruction problem, which are solved in the low rank TT representation using the block Alternating Linear Scheme. This bound is guaranteed to be above the energy norm of the total error, including the discretization error, the tensor approximation error, and the error in the solver of linear algebraic equations. Numerical examples with the Poisson equation and the Schroedinger equation with the Henon-Heiles potential in up to 40 dimensions will be presented to illustrate the efficiency of this approach.
We consider the rank-structured tensor approach for numerical modeling of long-range potentials in many-particle systems. The method of grid-based assembled tensor summation of the electrostatic potentials on 3D finite lattices [3] exhibits the computational complexity of the order of $O(L)$ which is much less than $O(L^3)$ in traditional Ewald-type summation.
Novel range-separated (RS) tensor format [4] applies to many-particle systems of general type. These can be the free space electrostatic potentials of large biomolecules or the multidimensional scattered data modeled by radial basis functions. The main advantage of the RS tensor format is that the rank of the canonical/Tucker tensor representing the sum of long range contributions from all particles in the collective potential depends only logarithmically on the number of particles $N$. The basic tool for calculation of the RS tensor representation is the reduced higher order SVD (RHOSVD) introduced in [5]. The representation complexity of the short range part is $O(N)$ with a small prefactor independent on the number of particles. The interaction energies and forces of the many-particle system can be computed by using only the long-range part of the collective potential, with representation complexity $O(n \mathrm{log} N )$, where n is the univariate grid size.
The new regularization scheme for the Poisson-Boltzmann equation (PBE) describing the electrostatic potential in biomolecules is based on the RS tensor representation to the discretized Dirac delta [2]. It leads to solving a single system of FDM/FEM equations only for the smooth long-range part of the initially singular right-hand side of the PBE [1]. The resulting solution of PBE is the sum of the long- and short-range parts, where the latter is precomputed by the direct tensor summations, without solving the PBE. The numerical examples are presented.
[1] P. Benner, V. Khoromskaia, B. N. Khoromskij, C. Kweyu and M. Stein. Computing Electrostatic Potentials of Biomolecules using Regularization based on the Range-separated Tensor Format. arXiv:1901.09864, 2019.
[2] Boris N. Khoromskij. Tensor Numerical Methods in Scientific Computing. De Gruyter, Berlin, 2018.
[3] Venera Khoromskaia and Boris N. Khoromskij. Tensor Numerical Methods in Quantum Chemistry. De Gruyter, Berlin, 2018.
[4] P. Benner, V. Khoromskaia and B. N. Khoromskij. Range-Separated Tensor Format for Many-particle Modeling. SIAM J Sci. Comput., 40 (2), A1034–A1062, 2018.
[5] B. N. Khoromskij and V. Khoromskaia. Multigrid Tensor Approximation of Function Related Arrays. SIAM J Sci. Comput., 31(4), 3002-3026, 2009.
We prove the convergence of the multigrid method for parameter-dependent symmetric positive definite linear systems. We present the smoothing and approximation properties for such problems and prove the smoothing property for the parameter-dependent damped Richardson as well as the parameter-dependent damped Jacobi method.
Our theoretical results require a parameter-dependent representation of the operator, the solution, the right-hand side, the smoother as well as the prolongation and restriction. We use tensor formats to derive such a representation for linear systems with affine-linear operators. For the parameter-dependent damped Jacobi method, we derive an approximation using exponential sums and prove the smoothing property for this representation.
In numerical experiments, we observe that the convergence behavior of the multigrid method with the damped Jacobi smoother in the parameter-dependent case seems to be grid size independent.
We discuss the problem of computing roots of systems of multivariate polynomials by approaches from numerical linear and multilinear algebra.The key concept in our approach is the Macaulay matrix, a large and highly structured matrix that contains the coefficients of the polynomials systems. The roots can be retrieved from the nullspace of this matrix. We then show how this root retrieval can be carried out using tensors and, in particular, tensor decompositions.
Automatic Face Recognition has become increasingly important in the past few years due to its several applications in daily life. Numerical linear algebra tools have been extensively used for classification purposes. However, since several factors can affect the image, multilinear algebra tools seem to be a natural choice to deal with image classification.
We propose a new algorithm based on Tensor-Train Decomposition and we compare it with other related methods, such as the ones based on SVD and HOSVD (High-Order Singular Value Decomposition).
get together on self-paid basis
We consider the nonsymmetric T-Riccati equation
$$
0 = \mathcal{R}_T(X):=DX+X^TA-X^TBX+C,\qquad (1)
$$
where $A,B,C,D\in\mathbb{R}^{n\times n}$ and sufficient conditions for the existence and uniqueness of a minimal (w.r.t. entry-wise comparison) solution $X_{\min}\in\mathbb{R}^{n\times n}$ are provided. To date, the nonlinear matrix equation (1) is still an unexplored problem in numerical analysis and both theoretical results and computational methods are lacking in the literature. We provide some sufficient conditions for the existence and uniqueness of a nonnegative minimal solution and discuss its efficient computation. Both the small-scale and the large-scale settings 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.
The numerical solution of the algebraic Riccati matrix equation
$$(1)\quad A^T X + XA − XSX + Q = 0,$$
where $A$, $S$, $Q\in\mathbb{R}^{n\times n}$ , is an interesting and still challenging task especially when the problem dimension is very large, say $n > 10^4$ , as the dense solution $X$ cannot be store and a memory-saving approximation has to be sought.
A vast portion of the literature focuses on the case of low-rank $S$ and $Q$ and many, diverse solution schemes have been developed to tackle this kind of equations. See, e.g., [1] and the references therein. In particular, the so-called low-rank methods compute a low-rank approximation to the exact solution $X$ in order to moderate the memory requirements of the adopted algorithm.
By exploiting novel results about the solution of Lyapunov equations with non low-rank right-hand side [2,3], a Newton iteration for (1) with a rank-structured $Q$ has been recently proposed [4]. However, also in such a scheme the matrix $S$ is still supposed to be low-rank.
In this talk we consider Riccati equations with banded, full-rank coefficent matrices $A$, $S$ and $Q$ and, by taking inspiration from some early results by Incertis [5], a fresh solution procedure that efficiently computes memory-saving approximate solutions is proposed. In particular, the structure of the computed solution $\tilde X$ depends on some properties of the matrices $A$ and $Q$ and we can represent $\tilde X$ in terms of either a banded matrix or a banded plus low-rank matrix maintaining a very low storage demand of the overall solution process.
Several numerical results are reported to illustrate the potential of the discussed method.
References
[1] P. Benner, Z. Bujanovic, P. Kürschner and J. Saak. A numerical comparison of solvers for large-scale, continuous-time algebraic Riccati equations. ArXiv Preprint: 1811.00850, 2018.
[2] D. Palitta and V. Simoncini. Numerical methods for large-scale Lyapunov equations with symmetric banded data. SIAM Journal on Scientific Computing, 40 (5): A3581–A3608, 2018.
[3] S. Massei, D. Palitta and L. Robol. Solving rank structured Sylvester and Lyapunov equations. SIAM Journal on Matrix Analysis and Applications, 39 (4): 1564–1590, 2018.
[4] D. Kressner, S. Massei and L. Robol. Low-rank updates and a divide-and-conquer method for linear matrix equations. SIAM Journal on Scientific Computing, 41 (2): A848–A876, 2019.
[5] F. C. Incertis. A new formulation of the algebraic Riccati equation problem. IEEE Transactions on Automatic Control, 26 (3): 768–770, 1981.
Our goal is the feedback stabilization of a two-dimensional two-phase Stefan problem coupled with (Navier-)Stokes equations. The Stefan problem can model solidification and melting of pure materials and gets its name from the purely algebraic Stefan condition which describes the coupling between the temperature of the material and its melting process.
After linearization and discretization, the stabilization problem results in a non-autonomous differential Riccati equation (DRE) with differential-algebraic structure. While the actual index of the resulting DAE needs further investigation, the problem features a structure that combines that of semi-explicit index-1 DAEs (resulting from the Stefan condition) and Stokes-type index-2 DAEs (resulting from the dynamics of the convection).
Moreover, the interface between the two phases in the domain evolves over time, which causes all coefficients of the resulting DRE to be time-varying. Thus, existing DRE solvers have to be adapted to this highly non-autonomous case which has significantly increased computational costs and memory requirements.
We present techniques to treat the specific structure of this coupled problem
in the solver for the DREs and show first results of the application of our feedback stabilization.
We consider the numerical solution of large-scale, differential matrix Riccati equations (DRE). Projection methods have recently arisen as a promising class of solution strategies. Existing approaches in this class focus on polynomial or extended Krylov subspaces as approximation space. We show that great computational and memory advantages are obtained with fully rational Krylov subspaces and we discuss several crucial issues such as efficient time stepping and stopping criteria. Numerical experiments illustrate the procedure's effectiveness.
We consider the differential Riccati equation,
$$\dot{X} = A^T X + X A - X BB^T X + C^T C.$$
The differential Riccati equation as well as the algebraic Riccati equation play important roles in applied mathematics like control theory and system theory. In our talk, we focus on the large-scale case. The numerical solution of these equation is challenging, in particular, because of the enormous amount of storage. A general approach, that has led to several algorithms, bases on an invariant subspace $Q \subseteq \mathbb R^{n\times n}$ such that $X(t) \in Q$ for all $t$. After identifying a suitable invariant subspace, we develop a Galerkin approach for the numerical solution of the differential Riccati equation. We review Davison-Maki methods for the numerical solution of the resulting equation.
Matrix equations have arisen as the natural setting for various PDE discretization methods such as finite differences, isogeometric analysis, spectral and finite elements.
Thanks to major recent computational advances, solving certain classes of linear matrix equations is a competitive alternative to dealing with the large (vector) linear systems classically stemming from the aforementioned discretizations. In this talk we support these considerations with examples from the numerical treatment of possibly time-dependent PDE problems.
In this paper, we present an efficient algorithm for model reduction of discrete-time index-2 descriptor systems arising in context of linear time invariant (LTI) control and stability of descriptor systems. We propose a balanced truncation based model order reduction (MOR) strategy which consists of two stages. In the first stage, we reformulate the descriptor system into a generalized system by manipulating its system structure. Once the reformulated generalized system is obtained, it is applicable for a balanced truncation-based model order reduction strategy. The second stage of our work focuses is on a MOR method where the descriptor system does not require to be transformed into the generalized system. Rather, the balanced truncation MOR method can directly be applied on the index-2 descriptor system. We also implement a Smith based iterative method to compute the solutions of the Lyapunov equations corresponding to the transformed system without having to compute the projection matrices. Results from various numerical simulations are included to verify the high performance and accuracy of our proposed algorithm.
The treatment of two-dimensional random walks in the quarter plane leads to Markov processes which involve semi-infinite matrices having Toeplitz or block Toeplitz structure plus a low-rank correction. In particular, finding the steady state probability distribution requires computing the minimal nonnegative solution of the quadratic matrix equation $A_1X^2+A_0X+A_{-1}=X$, where the coefficients $A_i$, $i=-1,0,1$, are semi-infinite nonnegative matrices, which are Toeplitz except for a low-rank correction.
We introduce the sets $\mathcal{QT}^d$ and $\mathcal{EQT}$ of semi-infinite matrices with bounded infinity norm [2]. The former is made by matrices representable as a sum of a Toeplitz matrix and a compact correction with columns having entries which decay to zero. The latter is formed by matrices in $\mathcal{QT}^d$ plus a further correction of the kind $ev^T$ for $e^T=(1,1,\ldots)$ and $v=(v_i)_{i\in Z^+}$ such that $\sum_{i=1}^\infty |v_i|<\infty$. We prove that $\mathcal{QT}^d$ and $\mathcal{EQT}$ are Banach algebras, i.e., they are Banach spaces with the infinity norm, closed under matrix multiplication. Moreover, we provide mild conditions under which the minimal nonnegative solution belongs either to $\mathcal{QT}^d$ or to $\mathcal{EQT}$.
These new classes allow to treat random walks with reset, which model queueing systems with breakdowns or restarts. These models cannot be handled by relying on the currently available Banach algebras [1]. Moreover, matrices in both classes can be approximated up to any precision by a finite number of parameters. This allows to introduce a matrix arithmetic and to extend the cqt-toolbox of [3], so that numerical algorithms valid for finite matrices can be applied to solve the quadratic matrix equation. The convergence of various fixed point iterations is analyzed [4].
[1] Bini, D.A., Massei, S., Meini, B.: Semi-infinite quasi-Toeplitz matrices with applications to QBD stochastic processes. Math. Comp. 87(314), 2811-2830 (2018)
[2] Bini, D.A., Massei, S., Meini, B., Robol, L.: A computational framework for two-dimensional random walks with restarts. Submitted for publication. arXiv:1909.11372v1. (2019)
[3] Bini, D.A., Massei, S., Robol, L.: Quasi-Toeplitz matrix arithmetic: a MATLAB toolbox. Numer. Algorithms 81(2), 741-769 (2019)
[4] Bini, D.A., Meini, B., Meng, J.: Solving quadratic matrix equations arising in random walks in the quarter plane. Submitted for publication. arXiv:1907.09796 (2019)
PDE-constrained optimization problems arise in a broad number of applications. The resulting large-scale saddle-point systems are challenging to solve and acquiring a full solution is often infeasible. We present a new framework to find a low-rank approximation to the solution by reformulating the system into a system of Sylvester-like matrix equations. These matrix equations are subsequently projected onto a small subspace via rational Krlyov-subspace iterations and we obtain a reduced problem by imposing a Galerkin condition on its residual. In our presentation we discuss implementation details and dependence on the problem parameters. Numerical experiments will illustrate the performance of the new strategy.
Consider the algebraic Riccati equation
$$A^{*}X + XA + C^*C - XBB^{*}X = 0$$
with large and sparse $A\in \mathbb{C}^{n\times n}$ and $B\in \mathbb{C}^{n\times m}$, $C\in \mathbb{C}^{p\times n}$. The goal is to find an approximate solution $X$ such that the rank of the Riccati residual remains equal to $p$. A method to determine an approximate solution to Riccati equations with this rank property is given by the RADI method [1], a generalization of the ADI method for linear matrix equations.
We present a new approach in which solves with $(A^*-\mu I_n)$ and the explicitly known residual factors are sufficient to obtain the unique low-rank residual solution.
[1] Benner, P., Bujanović, Z., Kürschner, P. et al. RADI: a low-rank ADI-type algorithm for large scale algebraic Riccati equations. Numer. Math. (2018) 138: 301. https://doi.org/10.1007/s00211-017-0907-5
A stochastic cubical tensor is a three-mode array (or hypermatrix) with nonnegative entries, whose $1$-mode fibers (i.e., columns) sum up to $1$. Such tensors appear in certain higher-order Markov chains and random walks with memory, exactly as stochastic matrices describe classical discrete Markov chains and random walks.
The interest in higher-order stochastic processes is significantly growing in recent years as they are much better at modeling high dimensional data and nonlinear dynamics in numerous applications. However, fundamental questions such as the convergence of the process towards a limiting distribution and the uniqueness of such a limit are still not well understood and are the subject of a rich recent literature.
We introduce a set of ergodicity coefficients for stochastic cubical tensors, by extending certain definitions known in the matrix case to the tensor setting. The proposed coefficients yield new explicit formulas that
(a) guarantee the uniqueness of nonlinear Perron eigenvectors of stochastic tensors,
(b) provide bounds on the sensitivity of such eigenvectors with respect to changes in the tensor and
(c) ensure the convergence of various higher-order Markov chains to the stationary distribution.
Joint work with Francesco Tudisco.
Graph (or network) analysis is concerned with knowledge discovery from graph data. Corresponding analysis methods usually have their roots in graph theory, matrix computations, and statistics.
A growing class of analysis algorithms consists of methods based on the graph's Laplacian matrix. These include traditional methods using Laplacian eigenvectors such as spectral partitioning, spectral clustering, or spectral embedding. In the past 10-15 years, Laplacian linear systems received significant interest in the context of graph analysis as well. They have applications in graph partitioning, graph sparsification, electrical networks, and graph robustness, to name just a few.
Our poster focuses on Laplacian linear systems in the context of centrality measures. Centrality measures indicate the importance of a node (or edge) in the graph based on its position and belong to the most important analysis methods. Typically, they assign a real number to each node (or edge), which yields a ranking. On the poster we review (i) the state of the art concerning efficient algorithms for computing Laplacian-based centrality measures and (ii) describe recent and ongoing work on accelerating such algorithms. We also portray briefly our implementation of the Lean Algebraic Multigrid (LAMG) solver in NetworKit, our growing open-source software for large-scale graph/network analysis.
Among other things, M.E.S.S. can solve Lyapunov and Riccati equations, and perform model reduction of systems in state space and structured differential algebraic form. While until version 1.0.1 the work focused on algebraic equations and autonomous differential Riccati equations (DREs), the latest release adds support also for non-autonomous DREs.
In the first place, the poster presented will give an overview on the matrix equations that appear in the low-rank parameter-dependent fluid-structure interaction (FSI) framework. In contrast to linear FSI problems, the equations that result after finite element discretization of parameter-dependent nonlinear FSI problems are not translatable to a single matrix equation. Such a discretization with respect to $m \in N$ shear moduli given by the set $\{\mu_i\}_{i \in \{1,...,m\}} \subset R$ and a number of $n \in N$ degrees of freedom yields equations of the form
$(A_0+\mu_iA_1+A_2(x_i) )x_i=b \quad \text{for} \quad i \in \{1,...,m\}$,
with $A_0$, $A_1$, $A_2(x_i) \in R^{n \times n}$, $b \in R^n$ and the finite element solution $x_i \in R^n$. $A_2(\cdot)$ is a discretization matrix of nonlinear operators and depends on the unknown. An equivalent translation of the $m$ equations above to a matrix equation, similar to the linear case, is not possible anymore. The poster will, in addition, illustrate how such parameter-dependent nonlinear FSI discretizations can be tackled on the basis of the Newton iteration.
We deduce a procedure to use balanced truncation for parameter-dependent differential-algebraic systems. For this, we have to solve parameter-dependent Lyapunov equations. Since solving large-scale Lyapunov equations for every parameter leads to very high computational costs we utilize the reduced-basis method to get a reduced Lyapunov equation that we can solve instead. In order to use this method, we first make the matrix pencil of the system strictly dissipative to apply error estimators in the reduced basis method. Additionally, we deal with the algebraic parts of the system by using projections that lead to ordinary differential equations.
Matrix equations and their solutions play a key role in several problems arising in science and engineering. Our primary focus lies on matrix equations appearing in model-order reduction problems, where several well-known procedures, such as balanced truncation and interpolation-based methods, can be set up in the matrix equations form. In the last five years, several model reduction schemes for various classes of systems have been proposed. These classes include quadratic-bilinear, polynomial, switched systems and quadratic output systems. The newly proposed MOR schemes for these classes require to solve particular types of matrix equations. In this poster, we provide an overview of the MOR methods and the newly arising matrix equations. Furthermore, we discuss the current methods to solve these matrix equations and the challenges in a large-scale setting.
In system theory, the so-called system Gramian matrices are operators encoding certain properties of an underlying input-output system. Usually, these system Gramians are computed as solutions to matrix equations, such as the Lyapunov equation and Sylvester equation. This means, the solution to certain matrix equations coincides with these system Gramians. Now, if the system Gramians are computable by other means than matrix equations, they still represent solutions to matrix equations. Empirical Gramians are such an alternative for system Gramian computation, which are based on their system-theoretic operator definition, and practically obtained via quadrature.
This contribution explores the connection between matrix equations, system Gramians and empirical Gramians, and proposes empirical Gramians as potential solver for matrix equations.
The MORLAB, Model Order Reduction LABoratoy, toolbox [1] is a free and open source software solution in MATLAB and GNU Octave for linear model reduction problems. Providing system theoretic model order reduction methods, the basis of the highly efficient implementation is built by a large variety of dense matrix equation solvers based on iterative methods like the matrix sign function, Newton or squared Smith iterations, and many more. Therefore, the toolbox can be used to solve continuous- or discrete-time Lyapunov, Sylvester and Riccati equations, as well as continuous-time Bernoulli equations in a lot of different formats and formulations.
References
[1] P. Benner and S. W. R. Werner.
MORLAB -- Model Order Reduction LABoratory (version 5.0), 2019.
See also http://www.mpi-magdeburg.mpg.de/projects/morlab.
doi:10.5281/zenodo.3332716
We consider $\mathcal{H}_2 \otimes \mathcal{L}_2$-optimal model order reduction of parametric linear time-invariant dynamical systems, where coupled matrix integral equations arise in the first-order necessary conditions (FONC). The quality of the reduced-order model is measured using the $\mathcal{H}_2$ norm for the parametric system, which is averaged in the $\mathcal{L}_2$-norm over the parameter domain.
We motivate the FONC and aim to satisfy them using an optimization-based approach that involves solving sequences of small-scale Lyapunov equations and tall and skinny Sylvester equations.
Dynamical low-rank approximation is introduced and a numerical integrator that computes a symmetric or skew-symmetric low-rank approximation to large symmetric or skew-symmetric time-dependent matrices that are either given explicitly or are the unknown solution to a matrix differential equation is presented. We show that low-rank time-dependent matrices are reproduced exactly, and the error behaviour is robust to the presence of small singular values of the solution or the approximate solution. The proposed structure preserving integrator and his properties are then extended to (anti-)symmetric low-rank Tucker tensors.
Port-Hamiltonian (pH) systems are a very important modeling structure for large classes of technical systems. In the optimal control problem for pH systems, the resulting boundary value problems and Lure/Riccati/Lyapunov have a double structure.
It is a big challenge to make use of this double structure inan effective way, but also in the solution of the associated structured eigenvalue problems. We will present the theory, first approaches, and major open questions, in particular for the case of pH descritor systems.
Evaluating the action of a matrix function on a vector, that is $x=f(\mathcal M)v$, is an ubiquitous task in applications. When the matrix $\mathcal M$ is large, subspace projection method, such as the rational Krylov method, are usually employed.
In this work, we provide a quasi-optimal pole choice for rational Krylov methods applied to this task when $f(z)$ is either Cauchy-Stieltjes or Laplace-Stieltjes (or, which is equivalent, completely monotonic) for a positive definite matrix $\mathcal M$.
Then, we consider the case when the argument $\mathcal M$ has the Kronecker structure $\mathcal M=I \otimes A - B^T \otimes I$, and is applied to a vector obtained vectorizing a low-rank matrix. This finds application, for instance, in solving fractional diffusion equation on rectangular domains.
We introduce an error analysis for the numerical approximation of $x$. Pole choices and explicit convergence bounds are given also in this case.
We describe a way to implement the matrix sign iteration $H_{k+1} = \frac12 (H_k^{\mathstrut}+H_k^{-1})$ on a dense Hamiltonian matrix of the form
$$H_k = \begin{bmatrix}A_k & B_kB_k^T\\ C_k^TC_k & -A_k^T & \end{bmatrix}$$ in such a way that the blocks in positions $(1,2)$ and $(2,1)$ are kept in low-rank factored form. The algorithm operates on their generators $B_k$ and $C_k$ directly, and relies on principal pivot transforms (PPTs) as its main building block; more specifically it makes use of the framework for factored-form PPTs in [Poloni, Strabic 2016]. We discuss the stability properties of the resulting algorithm, as well as applications that make use of the low-rank factors directly.