Due Giorni di Algebra Lineare Numerica e Applicazioni

20 Jan 2025, 08:30 → 21 Jan 2025, 22:00 Europe/Rome

Building E (Aula Magna "Fratelli Pontecorvo", Building E, Polo Fibonacci. Pisa)

Fabio Durastante (Università di Pisa), Leonardo Robol (University of Pisa)

Due Giorni di Algebra Lineare Numerica e Applicazioni @ UNIPI

This is the 20th workshop in a series dedicated to Numerical Linear Algebra and Applications. It aims at gathering the Italian Numerical Linear Algebra scientific community to discuss recent advances in this area and to promote the exchange of novel ideas and collaboration among researchers. We strongly encourage contributions by junior scholars on theoretical and computational topics, as well as on applications. Contributions covering other areas requiring tools from Numerical Linear Algebra are very welcome. 

Topics of interest include (but are not limited to): iterative methods for large-scale linear systems, preconditioning, eigen and singular values, numerical optimization, numerical dynamical systems, mathematics of machine learning, structured matrices, tensors, matrix functions, Krylov methods, randomized linear algebra algorithms, and high-performance computing for matrix computations. Applications include optimization and control, image processing and inverse problems, networks and graphs, machine learning and artificial neural networks, uncertainty quantification, finance, material sciences, life sciences, and industrial processes.

---

Steering committee

Michele Benzi, Marco Donatelli, Nicola Guglielmi, Beatrice Meini, Valeria Simoncini

---

Organizing committee

Paola Boito, Gianna Maria Del Corso, Fabio Durastante, Silvia Gazzola, Luca Heltai, Stefano Massei, Cecilia Pagliantini, Federico Poloni, Leonardo Robol.

---

Previous editions

1. 1997, Università di Pavia
2. 2000, Università di Pisa
3. 2002, Università di Pisa
4. 2003, Università di Pisa
5. 2004, CISM, Udine
6. 2005, Università dell'Insubria, Varese
7. 2006, Università di Pisa
8. 2007, Università di Padova
9. 2008, Università di Bologna
10. 2009, Università di Perugia
11. 2011, Università dell'Insubria, Como
12. 2012, Università di Genova
13. 2013, Università di Roma "Tor Vergata"
14. 2017, Università dell'Insubria, Como
15. 2018, Università di Padova
16. 2019, Università di Roma Sapienza
17. 2022, Università di Napoli Federico II
18. 2022, ALAMA-ALN2gg, University of Alcalá, Madrid
19. 2023, Gran Sasso Science Institute, L'Aquila

---

Sponsors 

		PRIN2022: Low-rank Structures and Numerical Methods in Matrix and Tensor Computations and their Application

Monday, 20 January

08:30 Registration

Morning Session I

1 Lanczos with compression for symmetric matrix functions

In this work we present a low-memory method for the approximation of the action of a symmetric matrix function $f(A) \in \mathbb{R}^{n \times n}$ on a vector $\mathbf b \in \mathbb{R}^n$, where the matrix $A$ is large and sparse. A popular approach for approximating $f(A) \mathbf b$ is the Lanczos algorithm. Given an orthonormal basis $Q_M \in \mathbb{R}^{n \times M}$ of the Krylov subspace $\mathcal{K}_M(A, \mathbf b) = \text{span} \{ \mathbf b, A \mathbf b, \dots, A^{M-1} \mathbf b \}$, the Lanczos approximation of $f(A) \mathbf b$ is given by
$$\mathbf f_M = Q_M f(T_M) Q_M^T \mathbf b, \qquad \text{where} \qquad T_M = Q_M^T A Q_M.$$ Although the basis $Q_M$ can be constructed with a short-term recurrence in a setting with limited memory, the computation of $\mathbf f_M$ still requires the storage of the whole Krylov basis. A simple way to overcome this issue and obtain a low-memory algorithm is to use the two-pass Lanczos approach, which however requires double the number of matrix vector products with $A$. Several other techniques have been proposed in the literature. In this talk we propose a novel algorithm to approximate $f(A) \mathbf b$ in a limited memory setting, by combining the Lanczos method with a basis compression procedure involving rational Krylov subspaces, which is employed whenever the basis grows beyond a certain size. The method that we propose has essentially the same convergence behaviour as Lanczos, with an additional error term that depends on rational approximation of the function $f$ and is typically negligible. The cost of the basis compression procedure is also negligible with respect to the cost of the Lanczos algorithm. In particular, the rational Krylov subspaces used for the compression of the Lanczos basis are built using small projected matrices, so their construction is cheap and does not involve the solution of any linear systems with the matrix $A$. The key theoretical result used in the development of our algorithm is a formula for the computation of low-rank updates of matrix functions using rational Krylov subspaces, derived in [1]. This update formula is applied iteratively to the tridiagonal projected matrix $T_M$, by splitting it into a block diagonal matrix and a rank-two off-diagonal part. The resulting algorithm is composed of a sequence of cycles: in each cycle, we first perform a fixed number of Lanczos iterations to expand the Krylov basis, and then we construct an appropriate rational Krylov supspace associated with the current projected matrix in order to compress the basis back to a fixed size, while also updating the approximation to $f(A) \mathbf b$. Numerical experiments demonstrate that our algorithm exhibits competitive performance when compared against other low-memory methods for $f(A) \mathbf b$. Further details can be found in the preprint [2].

[1] B. Beckermann, A. Cortinovis, D. Kressner, and M. Schweitzer, Low-rank updates of matrix functions II: rational Krylov methods, SIAM J. Numer. Anal., 59 (2021), pp. 1325-1347.

[2] A. A. Casulli and I. Simunec, A low-memory Lanczos method with rational Krylov compression for matrix functions, arXiv:2403.04390, 2024.

Speaker: Igor Simunec (École Polytechnique Fédérale de Lausanne)

2 Estimation of spectral gaps for sparse symmetric matrices

We propose and analyze an algorithm for identifying spectral gaps of a real symmetric matrix $A$ by simultaneously approximating the traces of spectral projectors associated with multiple different spectral slices. Our method utilizes Hutchinson's stochastic trace estimator together with the Lanczos algorithm to approximate quadratic forms involving spectral projectors. Instead of focusing on determining the gap between two particular consecutive eigenvalues of $A$, we aim to find all gaps that are wider than a specified threshold. By examining the problem from this perspective, and thoroughly analyzing both the Hutchinson and the Lanczos components of the algorithm, we obtain error bounds that allow us to determine the numbers of Hutchinson's sample vectors and Lanczos iterations needed to ensure the detection of all gaps above the target width with high probability. In particular, we conclude that the most efficient strategy is to always use a single random sample vector for Hutchinson's estimator and concentrate all computational effort in the Lanczos algorithm. Our numerical experiments demonstrate the efficiency and reliability of this approach.

Speaker: Michele Rinelli (KU Leuven)

3 On generating Sobolev orthonormal rational functions

In this research, we generate Sobolev rational functions with prescribed poles that are orthonormal with respect to a discrete Sobolev inner product through a long recurrence relation. The recurrence coefficients of this relation can be stored in a Hessenberg pencil. Our approach to generate Sobolev orthonormal rational functions involves reformulating the main problem to a Hessenberg pencil inverse eigenvalue problem as a matrix problem. This reformulation is infact based on the connection between Sobolev orthonormal rational functions and the orthonormal bases for rational Krylov subspaces generated by a Jordan-like matrix. A new updating procedure is then implemented to solve the equivalent Hessenberg pencil inverse eigenvalue problem by exploiting plane rotations. The performance of the introduced method is also examined through some numerical examples.

Speaker: Amin Faghih (KU Leuven)

4 An Augmented Lagrangian preconditioner for the control of the Navier--Stokes equations

Optimal control problems with PDEs as constraints arise very often in scientific and industrial applications. Due to the difficulties arising in their numerical solution, researchers have put a great effort into devising robust solvers for this class of problems. An example of a highly challenging problem attracting significant attention is the distributed control of incompressible viscous fluid flow problems. In this case, the physics is described by the incompressible Navier--Stokes equations. Since the PDEs given in the constraints are non-linear, in order to obtain a solution of Navier--Stokes control problems one has to iteratively solve linearizations of the problems until a prescribed tolerance on the non-linear residual is achieved.

In this talk, we present efficient and robust preconditioned iterative methods for the solution of the stationary incompressible Navier--Stokes control problem, when employing a Gauss--Newton linearization of the first-order optimality conditions. The iterative solver is based on an augmented Lagrangian preconditioner. By employing saddle-point theory, we derive suitable approximations of the $(1,1)$-block and the Schur complement. Numerical experiments show the effectiveness and robustness of our approach, for a range of problem parameters.

Speaker: Santolo Leveque (University of Houston)

5 Preconditioning Strategies for Generalized Krylov Subspace Methods for $\ell_p-\ell_q$ Minimization

In this talk I will present an extension of the Maximization-Minimization Generalized Krylov Subspace (MM-GKS) method for solving $\ell^p-\ell^q$ minimization problems, as proposed in [1], by introducing a right preconditioner aimed at accelerating convergence without compromising the quality of the computed solution. The original MM-GKS approach relies on iterative reweighting and projection onto subspaces of increasing dimensions, enabling efficient resolution of minimization problems. Our enhanced method leverages a carefully designed regularizing preconditioner, inspired by Iterated Tikhonov regularization, to address the inherent ill-conditioning of the problem. We demonstrate that our preconditioned MM-GKS method preserves the stability and accuracy of the original MM-GKS method, as validated by numerical results in image deblurring, showing significant reductions in CPU time.

This work is a collaboration with A. Buccini (University of Cagliari), M. Donatelli (University of Insubria), and L. Reichel (Kent State University).

Reference:

[1] A. Lanza, S. Morigi, L. Reichel, F. Sgallari, A generalized Krylov subspace method for $\ell_p-\ell_q$ minimization. SIAM Journal on Scientific Computing, (2015), 37(5), S30-S50.

Speaker: Marco Ratto (Università degli studi dell'Insubria)

10:40 Coffee Break

Morning Session II

6 Cross-eyed preconditioning for normal equations arising from PDEs

We explore preconditioning strategies for solving non-symmetric linear systems arising from PDE discretizations. In particular, we analyze the performance of Krylov subspace methods when applied to normal equations. We introduce the idea of cross-eyed preconditioning and discuss its advantages in cases involving convection-diffusion equations. By leveraging physical insights from PDEs, we propose novel preconditioners and evaluate their effectiveness through numerical experiments. Our results demonstrate that the physical nature of the normal equations can be leveraged to construct novel preconditioners, whose convergence properties are better understood compared to GMRES.

Speaker: Umberto Zerbinati (University of Oxford)

7 Matrix-oriented discretisations of phase-field and high-order diffusion problems

We consider a large class of arbitrarily high-order evolutionary PDEs, including phase-field models and polyharmonic reaction-diffusion problems on rectangular domains with Neumann boundary conditions. We propose a matrix-oriented approach that, after full discretisation, requires the solution of a sequence of algebraic matrix equations of Lancaster type. To this end, we apply lumped finite elements in space and a selection of time solvers. Hence, we solve the Allen-Cahn (2nd order), Cahn-Hilliard (4th order), polyharmonic diffusion problems (2m-th order, m $\in \mathbb{N}$), and polyharmonic reaction-diffusion systems. Since these models exhibit high-rank spatial solutions, we employ the so called reduced method in the spectral space. Since the Neumann matrices arising from lumped FEM are slightly different from the corresponding well-known matrices arising from central finite differences, we derive ex-novo the needed eigendecompositions in closed form. Moreover, we provide analytical bounds for the condition number of the vectorised counterpart to show that ill-conditioning can be avoided by our matrix approach. Numerical experiments illustrate the computational performances of the matrix-oriented discretisation in terms of execution time, accurate simulation of high-rank solutions such as Turing patterns and good conditioning.

Speaker: Massimo Frittelli (University of Salento, Department of Mathematics and Physics "E. De Giorgi")

8 Scalable approximation and solvers for ionic electrodiffusion in cellular geometries

The activity and dynamics of excitable cells are fundamentally regulated and moderated by extracellular and intracellular ion concentrations and their electric potentials. The increasing availability of dense reconstructions of excitable tissue at extreme geometric detail pose a new and clear scientific computing challenge for computational modelling of ion dynamics and transport. We present a scalable numerical algorithm [1] for solving the time-dependent and nonlinear KNP-EMI (Kirchhoff-Nernst-Planck Extracellular-Membrane-Intracellular) equations describing ionic electrodiffusion for excitable cells with an explicit geometric representation of intracellular and extracellular compartments and interior interfaces. We also introduce and specify a set of model scenarios of increasing complexity suitable for benchmarking. Our solution strategy is based on an implicit-explicit discretization and linearization in time, a mixed finite element discretization of ion concentrations and electric potentials in intracellular and extracellular domains, and an algebraic multigrid-based, inexact block-diagonal preconditioner for GMRES. Numerical experiments with up to $10^8$ unknowns per time step and up to 256 cores demonstrate that this solution strategy is robust and scalable with respect to the problem size, time discretization and number of cores.

1. Benedusi, Pietro, et al. "Scalable approximation and solvers for ionic electrodiffusion in cellular geometries." SIAM Journal on Scientific Computing 46.5 (2024): B725-B751.

Speaker: Pietro Benedusi

9 Blocking structures, approximation, and preconditioning

Large block-structured matrices with Toeplitz-type blocks of different sizes frequently arise in various applications, but pose computational issues when solving the associated linear systems. In our setting, the matrices $A_n$ are composed of (block rectangular) unilevel Toeplitz blocks defined by rectangular $s \times t$ matrix-valued generating functions. Under mild assumptions on the block dimensions, the asymptotic distribution of the singular values of the associated matrix sequences is recently known. Moreover, when the singular value symbol is Hermitian, the spectral symbol coincides with the singular value symbol.
Based on these theoretical findings, we develop a method to construct simplified block matrices that approximate the original matrices. These simplified matrices offer two key advantages: (a) they maintain the same singular value distributions as $\{A_n\}_n$; (b) they enable the solution of linear systems in $O(n \log n)$ arithmetic operations.
In this way, we propose a natural preconditioning strategy for linear systems with coefficient matrix $A_n$. We provide detailed singular value and spectral analyses of the preconditioned matrix sequences, and validate our approach through numerical experiments concerning the convergence of various (preconditioned) Krylov solvers

Speaker: Valerio Loi (Università degli studi dell'Insubria)

10 Eigenvalues and singular values of variable Toeplitz matrices and matrix-sequences, with application to variable two-step BDF approximations to parabolic equations

In the present talk we consider a general class of varying Toeplitz matrix sequences and we prove that they belong to the maximal *-algebra of GLT matrix sequences. We then examine specific structures within the class stemming from variable two-step backward difference formulae (BDF) approximations of parabolic equations as treated in previous works for studying the stability.

For both the general and the specific case we provide the GLT symbol and based on this, a study on the singular values and the eigenvalues of the related matrix sequences. In addition we provide a decomposition in at most rank-2 nonnegative definite matrices for the symmetrized case.

Speaker: NIKOS BARAKITIS (University of Insubria)

12:50 Lunch Break

The lunch break is offered by the conference and will be held at

Poster Blitz

Conveners: Aaron Melman (Santa Clara University), Alberto Bucci (Charles University), Alessandro Buccini (University of Cagliari), Alice Cortinovis (Università di Pisa), Andrea Adriani (University of Roma Tor Vergata), Andrea Benvenuti (Università di Pisa), Angelo Alberto Casulli (Gran Sasso Science Institute), Antonio Cicone (Università degli Studi dell'Aquila), Arturo De Marinis (Gran Sasso Science Institute), Fabio Durastante (Università di Pisa), Francesco Hrobat (Scuola Normale Superiore), Luisa Fermo (University of Cagliari), Mehdi Najafi Kalyani, Milo Viviani (Scuola Normale Superiore), Muhammad Faisal Khan, Nikita Deniskin (Scuola Normale Superiore), Patricia Diaz de Alba (University of Salerno), Piero Deidda (Scuola Normale Superiore, Pisa - Gran Sasso Science Institute, L'Aquila), Stefano Sicilia (Gran Sasso Science Institute)

Afternoon Session I

11 Preconditioned Low-Rank Riemannian Optimization for Symmetric Positive Definite Linear Matrix Equations

In this talk we focus on the numerical solution of large-scale symmetric positive definite matrix equations of the form $A_1XB_1^\top + A_2XB_2^\top + \dots + A_\ell X B_\ell^\top = F$, which arise from discretized partial differential equations and control problems. These equations frequently admit low-rank approximations of the solution $X$, particularly when the right-hand side matrix $F$ has low rank. For cases where $\ell \leq 2$, effective low-rank solvers have been developed, including Alternating Direction Implicit (ADI) methods for Lyapunov and Sylvester equations. For $\ell > 2$, several existing methods try to approach $X$ through combining a classical iterative method, such as the conjugate gradient (CG) method, with low-rank truncation. In this talk, we consider a more direct approach that approximates $X$ on manifolds of fixed-rank matrices through Riemannian nonlinear CG. A significant challenge is the integration of effective preconditioners into this first-order Riemannian optimization method. We propose novel preconditioning strategies, including a change of metric in the ambient space, preconditioning the Riemannian gradient, and a variant of ADI on the tangent space. Along with a rank adaptation strategy, the proposed method demonstrates competitive performance on a range of representative examples.

References
[1] Ivan Bioli, Daniel Kressner, and Leonardo Robol. Preconditioned Low-Rank Riemannian Optimization for Symmetric Positive Definite Linear Matrix Equations. Aug. 29, 2024. arXiv:2408.16416[math.NA].

Speaker: Ivan Bioli (Università di Pavia)

12 Tensor-Oriented LSQR for Tensor Least Squares Problems

We are interested in the numerical solution of the multiterm tensor least squares problem
$$ \min_{\mathcal{X}} \| \mathcal{F} - \sum_{i =0}^{\ell} \mathcal{X} \times_1 A_1^{(i)} \times_2 A_2^{(i)} \cdots \times_d A_d^{(i)} \|_F, $$ where $\mathcal{X}\in\mathbb{R}^{m_1 \times m_2 \times \cdots \times m_d}$, $\mathcal{F}\in\mathbb{R}^{n_1\times n_2 \times \cdots \times n_d}$ are tensors with $d$ dimensions (or modes), while $A_j^{(i)}\in\mathbb{R}^{n_j \times m_j}$ for every $i=1,\ldots,\ell$ and $j=1,\ldots,d$. The symbol $\times_j$ with $j=1,\ldots,d$ denotes the $j$-mode product of a tensor times a matrix. We are interested in the Tucker and Tensor-Train formats. Least squares tensor (and matrix) formulations have emerged in recent literature from different applications, including the numerical solution of PDEs, data science problems such as dictionary learning [1,2], control systems. The problem is challenging for the absence of direct methods that can efficiently handle multiple addends in tensorial form. In particular, the problem is extremely memory demanding even for modest values of each mode size, for large $\ell$.

In our presentation we will propose an implementation of Truncated Tensor-oriented LSQR, first introduced in vector form in [3] and studied in the matrix version in [6], and illustrate the potential of these new approaches on problems stemming from the
discretization of multidimensional PDEs and from tensor Dictionary Learning.

1. C. F. Dantas, J. E. Cohen, and R. Gribonval. "Learning Tensor-structured Dictionaries with Application to Hyperspectral Image Denoising". In: 27th European Signal Processing Conference, EUSPICO 2019, Coruna, Spain, September 2-6, 2019. IEEE, 2019, pp1-5.

2. C. F. Dantas, M. N. Da Costa, and R. da Rocha Lopes. "Learning dictionaries as sum of Kronecker products". In: IEEE Signal Processing Letters 24.5 (2017), pp. 559-563.

3. Christopher C. Paige, and Micheal A. Saunders. "LSQR: An algorithm for sparse linear equation and sparse least squares". In: ACM Transactions on Mathematical Software (TOMS) 8.1 (1982), pp. 43-71.

4. I. Oseledets. "Tensor-Train Decomposition". In: SIAM J. Scientific Computing 33 (Jan. 2011), pp. 2295-2317.

5. H. Al Daas et al. "Randomized Algorithms for Rounding in the Tensor-Train Format". In: \textit{SIAM Journal on Scientific Computing} 45.1 (2023), A47-A95.

6. V. Simoncini, L. Piccinini. "Truncated LSQR for matrix least squares problems". pp. 1-22, February 2024. To appear in Computational Optimization and Applications (COAP).

Speaker: Lorenzo Piccinini (Università di Bologna)

13 Tree Tensor Network Nyström

In this work, we present the tree tensor network Nyström (TTNN), an algorithm that extends recent research on streamable tensor approximation, such as for Tucker or tensor-train formats, to the more general tree tensor network format, enabling a unified treatment of various existing methods. Our method retains the key features of the generalized Nyström approximation for matrices, i.e. it is randomized, single-pass, streamable, and cost-effective. Additionally, the structure of the sketching allows for parallel implementation. We provide a deterministic error bound for the algorithm and, in the specific case of Gaussian dimension reduction maps, also a probabilistic one. We also introduce a sequential variant of the algorithm, referred to as sequential tree tensor network Nyström (STTNN), which offers better performance for dense tensors. Furthermore, both algorithms are well-suited for the recompression or rounding of tensors in the tree tensor network format. Numerical experiments highlight the efficiency and effectiveness of the proposed methods.

Speaker: gianfranco verzella (University of Geneva)

14 Riemann-Oracle: a general-purpose Riemannian optimizer to solve nearness problems in matrix theory

We propose an extremely versatile approach to address a large family of matrix nearness problems, possibly with additional linear constraints. Our method is based on splitting a matrix nearness problem into two nested optimization problems, of which the inner one can be solved either exactly or cheaply, while the outer one can be recast as an unconstrained optimization task over a smooth real Riemannian manifold. We observe that this paradigm applies to many matrix nearness problems of practical interest appearing in the literature, thus revealing that they are equivalent in this sense to a Riemannian optimization problem. We also show that the objective function to be minimized on the Riemannian manifold can be discontinuous, thus requiring regularization techniques, and we give conditions for this to happen. Time permitting, we demonstrate the practical applicability of our method by implementing it for a number of matrix nearness problems that are relevant for applications and are currently considered very demanding in practice. Among them, we recall the problem of finding the distance to singularity for matrix polynomials, the computation of the approximate gcd between to scalar polynomials and the nearest matrix with at least one eigenvalue
in a prescribed closed set, also known as distance to instability.

Speaker: Miryam Gnazzo (Gran Sasso Science Institute)

15 Inexact alternating projection methods with application to matrix completion

Consider the affine rank minimization problem [6] defined as
\begin{equation}
\begin{aligned}
\min & \quad \text{rank}(X) \
\textrm{s.t.} & \quad \mathcal{C}_\mathcal{A} (X) = b,
\end{aligned}
\tag{1}
\end{equation}
where $X \in \mathbb{R}^{n_1 \times n_2}$ is a matrix, and linear map $\mathcal{C}_\mathcal{A} : \mathbb{R}^{n_1 \times n_2} \rightarrow \mathbb{R}^m$ and vector $b \in \mathbb{R}^m$ are given.
This problem appears in diverse fields including system identification and control, Euclidean embedding, collaborative filtering, and image reconstruction.

Given two sets $A$ and $B$ with $A \cap B \neq \varnothing$, a feasibility problem consists in
\begin{equation}
\textrm{find} \ x \in A \cap B.
\tag{2}
\end{equation}

Problem (1) is a nonconvex optimization problem due to rank function, and it can be reformulated as (2) considering:
\begin{align}
A &= { X \ | \ \mathcal{C}_\mathcal{A} (X) = b } \
B &= { X \ | \ \text{rank}(X) \le r },
\end{align}
where the rank $r$ is known in advance.

In the 50s, the Alternating Projection Method (APM) was proposed by von Neumann [4] for solving (2) when $A,B$ are linear subspaces. It consists in the simple idea of projecting alternately onto $A$ and $B$ until the generated sequence converges to a point in the intersection. Formally, given $y^0 \in B$, we can define APM's iteration as follows:
\begin{equation}
\begin{cases}
x^k &\in P_A (y^{k-1}) \
y^k &\in P_B (x^{k})
\end{cases}
\quad k = 1, \dots
\tag{3}
\end{equation}
where $P_A, P_B$ are projection operators onto $A$ and $B$, respectively.

The method converges when $A,B$ are convex [1]. When one of the sets is nonconvex, APM converges only locally and further hypothesis on the sets and their intersection are required [5]. The difficulty comes from the projection onto the nonconvex set, which usually is unknown in closed form. To overcome this issue, inexact AP methods have been also proposed [2].

Many works in this field of literature, including those cited above, agree in considering the particular case of rank set an ``easy'' problem, despite the nonconvexity of $\mathcal{C}_r$. Indeed, from the Eckart-Young theorem, the projection of a matrix $X$ onto $\mathcal{C}_r$ is the SVD of $X$ truncated to the first $r$ singular values. Then, $P_{\mathcal{C}_r}$ is considered to be known in exact and closed form [2].

From a numerical perspective, however, it is well known that the SVD of a matrix, as well as its truncated version, is computed through iterative procedures that can be computationally expensive, especially when dealing with large matrices and that unavoidably truncation errors that, to the best of our knowledge, have not been considered in previous studies of methods dealing with the rank set, where the projection is assumed to be exact.

Additionally, in a more general setting, methods employing inexact projections adopt criteria that often involve the geometry of the sets without providing a thorough analysis of their implementation and efficiency from a numerical standpoint.

In this work, we introduce the inexact Regularized APM (iRAPM) that accounts for the inaccuracy of the numerical approximation of the truncated SVD and which is based on inexactness criteria that can be efficiently implemented.

Our inexactness criteria rely on a (problem dependent) subprocedure to replace projection onto the ``hard'' set. In the particular case of (2) with the rank set, we show that the procedure can be defined starting from the Lanczos algorithm for truncated SVD. Indeed, we show that the new criteria can be checked using byproduct quantities of the Lanczos method with negligible additional cost.

For our experiments, we used Lanczos implementation in PROPACK by [3] and we tested iRAPM on the Matrix Completion problem. Results show that iRAPM can be faster that classical APM retaining same reconstruction error.

Reference

1. Heinz H Bauschke and Jonathan M Borwein. On projection algorithms for solving convex feasibility
   problems. SIAM review, 38(3):367–426, 1996.
2. Dmitriy Drusvyatskiy and Adrian S Lewis. Local linear convergence for inexact alternating pro-
   jections on nonconvex sets. Vietnam Journal of Mathematics, 47:669–681, 2019.
3. Rasmis Munk Larsen. Lanczos bidiagonalization with partial reorthogonalization. Dept. Computer
   Science, Aarhus University, Technical report, DAIMI, PB(537), 1998.
4. John Von Neumann. Functional Operators (AM-22), Volume 2: The Geometry of Orthogonal
   Spaces. (AM-22). Princeton University Press, 1950.
5. D. Noll and A. Rondepierre. On Local Convergence of the Method of Alternating Projections.
   Found. Comput. Math., 16:425–455, 2016.
6. Benjamin Recht, Maryam Fazel, and Pablo A. Parrilo. Guaranteed minimum-rank solutions of
   linear matrix equations via nuclear norm minimization. SIAM Review, 52(3):471–501, 2010.

Speaker: Mr Mattia Silei (University of Florence)

Poster Session: Poster

16 Anti-Gauss cubature rules with applications to Fredholm integral equations on the square

We introduce the anti-Gauss cubature rule for approximating integrals defined on the square whose integrand function may have algebraic singularities at the boundaries. An application of such a rule to the numerical solution of Fredholm integral equations of the second-kind is also explored. The stability, convergence, and conditioning of the proposed Nystr\"om-type method are studied. The numerical solution of the resulting dense linear system is also investigated and several numerical tests are presented.

Speaker: Luisa Fermo (University of Cagliari)

17 Approximation properties of neural ODEs

We study the universal approximation property (UAP) of shallow neural networks whose activation function is defined as the flow of a neural ODE. We prove the UAP for the space of such shallow neural networks in the space of continuous functions. In particular, we also prove the UAP with the weight matrices constrained to have unit norm.

Furthermore, in [1] we are able to bound from above the Lipschitz constant of the flow of the neural ODE, that tells us how much a perturbation in input is amplified or shrunk in output. If the upper bound is large, then so it may be the Lipschitz constant, leading to the undesirable situation where certain small perturbations in input cause large changes in output. Therefore, in [2] we compute a perturbation to the weight matrix of the neural ODE such that the flow of the perturbed neural ODE has Lipschitz constant bounded from above as we desire. This leads to a stable flow and so to a stable shallow neural network.

However, the stabilized shallow neural network with unit norm weight matrices does not satisfy the universal approximation property anymore. Nevertheless, we are able to prove approximation bounds that tell us how poorly and how accurately a continuous target function can be approximated by the stabilized shallow neural network.

The results presented during this talk are being collected in [3].

1. N. Guglielmi, A. De Marinis, A. Savostianov, and F. Tudisco, Contractivity of neural ODEs: an eigenvalue optimization problem, arXiv preprint arXiv:2402.13092, 2024.
2. A. De Marinis, N. Guglielmi, S. Sicilia, and F. Tudisco, Stability of neural ODEs by a control over the expansivity of their flows, work in progress.
3. A. De Marinis, D. Murari, E. Celledoni, N. Guglielmi, B. Owren, and F. Tudisco, Approximation properties of neural ODEs, work in progress.

Speaker: Arturo De Marinis (Gran Sasso Science Institute)

18 DEIM-based hyper-reduction of metriplectic systems

Metriplectic or GENERIC systems describe complex dynamical systems that blend Hamiltonian and dissipative dynamics, i.e. systems that conserve energy and produce entropy. Maintaining the separation between Hamiltonian and dissipative part in the construction of approximate models is a challenging task. In this work, we employ a matrix-based implementation of the Discrete Empirical Interpolation Method (DEIM) to perform hyper-reduction of the Hamiltonian term in metriplectic systems. This approach is integrated with a structure-preserving reduced-order algorithm developed in [A. Gruber, M. Gunzburger, L. Ju, Z. Wang, CMAME 2023] that guarantees the retention of the properties needed for energy conservation and entropy production. We present some numerical results to test the performances of the proposed method. Our results indicate that, while maintaining similar accuracy and structure-preservation, the evaluation of the reduced model is computationally cheaper when using the DEIM-enhanced method.

Speaker: Andrea Benvenuti (Università di Pisa)

19 Geometric means of more than two matrix sequences in the case of hidden (asymptotic) structures

We consider the spectral distribution of the geometric mean of two or more Hermitian positive definite (HPD) matrix-sequences, under the assumption that all input matrix-sequences belong to the same Generalized Locally Toeplitz (GLT) $*$-algebra. As expected, the numerical experiments show that the geometric mean of $k$ positive definite GLT matrix-sequences forms a new GLT matrix-sequence, with the GLT symbol given by the geometric mean of the individual symbols. While the result is plain for $k=2$, it is highly non trivial for $k>2$, due to the limit process for defining the geometrc mean and due to the lack of a closed form expression. Theoretical tools for handling the difficult case are discussed.

1. G. Barbarino, C. Garoni, S. Serra-Capizzano, Block generalized locally Toeplitz sequences: theory and applications in the unidimensional case, Electr. Trans. Numer. Anal. 53 (2020), pp. 28–112.
2. G. Barbarino, C. Garoni, S. Serra-Capizzano, Block generalized locally Toeplitz sequences: theory and applications in the multidimensional case, Electr. Trans. Numer. Anal. 53 (2020), pp. 113–216.
3. D.A. Bini, B. Iannazzo, Computing the Karcher Mean of Symmetric Positive Definite Matrices, Linear Algebra and its Applications, 438 (2013), pp. 1700–1710.
4. C. Garoni, S. Serra-Capizzano, Generalized locally Toeplitz sequences: theory and applications. Vol. I, Springer, Cham, 2017.
5. C. Garoni, S. Serra-Capizzano, Generalized locally Toeplitz sequences: theory and applications. Vol. II, Springer, Cham, 2018.
6. M.F. Khan, S. Serra-Capizzano, Geometric means of more than two matrix-sequences in the case of hidden (asymptotic) structures, preprint 2024.

Speaker: Muhammad Faisal Khan

20 Higher Order Data Completion for Undersampled Raster Scans in Spectromicroscopy

In this poster we present a novel iterative algorithm for low rank tensor completion, tailored to the recovering the missing entries in undersampled X-ray spectromicroscopy data, which are used to study material distributions. Compared to established techniques that rely on data matricizations and low-rank matrix completion, the new method allows the selection of robust sampling patterns, tensor multi-rank and undersampling ratio, while minimising the impact of undersampling on the data analysis. Results obtained on real data will be illustrated.

Speaker: Silvia Gazzola

21 Iterated $\ell^2-\ell^q$ regularization

In numerous fields of science and engineering we are faced with problems
of the form $$A{x}+{\eta}=b^\delta,$$ where $A\in\mathbb{R}^{m\times n}$ is a large and severely ill-conditioned matrix, i.e., such that its singular values decay rapidly to zero without any significant gap between the smallest ones, ${x}\in\mathbb{R}^n$ is the unknown quantity we wish to recover, ${\eta}\in\mathbb{R}^m$ collects the unknown measurement and discretization errors, and is often referred to as noise, and ${b}^\delta\in\mathbb{R}^m$ represents the measured and perturbed data. We will assume that the noise that corrupts the data is white Gaussian. Problems of this form are referred to as discrete ill-posed inverse problems. Due to the ill-conditioning of $A$ and the presence of the noise ${\eta}$, the naive solution of the system of equations above is a poor approximation of the desired one. To reduce the sensitivity of the problem to the presence of noise, we consider the $\ell^2-\ell^q$ regularization $$\arg\min_{{x}}\frac{1}{2}\left\|A{x}-{b}^\delta\right\|_2^2+\frac{\mu}{q}\left\|L{x}\right\|_q^q,$$ where $0<q\leq2$, $\mu>0$ is the so-called regularization parameter, $L\in\mathbb{R}^{s\times n}$ is referred to as regularization operator, and $||z||_q^q=\sum_{i=1}^{n}\lvert z_i \rvert^q.$ We refer to latter quantity as $\ell^q$-norm for any $0<q\leq 2$, even though for $q<1$ it is not a norm as it does not satisfy the triangular inequality; see, e.g., [2] for a discussion on $\ell^p-\ell^q$ regularization and [1] for a software implementation. To further improve the quality of the computed solutions and to provide a theoretically sound convergence analysis, we consider the refinement strategy $$\begin{cases} {r}^{(k)}={b}^\delta-A{x}^{(k)},\\ {h}^{(k)}= \arg\min_{{h}}\frac{1}{2}\left\|A{h}-{r}^{(k)}\right\|_2^2+\frac{\mu}{q}\left\|L{h}\right\|_q^q,\\ {x}^{(k+1)}={x}^{(k)}+{h}^{(k)}. \end{cases}$$ We show that, under reasonable assumptions, the iterations are stopped after finitely many steps, if the Discrepancy Principle is employed as stopping criterion, and that the resulting method is a regularization method. Selected numerical results show the good performances of the proposed approach.
References

1. A. Buccini and L. Reichel, Software for limited memory restarted $\ell^p-\ell^q$ minimization methods using generalized Krylov subspaces, Electronic Transactions on Numerical Analysis, 61, 66-91 (2024).
2. A. Buccini, M. Pragliola, L. Reichel, and F. Sgallari, A comparison of parameter selection rules for $\ell^p-\ell^q$ minimization, Annali dell'Università di Ferrara, 68, 441-463 (2022).
   

Speaker: Alessandro Buccini (University of Cagliari)

22 Lanczos with compression for symmetric eigenvalue problems

This work focuses on computing a few smallest or largest eigenvalues and their corresponding eigenvectors of a large symmetric matrix $A$.

Krylov subspace methods are among the most effective approaches for solving large-scale symmetric eigenvalue problems. Given a starting vector $\mathbf{q}_1$, the $M$-th Krylov subspace is generated by repeatedly multiplying a general square matrix $A$ with $\mathbf{q}_1$:

$$\mathcal{K}_{M}(A,\mathbf{q}_1) := \mathrm{span} \{\mathbf{q}_1, A\mathbf{q}_1, \dotsc, A^{M-1}\mathbf{q}_1\}.$$ An orthonormal basis $Q_{M}$ for the Krylov subspace $\mathcal{K}_{M}(A,\mathbf{q}_1)$ is constructed using the Lanczos process. The Rayleigh--Ritz process is then applied to $Q_{M}$ to extract eigenvalue and eigenvector approximations, known as Ritz values and Ritz vectors. A significant challenge of Krylov subspace methods is the need to store $Q_{M}$. For slow convergence (with respect to $M$), available memory may be exhausted before achieving satisfactory approximations. Popular algorithms for large-scale eigenvalue problems address this issue by combining the Lanczos process with restarting. As an alternative to restarting, this work proposes a novel compression approach based on Rational Krylov subspaces associated with small matrices, limiting the Lanczos method's memory requirements. To provide intuition for our approach, suppose the spectrum of $A$ is ordered such that: $$\lambda_1 \le \cdots \le \lambda_m < \tau < \lambda_{m+1} \le \cdots \le \lambda_n,$$ where a shift $\tau$ separates the smallest $m \ll n$ eigenvalues to be computed from the rest. Let $\chi_\tau(x)$ denote the step function that equals $1$ for $x < \tau$ and $0$ otherwise. In theory, $\chi_\tau(A)\mathbf{q}_1$ serves as a perfect filter: the Lanczos method with this starting vector would yield the exact eigenvalues $\lambda_1, \ldots, \lambda_m$ after $m$ steps for a generic $\mathbf{q}_1$. In practice, however, evaluating $\chi_\tau(A)\mathbf{q}_1$ exactly is computationally infeasible. Instead, a rational approximation of $\chi_\tau$ can be used to compress $\mathcal{K}_{M}(A,\mathbf{q}_1)$. This approach parallels the developments in [2], which employs a rational approximation of a general function $f$ to compress the Lanczos method for approximating $f(A)\mathbf{b}$. To make this compression practical for eigenvalue problems, several additional components are required. Among these, a reorthogonalization procedure is essential, as numerical loss of orthogonality can lead to a failure in convergence.

We present a series of numerical experiments involving matrices from various applications. These experiments demonstrate that, in terms of matrix-vector products, our new method is consistently competitive with or superior to the Krylov--Schur method, often achieving significant improvements.

References

1. A. C., D. Kressner, N. Shao. Lanczos with compression for symmetric
   eigenvalue problems. In preparation.

2. A. C., I. Simunec. A low-memory Lanczos method with rational Krylov
   compression for matrix functions. arXiv, 2024.

Speaker: Angelo Alberto Casulli (Gran Sasso Science Institute)

23 Lanczos with compression for symmetric Lyapunov equations

In this work, we present a low-memory variant of the Lanczos algorithm for the solution of the Lyapunov equation
\begin{equation}
AX + XA = \mathbf{c}\mathbf{c}^T, \tag{1}
\end{equation}
where $A$ is a large-scale symmetric positive-definite matrix and $\mathbf{c}$ is a vector.

The classical Lanczos method consists in building an orthonormal basis $Q_M$ for the polynomial Krylov subspace
$$\mathcal{K}_M(A,\mathbf{c}) = span(\mathbf{c},A\mathbf{c}, \dots, A^{M-1}\mathbf{c})$$ and in approximating the solution $X$ with $Q_MY_MQ_M^T$, where $Y_M$ solves the projected equation $$ Q_M^TAQ_M Y_M + Y_MQ_M^TAQ_M = (Q_M^T \mathbf{c}) (Q_M^T \mathbf{c})^T.$$ The Lanczos algorithm often requires a relatively large $M$ to obtain a good approximation of the solution, which can lead to memory issues due to the storage demands of $Q_M$. Furthermore, the solution $X$ can be well approximated by a low-rank matrix, whose rank is significantly smaller than $M$, i.e. the dimension of the polynomial Krylov subspace. An alternative approach is to use a rational Krylov subspace instead of a polynomial one. Using the Zolotarev poles as the poles of the rational Krylov subspace, it is possible to approximate the solution $X$ by a low-rank matrix with the guarantee that the residual has norm smaller than a prescribed quantity ([$2$]). The rank of the computed approximate solution is usually close to the numerical rank of the real solution. The main drawback is that this method requires solving multiple shifted linear systems involving the matrix $A$, which is prohibitive if $A$ is large. Mimicking the approach in [$3$], our method employs a polynomial Krylov subspace to approximate the solution of (1) while leveraging rational Krylov subspaces associated with small matrices to compress the Lanczos basis. This method accesses $A$ only through matrix-vector products and requires the storage of only a few vectors from the polynomial Krylov subspace instead of the entire Lanczos basis, producing an approximate solution whose rank is independent of the dimension of the involved polynomial Krylov subspace. The computational cost of the proposed algorithm is dominated by the construction of the Lanczos basis, and the compression steps do not require additional matrix-vector products involving $A$. Furthermore, theoretical results demonstrate that the algorithm achieves an approximation error comparable to that of the standard Lanczos algorithm, with an additional error term that can be bounded a priori using Zolotarev numbers. In practice, this additional error is negligible compared to the Lanczos error.

Numerical experiments show that the behavior of the proposed algorithm is comparable to that of the Lanczos algorithm without reorthogonalization, both in terms of matrix-vector products and quality of the approximated solution. Comparisons with existing low-memory variants of the Lanczos method demonstrate competitive performance in terms of accuracy, computational cost, and runtime.

1. A. A. Casulli, F. H., D. Kressner, Lanczos with compression for symmetric Lyapunov equations, In preparation.
2. B. Beckermann, An error analysis for rational Galerkin projection applied to the Sylvester
   equation, SIAM J. Numer. Anal., 49 (2011), pp. 2430--2450.
3. A. A. Casulli and I. Simunec. A low-memory Lanczos method with rational
   Krylov compression for matrix functions, arXiv, 2024.

Speaker: Francesco Hrobat (Scuola Normale Superiore)

24 Linear FDEM data inversion via Landweber interations in Banach spaces

The practical motivation for this work is to reconstruct electromagnetic properties of the earth superficial layer using measurements taken above the ground. We approach this problem by inverting frequency-domain electromagnetic data through a well-established linear integral model, employing three different collocation methods to approximate the solution as a linear combination of linearly independent basis functions. This discretization results in a linear system that is highly ill-conditioned. To address this challenge, we apply an iterative regularization technique based on Landweber iterations in Banach spaces, which is particularly effective for reconstructing solutions that exhibit discontinuities or low smoothness—common characteristics in many imaging applications. Several numerical experiments demonstrate the superior performance of our method compared to other regularization techniques.

References:

1. P. Diaz de Alba, C. Estatico, M. Lazzaretti, G. Rodriguez. Linear FDEM subsoil data inversion in Banach spaces. Submitted to Electronic Transactions on Numerical Analysis (2024).

Speaker: Patricia Diaz de Alba (University of Salerno)

25 Minimum-norm solutions of the non-symmetric semidefinite Procrustes problem

Given two matrices $X,B\in \mathbb{R}^{n\times m}$ and a set $\mathcal{A}\subseteq \mathbb{R}^{n\times n}$, a Procrustes problem consists in finding a matrix $A \in \mathcal{A}$ such that the Frobenius norm of $AX-B$ is minimized. When $\mathcal{A}$ is the set of the matrices whose symmetric part is positive semidefinite, we obtain the so-called non-symmetric positive semidefinite Procrustes (NSPSDP) problem. The NSPSDP problem arises in the estimation of compliance or stiffness matrix in solid and elastic structures.
If $X$ has rank $r$, Baghel et al. in 2022 proposed a three-step semi-analytical approach:
(1) construct a reduced NSPSDP problem in dimension $r\times r$,
(2) solve the reduced problem by means of a fast gradient method with a linear rate of convergence, and
(3) post-process the solution of the reduced problem to construct a solution of the larger original NSPSDP problem.
In this talk, we revisit this approach by Baghel et al. and we identify an unnecessary assumption used by the authors leading to cases where their algorithm cannot attain a minimum and produces solutions with unbounded norm. In fact, revising the post-processing phase of their semi-analytical approach, we show that the infimum of the NSPSDP problem is always attained, and we show how to compute a minimum-norm solution. We also prove that the symmetric part of the computed solution has minimum rank bounded by $r$, and that the skew-symmetric part has rank bounded by $2r$. Several numerical examples show the efficiency of this algorithm, both in terms of computational speed and of finding optimal minimum-norm solutions.

Speaker: Stefano Sicilia (Gran Sasso Science Institute)

26 Mittag-Leffler Communicability Distance and shortest geodesics on networks

The Mittag-Leffler functions interpolate between the exponential and the resolvent functions, and they can be used for walk-based indices in networks, allowing to tune the weight of longer closed walks.
Here we study the Communicability based on the Mittag-Leffler function, and define Mittag-Leffler communicability distance functions, which are proved to be Euclidean and spherical. We then geometrize the graphs based on these distances and obtain the corresponding geodesics based on them. We finally analyze these geodesics for different kinds of real-world networks.

Speaker: Nikita Deniskin (Scuola Normale Superiore)

27 On the highly entrywise accurate methods for solution of multilinear PageRank problem

In this talk, we consider the solution of the multilinear PageRank problem, which can be reformulated as a nonlinear system of equations in which the coefficient matrix is an M-matrix. In [1], Alfa et al., demonstrated that the inverse of a nonsingular M-matrix can be determined with a high relative entrywise accuracy by using a triplet representation of the M-matrix through a method known as the GTH algorithm. This variant of Gaussian elimination allows for the computation of a numerical inverse with comparable entrywise relative accuracy. Our approach combines the GTH algorithm with either the Newton method or block iterative schemes (such as block Jacobi, block Gauss-Seidel, and block Successive Over-Relaxation iterations). Specifically, in each iteration of the Newton method, we use the GTH algorithm to compute the inverse. All sub-matrices of the coefficient matrix are nonsingular M-matrices and can also be represented in a triplet format; So, at each step in the block iterative schemes, the GTH algorithm is used for inverse computation, which ensures high relative entrywise accuracy. We focus on the cases where the problem approaches a critical case, potentially causing the coefficient M-matrix become singular. In cases where the M-matrix becomes singular, we consider shifted versions of the algorithms. Our numerical results show the efficiency and reliability of our approach.

References:

1. Alfa, A., Xue, J., & Ye, Q. (2002). Accurate computation of the smallest eigenvalue of a diagonally dominant M-matrix. Mathematics of computation, 71(237), 217-236.

2. Alfa, A. S., Xue, J., & Ye, Q. (2002). Entrywise perturbation theory for diagonally dominant M-matrices with applications. Numerische Mathematik, 90, 401-414.

3. O'Cinneide, C. A. (1996). Relative-error bounds for the LU decomposition via the GTH algorithm. Numerische Mathematik, 73, 507-519.

Speaker: Mehdi Najafi Kalyani (Dipartimento di Informatica, Università di Pisa, Pisa, Italy)

28 On the time-frequency representation of nonstationary signals

The time-frequency analysis of nonlinear and nonstationary processes is, in general, a challenging task. Standard techniques, like short time Fourier transform, and wavelet transform are limited in addressing the problem. An alternative way to tackle the problem is to first decompose the signal into simpler components and then analyze them separately. This is the idea behind what is called the Hilbert Huang Transform (HHT),
published originally in 1998. HHT has had and is still having a big impact on many fields of research.

HHT is made of two stages: the decomposition of the signal into simple oscillatory components, called Intrinsic Mode Functions (IMFs), obtained via the Empirical Mode Decomposition (EMD) method, and the time-frequency representation of these components obtained via the Hilbert Transform. However, the mathematical properties of EMD and its generalizations, like the Ensemble EMD, are still under investigation. For this reason, an alternative technique, called Fast Iterative Filtering (FIF), was proposed a few years ago [1]. Furthermore, the Hilbert Transform, even though mathematically sound, has limitations when applied to real-life data, especially strongly nonstationary or noisy ones.

In this talk, we review the mathematical properties of FIF and the decomposition it produces. In particular, we will talk about the energy preservation and orthogonality of this decomposition. Furthermore, we introduce an alternative method for the time-frequency representation of the IMFs called IMFogram, showing its robustness to noise, locality, and ability to produce, in a fast manner, crisp and artifact-free time-frequency plots [2]. We conclude by presenting applications of this approach to real-life signals and future directions of research.

1. A. Cicone, H. Zhou. Numerical Analysis for Iterative Filtering with New Efficient Implementations Based on FFT. Numerische Mathematik, 147 (1), pages 1-28, 2021.
2. A. Cicone, W. S. Li, H. Zhou. New theoretical insights in the decomposition and time-frequency representation of nonstationary signals: the IMFogram algorithm. Applied and Computational Harmonic Analysis, 71, 101634, 2024.

Speaker: Antonio Cicone (Università degli Studi dell'Aquila)

29 PSCToolkit: Bringing Linear Algebra to the Exascale

The growth of computational capabilities is ushering the era of exascale computing, presenting unprecedented opportunities and challenges for scientific applications. Sparse Linear algebra underpins a wide range of these applications: from simulating the transition to decarbonized energy in Europe to the construction of digital twins of the human body. The demands of exascale systems necessitate innovative approaches to algorithm design, scalability, and performance optimization.

PSCToolkit is a library designed to address these challenges by providing a scalable, portable, and efficient suite of linear algebra tools optimized for exascale architectures in the MPI+X framework. By leveraging hybrid parallelism and advanced GPU acceleration the PSCToolkit suite of libraries strives for optimal performance on heterogeneous systems. The toolkit integrates parallel sparse BLAS with Krylov methods and preconditioners.

In this poster, we present the architecture of PSCToolkit, highlight key design decisions [1]. We describe the most recent We describe the most recent additions from an algorithmic point of view [3] and demonstrate its performance through benchmarks on exascale-class systems [2,3].

1. D’Ambra, Pasqua, Fabio Durastante, and Salvatore Filippone. "Parallel Sparse Computation Toolkit." Software Impacts 15 (2023): 100463.
2. D'Ambra, Pasqua, Fabio Durastante, and Salvatore Filippone. "PSCToolkit: solving sparse linear systems with a large number of GPUs." arXiv preprint arXiv:2406.19754 (2024).
3. D'Ambra, Pasqua, Fabio Durastante, Salvatore Filippone, Stefano Massei, and Stephen Thomas. "Optimal Polynomial Smoothers for Parallel AMG." arXiv preprint arXiv:2407.09848 (2024).

Speaker: Fabio Durastante (Università di Pisa)

30 Quadratic projectable Runge--Kutta methods

Runge-Kutta methods are affine equivariant: applying a method before or after an affine change of variables yields the same numerical trajectory. However, for some applications, one would like to perform numerical integration after a quadratic change of variables. For example, in Lie-Poisson reduction, a quadratic transformation reduces the number of variables in a Hamiltonian system, yielding a more efficient representation of the dynamics. Unfortunately, directly applying a symplectic Runge-Kutta method to the reduced system generally does not preserve its Hamiltonian structure, so many proposed techniques require computing numerical trajectories of the original, unreduced system.

In this talk, we show when a Runge-Kutta method in the original variables descends to a numerical integrator expressible entirely in terms of the quadratically transformed variables. In particular, we show that symplectic diagonally implicit Runge-Kutta (SyDIRK) methods, applied to a quadratic projectable vector field, are precisely the Runge-Kutta methods that descend to a method (generally not of Runge-Kutta type) in the projected variables. We illustrate our results with several examples in both conservative and non-conservative dynamics.

Speaker: Milo Viviani (Scuola Normale Superiore)

31 Randomized Schatten norm estimation

In this talk, we analyze the variance of a stochastic estimator for computing Schatten norms of matrices. The estimator extracts information from a single sketch of the matrix, that is, the product of the matrix with a few standard Gaussian random vectors. While this estimator has been proposed and used in the literature before, the existing variance bounds are often pessimistic. Our work provides a sharper upper bound on the variance and we also give estimates of the variance that work well for numerically low-rank matrices. Our theoretical findings are supported by numerical experiments, demonstrating that the new bounds are significantly tighter than the existing ones.

Speaker: Alice Cortinovis (Università di Pisa)

32 Randomized sketched TT-GMRES for linear systems with tensor structure

We consider the use of sketching to improve the performances and usability of Krylov tensor methods; in this talk, we will consider the case study of TT-GMRES, and describe the many possibilities that arise where sketching is helpful to improve performances and accuracy, often in somewhat unexpected ways. Our ideas build upon sketching in GMRES and streaming low-rank approximation algorithms that were recently introduced. This is joint work with Leonardo Robol (University of Pisa) and Davide Palitta (University of Bologna).

1. A. B., D. Palitta and L. Robol, 2024. Randomized sketched TT-GMRES for linear systems with tensor structure. arXiv preprint, arXiv:2409.09471.

Speaker: Dr Alberto Bucci (Charles University)

33 Revisiting the notion of approximating class of sequences for handling approximated PDEs on moving or unbounded domains

In the current talk we consider sequences of matrix-sequences $\{\{B_{n,t}\}_n\}_t$ of increasing sizes depending on a parameter $n$ and equipped with an additional parameter $t>0$. For every fixed $t>0$, we assume that each $\{B_{n,t}\}_n$ possesses a canonical spectral/singular values symbol $f_t$ defined on a measurable set $D_t\subset \mathbb{R}^{d}$ of finite measure, $d\ge 1$. Furthermore, we assume that $ \{ \{ B_{n,t}\}_n\}_t $ is an approximating class of sequences (a.c.s.) for a matrix-sequence $ \{ A_n \}_n $ and that $ \bigcup_{t > 0} D_t = D $ with $ D_{t + 1} \supset D_t $. Under such assumptions and via the notion of a.c.s, we prove results on the canonical distributions of $ \{ A_n \}_n $, whose symbol, when it exists, can be defined on the, possibly unbounded, domain $D$ of finite, or even infinite, measure.

We then extend the concept of a.c.s. to the case where the approximating sequence $ \{\{ B_{n,t}\}_n\}_t $ has possibly a different dimension than the one of $ \{ A_n\}_n $. This concept seems to be particularly natural when dealing, e.g., with the approximation both of a partial differential equation (PDE) and of its (possibly unbounded, or moving) domain $D$, using an exhausting sequence of domains $\{ D_t \}_t$.

The theory in the present talk seems to have promising developments both in terms of theoretical results and, particularly, in terms of numerical applications. Some of these possible developments will be discussed in the talk.

Few numerical tests and visualizations are presented in order to corroborate the theoretical aspects.

1. Adriani, A., Furci, I., Garoni, C., Serra-Capizzano, S.: Block multilevel structured matrix-sequences and their spectral and singular value canonical distributions: general theory and few emblematic applications. Preprint (2024).

2. Adriani, A., Schiavoni-Piazza, A.J.A., Serra-Capizzano, S.: Block structures, g.a.c.s. approximation, and distributions. Special Volume in Memory of Prof. Nikolai Vasilevski (2024), to appear.

3. Adriani, A., Schiavoni-Piazza, A.J.A., Serra-Capizzano, S., Tablino-Possio, C.: Revisiting the notion of approximating class of sequences for handling approximated PDEs on moving or unbounded domains. Electronic Transactions on Numerical Analysis (2024), to appear.

4. Barakitis, N., Ferrari, P., Furci, I., Serra-Capizzano, S.:
   An extradimensional approach for distributional results: the case of $2\times 2$ block Toeplitz structures, {\it Springer Proceedings on Mathematics and Statistics} (2024), to appear.

5. Furci, I., Adriani, A., Serra-Capizzano, S.: Block structured matrix-sequences and their spectral and singular value canonical distributions: a general theory. In arXiv:2409.06465 (2024).

Speaker: Andrea Adriani (University of Roma Tor Vergata)

34 SCALAR POLYNOMIAL ZERO BOUNDS FROM MATRIX POLYNOMIALS

Useful properties of scalar polynomials can be obtained by embedding them into the larger framework of their generalization to matrix polynomials. As an example of this, we derive explicit simple upper bounds on the magnitudes of scalar polynomial zeros that are close approximations of the Cauchy radius. Their complexity is of the order of a single polynomial evaluation.
The Cauchy radius is an upper bound that is optimal among all bounds depending only on the moduli of the coefficients, but it has the disadvantage of requiring the solution of a nonlinear equation, so that the close explicit approximation we obtain is a useful substitute.

Speaker: Aaron Melman (Santa Clara University)

35 The Joint Spectral Radius of Neural Networks

Feed-forward neural networks can be interpreted as dyscrete switched dynamical systems where the switching rule is determined by the training process. The stability of switched systems has been extensively investigated in the linear case, where the key ingredient is given by the joint spectral radius of the family of matrices that determines the dynamics of the system. On the other hand, neural networks used in modern deep learning typically require nonlinear subhomogeneous activation functions for which new theoretical investigation is needed. Thus, we investigate switched systems that alternate maps from a (possibly infinite) class $\mathcal{F}$ of nonlinear functions, following an unknown switching rule. In particular we focus on families of sub-homogeneous nonlinear functions that preserve the ordering induced by a cone. We introduce the notion of nonlinear joint spectral radius (JSR) of a family of such functions. Then we show that the JSR provides information about the stability of the system. After that, we investigate properties of the nonlinear JSR, including a dual formulation in terms of monotone prenorms of the functions in $\mathcal{F}$, and a formulation in terms of the spectral radii of the functions in the semigroup generated by $\mathcal{F}$. Finally, making use of the different formulations of the JSR, we present an algorithm devoted to the computation of the nonlinear JSR. Our algorithm, inspired by the polytopal algorithm used in the linear case, iteratively builds a monotone extremal prenorm for the system in terms of the Minkowski functional of a finitely generated subset of the cone.

Speaker: Piero Deidda (Scuola Normale Superiore, Pisa - Gran Sasso Science Institute, L'Aquila)

20:00 Social Dinner

Tuesday, 21 January

Morning Session III

36 Walk-induced strategies for updating Katz centrality in simple graphs

Katz centrality is one of the most widely used indices in network science to determine the relative influence of a node within a network. Its unique feature is that it can be expressed as the solution of the (typically sparse) linear system $(I- \alpha A) \mathbf{x} = \mathbf{1}$, where $A$ is the adjacency matrix of the underlying graph and $\alpha>0$ is a dumping parameter such that $\alpha \rho(A)<1$, making it computable using broadly available and efficient linear algebra methods. However, repeated calculation of the measure can still be challenging in practical applications, in which simulating the sequential removal of a set of nodes/edges is often used to study issues related to network robustness and communication bottlenecks. For large networks, it may be unfeasible to recompute the vector of centrality scores after each elimination, so it is rather useful to devise updating strategies instead.

In this talk, we will provide explicit formulas that characterise the loss of walks that the network faces and, consequently, how Katz centrality is affected when a set of edges (resp., nodes) is removed from the network. Then, we will use these formulas to derive efficient and effective strategies for approximately updating the vector of Katz centrality, significantly reducing the computational burden of recomputing it repeatedly from scratch. We will also show bounds on the modifications of the induced total network communicability.

1. F. Arrigo, D. Bertaccini, A. Filippo, Updating Katz centrality by counting walks. Submitted.
2. M. Benzi, C. Klymko. Total communicability as a centrality measure. Journal of Complex Networks (2013), no. 1, 124–149.
3. L. Katz, A new status index derived from sociometric analysis, Psychometrika, 18 (1953), 765 pp. 39–43.

Speaker: Alessandro Filippo (University of Rome Tor Vergata)

37 Explicit Characterization of Consensus in a Distributed Estimation Problem on Chain Graphs

This Abstract briefly describes the original results that are to be submitted to the international journal IEEE Transactions on Automatic Control. In particular, the problem of looking for time-varying vectors ${\hat{\Theta}}(t)$ - named $\Theta$-estimates - which exponentially converge to the unknown constant parameter vector $\Theta \in \mathbb{R}^{m}$ defined by the set of linear time-varying equations:
\begin{eqnarray}
y_{i}(t) = \phi_{i}^{\rm T}(t) \Theta, \ \ \ i=1,\ldots,p \ \ \ \ \ \ \ \ \ (1)
\end{eqnarray}
is addressed, where $y_{i}$ are the locally measured outputs and $\phi_{i}(\cdot): \mathbb{R}_{0}^{+} \to \mathbb{R}^{m}$ are the local regressor vectors, $i=1,\ldots,p$, each of them assumed to be available at the running time at each node of the graph. This general problem is referred to as the estimation problem for (1) on a graph. Now, when there are no edges connecting nodes of the graph that are far from each other, the problem passes from being full-graph-knowledge based to being partial-graph-knowledge based. This way, the burden of information that has to be communicated to the various measurement/estimation nodes of the graph might be (even largely) reduced. Specifically, the original contribution consists in showing that, under the weakest $\Theta$-identification condition on the graph [CHEN]:

(Cooperative PE Condition) The corresponding regressor matrix $\Phi^{\rm T}(\cdot)\in \mathbb{R}^{p \times m}$
\begin{eqnarray}
\Phi^{\rm T}(\cdot) &=& \left[ \begin{array}{cccc}
\phi_{1}^{\rm T}(\cdot) \
\ldots \
\phi_{p}^{\rm T}(\cdot)
\end{array} \right]
\end{eqnarray}
is persistently exciting (PE), i.e., there exist (known) positive reals $c_{p}$ and $T_{p}$ such that the following condition [$\mathbb{I} \in \mathbb{R}^{m \times m}$] holds:
\begin{eqnarray}
\int_{t}^{t+T_{p}}\Phi(\tau)\Phi^{\rm T}(\tau) {\rm d}\tau \geqq c_{p}\mathbb{I}, \ \ \ \forall \ t \geq 0,
\end{eqnarray}

a set of suitably tailored differential equations for the time-dependent vectors ${\hat{\Theta}}^{[i]}(t)$, all of them converging to the unknown $\Theta$ and by thus achieving consensus, can be (redundantly) designed at each node $i=1,\ldots,p$. This is proved under the condition that nodes undirectedly connected in series (undirected chain graph [CONCAS]) are considered, so that each estimation scheme at the node can share information - namely, its own $\Theta$-estimate - with the neighbours only, one for node $1$ and $p$, two for the remaining nodes. Thus the aforementioned problem belongs to the class of problems known as {\em multi-agent estimation with consensus}.
The general scenario is the one given by a multisensor network in which a distributed parameter estimator is to be designed on the basis of distributed sensing. This situation typically arises in the case in which mobile sensors measure the distribution of an unknown quantity over a field and, due to the size of the region, it is computationally inefficient or just infeasible to visit every point in the space to collect data [MATVEEV]. Now, agents at the nodes face a local identification problem, in which they cannot consistently estimate the parameter vector in isolation, so they have to engage in communication with their neighbours. In particular, estimate-consensus has to be achieved through a sort of penalization of the mismatch between the parameter estimates.
Indeed, in this Abstract, by assessing the positive definite nature of a quadratic form associated with a tridiagonal block structure, we prove the exponential achievement - with an explicit characterization of the exponential convergence - of the distributed parameter estimation task on an undirected chain graph, in which an original neighbourhood-based decentralized parallel architecture (aiming at reducing the computational burden of a centralized estimation scheme on the graph) is adopted. It is not necessary for $\Theta$ to be fully identifiable at each node. Instead, it is sufficient for the nodes in the graph to cooperatively provide information at the local stage. Indeed, the derivations here are essentially in accordance with [CHEN], however, in contrast to [CHEN], which - for general graph topologies - uses weaker contradiction arguments to prove exponential consensus, here original proofs of convergence are able to provide an explicit characterization of the exponentially achieved consensus.
Nevertheless, the problem of identifying time-varying parameters that are periodic with known periods can be innovatively solved as well. This constitutes an additional original contribution of the paper: adaptive tools can be directly replaced by repetitive learning ones [VERRELLI] within the same theoretical framework, where the asymptotic consensus is successfully ensured under identification mechanisms based on the information exchange between neighbours.

• W. Chen, C. Wen, S. Hua, and C. Sun, Distributed cooperative adaptive identification and control for a group of continuous‑time systems with a cooperative PE condition via consensus, IEEE Transactions on Automatic Control, vol. 59(1), pp. 91‑106, 2014.

• A. Concas, L. Reichel, G. Rodriguez, and Y. Zhang, Chained graphs and some applications, Applied Network Science, vol. 6(39), 2021.

• A.S. Matveev, M. Almodarresi, R. Ortega, A. Pyrkin, and S. Xie, Diffusion-based distributed parameter estimation through directed graphs with switching topology: Application of dynamic regressor extension and mixing, IEEE Transactions on Automatic Control, vol. 67(8), pp. 4256-4263, 2022.

• C.M. Verrelli, Exponential output tracking via a saturated RLC for a class of nonlinear systems under PSD of periodic references, Systems and Control Letters, vol. 192, 2024.

Speaker: Prof. Cristiano Maria Verrelli (University of Rome Tor Vergata, Electronic Engineering Department)

38 Modelling advection on distance-weighted directed networks

We introduce a mathematical model to capture and describe the dynamics of the advection differential equation on distance-weighted directed graphs, with applications to various networked systems. The primary objective of our model is to generalize advection processes—which traditionally describe phenomena like fluid flow or traffic movement—by formulating them within the framework of discrete network structures.

Our approach begins with defining a set of essential properties, or axioms, that any discrete advection operator must satisfy. These axioms ensure that the operator reflects the fundamental characteristics of advection processes in continuous spaces, such as directional flow, conservation properties, and respect for edge weights based on distance. We demonstrate that, under these conditions, there exists an essentially unique operator that fulfills all prescribed axioms, providing a robust and mathematically sound foundation for discrete advection on directed graphs.

To validate our theoretical results, we present both analytical and numerical examples that showcase the operator's behavior under various configurations. As a practical demonstration, we apply the model to simulate traffic flow in a transportation network.

Speaker: Francesco Zigliotto (Scuola Normale di Pisa)

39 Generalized approximating class of sequences and asymptotic block structures with rectangular Toeplitz blocks.

In the talk, the problem of identifying the distribution in the Weyl sense (singular values and eigenvalues) is treated, when block matrix-sequences with rectangular Toeplitz blocks are considered. The case of sequences whose block sizes tend to rational numbers is of more classical flavour and is concisely introduced.
Using recent results regarding the notion of generalized approximating class of sequences (g.a.c.s), a symbol for such sequences is computed when the proportions between the different block sizes tend to irrational numbers. In doing so, the flexibility of the idea of g.a.c.s. is emphasized, when dealing with integro-differential equations on unbounded or moving domains.
Theoretical results also suggest practical preconditioning strategies for designing fast preconditioned Krylov solvers for the related large linear systems.

References:

1. Adriani A.,Schiavoni-Piazza A.J.A.,Serra-Capizzano S.,"Blocking structures, g.a.c.s. approximation, and distributions." Special Volume in Memory of Prof.Nikolai Vasilevski (2024), to appear.
2. Adriani A.,Schiavoni-Piazza A.J.A.,Serra-Capizzano S.,"Revisiting the notion of approximating class of sequences for handling approximated PDEs on moving or unbounded domains." Electronic Transactions on Numerical Analysis (2024), to appear.
3. Barakitis N., Ferrari P., Furci I., Serra-Capizzano S., "An extradimensional approach for distributional results: the case of 2x2 block Toeplitz Structures." Springer Proceedings on Mathematics and Statistics (2024), to appear.
4. Furci I., Adriani A., Serra-Capizzano S., "Block structured matrix-sequences and their spectral and singular value distributions: a general theory." In arXiv:2409.06465 (2024)

Speaker: Alec Jacopo Almo Schiavoni Piazza (SISSA)

40 SVD-based Line Integral Methods for preserving multiple invariants of Hamiltonian problems

In this talk, we generalize the class of energy-conserving Runge-Kutta methods, named Hamiltonian Boundary Value Methods, to handle the numerical solution of Hamiltonian problems with additional independent invariants besides the Hamiltonian. The proposed strategy relies on the solution of a perturbed problem, where a minimum norm perturbation is computed by resorting to the singular value decomposition. The analysis of the approach is given and numerical tests are reported, to make evidence of the theoretical findings and assess their effectiveness.

1. L. Brugnano, F. Iavernaro. Line Integral Methods for Conservative Problems. Chapman et Hall/CRC, Boca~Raton, FL, USA, (2016).
2. L. Brugnano, F. Iavernaro. Modified line integral methods for conservative problems with multiple invariants. AIP Conference Proceedings 1648 150010 (2015).
3. L. Brugnano, F. Iavernaro. Line Integral Methods which preserve all invariants of conservative problems. J. Comput. Appl. Math. 236 3905–3919 (2012).
4. L. Brugnano, Y. Sun. Multiple invariants conserving Runge-Kutta type methods for Hamiltonian problems. Numer. Algorithms 65 611–632 (2014).

Speaker: Dr Gianmarco Gurioli (University of Florence)

10:40 Coffee Break

Morning Session IV

41 Randomized Gauss-Newton methods for large scale nonlinear least-squares

We address the solution of large-scale nonlinear least-squares problems by stochastic Gauss-Newton methods combined with a line-search strategy. The algorithms proposed have per-iteration computational complexity lower than classical deterministic methods, due to the employment of random models inspired by randomized linear algebra tools. Under suitable assumptions, the stochastic optimization procedures can achieve a desired level of accuracy in the first-order optimality condition. We discuss the construction of the random models and the iteration complexity results to drive the gradient below a prescribed accuracy.

Speaker: Benedetta Morini (Università di Firenze)

42 Row-aware Randomized SVD with applications

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

1. N. Halko, P. G. Martinsson, and J. A. Tropp. Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions. SIAM Review \textbf{53.2}, 217–288 (2011). DOI: 10.1137/090771806
2. Davide Palitta and Sascha Portaro. Row-aware Randomized SVD with applications. arXiv: 2408.04503 [math.NA] (2024). URL: https://arxiv.org/abs/2408.04503.

Speaker: Sascha Portaro (University of Bologna)

43 A Deep-QLP Decomposition Algorithm and Applications

Singular value decomposition (SVD) is a fundamental tool in data analysis and machine learning. Starting from the Stewart’s QLP decomposition [1], we propose an innovative Deep-QLP decomposition algorithm for efficiently computing an approximate Singular Value Decomposition (SVD) based on the preliminary work in [2]. Given a specified tolerance $\tau$, the algorithm automatically computes a positive integer $f$ and a factorization $\mathcal{U}_f \mathcal{L}_f^D \mathcal{V}_f^T$, with $\mathcal{L}_f^{D}$ diagonal matrix, $\mathcal{U}_f,\mathcal{V}_f$ matrices of rank $f$ with orthonormal columns such that $ \|A-\mathcal{U}_{f}\mathcal{L}^{D}_{f}\mathcal{V}_{f}^{T}\|_2 \leq 3\tau \|A\|_2.$
The Deep-QLP algorithm stands out for its ability to return an approximation of the largest singular values, based on a fixed tolerance, to achieve significant dimensionality reduction while simultaneously preserving essential information in the data. In addition, it can also be used to return an approximation of the smallest singular values that can be used in some applications.
The algorithm has been successfully integrated with the randomized SVD [3], making
the Deep-QLP algorithm particularly effective for sparse matrices, which are prevalent
in numerous applications such as text mining.
Several numerical experiments have been conducted, demonstrating the effectiveness of the proposed method.

References

1. Gilbert W. Stewart. The QLP approximation to the singular value decomposition. SIAM J. Sci. Comput., 20:1336–1348, 1999.
2. Antonella Falini and Francesca Mazzia. Approximated iterative QLP for change detection in hyperspectral images. AIP Conference Proceedings, 3094(1):370003, 06 2024.
3. Nathan Halko, Per-Gunnar Martinsson, and Joel A Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM Review, 53(2):217–288, 2011.

Speaker: Cristiano Tamborrino (Università degli Studi di Bari)

44 Extracting accurate singular values from approximate subspaces

We present a comprehensive matrix perturbation analysis of methods for extracting singular value from approximate subspaces. Given (orthonormal) approximations $\tilde{U}$ and $\tilde{V}$ to the left and right subspaces spanned by the leading singular vectors of a matrix A, we discuss methods to approximate the leading singular values of A and study their accuracy. In particular, we focus our analysis on the generalized Nyström approximation, as surprisingly, it is able to obtain significantly better accuracy than classical methods, namely Rayleigh-Ritz and (one-sided) projected SVD. A key idea of the analysis is to view
the methods as finding the exact singular values of a perturbation of A. In this context, we derive a matrix perturbation result that exploits the structure of such $2\times 2$ block matrix perturbation. We then obtain bounds on the accuracy of the extracted singular values. This leads to sharp bounds that predict well the approximation error trends and explain the difference in the behaviour of these methods. Our theoretical results translate into practical insights on how best to approximate singular values given approximate subspaces. Specifically, when approximations to both left and right subspaces are available, the generalized Nyström approach consistently outperforms other single-pass methods by achieving higher accuracy in the leading singular values, as indicated by our bound analysis.

To validate our theoretical predictions, we conducted numerical experiments. These tests highlight how the derived bounds exhibit tighter estimates than classical bounds, effectively capturing the trends observed empirically. For instance, the bounds reflect a more gradual change in error across the leading singular values.

1. Matrix perturbation analysis of methods for extracting singular values from approximate singular subspaces, L.Lazzarino, H. Al Daas, Y. Nakatsukasa, 2024, arXiv:2409.09187

Speaker: Lorenzo Lazzarino (University of Oxford)

45 The Perron-Frobenius theory for matrix semigroups: a geometrical approach

The Perron-Frobenius theory was extended to multiplicative matrix semigroups relatively recently and has found many applications to synchronizing automata, nonhomogeneous Markov chains, linear dynamical systems, graphs, etc. This theory connects the spectral and combinatiorial properties of semigroups of nonnegative matrices. In particular, the concept of primitivity is especially important. A set of nonnegative matrices is called primitive if it possesses at least one strictly positive product. There are rather surprising links between primitive semigroups and the geometrical problem of self-similar partitions of a convex body. They can produce short and clear proofs of some known results on primitivity, for example, the characterization of primitive families, the existence of a common invariant affine subspace for matrices of non-synchronizing automata, etc. A new approach is also useful to study Hurwitz primitive (or k-primitive) semigroups. We derive possible generalizations of those results to sets of matrices with a common invariant cone.

Speaker: Vladimir Protasov (University of L’Aquila)

12:50 Lunch Break

Afternoon Session II

46 Provably convergent, extrapolated algorithms for Nonlinear Matrix Decomposition problems

In this contribution I propose a new problem in low-rank matrix factorization, that is the Nonlinear Matrix Decomposition (NMD): given a sparse nonnegative matrix~$X$, find a low-rank matrix $\Theta$ such that $X \approx f(\Theta)$, where $f$ is element-wise nonlinear. I will focus on the so-called ReLU-NMD, where $f(\cdot) = \max(0, \cdot)$, the rectified unit (ReLU) non-linear activation.

At first, I will provide a brief overview of the motivations, possible interpretations of the model and its connection with neural networks, explaining the idea that stands behind ReLU-NMD and how nonlinearity can be exploited to get low-rank approximation of given data. Secondly, I will talk about standard approaches to model ReLU-NMD as an optimization problem.

Then, I will introduce the Block-Coordinate-Descent (BCD) method along with some convergence guarantees. Moreover, I will show how the BCD scheme applies to ReLU-NMD (BCD-NMD) and how it can be accelerated via extrapolation, maintaining its convergence properities also in the non-convex setting (eBCD-NMD).

Finally, I will illustrate the effectiveness of the proposed algorithms on synthetic and real-world data sets, providing some possible applications.

Speaker: Giovanni Seraghiti (Unifi/Umons)

47 DeepONet-based approximation of inverse operators: A data-driven approach to address the low-frequency part of the eigenspectrum

Preconditioners are essential tools for efficiently solving linear systems arising from the discretization of PDEs. Traditional single-level approaches like Jacobi, Incomplete LU factorization (ILU), and Factorized Sparse Approximate Inverse (FSAI)$^{1}$ are effective in reducing high frequency error components but struggle with low-frequency components. The basic idea of multigrid and two-level domain decomposition is to achieve efficiency by balancing coarse-grid correction with smoothing.

In this work, we propose a novel method for approximating the inverse of discrete operators by leveraging DeepONet, a supervised learning framework for nonlinear operators$^2$. This Deep Learning based method is well-suited to focus on low-frequency features$^3$. Therefore, instead of employing a conventional coarse grid in the multigrid sense, we construct an analogous structure through a DeepONet trained on vectors representing low frequencies.

Alternated to classical single-level preconditioners, the DeepONet approximation compensates the action on the low-frequency part of the eigenspectrum, thus accelerating the convergence. Preliminary test cases are investigated and presented in order to analyze the potential of the proposed approach, in view of its possible application in matrix-free multi-level methods.

1. C. Janna, M. Ferronato, F. Sartoretto, G. Gambolati. FSAIPACK: A Software Package for High-Performance Factored Sparse Approximate Inverse Preconditioning. ACM Trans. Math. Softw. 41.2, 1-26 (2015).
2. L. Lu, P. Jin, G. Pang, Z. Zhang, G. E. Karniadakis. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature machine intelligence 3.3, 218-229 (2021).
3. E. Zhang, A. Kahana, A. Kopaničáková, E. Turkel et al. Blending neural operators and relaxation methods in PDE numerical solvers. Nature machine intelligence 6, 1303–1313 (2024).

Speaker: Caterina Millevoi (University of Padova)

48 On the ideality of real-world datasets for photometric stereo under unknown lights position

Photometric stereo is a computer vision technique for reconstructing the shape of a three-dimensional object starting from digital images. Several assumptions are required but they are rarely verified in experimental datasets. Specifically, the object under observation should behave as a Lambertian reflector, with light sources positioned at an infinite distance, along known directions. In real-world scenarios these assumptions can not be satisfied, indeed is impossible to accurately measure the relative position between the light sources and the target. This situation is common in archaeological applications. Although the Hayakawa method determines an estimation of the light source position starting from the data images, in some cases it breaks down when some images deviate from ideality. To identify which images from a given dataset should be selected to produce a better reconstruction, we introduce and discuss two measures of data ideality. Then, we investigate the effectiveness of these indicators using synthetic and experimental datasets.

1. E. Crabu, F. Pes, G. Rodriguez, and G. Tanda. Ascertaining the Ideality of Photometric Stereo Datasets under Unknown Lighting. Algorithms, 16(8), 375 (2023).
2. E. Crabu, F. Pes, G. Rodriguez, and G. Tanda. Color reconstruction by photometric stereo with unknown lighting of rock art carvings found in two Sardinian Domus de Janas. Submitted (2024).

Speaker: Federica Pes (Università di Cagliari)

49 Fractional Laplacian and ADMM for glyph extraction

In archaeology it is a common task to extract incisions or glyphs from a surface. This procedure is usually done manually and, therefore, it is prone to errors and it can be extremely time consuming. In this talk we present a variational model to automatically extract these incisions from a smooth surface.

We model this problem in the following way. Let $\mathbf{x}\in\mathbb{R}^n$ be a vector containing a sampling of the archaeological surface, we wish to find two vectors $\mathbf{x}_s^*$ and $\mathbf{x}_g^*$ such that $\mathbf{x}=\mathbf{x}_s^*+\mathbf{x}_g^*$, where $\mathbf{x}_s^*$ is smooth and contains the background and $\mathbf{x}_g^*$ is sparse and contains the glyph. To this aim we consider the model
$$ \begin{split} \left(\mathbf{x}_s^*,\mathbf{x}_g^*\right)=&\arg\min_{\mathbf{x}_s,\mathbf{x}_g\in\mathbb{R}^{n\times n}} \frac{1}{2}\left\|L^\alpha\textbf{x}_s\right\|_2^2+\mu\left\|\mathbf{x}_g\right\|_1,\\&\quad \mbox{s.t. }\mathbf{x}_s+\mathbf{x}_g=\mathbf{x}, \end{split} $$ where $\mu>0$, $\alpha\in[1,2]$, $$\|\mathbf{x}\|_p^p=\sum_{i=1}^{n}|\mathbf{x}_i|^p,$$ and $L\in\mathbb{R}^{n\times n}$ denotes the Laplacian operator. To perform the minimization,
we employ the Alternating Direction Multiplier Method (ADMM). We provide a procedure to generate
realistic synthetic data and we show the performances of the proposed method on this kind of data.

1. S. Boyd, N. Parikh, E.Chu, B. Peleato, J. Eckstein, et al, {Distributed optimization and statistical learning via the alternating direction method of multipliers}, Foundations and Trends® in Machine learning 3 (2011) 1–122.
2. E. Di Nezza, G. Palatucci, E. Valdinoci, Hitchhiker’s guide to the fractional Sobolev spaces, Bulletin des sciences mathématiques 136 (2012) 521–573.
3. A. Gholami, S. Gazzola, Automatic balancing parameter selection for Tikhonov-TV regularization, BIT Numerical Mathematics 62 (2022) 1873–1898.

Speaker: Mr Andrea Azzarelli (University of Cagliari)

50 RandONets: Shallow networks with random projections for learning nonlinear operators

Neural Operators such as DeepONets have been recently introduced to approximate nonlinear operators with a focus on solution operators of PDEs. However, their implementation requires the use of deep neural networks whose training is performed in a high-dimensional space of parameters and hyperparameters. This, coupled with the need for significant computational resources, creates challenges for achieving high numerical accuracy at limited resources. Very recently, inspired by DeepONets, we have introduced RandONets [1] shallow networks embedding the input space with random projections, thus using specialized numerical analysis techniques from linear algebra for their training to efficiently and accurately learn linear and nonlinear operators. We prove that RandONets are universal approximators of nonlinear operators. Furthermore, we assess their performance with a focus on the approximation of evolution operators (right-hand-sides (RHS)) of PDEs. We demonstrate that, RandONets outperform DeepONets by several orders of magnitude both in terms of numerical approximation accuracy and computational cost. In summary, our work shows that carefully designed “light” neural networks, aided by tailor-made numerical analysis methods, can provide significantly faster and more accurate approximations of nonlinear operators compared to deep neural networks.

[1] Fabiani, G., Kevrekidis, I. G., Siettos, C., & Yannacopoulos, A. N. (2025). RandONets: Shallow networks with random projections for learning linear and nonlinear operators. Journal of Computational Physics, 520, 113433.

Speaker: Prof. Constantinos Siettos (Università degli Studi di Napoli Federico II)

51 Efficient Regularization and Numerical Reconstruction Methods for Inverse Source Problems in Generalized Diffusion Equation

We address two inverse source problems when determining a space-dependent source term and a time-dependent coefficient for a two-dimensional generalized diffusion equation. These problems are ill-posed in the Hadamard sense, where small perturbations in the data can lead to uncontrolled variations in the solution. From a analytic viewpoint we provide existence and uniqueness results for the solutions of these problems under appropriate over-specified and regularity conditions. From a computational viewpoint, we use a quasi-boundary value regularization to reconstruct the underlying approximate source term. In the reconstruction, by employing a finite difference discretization of the regularized problem, we encounter a two-by-two (structured) block linear system. We propose efficient (preconditioned) Krylov solvers and we present numerical experiments to demonstrate the effectiveness of our method.

1. A. Ilyas, S. A. Malik, Direct and some inverse problems for a generalized diffusion equation with variable coefficients, Computational and Applied Mathematics, 43(6), (2024), 364 pp.
2. H.-K. Pang, H.-H. Qin, S. Ni, Sine transform based preconditioning for an inverse source problem of time-space fractional diffusion equations. Journal of Scientific Computing, 100(3), Paper No. 74,(2024), 32 pp.
3. C. Garoni, S. Serra-Capizzano, Generalized locally Toeplitz sequenecs: theory and applications, Vol. II, Springer, Cham, 2018.
4. A. Ilyas, S. Serra-Capizzano, Inverse Source Problems for Identifying Time and Space-Dependent Coefficients in a 2D Generalized Diffusion Equation, preprint 2024.

Speaker: Asim Ilyas (University of Insubria)

16:20 Closing Remarks