Projekt: NA-KEX VT2025
| Handledare | Projektets titel | Beskrivning av projektet |
| Krantz/Tornberg | Error estimation for numerical integration in boundary integral methods | This project aims to enhance the efficiency of evaluating so-called layer potentials in the MATLAB package chunkIE (https://github.com/fastalgorithms/chunkie), which is used to solve partial differential equations (PDEs) via boundary integral methods. Layer potentials are integrals over the boundary of a domain that represent solutions to PDEs when using boundary integral methods. For points far from the boundary, standard numerical integration (quadrature) methods, such as Gaussian or trapezoidal rules, suffice for accurate approximation of the layer potentials. However, for points near the boundary, where the integrand exhibits rapid variation, specialized quadrature methods are necessary to maintain accuracy. The goal of this thesis is to derive and implement error estimates for regular quadrature in chunkIE, enabling it to automatically determine when standard quadrature methods are sufficient based on a user-defined error tolerance. |
| Krantz/Tornberg | Rational vs. polynomial approximations of challenging functions | This project explores the approximation of functions with singularities and poles using both polynomial and rational approximation techniques. While polynomial approximation works well for smooth functions, it often struggles with functions that exhibit singularities or poles, such as f(x) = |x| on [-1,1]. In contrast, rational approximation, particularly through the AAA algorithm, efficiently handles these challenges by incorporating poles into the approximating function. The project begins with a study of polynomial and rational approximations, aiming to compare their performance in terms of accuracy, efficiency, and convergence when applied to functions with challenging features like singularities and poles. The MATLAB package Chebfun (https://www.chebfun.org/) offers powerful tools that can assist in the implementation and testing of these approximation techniques. |
| Lorentzon/Jarlebring | Effektiv Beräkning av Matrisfunktioner | I detta projekt utforskar vi hur klassiska metoder för att beräkna matrisexponetialfunktionen kan förbättras och anspassas genom användning av andra tekniker och tillämpning på andra funktioner, såsom sin(A) och cos(A). Som utgångspunkt använder vi artikeln "Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later". Projektets mål är att implementera och anpassa alternativa tekniker för polynomevaluering, och undersöka deras precision och effektivitet vid beräkning av olika matrisfunktioner, inklusive förbättring av evaluering av sin(A) och cos(A) med hjälp av Eulers formel. |
| Lorentzon/Jarlebring | Analys av Nätverk genom Matrisfunktioner | I detta projekt utforskar vi hur matrisfunktioner, såsom matrisexponetialfunktionen, kan användas för att analysera egenskaper hos nätverk. Genom att utnyttja tekniker från grafteori och linjär algebra fokuserar vi på viktiga aspekter som centralitet och kommunikabilitet i nätverket. Som utgångspunkt använder vi artikeln "Network Properties Revealed through Matrix Functions". Projektets mål är att implementera tekniker för att beräkna matrisfunktioner och tillämpa dem på analys av olika nätverk, med fokus på att undersöka deras tillämpbarhet och effektivitet. |
| Justus Sagemüller | Generative AI through non-Gaussian “diffusion”? | Diffusion models are the backbone of many “AI” image generator systems. The idea behind them is to sample from a distribution (e.g. of photographic images) without explicitly knowing the distribution. Instead, they combine sampling from a simpler distribution (in practice, pure Gaussian noise) with a score function trained from example data, which corrects the generated samples towards realistic images. This procedure has a good theoretical backing from the field of stochastic differential equations, which is however reliant on the Gaussian prior. In many scientific applications such as electron microscopy or astronomy, very similar denoising models are used, but the noise is often not Gaussian but e.g. Poisson distributed. This bachelor project will investigate empirically how much the practical realisation of the popular Stable Diffusion image generator depends on Gaussian noise statistics, respectively how the results change when using various stages of discrete Poisson noise instead. The outcome is of interest for both refinement of theory behind diffusion models, and for applications where it is not sufficient to just generate visually convincing images, but also to fulfill particular statistical properties. |
| Justus Sagemüller | Fast Fourier Transforms and spatial rotations | Fourier series and -transforms (FTs) have many applications. One class of use cases is based on the fact that they make cyclical shifts – in other words, rotations in a suitable coordinate system – a matter of only complex-scalar multiplication. This underlies their utility for uses ranging from audio signal processing to tomographic ray transforms. In practice, the (discrete/discretized) FTs are almost always performed numerically with a Fast Fourier Transform (FFT) algorithm. Usually, this is thought of as an implementation detail for the FT, and FFT algorithms are developed separately from the algorithms that use the Fourier decompositions. However, it is also interesting to consider the recursive structure of the FFT as a processing pipeline of its own right, among other reasons because it has certain similarity to Deep Learning architectures such as U-Nets. This bachelor thesis will investigate how geometric operations on image data manifest themselves in the intermediate representations which an FFT (as it might be used in a ray transform) produces from them. This is of interest both for performance optimisation of such algorithms, as well as possibly generalizations of these algorithms to nonlinear and/or higher-dimensional use cases involving machine learning. |
| Nissan | Numerical methods for wave propagation (two projects) | Everywhere around us we can observe or experience phenomena described by propagating waves. Music and sound are composed by mechanical waves. Electromagnetic waves are used in health care both for diagnostics and therapeutic purposes. Seismic waves are electromagnetic waves that travel through the Earth’s layers after an earthquake. Waves can travel for long times and over distances that far exceeds their wavelengths. High order accurate finite difference (FD) methods are often efficient for propagating waves: high order methods tend to better capture dispersion relations (wave speeds) compared to low order methods and FD methods can be implemented to run efficiently in parallel codes. In order for the methods to be efficient, it is advantageous to use adaptive grids, i.e. numerical grids that are more refined in regions where the solution varies a lot and coarse grids in regions with little activity. The purpose of this project is to improve properties of FD methods applied to partial differential equations (PDE) on adaptive grids. There are two available projects outlined broadly below: Project 1: Improving interpolation. The numerical solution in regions with different grid spacings is coupled using interpolation. You will investigate how to improve the interpolation with respect to improving accuracy of the method, minimizing dispersion errors, minimizing numerical reflections between grids of different grid spacings etc. The findings can be tested and visualized in 2D codes for wave propagation. Project 2: Moving grids. By letting a refined grid move along with a propagating wave, higher accuracy can ideally be obtained at a lower cost compared to when using a uniform grid or an adaptive grid that is stationary in time. You will investigate how to discretize equations on an adaptive grid and how to move the grid and the numerical solution in a robust way. |
| Öktem | High performance numerical linear algebra in medical imaging | Background ---------------- Tomography is a technique to non-invasively image the interior structure of an object. It is based on repeatedly probing the object with particles/waves that interact with the interior structure of the object. from different directions surrounding the object. This data (sinogram) represents noisy indirect observations of the interior structure one seeks to determine. Tomographic imaging therefore involves a mathematical procedure (reconstruction method) for computationally recovering the interior structure from a given sinogram. The most common approach is x-ray computed tomography (CT), which probes the object/patient with x-ray photons. The interior structure (=image) of the object is here represented by the (spatially varying) linear attenuation coefficient, which mathematically is a real-valued function define on 3- or 2-dimensional space. Data (sinogram) consists x-ray intensity measurements along a family of lines passing through the object. The "acquisition geometry" is a precise description of the family of lines one has. One can now adopt a simplified physics model for the interaction of x-rays with matter (Lambert--Beer law). Combining this with appropriate processing on the measured data allows one to recast reconstruction as solving a system of linear equations (each equation corresponds to one measurement, i.e. intensity recorded at the detector of x-ray travelling along a single line through the object). The corresponding matrix (system matrix) is large, but it is also sparse since a line only passes through a few pixels/voxels in the object. Taking the above viewpoint that tomographic image reconstruction amounts to solving a linear system of equations, which opens up for using methods from linear algebra. A difficulty is that the sparse system matrix is very large for realistic data acquisition settings. Furthermore, in many applications, like clinical CT, the linear system is overdetermined (more equations than unknowns), so there are infinitely many possible solutions. The normal setup is then to consider the least-squares solution, but the matrix for the normal equations (which is not sparse) is very large and it often has a high condition number. Hence, it is not desirable to compute the least-squares solution. In some applications, like sparse view industrial CT, one can also end-up with a system that is underdetermined (more unknowns than equations). That setting is even more challenging than the overdetermined setting. In conclusion, one seeks an approximate solution to the linear system that can be computed in a stable way from noisy data. Project ---------- The BSc project aims to implement a number of methods from linear algebra that are designed to approximately solve least squares problems related to large overdetermined linear systems. Implementations of the various methods will be in python/NumPy/SciPy. The specific tasks are as follows: (a) Implement code for generating the system matrix for 2D parallell and fan beam acquisition with full, full but sparse and limited angle acquisition geometries. The system matrix should be stored using sparse matrix routines in the SciPy package. (b) Implement code for reading a 2D grey-scale image, then using (a) to generate noise-free and noisy sinograms for various parallell and fan-beam acquisition geometries. Consider full data and some limited data settings, like sparse view and limited angle acquisition. (c) Implement least squares and truncated singular value decomposition methods for reconstruction, i.e., solving the linear systems whose system matrices are given in (a) with data from (b). (d) Implement Kacmarz and Landweber methods for solving the linear systems whose system matrices are given in (a) and data in (b). The implementations should avoid storing the system matrix in memory. (e) Import system matrix and experimental data from the Helsinki walnut data repository, apply (c)-(d) for reconstruction. (f) Depending on the time/progress, consider implementing variants of Landweber (Cimmino, CAV, DROP, SART, ...) for reconstruction on problems in (a)-(b) and (e). (g) Depending on the time/progress, the project may also involve summarising the theory underlying the reconstruction methods. |
| Öktem | PDE based image denoising & reconstruction | Background ---------------- Imaging is a key technology in society. However, images are unavoidably corrupted by noise, so much effort has been directed towards developing denoising schemes. A specific challenge in denoising is to reduce noise without losing image features, like edges, corners, and other sharp structures. In denoising, the observed data is the noisy image and the goal is to recover the original noise-free image. In a more general setting, observed data is an indirect noisy observation of the original noise-free image. Image reconstruction refers to reverse engineering the clean image from such noisy indirect observations. Project ---------- The BSc project aims to implement a number of methods from PDE theory for denoising/reconstruction of images. Implementations of the various methods will be in the python-package ODL (https://github.com/odlgroup/odl). The specific tasks are as follows: (a) Theory: Derive the PDE given by the Euler-Lagrange equations that are obtained from Tikhonov variational regularisation of the denoising problem, see [1,2]. (b) Theory: Repeat (a) in the setting of total variation regularisation for denoising/reconstruction. (d) Implement numerical solvers for the PDEs in (a)-(b). For finite difference methods, use in Scikit-fdiff, for finite elements methods use FEniCS. Test performance of the various methods on image denoising. (e) Theory: Extend the formulations in (a)-(b) to image reconstruction where observed data is a noisy realisation of a linear transformation (forward operator) of the image one seeks to recover. (d) If time permits, extend solvers in (d) to handle image reconstruction by PDE in (e). Use the python-package ODL to encode the forward operator, which can be tomography or convolution. Test performance of the various PDE methods for image reconstruction on tomographic image reconstruction and image deconvolution. References --------------- [1] https://tristanvanleeuwen.github.io/IP_and_Im_Lectures/variational_formulations.html |
| Ström/Tornberg/Öktem | Can Statistical PDE-solvers overcome the Curse of Dimensionality? | A major difficulty in computational mathematics is problems that involve searching high- or even infinite-dimensional spaces, such as partial differential equations (PDE). Most numerical solvers that exist for PDE can handle up to 3 spatial dimensions, but struggle beyond that since the computational time increases exponentially with dimension. Luckily, many physical problems exhibit some low-dimensional structure (see the "manifold hypothesis"). This is partly because, even if the physical models themselves are formulated very generally by mathematicians, the conditions that they operate under are limited in practice. Recently, machine learning has been shown to scale favourably with dimension on some major high-dimensional problems like image generation (Dall-E, ~ 1million dimensions), text synthesis (GPT-4, ~3 million dimensions) as well as in the hard sciences, with protein folding (Apha-fold, ~1 million dimensions). One of the explanations of why machine learning has been so successful at problems in high dimensions, is that, since machine learning models are specifically constructed from actual data, it can efficiently find the low-dimensional structure of the problems. This project will examine the capacity of machine learning to solve high-dimensional PDE-problems, both from a theoretical, and a practical perspective. We will use a simple yet powerful technique known as random Fourier features, since it is not so complicated to analyse, yet still complicated enough to model complex behaviour. The end goal of the project is to answer the question: Can statistical PDE-solvers overcome the curse of dimensionality, in practice and in theory? On the theoretical side, we will do a literature search on random fourier features and its theoretical guarantees, Physics-Informed Neural Nets (PINN) and attempt to combine the to for solving a PDE problem. On the experimental side, we will look at some toy problems like Poisson, Advection and, if time allows, some more complicated PDE. |
| Ström/Tornberg/Öktem | Statistically Optimal Bases for Linear PDE Problems | A major difficulty in computational mathematics is problems that involve searching high- or even infinite-dimensional spaces, such as partial differential equations (PDE). A common approach to solving such problems numerically utilize a so-called Petrov-Galerkin method, in which the search space is projected to a lower-dimensional subspace, characterized by a finite number of basis vectors. A key difficulty in the Petrov-Galerkin method is how one should choose such a basis. Ideally, a good basis should make the problem easier to solve, which usually means that the matrix that comes from writing the problem as a system of equations should have the following properties: 1) Have a low condition number 2) Be sparse (contain few nonzeros) 3) Be well approximated with only a few basis functions. Typically, the optimal basis that satisfies 1-3, depends on the type of data that is given as input to the problem, and the problem itself. The core question in this project is: If we have access to random Data samples that are interesting for the problem, is it then possible to design a basis that satisfies (1)-(3) in the statistical sense? I.e. that with high probability, (1)-(3) will hold. The project will start simple by examining some existing bases, and specifically how they affect the points (1)-(3) on a simple PDE problem. We will investigate the Fourier basis, as well as some classical FEM-based bases, and the so-called Krylov-basis. We will also touch on some seminal works. We will then construct a data-informed Galerkin projection, and investigate its performance on a set of known problems. Some key problems in the project will be a) How can we quantify to what extent (1)-(3) hold ? b) How does the presence of randomness affect these quantities, and c) Can we control these quantities if we know something about the data? |
| Lithell/Jarlebring |
Projektförslag: Utforskning av AAA-metoden för Rationell Approximation |
AAA-metoden (Adaptive Antoulas–Anderson) är en effektiv algoritm för rationell approximation, särskilt användbar för analytiska funktioner. I detta projekt kommer vi att reproducera resultat från artikeln "The AAA Algorithm for Rational Approximation" och undersöka metodens prestanda för mer utmanande funktioner. Projektmålen är att: Implementera AAA-metoden i MATLAB och reproducera figurer från artikeln. Testa metoden på funktioner med flera singulariteter eller komplexa domäner. Jämföra AAA-metoden med andra rationella approximationsmetoder, som Vector Fitting. Projektet erbjuder en inblick i modern numerisk analys och praktisk tillämpning av rationell approximation. |
