Seminari MAT/08 Team

This seminar addresses the problem of image segmentation through an accurate high-order scheme based on the level-set method. In this approach, the curve evolution is described as the 0-level set of a representation function, but the velocity that drives the curve to the boundary of the object has been modified in order to obtain a new velocity with additional properties that are extremely useful to develop a more stable high-order approximation with a small additional cost. The approximation scheme proposed here is the 2D version of an adaptive “filtered” scheme, which combines two building blocks (a monotone scheme and a high-order scheme) via a filter function and smoothness indicators that allow one to detect the regularity of the approximate solution adapting the scheme in an automatic way. Some numerical tests on synthetic and real images confirm the accuracy of the proposed method and the advantages given by the new velocity.

Multilayer networks are a type of complex network that consist of multiple layers, where each layer represents a different type of connection or interaction between the same set of nodes. These networks are used to model systems where entities are connected in multiple ways simultaneously, capturing the complexity of real-world relationships better than traditional single-layer networks. Through a particular interlayer structure, the dynamical evolution of a complex system over time can be represented. Computing the centrality of a temporal network can improve our understanding of how the most important nodes in a network change over time. Our focus is centered on the computation of the centralities of a multilayer temporal network whose modifications over time consist of low-rank updates of the edge adjacency matrix of a transport network. Using Krylov subspace methods for matrix function approximations, we will exploit the particular structure of the problem to gain some computational advantages and modeling insights.

We introduce a novel procedure for computing an SVD-type approximation of a tall matrix A. Specifically, we propose a randomization-based algorithm that improves the standard Randomized Singular Value Decomposition (RSVD). Most significantly, our approach, the Row-aware RSVD (R-RSVD), explicitly constructs information from the row space of A. This leads to better approximations to Range(A) while maintaining the same computational cost. The efficacy of the R-RSVD is supported by both robust theoretical results and extensive numerical experiments. Furthermore, we present an alternative algorithm inspired by the R-RSVD, capable of achieving comparable accuracy despite utilizing only a subsample of the rows of A, resulting in a significantly reduced computational cost. This method, that we name the Subsample Row-aware RSVD (Rsub-RSVD), is supported by a weaker error bound compared to the ones we derived for the R-RSVD, but still meaningful as it ensures that the error remains under control. Additionally, numerous experiments demonstrate that the Rsub-RSVD trend is akin to the one attained by the R-RSVD when the subsampling parameter is on the order of n, for a m×n A, with m >> n. Finally, we consider the application of our schemes in two very diverse settings which share the need for the computation of singular vectors as an intermediate step: the computation of CUR decompositions by the discrete empirical interpolation method (DEIM) and the construction of reduced-order models in the Loewner framework, a data-driven technique for model reduction of dynamical systems.

An inverse problem is the task of retrieving an unknown quantity from indirect observations. When the model describing the measurement acquisition is linear, this results in the inversion of a linear operator (a matrix, in a discrete formulation) which is usually ill-posed or ill-conditioned. A common strategy to tackle ill-posedness in inverse problems is to use regularizers, which are (families of) operators providing a stable approximation of the inverse map. Model-based regularization techniques often leverage prior knowledge of the exact solution, such as smoothness or sparsity with respect to a suitable representation; on the other side, in recent years many data-driven methods have been developed in the context of machine learning. Those techniques tackle the approximation of the inverse operator in suitable spaces of parametric functions (i.e., neural networks) and rely on large datasets of paired measurements and ground-truth objects. In this talk, I will focus on hybrid strategies, which aim at blending model-based and data-driven approaches, providing both satisfying numerical results and sound theoretical guarantees. I will describe a general framework to comprise many existing techniques in the theory of statistical learning, also reporting some recent theoretical advances (in the direction of generalization guarantees). I will help the discussion by presenting some relevant examples in the context of medical imaging and, specifically, in computed tomography.

We are interested in the numerical solution of the matrix least squares problem min_X ∥AXB + CXD-F ∥_F , with A and C full column rank, B, D full row rank, F an n×n matrix of low rank, and ∥•∥_F the Frobenius norm. We derive a matrix-oriented implementation of LSQR, and devise an implementation of the truncation step that exploits the properties of the method. Experimental comparisons with the Conjugate Gradient method applied to the normal matrix equation and with a (new) sketched implementation of matrix LSQR illustrate the competitiveness of the proposed algorithm. We also explore the applicability of our method in the context of Kronecker-based Dictionary Learning, and devise a representation of the data that seems to be promising for classification purposes.

Tensors have become widely used in various domains due to their practicality. Tensor factorization techniques are used to solve computationally demanding problems, analyze large datasets, and refine descriptions of complex phenomena. This presentation outlines the development of my research on tensors, including an overview of commonly used tensor methods and their applications in various fields such as remote sensing, multilinear algebra, numerical simulation, and signal processing. Criteria for selecting the most appropriate tensor technique depending on the problem under consideration will be emphasized. The presentation aims to outline the advantages and limitations inherent in these techniques. It explores the challenges and offers insights into current research directions driven by real-world, computational, and applied problems.

Sunlight constitutes an abundant and endless natural fuel, available worldwide. In a society where a substantial part of the global energy yield is being directly expended at the city scale, urban areas appear as serious candidates for the production of solar energy. Their intrinsic complexity yet makes it challenging. The morphological heterogeneity between urban geometries and intricacy of their materials optical properties especially contribute together to causing important spatiotemporal variations in the distribution of incident solar radiations. The field of irradiance received by a specific urban region (e.g. façade, building, district) may thus rapidely become the result of complex miscellaneous interactions between many degrees of freedom. Besides, Principal Component Analysis (PCA) has been widely validated as an efficient algorithm to identify the principal behavioural features, or modes of variability, of a high-dimensional phenomenon. An approach is proposed here for analysing the variations in space and time of the solar resource within an urban context by means of PCA. A parametric investigation is conducted on a set of theoretical 100×100 m² urban districts, defined as arrangements of cuboid-like buildings, with various typological indicators (Total Site Coverage, Average Building Height) and surface materials (Lambertian, highly-specular) at three different latitudes. For each configuration, the distribution of irradiance incident on the facets of the central building is modelled via backwards Monte-Carlo ray tracing over a full year and under clear sky conditions, with a 15 min timestep and 1 m spatial resolution. PCA is subsequently applied to the simulated radiative fields to extract dominant modes of variation. First results validate energy-based orthogonal decompositions like PCA as efficient tools for characterising the variability distribution of multivariate phenomena in this context, allowing for the identification of district areas subjected to important spatial and temporal variations of the solar resource. Characteristic time scales are clearly represented across successive orders of decomposition. Information about the district morphology is also obtained, with the contribution of surrounding geometries being portrayed by specific spatial modes. Similar prevalent variables are further repetitively encountered across multiple evaluated surfaces, but at different modal ranks.

CholeskyQR is a popular algorithm for QR factorization in both academia and industry. In order to have good orthogonality, CholeskyQR2 is developed by repeating CholeskyQR twice. Shifted CholeskyQR3 introduces a shifted item in order to deal with ill-conditioned matrices with good orthogonality. This talk primarily focuses on deterministric methods. We define a new matrix norm and make improvements to the shifted item and error estimations in CholeskyQR algorithms. We use such a technique and provide an analysis to some sparse matrices in the industry for CholeskyQR. Moreover, we combine CholeskyQR and our new matrix norm with randomized models for probabilistic error analysis and make amelioration to CholeskyQR. A new 3-step algorithm without CholeskyQR2 is also developed with good orthogonality.

We introduce the definition of tensorized block rational Krylov subspaces and their relation with multivariate rational functions, extending the formulation of tensorized Krylov subspaces introduced in [2]. Moreover, we develop methods for the solution of tensor Sylvester equations with low multilinear or Tensor Train rank, based on projection onto a tensor block rational Krylov subspace. We provide a convergence analysis and some strategies for poles selection based on the techniques developed in [1]. References [1] A. A. Casulli and L. Robol. “An effcient block rational Krylov solver for Sylvester equations with adaptive pole selection”. In: arXiv preprint arXiv:2301.08103 (2023). [2] D. Kressner and C. Tobler. “Krylov subspace methods for linear systems with tensor product structure”. In: SIAM Journal on Matrix Analysis and Applications 31.4 (2010), pp. 1688–1714

The widespread interest in quantum computation has motivated (among others) applications to network analysis, where quantum advantage may turn out to be especially beneficial for the treatment of large-scale problems. In the past years, several authors have proposed the use of quantum walks - as opposed to classical random walks - in the definition and analysis of centrality measures for graphs, which in turn are the basis for ranking algorithms. Here we focus on unitary continuous-time quantum walks (CTQW) applied to directed graphs and propose new quantum algorithms for hub and authority ranking of the nodes. In particular we explore - the choice of Hamiltonian operator that defines the time evolution of a CTQW, - the choice of the initial state of the system (which turns out to have a non-negligible effect on the final ranking), and perform numerical comparisons with well-known classical ranking algorithms such as HITS and PageRank.

As part of the Una Europa funded grant TOC4Deep (Tensor-based Optimal Control Approaches for Deep Learning) a series of half-day workshops will be held over the next 6 months to encourage scientific networking and discussions between the project universities of Edinburgh, Bologna and KU Leuven. The second of these workshops will take place on Wednesday 06th July, with a main focus on computer scientists' perspective on machine and deep learning. The workshop will be hybrid, with in person attendance in Seminario II, Dipartimento di Matematica, University of Bologna, or via Zoom.

As part of the Una Europa funded grant TOC4Deep (Tensor-based Optimal Control Approaches for Deep Learning) a series of half-day workshops will be held over the next 6 months to encourage scientific networking and discussions between the project universities of Edinburgh, Bologna and KU Leuven. The first of these workshops will take place on Thursday 21st April, with a focus on research relating to optimization and optimal control. The workshop will be hybrid, with in person attendance in JCMB 5323 or via Zoom (link to follow). The timetable and speaker information is below: Time (BST) - please mind the time zone 9-9.30 - Introduction/TOC4Deep presentation 9.30-10.15 - John Pearson, University of Edinburgh "Preconditioned Iterative Methods for Multiple Saddle-point Systems Arising from PDE-constrained Optimization" 10.15-10.30 - Break 10.30-11.15 - Wim Michiels, KU Leuven "Stability, Robustness Analysis and Model Order Reduction of Periodic Control Systems with Delay" 11.15-12 - Margherita Porcelli, Università di Bologna "A spectral PALM algorithm for Dictionary Learning"

In this seminar, I overview the research work carried out by IMATI-CNR on the extension of the Hough transform (HT) to recognize families of planar and spatial curves and surface primitives on 3D objects and point clouds. Recent developments are transforming the HT into a tool computationally affordable even outside the classical context of recognition of lines, circles and ellipses in the plane and planes and spheres in the space. In particular, I will focus on some applications such as the characterization of curves and complex patterns on artefacts and the creation of geometric models with curvilinear elements.

We consider convex optimisation problems defined in the variable exponent Lebesgue space L^p(·)(Ω), where the functional to minimise is the sum of a smooth and a non-smooth term. Compared to the standard Hilbert setting traditionally considered in the framework of continuous optimisation, the space L^p(·) (Ω) has only a Banach structure which does not allow for an identification with its dual space, as the Riesz representation theorem does not hold in this setting. This affects the applicability of well-known proximal (a.k.a. forward-backward) algorithms, since the gradient of the smooth component here lives in a different space than the one of the iterates. To circumvent this issue, the use of duality mappings is required; they link primal and dual spaces in a nonlinear fashion, thus allowing a sensible definition of the algorithmic iterates. However, such nonlinearity introduces further difficulties in the definition of the proximal (backward) step and, overall, in the convergence analysis of the algorithm. To overcome the non-separability of the natural induced norm on L^p(·)(Ω), we consider modular functions allowing for a an appropriate definition of proximal algorithms in this setting for which convergence properties in function values can be proved. Some numerical examples showing the flexibility of our approach in comparison with standard (Hilbert, L^p with constant p) algorithms on some exemplar inverse problems (deconvolution, denoising) are showed.

We put established Krylov subspace methods and Gauss quadrature rules to new use by generalizing the class of matrix function-based centrality measures from single-layered to multiplex networks. Our approach relies on the supra-adjacency matrix as the network representation, which has already been used to generalize eigenvector centrality to temporal and multiplex networks. We discuss the cases of unweighted and weighted as well as undirected and directed multiplex networks and present numerical studies on the convergence of the respective methods, which typically requires only few Krylov subspace iterations. The focus of the numerical experiments is put on urban public transport networks.

Optimization problems subject to PDE constraints form a mathematical tool that can be applied to a wide range of scientific processes, including fluid flow control, medical imaging, biological and chemical processes, and many others. These problems involve minimizing some function arising from a physical objective, while obeying a system of PDEs which describe the process. It is necessary to obtain accurate solutions to such problems within a reasonable CPU time, in particular for time-dependent problems, for which the “all-at-once” solution can lead to extremely large linear systems. In this talk we consider Krylov subspace methods to solve such systems, accelerated by fast and robust preconditioning strategies. A key consideration is which time-stepping scheme to apply — much work to date has focused on the backward Euler scheme, as this method is stable and the resulting systems are amenable to existing preconditioners, however this leads to linear systems of even larger dimension than those obtained when using other (higher-order) methods. We will summarise some recent advances in addressing this challenge, including a new preconditioner for the more difficult linear systems obtained from a Crank-Nicolson discretization, and a Newton-Krylov method for nonlinear PDE-constrained optimization. At the end of the talk we plan to discuss some recent developments in the preconditioning of multiple saddle-point systems, specifically positive definite preconditioners which may be applied within MINRES, which may find considerable utility for solving optimization problems as well as other applications. This talk is based on work with Stefan Güttel (University of Manchester), Santolo Leveque (University of Edinburgh), and Andreas Potschka (TU Clausthal).

The tensor rank decomposition or canonical polyadic decomposition (CPD) is a generalization of a low-rank matrix factorization from matrices to higher-order tensors. In many applications, multi-dimensional data can be meaningfully approximated by a low-rank CPD. In this talk, I will describe a Riemannian optimization method for approximating a tensor by a low-rank CPD. This is a type of optimization method in which the domain is a smooth manifold, i.e. a curved geometric object. The presented method achieved up to two orders of magnitude improvements in execution time for challenging small-scale dense tensors when compared to state-of-the-art nonlinear least squares solvers.