Cicon Structured Matrices in Numerical Linear Algebra. Cortona 2017.
Topics of the workshop are:
  • Structured matrix analysis including (but not limited to) Toeplitz, Hankel, Vandermonde, banded, semiseparable, Cauchy, Hessenberg, mosaic, block, multilevel matrices and the theoretical and applicative problems from which they are originated (structured problems);
  • Applications involving structured matrices including (but not limited to) interpolation, integral and differential equations, least squares and regularization, polynomial computations, matrix equations, control theory, queueing theory and Markov chains, image and signal processing;
  • Design and analysis of algorithms for the solution of structured problems.
The aims of the meeting are:
  • presenting recent results on theory, algorithms and applications concerning structured problems in numerical linear algebra and matrix theory;
  • reviewing and discussing methodologies and the related algorithmic analysis;
  • improving collaborations between theoretical and applied research;
  • tracing the state-of-the-art together with the main directions of future research;
  • fostering the contacts between PhD students, postdocs and young scholars with the most advanced research groups;
  • increasing the collaboration between European research groups.
The workshop will be held in Cortona, Italy, at Il Palazzone, a XVI century monumental palace, also known as Villa Passerini, at walking distance from the center (see map below).

Cortona can be reached:
  • by plane. Landing to the airports of Rome (ROM), Florence (FLR) or Perugia (PEG) and continuing the travel by train.
  • by train. Stopping to one of the two railway stations: Terontola-Cortona (about 15 km from the village) and Camucia-Cortona (about 5 km). Duration of the trip: from Rome 2h30'; from Florence, 1h30'; from Perugia, 0h35'; from Milan, 4h30'. Check timetables.
    The transportation from the railway station to Cortona will be organized locally.
  • by car. Get directions from Rome Fiumicino, from Firenze, from Perugia.

Lidia Aceto, Università di Pisa, Italy
Alessandra Aimi, Università di Parma, Italy
Francesca Arrigo, University of Strathclyde, Glasgow, UK
Giovanni Barbarino, Scuola Normale Superiore, Pisa, Italy
Michele Benzi, Emory University, Atlanta, USA
Luca Bergamaschi, Università di Padova, Italy
Daniele Bertaccini, Università di Roma Tor Vergata, Italy
Davide Bianchi, Università dell'Insubria, Como, Italy
Dario A. Bini, Università di Pisa, Italy
Paola Boito, Université de Limoges and LIP-ENS de Lyon, France
Matthias Bolten, University of Wuppertal, Germany
Claude Brezinski, Université de Lille, France
Raymond Chan, The Chinese University of Hong Kong, Hong Kong
Antonio Cicone, Università dell'Aquila, Italy
Stefano Cipolla, Università di Roma Tor Vergata, Italy
Anna Concas, Università di Cagliari, Italy
Francisco De Terán, Universidad Carlos III de Madrid, Spain
Gianna M. Del Corso, Università di Pisa, Italy
Pietro Dell'Acqua, Università dell'Aquila, Italy
Fabio Di Benedetto, Università di Genova, Italy
Carmine Di Fiore, Università di Roma Tor Vergata, Italy
Marco Donatelli, Università dell'Insubria, Como, Italy
Froilán M. Dopico, Universidad Carlos III de Madrid, Spain
Fabio Durastante, Università di Roma Tor Vergata, Italy
Yuli Eidelman, Tel Aviv University, Israel
Claudio Estatico, Università di Genova, Italy
Massimiliano Fasi, The University of Manchester, UK
Dario Fasino, Università di Udine, Italy
Caterina Fenu, Università di Cagliari, Italy
Carlo Garoni, Università della Svizzera italiana, Lugano, Switzerland
Luca Gemignani, Università di Pisa, Italy
Nicola Guglielmi, Università dell'Aquila, Italy
Thomas Huckle, Technische Universität München, Germany
Bruno Iannazzo, Università di Perugia, Italy
Carlo Janna, Università di Padova, Italy
Thomas Kailath, Stanford University, USA
Daniel Kressner, EPFL, Lausanne, Switzerland
Carla Manni, Università di Roma Tor Vergata, Italy
Stefano Massei, EPFL, Lausanne, Switzerland
Nicola Mastronardi, CNR Bari, Italy
Mariarosa Mazza, Max Planck Institute, München, Germany
Beatrice Meini, Università di Pisa, Italy
Marilena Mitrouli, National and Kapodistrian University of Athens, Greece
Vanni Noferini, University of Essex, UK
Silvia Noschese, Università di Roma La Sapienza, Italy
Dimitrios Noutsos, University of Ioannina, Greece
Davide Palitta, Università di Bologna, Italy
Victor Y. Pan, City University of New York, USA
Bor Plestenjak, University of Ljubljana, Slovenia
Federico Poloni, Università di Pisa, Italy
Daniel Potts, University of Chemnitz, Germany
Stefano Pozza, Charles University, Prague, Czech Republic
Michela Redivo Zaglia, Università di Padova, Italy
Lothar Reichel, Kent State University, Kent, USA
Leonardo Robol, ISTI-CNR, Pisa, Italy
Paraskevi Roupa, University of Athens Panepistimiopolis, Greece
Stefano Serra Capizzano, Università dell'Insubria, Como, Italy
Debora Sesana, Università dell'Insubria, Como, Italy
Valeria Simoncini, Università di Bologna, Italy
Hendrik Speleers, Università di Roma Tor Vergata, Italy
Daniel Szyld, Temple University, Philadelphia, USA
Francoise Tisseur, University of Manchester, UK
Francesco Tudisco, University of Strathclyde, Glasgow, UK
Eugene Tyrtyshnikov, Russian Academy of Sciences, Russia
Marc Van Barel, KU Leuven, Leuven, Belgium
Cornelis Van Der Mee, Università di Cagliari, Italy
Raf Vandebril, KU Leuven, Leuven, Belgium
Paris Vassalos, Athens University of Economics and Business, Greece
Joab Winkler, University of Sheffield, UK
The book of abstracts is avalaible here

Move the mouse over the names to read titles and abstracts of the talks and posters. Use the arrow keys for long abstracts.

Monday Tuesday Wednesday Tuesday Friday
08.45-09.00 Opening
Chan Fractional-order Total Variation Model for Image Restoration
Fractional-order derivative is attracting interest from researchers working in image processing because it helps to preserve more texture than total variation when noise is removed. There are several definitions of fractional-order derivative, but usually the Grunwald-Letnikov fractional-order derivative is used. To preserve structure, one will consider the Direchlet homogeneous boundary condition. The corresponding discrete partial fractional-order derivative operator is a full lower triangular Toeplitz matrix. This increases the complexity in solving the imaging models. In this talk, we will discuss ways to modify the fractional-order partial derivative operator to speed up the computation.
Research supported by HKRGC, XJTLU and JSNSF
Benzi Iterative Methods for Linear Systems with Double Saddle Point Structure
We consider several iterative methods for solving a class of linear systems with double saddle point structure. Both Uzawa-type stationary methods and block preconditioned Krylov subspace methods are discussed. We present convergence results and eigenvalue bounds together with illustrative numerical experiments using test problems from two different applications: a mixed-hybrid discretization of the potential fluid flow problem, and finite element modeling of liquid crystal directors.
Reichel Some structured matrix problems in Gauss-type quadrature
It is well-known that Gauss quadrature rules can be computed by evaluating the integrand at a symmetric tridiagonal matrix. This talk describes several generalizations that can be applied to approximate certain matrix functions and estimate the error in the approximation so determined. These generalizations involve symmetric and nonsymmetric tridiagonal and block tridiagonal matrices.
Tisseur The Structured Condition Number of a Differentiable Map Between Matrix Manifolds
We study the structured condition number of differentiable maps between smooth matrix manifolds and present algorithms to compute it. We analyze automorphism groups, and Lie and Jordan algebras associated with a scalar product as a special case of smooth matrix manifolds (these include orthogonal matrices, symplectic matrices, Hamiltonian matrices, ...). For such manifolds, we derive a lower bound on the structured condition number that is cheaper to compute than the structured condition number. We provide numerical comparisons between the structured and unstructured condition numbers for the principal matrix logarithm and principal matrix square root of matrices in automorphism groups as well as for maps between matrices in automorphism groups and their structured polar decomposition and structured matrix sign decomposition. We show that our lower bound can be used as a good estimate for the structured condition number when the matrix argument is well-conditioned. We show that the structured and unstructured condition numbers can differ by several orders of magnitude, thus motivating the development of algorithms preserving structure.
This is joint work with Bahar Arslan and Vanni Noferini.
Tyrtyshnikov Multidimensional matrices, optimization and quadratic kernels
We discuss the links between multidimensional matrices and some optimization problems, especially the so called singular problems where the classical Newton method cannot be applied. In particular, we show how some algorithms with quadratic convergence can be constructed in the singular case and how the constructions are related with the triviality of the quadratic kernel of some matrix subspaces.
Estatico Numerical linear algebra and regularization in Banach spaces
We consider the functional linear equation $Af=g$ arising in image restoration, characterized by a structured and ill-posed linear operator $A:X \longrightarrow Y$ between two Banach spaces $X$ and $Y$. In this talk, the discretized linear system is solved in the framework of the regularization theory in Banach spaces, where the variational approach based on non-linear iterative minimization of the residual in the dual spaces $X^*$ and $Y^*$ is used. Within this framework, some relationships between classical algorithms of numerical linear algebra and the more recent non-linear iterative regularization algorithms in Banach spaces are analyzed. In particular, the well-known class of iterative projection algorithms (including Cimmino, SIRT and DROP methods), as well as the class of (conjugate) gradient methods will be investigated within the context of regularization in Banach spaces.
Tudisco Modularity matrices and community detection under the degree-corrected stochastic block model
A community in a network $G$ is roughly defined as a set of nodes being highly interconnected and poorly linked with the rest of $G$. Conversely, an anti-community is a node set being loosely connected internally but having many external connections. We propose a spectral method based on the extremal eigenspaces of generalized modularity matrices to simultaneously look for communities and anti-communities within a given graph or network [2]. We provide theoretical and experimental evidence in support of the proposed method. Moreover, we discuss the behavior of the proposed strategy under the degree-corrected stochastic block model (DC-SBM) [1]. The latter is a random graph model which generates graphs with prescribed clustering structure and heterogeneous degree distributions. We present a complete spectral analysis of the expected modularity matrix $M$ of a DC-SBM random graph. In particular, we show that $M$ can be written as the Khatri-Rao product of a small matrix associated to the clustering structure of the model and a rank-one, block structured matrix related to the degree distribution. [1] K. Chaudhuri, F. Chung, A. Tsiatas. Spectral clustering of graphs with general degrees in the extended planted partition model. COLT 2012, Journal of Machine Learning Research, 23 (2012), 35.1-35.23, [2] D. Fasino, F. Tudisco. A modularity based spectral method for simultaneous community and anti-community detection. Linear Algebra Appl., 2017, to appear.
Guglielmi On some matrix nearness problems in matrix theory
Matrix nearness problems arise naturally in matrix theory when measuring robustness of some property of a given matrix or when aiming to correct a matrix which does not present a desired property. We discuss ODE-based methods to deal with this kind of problems with application to stability, controllability, common divisibility of polynomials and other properties. In all cases we are particularly interested to the structure preservation of the matrix, which is usually a relevant issue. Numerical examples will be shown in order to illustrate the behavior of the proposed methodology. The talk is inspired by joint works with Christian Lubich (Universitaet Tuebingen).
Noferini Eigenvalue avoidance for structured matrices depending on a real parameter
Let $H(t)$ be a Hermitian matrix that depends continuously on the real parameter $t$. If the eigenvalues of $H(t)$ are plotted against $t$, it can typically be observed that the trajectories of the eigenvalue functions appear to collide; however, they undergo a last-minute repulsive effect, thus avoiding intersections. Although it is of course possible to construct examples that do not follow this paradigm, this phenomenon happens with probability 1 in the set of $n\times n$ Hermitian matrices whose entries are continuous functions. Eigenvalue avoidance of self-adjoint operators was first observed experimentally by quantum physicists in the early XX century. In 1929, J. von Neumann and E. P. Wigner were the first to explain eigenvalue avoidance rigorously, proving via geometric arguments that it is generic. I plan to first briefly expose the von Neumann-Wigner approach. Then, I will present a number of new results, extending the theory to various classes of structured matrices. I would like to discuss for which classes eigenvalue avoidance occurs, for which classes it does not, and for which classes the question is still open (to the date in which this abstract is being written).
Vandebril Eigenvalues and eigenvectors of matrix polynomials
In the last decade matrix polynomials have been investigated with the primary focus on adequate linearizations and good scaling techniques for computing their eigenvalues and eigenvectors. We propose a new method for computing a factored Schur form of the associated companion pencil. The algorithm has a quadratic cost in the degree of the polynomial and a cubic one in the size of the coefficient matrices. Also the eigenvectors can be computed at the same cost. The algorithm is a variant of Francis's implicitly shifted QR algorithm applied on the companion pencil. A preprocessing unitary equivalence is executed on the matrix polynomial to simultaneously bring the leading matrix coefficient and the constant matrix term to triangular form before forming the companion pencil. The resulting structure allows us to stably factor each matrix of the pencil as a product of $k$ matrices of unitary-plus-rank-one form, admitting cheap and numerically reliable storage. The problem is then solved as a product core chasing eigenvalue problem. The algorithm is normwise backward stable. Computing the eigenvectors via reordering the Schur form is discussed as well.
Bianchi Nonstationary preconditioned iterative methods with structure preserving property for image deblurring
Image deblurring is the process of reconstructing an approximation of an image from blurred and noisy measurements. By assuming that the point spread function (PSF) is spatially-invariant, the observed image $g(x, y)$ is related to the true image $f (x, y)$ via the integral equation \[ g(x,y)=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty} h(x-x',y-y')f(x',y')dx'dy'+\eta(x,y)\qquad\qquad\qquad (1) \] for every $(x,y)\in\Omega\subset \mathbb R^2$ and where $\eta(x,y)$ is the noise. By collocation of the previous integral equation on a uniform grid, we obtain the grayscale images of the observed image, of the true image, and of the PSF, denoted by $G$, $F$, and $H$, respectively. Since collected images are available only in a finite region, the field of view (FOV), the measured intensities near the boundary are affected by data outside the FOV. Given an $n \times n$ observed image $G$ (for the sake of simplicity we assume square images), and a $p \times p$ PSF with $p \le n$, then $F$ is $m \times m$ with $m = n + p - 1$. Denoting by $g$ and $f$ the stack ordered vectors corresponding to $G$ and $F$, the discretization of (1) leads to the under-determined linear system \[ \boldmath g=A\boldmath f+\boldmath\eta, \qquad\qquad\qquad\qquad\qquad\qquad \qquad\qquad\qquad (2) \] where the matrix $A$ is of size $n^2 \times m^2$. When Imposing proper Boundary Conditions (BCs), the image $A$ becomes square $n^2 \times n^2$ and in some cases, depending on the BCs and the symmetry of the PSF, it can be diagonalized by discrete trigonometric transforms. For example, the matrix $A$ is block circulant circulant block (BCCB) and it is diagonalizable by Discrete Fourier Transform (DFT), when periodic BCs are imposed. Due to the ill-posedness of (1), A is severely ill-conditioned and may be singular. In such case, linear systems of equations (2) are commonly referred to as linear discrete ill-posed problems. Therefore a good approximation of f cannot be obtained from the algebraic solution (e.g., the least-square solution) of (2), but regularization methods are required. The basic idea of regularization is to replace the original ill-conditioned problem with a nearby well-conditioned problem, whose solution approximates the true solution. Thresholding iterative methods are recently successfully applied to image deblurring problems. We investigate the modified linearized Bregman algorithm (MLBA) used in image deblurring problems, with a proper treatment of the boundary artifacts. The fast convergence of the MLBA depends on a regularizing pre-conditioner that could be computationally expensive and hence it is usually chosen as a block circulant circulant block (BCCB) matrix, diagonalized by discrete Fourier transform. We show that the standard approach based on the BCCB preconditioner may provide low quality restored images. Indeed, in order to get an effective preconditioner it is crucial to preserve the structure of the coefficient matrix that depends on the BCs, in the image deblurring problems. Motivated by a recent nonstationary preconditioned iteration [4], we propose a new algorithm that combines such method with the MLBA and it is structure preserving.
[1] Y. Cai, M. Donatelli, D. Bianchi and T.Z. Huang, Regularization preconditioners for frame-based image deblurring with reduced boundary artifacts, SIAM J. Sci. Comput., 38(1) (2016) 164--189.
[2] J. F. Cai, S. Osher, and Z. Shen, Linearized Bregman iterations for frame- based image deblurring, SIAM J. Imaging Sci., 2--1 (2009), pp. 226--252.
[3] P. Dell'Acqua, M. Donatelli, C. Estatico and M. Mazza, Structure preserving preconditioners for image deblurring, J. Sci. Comput., 72(1) (2017) 147--171.
[4] M. Donatelli and M. Hanke, Fast nonstationary preconditioned iterative methods for ill-posed problems, with application to image deblurring, Inverse Problems, 29 (2013) 095008.
Arrigo Non-backtracking walk centrality measures for networks
The task of community detection and node centrality measurement can both be addressed by quantifying traversals around a network: either random or deterministic. In particular, the concept of a walk, which allows nodes and edges to be revisited, forms the basis of the Katz centrality and the total communicability, whose limiting behaviour corresponds to eigenvector centrality. In this talk we describe how to extend certain walk-based centrality measures to their non-backtracking analogue, where walks that include at least one back-and-forth flip between a pair of nodes are not allowed. We further show that these new measures may be described in terms of standard walks taking place in certain multilayer networks, and can thus be computed by working on 3x3 or 2x2 block matrices.
Bergamaschi Spectral Low-rank Preconditioners for Large Linear Systems and Eigenvalue Problems
Fast solution of large and sparse SPD linear systems by Krylov subspace methods is usually prevented by the presence of eigenvalues near zero of the coefficient matrix $A$. This is particularly true when computing the smallest or interior eigenvalues where, using Lanczos' or the Jacobi-Davidson approach, a system like $(A - \sigma I)x = r$ has to be repeatedly solved, with $\sigma$ close to the wanted eigenvalue.
We propose and discuss how cost-effective spectral information on the coefficient matrix can be used to construct a spectral preconditioner, i.e., a low-rank modification of a given approximate inverse of $A$, $P_0 \approx A^{-1}$. The spectral preconditioner usually moves away from zero the smallest eigenvalues of the preconditioned matrices with a consequent, sometimes dramatic, reducing of the condition number and speeding up of the iterative process.
Given a set of very roughly approximated smallest eigenvalues of $A:$ $\lambda_1 , \ldots , \lambda_p$ and corresponding approximate eigenvectors $v_1 , \ldots , v_p$, and defining the matrices $\Lambda_p =\hbox{diag}(\lambda_1 , \ldots, \lambda_p)$ and $V_p = [v_1 , . . . , v_p ]$ we investigate the properties of two classes of preconditioners ([3, 4]), among the others, defined as \[\begin{split} &P= P_0 + V_p \Lambda_p^{-1}V_p^T \\ &P= V_p\Lambda_p^{-1}V_p^T + (I - V_p V_p^T )P_0 (I - V_p V_p T ) \end{split} \]
We will provide an extensive testing of such preconditioners which will be shown to provide an important acceleration of iterative eigensolvers (see e.g. [2]) as well as of very ill-conditioned sequences of linear systems [1].
[1] L. Bergamaschi, E. Facca, A. Martínez, and M. Putti, Spectral preconditioners for the efficient numerical solution of a continuous branched transport model, Proc. Cedya + CMA Conf. Cartagena, Spain, 2017.
[2] L. Bergamaschi and A. Martínez, Two-stage spectral preconditioners for iterative eigensolvers, Numer. Lin. Alg. Appl., 24 (2017), pp. 1--14.
[3] B. Carpentieri, I. S. Duff, and L. Giraud, A class of spectral two-level preconditioners, SIAM J. Sci. Comput., 25 (2003), pp. 749--765 (electronic).
[4] A. Martínez, Tuned preconditioners for iterative SPD eigensolvers, Numer. Lin. Alg. Appl., 23 (2016), 427--443.
Noschese Computing Structured Pseudospectrum Approximations
In many applications it is important to know location and sensitivity to perturbations of eigenvalues of matrices or polynomial matrices. The sensitivity commonly is described with the aid of condition numbers or by computing pseudospectra; see, e.g., [1, 2]. However, the computation of pseudospectra is very demanding computationally. We propose a new approach to computing approximations of pseudospectra of both matrices [3] and matrix polynomials [4] by using rank-one or projected rank-one perturbations that respect the given structure, if any. These perturbations are inspired by Wilkinson's analysis of eigenvalue sensitivity [5]. Numerical examples show this approach to perform much better than methods based on random rank-one perturbations both for the approximation of structured and unstructured (polynomial) pseudospectra.
[1] L. N. Trefethen and M. Embree, Spectra and Pseudospectra. Princeton University Press, Princeton, 2005.
[2] N. J. Higham and F. Tisseur, More on pseudospectra for polynomial eigenvalue problems and applications in control theory. Linear Algebra Appl., 351 (2002), pp. 435-453.
[3] S. Noschese and L. Reichel, Approximated structured pseudospectra. Numer. Linear Algebra Appl., 24 (2017) e2082.
[4] S. Noschese and L. Reichel, Computing Unstructured and Structured Polynomial Pseudospectrum Approximations. Submitted.
[5] J. H. Wilkinson, The Algebraic Eigenvalue Problem. Oxford University Press, 1965.
Robol Backward error analysis for core-chasing algorithms
Fast methods that exploit the structure of the underlying problem are typically attractive because of their performances and, sometimes, because of the reduced storage cost. However, the improved performance can come at the cost of stability, and it is often challenging to prove satisfying backward stability results.
In this talk, we focus on core-chasing algorithms, that allow the efficient development of methods that exploit unitary plus rank $1$ or unitary plus rank $k$ structure. In particular, we analyze the case of polynomial rootfinding, that relies on the computation of the eigenvalues of a matrix with unitary plus rank $1$ structure (or a similarly structured pencil), and we will show that this approach can be backward stable on the polynomial coefficients.
We discuss the differences between QZ and QR, and we prove that -- using structured core-chasing methods -- the backward error is also structured. This leads to interesting results, and shows that the QR iteration is backward stable on the polynomial coefficients in unexpected situations.
10.30-11.00 Coffee Break Coffee Break
Coffee Break
Coffee Break
Coffee Break
Aimi On the Energetic Galerkin BEM and its algebraic reformulation
The Energetic Galerkin Boundary Element Method (BEM) is a discretization technique for the numerical solution of wave propagation problems, introduced in [1] and applied in the last decade to scalar wave propagation inside bounded domains or outside bounded obstacles, in 1D, 2D and 3D space dimension.
The differential initial-boundary value problem at hand is converted into a space-time Boundary Integral Equation (BIE), which is then written in weak form through energy considerations and discretized by a Galerkin approach.
The talk will focus on the extension of this method in the context of 2D soft and hard scattering of damped waves, taking into account both viscous and material damping.
Details will be given on the algebraic reformulation of Energetic Galerkin BEM, i.e., on the so-called time-marching procedure that gives rise to a linear system whose matrix has a Toeplitz lower triangular block structure. Numerical results confirm accuracy and stability of the proposed technique, already proved for the numerical treatment of undamped wave propagation problems in several space dimensions [2-3] and for the 1D damped case [4-5].
[1] Aimi, A. and Diligenti, M., A new space-time energetic formulation for wave propagation analysis in layered media by BEMs, Int. J. Numer. Meth. Engrg. 75, 1102--1132 (2008).
[2] Aimi, A., Diligenti M. and Panizzi S., Energetic Galerkin BEM for wave propagation Neumann exterior problems, Computer Model. Engrg. Sciences 1(1), 1--33 (2009).
[3] Aimi, A., Diligenti, M., Frangi, A. and Guardasoni, C., Neumann exterior wave propagation problems: Computational aspects of 3D energetic Galerkin BEM, Comput. Mech. 51(4), 475--493 (2013).
[4] Aimi, A. and Panizzi, S., BEM-FEM coupling for the 1D Klein-Gordon equation, Numer. Methods Partial Diff. Equations 30(6), 2042--2082 (2014).
[5] Aimi, A., Diligenti, M. and Guardasoni, C., Energetic BEM-FEM coupling for the numerical solution of the damped wave equation, Adv. Comput. Math. 43, 627-651 (2017).
Van Barel Matrices in polynomial system solving
This talk is about the role of matrix computations in the problem of solving systems of polynomial equations. Let $k$ be an algebraically closed field and let $p_1 = p_2 = ... = p_n = 0$ define such a system in $k^n$: $p_i \in k[x_1,...,x_n]$. Let $I$ be the ideal generated by these polynomials. We are interested in the case where the system has finitely many isolated solutions in $k^n$. It is a well known fact that this happens if and only if the quotient ring $k[x_1,..., x_n]/I$ is finite dimensional as a k-vector space. The multiplication endomorphisms of the quotient algebra provide a natural linear algebra formulation of the root finding problem. Namely, the eigenstructure of the multiplication matrices reveals the solutions of the system. These multiplication matrices can be calculated from the coefficients of the $p_i$, for example by using Groebner bases. The computations make an implicit choice of basis for $k[x_1,...,x_n]/I$, which from a numerical point of view is not a very good choice. Significant improvement can be made by using elementary numerical linear algebra techniques on a Macaulay-type matrix. In this talk we will present thi technique and show how the resulting method can handle challenging systems.
Szyld Singular values of certain almost block Toeplitz matrices
We analyze a certain class of "almost" block Toeplitz matrices arising in the analysis of Optimized Schwarz methods for the numerical solution of certain PDEs. The aim is to show that these matrices are contracting.
Mastronardi Revisiting the perfect shift strategy in the Implicitly Shifted QR algorithm
In this talk we revisit the Implicit-Q Theorem and analyze the problem of performing a QR-step on an unreduced Hessenberg matrix $H$ when we know an "exact" eigenvalue $\lambda_0$ of $H$. Under exact arithmetic, this eigenvalue will appear on diagonal of the transformed Hessenberg matrix $H_1$ and will be decoupled from the remaining part of the Hessenberg matrix, thus resulting in a deflation. But it is well known that in finite precision arithmetic the so-called perfect shift could get blurred and the eigenvalue $\lambda_0$ can not be deflated and/or is perturbed significantly. In this talk we develop a new strategy for computing such a QR step so that the deflation is indeed successful. The method is based on the preliminary computation of the corresponding eigenvector $x$ such that the residual $(H-\lambda_0)/x$ is sufficiently small. The eigenvector is then transformed to a unit vector by a sequence of Givens transformations, which are also performed on the Hessenberg matrix. Such a QR step is the basic ingredient of the QR method to compute the Schur form, and hence the eigenvalues of an arbitrary matrix. But it also is a crucial step in the reduction of a general matrix $A$ to its Weyr form. It is in fact this last problem that lead to the development of this new technique.
Brezinski Shanks sequence transformations and their links to linear algebra methods
In this talk present a general framework for Shanks transformations of sequences of elements in a vector space. It is shown that the Minimal Polynomial Extrapolation (MPE), the Modified Minimal Polynomial Extrapolation (MMPE), the Reduced Rank Extrapolation (RRE), the Vector Epsilon Algorithm (VEA), the Topological Epsilon Algorithm (TEA), and Anderson Acceleration (AA), which are standard general techniques designed for accelerating arbitrary sequences and/or solving systems of linear and nonlinear equations, all fall into this framework. Their properties and their connections with QuasiNewton and Broyden methods are studied. Then, we exploit this framework to compare these methods. In the linear case, it is known that AA and GMRES are 'essentially' equivalent in a certain sense while GMRES and RRE are mathematically equivalent. The talk also discusses the connection between AA, the RRE, the MPE, and other methods in the nonlinear case.
Mazza Spectral analysis and spectral symbol for pure and stabilized 2D curl-curl operator with applications to the related iterative solutions
In this work, we focus on large and highly ill-conditioned linear systems of equations arising from various formulations of the Maxwell equations appearing, e.g., in Time Harmonic Maxwell as well as in the MagnetoHydroDynamics.
First, we consider a compatible B-Spline discretization based on a discrete De Rham sequence of the 2D curl-curl operator stabilized with zero-order term, and we show that the sequence of the coefficient matrices belongs to the Generalized Locally Toeplitz class. Moreover, looking at the entries of the coefficient matrix, we compute the symbol describing its asymptotic eigenvalue distribution, as the matrix size diverges. Thanks to this spectral information we show that the coefficient matrix is affected by three severe sources of ill-conditioning related to the relevant parameters: the matrix size, the spline degree and the stabilization parameter. As a consequence, when used for solving the associated linear systems, Conjugate Gradient type methods are extremely slow and their convergence rate is not robust with respect to the parameters.
On this basis, we replace the zero-order stabilization with a divergence-type one and we spectrally analyze the corresponding B-spline discretization matrix-sequence. The retrieved spectral information is then used to design a 2D vector extension of a multi-iterative approach already used in the literature for the scalar Laplacian operator. Finally, a variety of numerical tests and some open problems are discussed.
Plestenjak Minimal determinantal representations of bivariate polynomials
It is known since Dixon's 1902 paper that every bivariate polynomial $p$ of degree $n$ admits a determinantal representation $p(x,y)=\det(A+xB+yC)$ with $n\times n$ symmetric matrices.
However, the construction of such matrices is far from trivial and up to now there have been no efficient numerical algorithms, even if we do not insist on matrices being symmetric.
We present the first numerical construction that returns a determinantal representation with $n\times n$ matrices for a square-free bivariate polynomial of degree $n$, which, with the exception of the symmetry, agrees with Dixon's result.
For a non square-free polynomial one can combine it with a square-free factorization to obtain a representation of order $n$.
Our motivation is a numerical method for solving systems of bivariate polynomials as two-parameter eigenvalue problems. Symmetry is not important for this particular application.
The resulting numerical method for the roots of a system of two bivariate polynomials is competitive with some existing methods for polynomials of small degree.
Huckle Preconditioning for sparse and structured matrices
We consider linear systems of equations related to sparse matrices $A$. Incomplete LU decomposition usually leads to a good preconditioner. Here, we will present iterative methods for computing ILU and MILU that are also easy to parallelize. The convergence of these iterative methods depends strongly on the condition number of $A$. Furthermore, the condition number of L and U should be large to obtain an efficient preconditioner. So fast converging iterative ILU and MILU methods are necessary in this case. For solving the resulting sparse triangular systems we use the Jacobi method with an incomplete sparse approximate inverse preconditioner. Furthermore, the Jacobi iteration can be accelerated by using the Euler expansion. In view of the triangular structure of L and U, the accelerated Jacobi iteration is guaranteed to converge fast. Overall, this results in efficient, parallel methods for solving sparse linear systems. As special case we discuss the presented methods for structured matrices via the symbol.
Noutsos Extensions of M-matrices and their properties concerning the Schur Complement
It is well known that $M-$matrices are the square matrices of the form $A=sI-B$, where $B$ is entry-wise nonnegative and $s \geq \rho(B)$. This class of matrices has many applications in Ordinary and Partial Differential Equations, Integral equations, Economics, Linear complementarity problems, Population Dynamics, Markov chains, Theory of Games and Control Theory. Some extended classes of $M-$matrices, with many applications in the same areas of Science, have been proposed in the last decades. The class of ${M_v}-$matrices contains the matrices of the form $A=sI-B$, where $B$ is eventually nonnegative and $s \geq \rho(B) $. The class of Generalized $M-$matrices ($GM-$matrices) is defined following the same reasoning except that the matrix $B$ has the Perron-Frobenius property. Very useful properties concerning $M-$matrices, Inverse $M-$matrices and their Schur complements are known. In this work we discuss extensions of such properties concerning ${M_v}-$ and $GM-$matrices, Inverse ${M_v}-$ and Inverse $GM-$matrices with respect to their Schur complements. Numerical examples are shown to confirm the obtained theoretical results.
Redivo-Zaglia Shanks sequence transformations and the $\varepsilon$-algorithms. Theory and applications.
Let $(\mathbf S_n)$ be a sequence of elements of a vector space $E$ on a field $\mathbb K$ ($\mathbb R$ or $\mathbb C$) which converges to a limit $\mathbf S$. If the convergence is slow, it can be transformed, by a {\it sequence transformation}, into a new sequence or a set of new sequences which, under some assumptions, converges faster to the same limit. When $E$ is $\mathbb R$ or $\mathbb C$, a well known such transformation is due to Shanks, and it can be implemented by the scalar $\varepsilon$--algorithm of Wynn.
This transformation was generalized in several different ways to sequences of elements of a vector space $E$. When $E$ is $\mathbb R^p$ or $\mathbb C^p$, this generalization leads to the {\it Reduced Rank Extrapolation} (RRE) and to the {\it Minimal Polynomial Extrapolation} (MPE). For a general vector space $E$, the {\it Modified Minimal Polynomial Extrapolation} (MMPE) and the {\it topological Shanks transformations} are obtained. The interest of these last two generalizations is that they can treat sequences or matrices or even tensors, and that they can be recursively implemented, the first one by the $S\beta$ algorithm of Jbilou and the second ones by the topological $\varepsilon$--algorithms of Brezinski [1].
However, the topological $\varepsilon$--algorithms are quite complicated since they possess two rules, they require the storage of many elements of $E$, and the duality product with an element $\mathbf y$ is recursively used in their rules. Recently, simplified versions of these algorithms were obtained and called the {\it simplified topological $\varepsilon$--algorithms} [2]. They have only one recursive rule instead of two, they require less storage than the initial algorithms, elements of the dual vector space $E^*$ of $E$ no longer have to be used in the recursive rules but only in their initializations, the numerical stability is improved, and it was possible to prove theoretical results on them.
In this talk, we present the {\it simplified topological $\varepsilon$--algorithms} and the Matlab package {\tt EPSfun} [3], available in the public domain library {\it netlib}, for implementing and using them.
Then, we give applications to the solution of linear and nonlinear systems of vector and matrix equations, to the computation of matrix functions, to the solution of nonlinear Fredholm integral equations of the second kind, and to the computation of tensor $l^p$-eigenvalues and eigenvectors.
[1] C. Brezinski, Généralisation de la transformation de Shanks, de la table de Padé et de l'$\varepsilon$--algorithme, Calcolo, 12 (1975) 317--360.
[2] C. Brezinski, M. Redivo--Zaglia, The simplified topological $\varepsilon$--algorithms for accelerating sequences in a vector space, SIAM J. Sci. Comput., 36 (2014) A2227--A2247.
[3] C. Brezinski, M. Redivo--Zaglia, The simplified topological $\varepsilon$-algorithms: software and applications, Numer. Algorithms 74 (2017), 1237--1260.
Vassalos A general tool for determining asymptotic spectral distribution of hermitian matrix sequences
The approximation theory for sequences of matrices with increasing dimension is a topic having both theoretical and practical interest. In this talk, we consider sequences of Hermitian matrices with increasing dimension, and we provide a general tool for determining the asymptotic spectral distribution of a 'difficult' sequence ${A_n}_n$ from the one of 'simpler' sequences ${B_{n,m}}_n$ that approximate ${A_n}_n$ when $m \rightarrow \infty$. The tool is based on the notion of an approximating class of sequences (a.c.s.), and it is applied in a more general setting. As an application we illustrate how it can be used in order to derive the famous Szego theorem on the spectral distribution of Toeplitz matrices.
This is a joint work with C. Garoni and S. Serra Capizzano.
De Terán Polynomial root-finding using companion matrices
The use of companion matrices to compute the roots of scalar polynomials is a standard approach (it is for instance, the one followed by the command roots in MATLAB). It consists in computing the roots of a scalar polynomial as the eigenvalues of a companion matrix. In this talk, I will review several numerical and theoretical issues on this topic. I will pay special attention to the backward stability of solving the polynomial root-finding problem using companion matrices. More precisely, to this question: Even if the computed eigenvalues of the companion matrix are the eigenvalues of a nearby matrix, does this guarantee that they are the roots of a nearby polynomial? Usually, the companion matrix approach focuses on monic polynomials, since one can always divide by the leading coefficient, if necessary. But, is it enough for the backward stability issue to focus on monic polynomials? I will also pay attention to some other (more theoretical) questions like: How many companion matrices are there and what do they look like?
Potts High dimensional approximation with trigonometric polynomials
In this talk, we present algorithms for the approximation of multivariate functions by trigonometric polynomials. The approximation is based on sampling of multivariate functions on rank-1 lattices. To this end, we study the approximation of functions in periodic Sobolev spaces of dominating mixed smoothness. The proposed algorithm based mainly on a one-dimensional fast Fourier transform, and the arithmetic complexity of the algorithm depends only on the cardinality of the support of the trigonometric polynomial in the frequency domain. Therefore, we investigate trigonometric polynomials with frequencies supported on hyperbolic crosses and energy based hyperbolic crosses in more detail. Furthermore, we present algorithms where the support of the trigonometric polynomial is unknown.
Fasi Computing the action of the weighted geometric mean of two large-scale matrices on a vector
We consider two classes of methods for computing the product of the weighted geometric mean of two large-scale positive definite matrices and a vector. We derive algorithms based on quadrature formulae for the matrix $p$th root and on the Krylov subspace, and compare these approaches in terms of convergence speed and execution time. By exploiting an algebraic relation between the weighted geometric mean and its inverse, we show how these methods can be used to efficiently solve large and sparse linear systems whose coefficient matrix is a weighted geometric mean.
Del Corso An implicit QR method for unitary plus low rank matrices
In this talk we present an implicit QR method for the approximation of the eigenvalues of unitary plus low rank matrices. Block companion linearizations of matrix polynomials are, perhaps, the most relevant cases of matrices of this form. The problem of computing the roots of a scalar polynomial has received a lot of attention in the last decade [2, 3, 4, 5] and recently the block companion case has been considered by Aurentz and others [1]. The authors gave an implicit QZ method whose cost is $O(d^2 k^3 )$, $d$ being the degree of the matrix polynomial, and $k$ the size of the matrices of the polynomial. This bound is claimed by the authors to be the best one can achieve using the QR or QZ approach. In this talk we present a new implicit QR method for dealing with the eigenvalue problem in the more general case of unitary plus low rank matrices. The algorithm is based on the representation of the matrix as follows \[ A = V (H + EZ^T )W^T , \] where $V$ and $W$ are the product of $k$ upper Hessenberg matrices, $H$ is a Hessenberg matrix , $E = [I_k , O]^ T$ and $Z$ is full rank $n \times k$ matrix. A single QR step can be implemented using Givens transformations, with a cost of $O(nk)$ per iteration, giving a total cost of $O(n^2 k)$. In the case of the block companion linearization we obtain the same asymptotic cost of the method proposed in [1], since $n = kd$. Numerical experiments show that the method is backward stable, and we are able to recover, when the polynomial is not too badly ill conditioned, with high accuracy the eigenvalues of the matrix.
[1] J. Aurentz, T. Mach, L. Robol, R. Vandebril, and D. S. Watkins. Fast and backward stable computation of the eigenvalues and eigenvectors of matrix polynomials. {\em ArXiv e-prints}, November 2016.
[2] Jared L. Aurentz, Raf Vandebril, and David S. Watkins. Fast computation of roots of Companion, Comrade, and related matrices. {\em BIT}, 54(1):85--111, 2014.
[3] R. Bevilacqua, G. M. Del Corso, and L. Gemignani. Compression of unitary rank--structured matrices to CMV-like shape with an application to polynomial rootfinding. {\em J. Comput. Appl. Math.}, 278:326--335, 2015.
[4] D. A. Bini, P. Boito, Y. Eidelman, L. Gemignani, and I. Gohberg. A fast implicit QR eigenvalue algorithm for companion matrices. {\em Linear Algebra Appl.}, 432(8):2006--2031, 2010.
[5] S. Chandrasekaran, M. Gu, J. Xia, and J. Zhu. A fast QR algorithm for companion matrices. In {\em Recent advances in matrix and operator theory}, volume 179 of{\em Oper. Theory Adv. Appl.}, pages 111--143. Birkh\"auser, Basel, 2008.
12.30-14.30 Lunch Break Lunch Break Lunch Break Lunch Break Lunch Break
Barbarino Equivalence between measurable functions and GLT sequences
The theory of Generalized Locally Toeplitz (GLT) sequences established a bridge between the space of matrix sequences and the space of measurable functions.
Recently, it has been proved that the algebra of GLT sequences is actually isomorphic to a space of measurable function, and that the notion of Approximating Class of Sequences (acs) is induced by a pseudometric linked to the convergence in measure.
The latter is metrizable, and for a wide class of distances inducing this convergence it is possible to build a class of corresponding pseudometrics that provides the GLT algebra with a structure isometrically isomorphic to the measurable functions.
These results are used to prove the completeness of the space of matrix sequences, and to provide new tools to test the convergence of matrix sequences through the notion of acs.
Dopico Matrix polynomials with bounded rank and degree: generic eigenstructures and explicit descriptions
Low rank perturbations of matrices, matrix pencils, and matrix polynomials appear naturally in many applications where just a few degrees of freedom of a complicated system are modified. As a consequence many papers have been published in the last 15 years on this type of problems for matrices and pencils, but just a few for matrix polynomials. A possible reason of this lack of references on low rank perturbations of matrix polynomials is that the set of matrix polynomials with bounded (low) fixed rank and fixed degree is not easy to describe when the rank is larger than one. The purpose of this talk is to describe such set both in terms of its generic eigenstructures and in terms of products of two factors.
Meini Solving quadratic matrix equations with infinite quasi-Toeplitz coefficients
Quadratic matrix equations of the kind $A_1X^2+A_0X+A_{-1}=X$, where the coefficients are nonnegative semi-infinite quasi-Toeplitz matrices, are encountered in certain Quasi-Birth-Death (QBD) stochastic processes, as the tandem Jackson queue or the reflecting random walk in the quarter plane.
We provide a numerical framework for approximating the minimal nonnegative solution of these equations which relies on semi-infinite quasi-Toeplitz matrix arithmetic. This arithmetic does not perform any finite size truncation, instead it approximates the infinite matrices through the sum of a banded Toeplitz matrix and a compact correction.
In particular, we show that the algorithm of Cyclic Reduction can be effectively applied and can approximate the infinite dimensional solutions with quadratic convergence at a cost which is comparable to that of the finite case. This way, we may compute a finite approximation of the sought solution, as well as of the invariant probability measure of the associated QBD process, within a given accuracy. Numerical experiments, performed on a collection of benchmarks, confirm the theoretical analysis.
Donatelli Spectral analysis and multigrid preconditioners for space-fractional diffusion equations
Fractional partial diffusion equations (FDEs) are a generalization of classical partial differential equations, used to model anomalous diffusion phenomena. Several discretization schemes (finite differences, finite volumes, etc.) combined with (semi)-implicit methods leads to a Toeplitz-like matrix-sequence. In the constant diffusion coefficients case such a matrix-sequence reduces to a Toeplitz one, then exploiting well-known results on Toeplitz sequences, we are able to describe its asymptotic eigenvalue distribution. In the case of nonconstant diffusion coefficients, we show that the resulting matrix-sequence is a generalized locally Toeplitz (GLT) and then we use the GLT machinery to study its singular value/eigenvalue distribution as the matrix size diverges (see [3]).
The new spectral information is employed for analyzing preconditioned Krylov and multigrid methods recently appeared in the literature, with both positive and negative results. Moreover, such spectral analysis guides the design of new preconditioning and multigrid strategies. We propose new structure preserving preconditioners with minimal bandwidth (and so with efficient computational cost) and multigrid methods for 1D and 2D problems (see [1] for the 1D case and [2] for the 2D case). Some numerical results confirm the theoretical analysis and the effectiveness of the new proposals.
[1] M. Donatelli, M. Mazza, S. Serra-Capizzano: "Spectral analysis and structure preserving preconditioners for fractional diffusion equations", J. Comput. Phys., Vol. 307, pp. 262--279, 2016.
[2] H. Moghaderi, M. Dehghan, M. Donatelli, M. Mazza: "Spectral analysis and multigrid preconditioners for two-dimensional space-fractional diffusion equations", arXiv, 1926487.
[3] S. Serra-Capizzano: "The GLT class as a generalized Fourier Analysis and applications", Linear Algebra Appl. Vol. 419, pp. 180--233, 2006.
Speleers Design of fast multigrid solvers for isogeometric analysis: a symbol approach
In this talk we focus on the numerical solution of elliptic problems by means of isogeometric analysis. By exploiting specific spectral properties compactly described by a symbol, we are able to design efficient multigrid methods for the fast solution of the related linear systems.
Despite the theoretical optimality, the convergence rate of multigrid methods with classical stationary smoothers worsens exponentially as the involved spline degrees increase. With the aid of the symbol, we can give a theoretical interpretation of this exponential worsening. Moreover, a proper factorization of the symbol allows us to develop a preconditioned conjugate gradient smoother, in the spirit of the multi-iterative strategy, that results in a good multigrid convergence rate, independent of both matrix size and spline degree. Numerical experiments confirm the effectiveness of our proposal and the numerical optimality.
This is joint work with Marco Donatelli, Carlo Garoni, Carla Manni and Stefano Serra-Capizzano.
Boito An algorithm for efficient solution of shifted quasiseparable linear systems and applications
Solving sets of shifted linear systems $(A + \sigma_k I) x_k = b_k$ is a well-studied problem in numerical linear algebra. In many cases, the matrix $A$ is structured. This work focuses on the case where A is quasiseparable, that is, the off-diagonal blocks of A have low rank. This is often the case, for instance, when discretizing differential operators. We present a fast algorithm that relies on structured QR factorization and exploits both the quasiseparable and the shifted structure of the problem. The main application we consider is the computation of the product of matrix functions times a vector via rational approximations or contour integrals. The same approach can also be used for the solution of linear matrix equations.
Palitta Numerical methods for Lyapunov matrix equations with banded symmetric data
We are interested in the numerical solution of the large-scale Lyapunov equation $AX + XA^T = C$, where $A, C \in\mathbb R^{ n\times n}$ are both large and banded matrices. We suppose that $A$ is symmetric and positive definite and $C$ is symmetric and positive semidefinite.
While the case of low-rank $C$ has been successfully addressed in the literature, the more general banded setting has not received much attention, in spite of its possible occurrence in applications. In this talk we aim to fill this gap. It has been recently shown that if $A$ is well conditioned, the entries of the solution matrix $X$ decay in absolute value as their indexes move away from the sparsity pattern of $C$. This property can be used in a memory-saving matrix-oriented Conjugate Gradient method to obtain a banded approximate solution.
For $A$ not well conditioned, the entries of X do not sufficiently decay to derive a good banded approximation. Nonetheless, we show that it is possible to split $X$ as $X = Z_b + Z_r$ , where $Z_b$ is banded and $Z_r$ is numerically low rank. We thus propose a novel strategy that efficiently approximates both $Z_b$ and $Z_r$ with acceptable memory requirements. Numerical experiments are reported to illustrate the potential of the discussed method.
Bolten Analysis of parallel time integrators using structured matrices
With the advent of large scale high performance computers with 100.000s of cores scalability of numerical algorithms becomes one of the most important aspects. In many applications problems do not only depend on spatial degrees of freedom, but they evolve during time. Usually, the speedup of a pure spatial distribution of these problems saturates at some point, requiring to rethink the parallelization strategy. One option are parallel time integrators that decompose the domain in all dimensions.
A variety of parallel time integration schemes exists, e.g., the well-known parareal method by Lions, Maday, and Turinici, PFASST by Emmett and Minion, or space-time multigrid that has been considered already by Hackbusch. The analysis of these methods often leads to structured matrices. In many cases a bidiagonal block Toeplitz matrix has to be studied in order to assess the properties of the numerical method. Parallel-in-time integrators share many similarities with multigrid methods, so the theory and methodology developed in the context of multigrid methods for structured matrices can be applied straightforwardly in this case.
We will present the analysis of some time-parallel method using structured matrices. Using these techniques we are able to analyze the behavior of the studied methods. The results provide a better understanding of parallel time integrators and allow to design more efficient methods.
15.30-16.00 Poster blitz
Eidelman Rational approximations and fast quasiseparable algorithms for matrix and operator functions and differential equations
We discuss methods for computations of matrix and operator functions via rational approximations. Using quasiseparable representations of matrices we obtain fast algorithms to solve these problems. The model examples we use are various problems for differential equations.
Massei Solving large scale Sylvester and Lyapunov equations with quasiseparable coefficients
In this talk we address the problem of efficiently solving $AX+XB=C$ in the case where the coefficients $A,B$ and $C$ are $m\times m$ matrices with rank structured off-diagonal blocks. We prove that under reasonable assumptions the structure is present in the solution. We design and test different strategies that --- combined with the hierarchical matrices technology --- are able to solve the equation with linear-polylogarithmic complexity. Extension to the treatment of certain generalized Sylvester equation is discussed.
Sesana Multigrid methods for block-circulant linear systems
Many problems in physics, engineering and applied sciences are governed by functional equations (in differential or integral form) that do not admit closed form solution and, therefore, require numerical discretization techniques which often involve the solution of linear systems of large size. The main novelty contained in the literature treating structured matrices is the use of the \textit{symbol}, that is a function $f$ that provides a compact description of the asymptotic global behavior of the spectrum of the sequence of the coefficient matrices $\{A_n\}$, obtained from the discretization of a differential or integral problem when the fineness parameter tends to zero.
Among the iterative methods available to solve such systems, in the last 25 years Multigrid methods have gained a remarkable reputation as fast solvers for structured matrices associated to shift invariant operators, where the size $n$ of the problem is large. The convergence analysis of two particular Multigrid method, the Two-grid and V-cycle, has been handled in a compact and elegant manner by studying few analytical properties of the symbol $f$ associated with the sequence of coefficient matrices.
In the cases taken into account, especially concerning Toeplitz matrices, this symbol is a (multivariate) scalar-valued function $f$, while much remains to be done in the case of a matrix-valued symbols, which are obtained for example in the discretization of systems of equations.
In this talk our aim is to fill this gap, generalizing some of the existing proofs for the Two-grid and the V-cycle method for systems with matrix in algebra, such as circulant, Hartley and tao, to the case where the latter have a matrix-valued symbol. The next step is the extension for matrices non in algebra, such as Toeplitz matrices.
16.00-16.30 Coffee Break
Poster session
Coffee Break
Coffee Break
16.30-17.00 Poster session
Kressner Low-rank updates of matrix functions
This talk is concerned with the development and analysis of fast algorithms for updating a matrix function $f(A)$ if $A$ undergoes a low-rank change $A+L$. For example, when $A$ is the Laplacian of an undirected graph then removing one edge or vertex of the graph corresponds to a rank-one change of $A$. The evaluation and update of such matrix functions (or parts thereof) is of relevance in the analysis of network community structure of networks and signal graph processing. Our algorithms are based on the tensorization of polynomial or rational Krylov subspaces involving $A$ and $A^T$. The choice of a suitable element from such a tensorized subspace for approximating $f(A+L)-f(A)$ is straightforward in the symmetric case but turns out to be more intricate in the nonsymmetric case. We show that the usual convergence results for Krylov subspace methods for matrix functions can be extended to our algorithms. If time permits, we will also touch upon the closely related problem of applying the Fréchet derivative of a matrix function to a low-rank matrix.
This is joint work with Bernhard Beckermann and Marcel Schweitzer.
Mitrouli Estimating bilinear forms for some large scale computation problems
A spectrum of applications arising from Statistics, Machine Learning, Network Analysis require computation of bilinear forms $x^Tf (A)y$, where $A$ is a diagonalizable matrix and $x$, $y$ are given vectors. In this work we are interested in efficiently computing bilinear forms primarily due to their importance in several contexts. For large scale computation problems it is preferable to achieve approximations of bilinear forms without exploiting the whole matrix function. For this purpose an extrapolation procedure has been developed, attaining the approximation of bilinear forms with one, two or three term estimates in a complexity of square order. Furthermore, a prediction approach based on Aitken's acceleration scheme has also been developed computing alternative estimates again in a square complexity. Both schemes are characterized by easy applicable formulae of low complexity that can be implemented in vectorized form.
Cicone Spectral and convergence analysis of the discrete Adaptive Local Iterative Filtering method by means of Generalized Locally Toeplitz sequences
In this talk we introduce a newly developed algorithm, the Adaptive Local Iterative Filtering (ALIF). This is an alternative technique to the well known Empirical Mode Decomposition method for the decomposition of nonstationary real life signals. This last algorithm has been successfully applied in many fields of research in the last two decades.
The mathematical convergence analysis of all the aforementioned algorithms is still lacking. For this reason, focusing the attention on the discrete version of ALIF, we show how recent results about sampling matrices and, in particular, the theory of Generalized Locally Toeplitz sequences allow to perform a spectral analysis of the matrices involved in the iterations. In particular we are able to study the eigenvalue clustering and the eigenvalue distribution for such matrices; we provide a necessary condition for the convergence of the Discrete ALIF method and we derive a simple criterion to construct filters, needed in the ALIF algorithm, that guarantee the fulfillment of the aforementioned necessary condition. We show some numerical examples and discuss about important open problems that are waiting to be tackled.
Pan Low Rank Approximation: Failing but Accurate Superfast Algorithms, Pre-processing and Extensions
Low rank approximation (hereafter LRA) of a large matrix has routinely been computed by using much fewer memory units and arithmetic operations than the input matrix has entries, but so far no proof supports the consistently observed accuracy of such superfast algorithms.
We first explain why such a proof has been missing -- we specify a small family of matrices (we call them $\pm \delta$-{\em matrices}) for most of which LRA computed by any superfast algorithm is no better than just by the matrix filled with zeros.
Then we analyze superfast LRA algorithms of three kinds and prove that for each of them the class of such hard inputs is narrow: these algorithms output close LRAs to the average matrix allowing LRA, to any matrix allowing LRA unless it degenerates like $\pm \delta$-matrices, and with a high probability even to such a degenerate matrix if it is pre-processed with a random Gaussian, SRHT or SRFT multiplier; moreover empirical efficiency of these multipliers is matched by our much sparser multipliers.
We also provide some recipes and formal support for enhancing the accuracy of such superfast LRA and discuss further research directions.
Our progress relies on our new insights and techniques, should encourage application of simplified heuristics, the average case analysis and randomized multiplicative pre-processing in matrix computations, and can be immediately extended to various important computational areas, e.g., tensor decomposition, but we also show a novel surprising extension to dramatic acceleration of the Conjugate Gradient algorithms.
Dell'Acqua Taylor boundary conditions for accurate image restoration
In recent years, several efforts were made in order to introduce boundary conditions for deblurring problems that allow to get accurate reconstructions. This resulted in the birth of Reflective, Anti-Reflective and Mean boundary conditions, which are all based on the idea of guaranteeing the continuity of the signal/image outside the boundary. Here we propose new boundary conditions that are obtained by suitably combining Taylor series and finite difference approximations. Moreover, we show that also Anti-Reflective and Mean boundary conditions can be attributed to the same framework. Numerical results show that, in case of low levels of noise and blurs able to perform a suitable smoothing effect on the original image (e.g. Gaussian blur), the proposed boundary conditions lead to a significant improvement of the restoration accuracy with respect to those available in the literature.
Garoni Symbol Analysis of Differential Eigenproblems
Isogeometric Analysis (IgA) and Finite Element Analysis (FEA) are two distinguished numerical methods for the numerical solution of differential problems. While FEA is a very popular technique, which dates back to the 1950s, IgA has been introduced only recently, between 2005 and 2009. However, due to its capability to enhance the connection between numerical simulation and Computer-Aided Design (CAD) systems, IgA is gaining more and more attention over time.
In this presentation, we focus on the numerical solution of a simple eigenvalue problem by means of IgA and FEA. We show that IgA is superior to FEA in the spectral approximation, because, while the numerical eigenvalues computed by IgA reproduce (approximately) all the exact eigenvalues, only a small portion of the numerical eigenvalues computed by FEA can be considered as approximations of the exact eigenvalues.
Our analysis is based on the notion of "symbol", which will be introduced in the talk along with IgA and FEA.
Gemignani On the Design of Fast EigenSolvers for Diagonal plus Low-Rank Matrices
Most of researches on rank structured matrix technology for structured eigenproblems stem from the design of efficient eigensolvers for diagonal plus low-rank matrices. In this talk we present two fast algorithms for Hessenberg reduction of a structured matrix $A = D + UV^H$ where $D$ is a real or unitary $n \times n$ diagonal matrix and $U, V \in\mathbb{C}^{n \times k}$. The proposed algorithm for the real case exploits a two--stage approach by first reducing the matrix to a generalized Hessenberg form and then completing the reduction by annihilation of the unwanted sub-diagonals. It is shown that the novel method requires $O(n^2k)$ ops and it is significantly faster than other reduction algorithms for rank structured matrices. The method is then extended to the unitary plus low rank case by using a block analogue of the CMV form of unitary matrices. It is shown that a block Lanczos-type procedure for the block tridiagonalization of $\Re(D)$ induces a structured reduction on $A$ in a block staircase CMV--type shape. Then, we present a numerically stable method for performing this reduction using unitary transformations and we show how to generalize the sub-diagonal elimination to this shape, while still being able to provide a condensed representation for the reduced matrix with cost $O(n^2k)$ ops.
Winkler A structure-preserving matrix method for the computation of multiple roots of a Bernstein basis polynomial
The Bernstein basis is used for the representation of curves and surfaces in geometric modelling systems because of its elegant geometric properties and enhanced numerical stability with respect to the power basis. One of the most important computations in these systems is the calculation of the points of intersection of curves and surfaces, which reduces to the calculation of the roots of a polynomial. Particular interest is focussed on multiple roots because they define conditions of tangency, and therefore smooth intersections, which are required for the reduction of stress-concentration factors at sharp corners, and ease of handling -- it is easier to handle an object with rounded corners than an object with sharp corners.
This paper considers the application of a structure-preserving matrix method for the computation of multiple roots of a Bernstein polynomial $f (y)$. Extensive use is made of the Sylvester matrix and its subresultant matrices for the calculation of the multiplicities of the distinct roots of $f (y)$, which involves many greatest common divisor computations. It is shown that the development of this polynomial root solver is significantly more complicated than its equivalent for the power basis form of $f (y)$ because of the combinatorial terms in the Bernstein basis functions. In particular, even if the coefficients of $f (y)$ are of the same order of magnitude, the presence of the combinatorial factors implies that, even for polynomials of modest degree, the entries in the matrices that arise in the computations may span several orders of magnitude, which can cause numerical problems that do not arise when the power basis form of $f (y)$ is considered. Furthermore, the presence of these combinatorial factors in the Sylvester matrix destroys its concatenated Toeplitz matrix structure, which increases the computational cost of the method. It is shown that the adverse effects of the combinatorial factors can be mitigated by multiplying the Sylvester matrix by a diagonal matrix, and that this operation improves the accuracy of the computed roots.
Computational results obtained from this root solver, and root solvers that exploit the geometric properties of the Bernstein basis, will be presented and it will be shown that the root solver described above yields considerably better results.
Aceto, Rational approximations to the fractional Laplacian
Fractional-order in space mathematical models, in which an integer-order differential operator is replaced by a corresponding fractional one, are becoming increasingly used since they provide an adequate description of many processes that exhibit anomalous diffusion. In this talk, in particular, we focus on the numerical solution of fractional in space reaction-diffusion equations on bounded domains under homogeneous Dirichlet boundary conditions. Using the matrix transfer technique, the fractional Laplacian operator is replaced by a matrix which, in general, is dense [3, 4]. The approach proposed is based on the approximation of this matrix by the product of two suitable banded matrices [1]. Since these two matrices will be involved in the linear algebra tasks required by the differential solver, it is fundamental that their condition numbers are not responsible for inaccuracy, otherwise there will be an a-priori barrier for the choice of the bandwidth. In order to face this problem, we use a simple but reliable strategy that allows to keep the conditioning under control [2]. Work in collaboration with Paolo Novati, Department of Mathematics and Geosciences, University of Trieste.
[1] L. Aceto, P. Novati, Rational approximation to the fractional Laplacian operator in reaction-diffusion problems. SIAM J. Sci. Comput., 39 (2017), no. 1, A214--A228.
[2] L. Aceto, P. Novati, Efficient implementation of rational approximations to fractional differential operators, submitted.
[3] M. Ilic, F. Liu, I. Turner, V. Anh, Numerical approximation of a fractional-in-space diffusion equation I, Fract. Calc. Appl. Anal., 8 (2005) 323--341.
[4] M. Ilic, F. Liu, I. Turner, V. Anh, Numerical approximation of a fractional-in-space diffusion equation (II)-with nonhomogeneous boundary conditions, Fract. Calc. Appl. Anal., 9 (2006) 333--349.
Bertaccini, Fast Solution of Fractional Differential Equation
In our talk we propose an innovative algorithm for the large (full) linear systems of time-dependent partial fractional differential equations discretized in time with linear multistep formulas, both in classical [1] and in boundary value form [2]. We use the short--memory principle to ensure the decay of the entries of sparse approximations of the discretized operator and its inverse. FGMRES with preconditioners based on short--memory principle as well are then used to solve the underlying sequence of linear systems.
The ideas above are implemented on GPU devices by the techniques for sparse approximate inverse preconditioners proposed in [3].
[1] Bertaccini, D., Durastante, F. (2017). Solving mixed classical and fractional partial differential equations using short memory principle and approximate inverses. Numer. Algorithms, 74(4), 1061--1082.
[2] Bertaccini, D., Durastante, F. (2017). Limited memory block preconditioners for fast solution of fractional PDEs. (submitted )
[3] Bertaccini, D., Filippone, S. (2016). Sparse approximate inverse preconditioners on high performance GPU platforms. Comput. Math. Appl., 71(3), 693--711.
Cipolla, Adaptive matrix algebras in unconstrained minimization
In this communication we will introduce some recent techniques which involve structured matrix spaces in the reduction of time and space complexity of BFGS-type minimization algorithms [1], [2]. Some general results for the global convergence of algorithms for unconstrained optimization based on a BFGS-type Hessian approximation scheme are introduced and it is shown how the constructibility of convergent algorithms suitable for large scale problems can be tackled using projections onto low complexity matrix algebras.
[1] C.Di Fiore, S.Fanelli, F.Lepore, P.Zellini, Matrix algebras in Quasi-Newton methods for unconstrained minimization, Numerische Mathematik, 94, (2003).
[2] S.Cipolla, C.Di Fiore, F.Tudisco, P.Zellini, Adaptive matrix algebras in unconstrained minimization, Linear Algebra and its Application, 471, (2015).
Concas, Matlab implementation of a spectral algorithm for the seriation problem
The seriation problem is an important ordering issue related to finding all the possible ways to sequence the elements of a certain set.
It is frequently used in archeology and has applications in many fields such as genetics, bioinformatics and graph theory.
We will present a Matlab implementation of a spectral method, which was presented in [1], to solve the seriation problem.
This algorithm is based on the use of the Fiedler vector of the Laplacian matrix associated to the problem and it describes the results in terms of a particular data structure called PQ-tree used for storing the admissible orderings. We will discuss the case of the presence of a multiple Fiedler value and some numerical examples of graphs for which is not possible to find a precise solution [2].
[1] Jonathan E. Atkins, Erik G. Boman, and Bruce Hendrickson. A spectral algorithm for seriation and the consecutive ones problem. SIAM Journal on Computing, 28(1):297--310, 1998.
[2] Anna Concas, Caterina Fenu, and Giuseppe Rodriguez. PQ-ser: a Matlab package for spectral seriation. In preparation.
Durastante, Fast Solution of Fractional Differential Equation
In our talk we propose an innovative algorithm for the large (full) linear systems of time-dependent partial fractional differential equations discretized in time with linear multistep formulas, both in classical [1] and in boundary value form [2]. We use the short--memory principle to ensure the decay of the entries of sparse approximations of the discretized operator and its inverse. FGMRES with preconditioners based on short--memory principle as well are then used to solve the underlying sequence of linear systems.
The ideas above are implemented on GPU devices by the techniques for sparse approximate inverse preconditioners proposed in [3].
[1] Bertaccini, D., Durastante, F. (2017). Solving mixed classical and fractional partial differential equations using short memory principle and approximate inverses. Numer. Algorithms, 74(4), 1061--1082.
[2] Bertaccini, D., Durastante, F. (2017). Limited memory block preconditioners for fast solution of fractional PDEs. (submitted )
[3] Bertaccini, D., Filippone, S. (2016). Sparse approximate inverse preconditioners on high performance GPU platforms. Comput. Math. Appl., 71(3), 693--711.

Fasino, Error Analysis of TT-format Tensor Algorithms
The {\em tensor train decomposition} is a representation technique which allows compact storage and efficient computations with arbitrary tensors [2]. Basically, a tensor train (TT) decomposition of a $d$-dimensional tensor ${\bf A}$ with size $n_1\times n_2 \times \cdots \times n_d$ is a sequence $G_1,\ldots,G_d$ of 3D tensors (the {\em carriages}) such that the size of $G_i$ is $r_{i-1}\times n_i \times r_i$ with $r_0 = r_d = 1$ (that is, $G_1$ and $G_d$ are ordinary matrices) and \[ {\bf A}(i_1,i_2,\ldots,i_d) = \sum_{\alpha_1,\ldots,\alpha_{d-1}} G_1(i_1,\alpha_1)G_2(\alpha_1,i_2,\alpha_2)\cdots G_d(\alpha_{d-1},i_d). \] The index $\alpha_i$ runs from $1$ to $r_i$, for $i=1,\ldots,d-1$, and the numbers $r_1,\ldots,r_{d-1}$ are the {\em TT-ranks} of ${\bf A}$.
An alternative viewpoint on this decomposition is \[ {\bf A}(i_1,i_2,\ldots,i_d) = G_1'(i_1)G_2'(i_2)\cdots G_d'(i_d), \] where now $G'(i_k)$ is an $r_{k-1}\times r_k$ matrix depending on the integer parameter $i_k$.
We will use the notation ${\bf A} = \mathrm{TT}(G_1,\ldots,G_d)$.
We present a backward error analysis of two algorithms found in [3] which perform computations with tensors in TT-format. The first one produces an exact or approximate TT-decomposition $G_1,\ldots,G_d$ of a tensor $\bf A$ given in functional form, depending on a tolerance $\varepsilon$.
If $\varepsilon = 0$ then the output of the algorithm is an exact TT-decomposition, that is, ${\bf A} = \mathrm{TT}(G_1,\ldots,G_d)$. If $\varepsilon > 0$ then $\mathrm{TT}(G_1,\ldots,G_d)$ is an $\mathcal{O}(\varepsilon)$-approximation of $\bf A$ which can realize significant savings in memory space.
The computational core of the algorithm is a suitable (approximate) matrix factorization that, in the original paper, relies on SVD computations.
We prove that analogous performances can be obtained by means of QR factorizations.
Moreover, we obtain a stability estimate that, in the exact arithmetic case, reduces to the error estimate given in [3] for the SVD-based decomposition and, in the floating point arithmetic case, separates the effects of the tolerance $\varepsilon$, the lack of orthogonality of certain intermediate matrices, and the stability properties of the basic factorization on the backward error of the computed TT-decomposition.
The second algorithm computes the multilinear form (also called {\em contraction}) of a given $d$-dimensional tensor in TT-format and vectors $v_1,\ldots,v_d$, \[ a = \sum_{i_1=1}^{n_1} \cdots \sum_{i_d=1}^{n_d} \sum_{\alpha_1,\ldots,\alpha_{d-1}} G_1(i_1,\alpha_1)G_2(\alpha_1,i_2,\alpha_2) \cdots G_d(\alpha_{d-1},i_d) v_1(i_1) \cdots v_d(i_d). \] By means of known error bounds for inner products in floating point arithmetic [2], we prove backward stability of the proposed algorithm under very general hypotheses on the evaluation order of the innermost summations. More precisely, if ${\bf A} = \mathrm{TT}(G_1,\ldots,G_d)$ and no underflows or overflows are encountered then the output $\widehat a$ computed by the algorithm in floating point arithmetic is the exact contraction of ${\bf \widehat A} = \mathrm{TT} (G_1 + \Delta G_1,\ldots,G_d + \Delta G_d)$ and $v_1,\ldots,v_d$ where $|\Delta G_i| \leq (n_i+r_{i-1})u|G_i| + \mathcal{O}(u^2)$ and $u$ is the machine precision.
[1] C.-P.~Jeannerod, S.~M.~Rump. Improved error bounds for inner products in floating-point arithmetic. {\em SIMAX} 34 (2013), 338--344.
[2] I.~Oseledets. Tensor-train decomposition. {\em SIAM J. Sci. Comput.} 33 (2011), 2295--2317.
[3] I.~Oseledets, E.~Tyrtyshnikov. TT-cross approximation for multidimensional arrays. {\em Lin.~Alg.~Appl.} 432 (2010).
Fenu, On the computation of the GCV function for Tikhonov Regularization
Tikhonov regularization is commonly used for the solution of linear discrete ill-posed problems with error-contaminated data. A regularization parameter, that determines the quality of the computed solution, has to be chosen. One of the most popular approaches to choosing this parameter is to minimize the Generalized Cross Validation (GCV) function. The minimum can be determined quite inexpensively when the matrix A that defines the linear discrete ill-posed problem is small enough to rapidly compute its singular value decomposition (SVD). We are interested in the solution of linear discrete ill-posed problems with a matrix A that is too large to make the computation of its complete SVD feasible. We will present two fairly inexpensive ways to determine upper and lower bounds for the numerator and denominator of the GCV function for large matrices A [1, 2]. The first one is based on Gauss-type quadrature and the second one on a low-rank approximation of the matrix A. These bounds are used to determine a suitable value of the regularization parameter. Computed examples illustrate the performance of the proposed methods.
[1] C. Fenu, L. Reichel, and G. Rodriguez. GCV for Tikhonov regularization via global Golub-Kahan decomposition. Numer. Linear Algebra Appl., 23:467-484, 2016.
[2] C. Fenu, L. Reichel, G. Rodriguez and H. Sadok. GCV for Tikhonov regularization by partial singular value decomposition. BIT Numer. Math., DOI 10.1007/s10543-017-0662-0, 2017.
Iannazzo, Nonnegative factorization of tensor grids
The Nonnegative Matrix Factorization (NMF) approximately decomposes a (large) matrix with nonnegative entries into the product of two matrices with nonnegative entries and with smaller dimensions.
A similar decomposition is desired when dealing with tensor grids, matrices whose entries are positive definite matrices and not just numbers. In this case, NMF is not suited because of the presence of possibly negative entries.
We present two generalizations of the NMF preserving the positive definite structure of tensor grids. They are based on the use of non-Euclidean geometries on the set of positive definite matrices.
Poloni, Counting Fiedler pencils with diagrams
Fiedler pencils (with repetitions) are a family of matrix pencils which generalizes the well-known companion form: they are linearizations of matrix polynomials, i.e., they provide a template to construct, given a matrix polynomial, a linear eigenvalue problem with the same eigenvalues and multiplicity. Fiedler pencils are constructed as product of special block matrices that act nontrivially only on two contiguous blocks. They have a rich structure that gives rise to many combinatorial properties.
We introduce a notation that associates to each pencil a diagram that depicts its action on the blocks. Using it, we can obtain visual proofs of several results in the theory, and we can solve several counting problems (such as "how many distinct Fiedler pencils with repetitions of a given dimension exist"). Among them, in particular, we are interested in counting Fiedler pencils associated to symmetric and palindromic matrix polynomials which preserve the same structure.
Pozza, Lanczos algorithm and the Gauss quadrature for linear functionals
Gauss quadrature can be naturally generalized to approximate quasi-definite linear functionals [1], where the interconnections with formal orthogonal polynomials, Padé approximants, complex Jacobi matrices and Lanczos algorithm are analogous to those in the positive definite case. As presented in [2], the existence of the n-weight complex Gauss quadrature corresponds to performing successfully the first $n$ steps of the non-Hermitian Lanczos algorithm. Some further results on the relation between the non-definite case and the look-ahead Lanczos algorithm will be shown.
[1] S. Pozza, M. Prani, Z. Strakos, Gauss quadrature for quasi-definite linear functionals, to appear in: IMA Journal of Numerical Analysis.
[2] S. Pozza, M. Prani, Z. Strakos, Lanzcos algorithm and the complex Gauss quadrature, submitted, May 2017.
Roupa Approximation of matrix-vector products and applications
The estimation of matrix-vector products such as $A^{1/2}b$, $\exp(A)b$, $A^{-1}b$, where $A\in\mathbb R^{p\times p}$ is a given diagonalizable matrix and $b\in\mathbb R^p$, appears in many applications arising from the fields of mathematics, statistics, mechanics, networks, machine learning and physics. In particular, the matrix-vector product $A^{1/2}b$ is used in sampling from a Gaussian process distribution. In network analysis the vector form $\exp(A)b$ determines the total communicability as a global measure of how well the nodes in a graph can exchange information [1].
In this work, methods for approximating matrix-vector products based on extrapolation [3], Krylov subspaces, [4] and polynomial approximations [2] are presented. The above methods for estimating vector forms are compared concerning the accuracy, the computational complexity and the execution time.
Several numerical examples will be presented to illustrate the effectiveness of these methods.
[1] F. Arrigo, M. Benzi, Updating and downdating techniques for optimizing network communicability, SIAM J. Sci. Comput., 38 (2016), pp. B25-B49.
[2] J. Chen, M. Anitescu, Y. Saad, Computing f (A)b via least squares polynomial approximations, SIAM J. Sci. Comput., 33 (2011), pp. 195-222.
[3] P. Fika, M. Mitrouli, P. Roupa, Estimates for the bilinear form $x^TA^{-1}y$ with applications to linear algebra problems, Elec. Trans. Numer. Anal., 43 (2014), pp. 70-89.
[4] N. J. Higham, Functions of Matrices: Theory and Computation, SIAM, Philadelphia, PA, USA, 2008.
A selection of pictures of the meeting will be posted after the workshop.