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...
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...
Multiple orthogonal polynomials (MOPs) arise in various applications, including approximation theory, random matrix theory, and numerical integration. To define MOPs, one needs multiple inner products. In this talk, we restrict our attention to the case of two inner products. These MOPs satisfy recurrence relations, and we focus specifically on the stepline recurrence relation.
We derive an...
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...
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...
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...
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...
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...
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$...
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...
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...
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...
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$...
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...
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...
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...
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,...
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,...
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$...
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...
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...
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...
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...
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 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...
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...
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...
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...
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$....
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...
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,...
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...
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...
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...
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...
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...
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...
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...
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...
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...
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....
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...
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...
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...
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...
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...
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...
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...
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...
