Search

We leverage the proximal Galerkin algorithm (Keith and Surowiec, Foundations of Computational Mathematics, 2024, DOI: 10.1007/s10208-024-09681-8), a recently introduced mesh-independent algorithm, to obtain a high-order finite element solver for variational problems with pointwise inequality constraints. This is achieved by discretizing the saddle point systems, arising from the latent variable proximal point method, with the hierarchical $p$-finite element basis. This results in discretized sparse Newton systems that admit a simple and effective block preconditioner. The solver can handle both obstacle-type, $u \leq \varphi$, and gradient-type, $|\nabla u| \leq \varphi$, constraints. We apply the resulting algorithm to solve obstacle problems with $hp$-adaptivity, a gradient-type constrained problem, and the thermoforming problem, an example of an obstacle-type quasi-variational inequality. We observe $hp$-robustness in the number of Newton iterations and only mild growth in the number of inner Krylov iterations to solve the Newton systems. Crucially we also provide wall-clock timings that are faster than low-order discretization counterparts.

We develop a semismooth Newton framework for the numerical solution of fixed-point equations that are posed in Banach spaces. The framework is motivated by applications in the field of obstacle-type quasi-variational inequalities and implicit obstacle problems. It is discussed in a general functional analytic setting and allows for inexact function evaluations and Newton steps. Moreover, if a certain contraction assumption holds, we show that it is possible to globalize the algorithm by means of the Banach fixed-point theorem and to ensure $q$-superlinear convergence to the problem solution for arbitrary starting values. By means of a localization technique, our Newton method can also be used to determine solutions of fixed-point equations that are only locally contractive and not uniquely solvable. We apply our algorithm to a quasi-variational inequality which arises in thermoforming and which not only involves the obstacle problem as a source of nonsmoothness but also a semilinear PDE containing a nondifferentiable Nemytskii operator. Our analysis is accompanied by numerical experiments that illustrate the mesh-independence and $q$-superlinear convergence of the developed solution algorithm.

We develop a sparse hierarchical $hp$-finite element method ($hp$-FEM) for the Helmholtz equation with rotationally invariant variable coefficients posed on a two-dimensional disk or annulus. The mesh is an inner disk cell (omitted if on an annulus domain) and concentric annuli cells. The discretization preserves the Fourier mode decoupling of rotationally invariant operators, such as the Laplacian, which manifests as block diagonal mass and stiffness matrices. Moreover, the matrices have a sparsity pattern independent of the order of the discretization and admit an optimal complexity factorization. The sparse $hp$-FEM can handle radial discontinuities in the right-hand side and in rotationally invariant Helmholtz coefficients. We consider examples such as a high-frequency Helmholtz equation with radial discontinuities, the time-dependent Schr\"odinger equation, and an extension to a three-dimensional cylinder domain, with a quasi-optimal solve, via the Alternating Direction Implicit (ADI) algorithm.

We show, in one dimension, that an $hp$-Finite Element Method ($hp$-FEM) discretisation can be solved in optimal complexity because the discretisation has a special sparsity structure that ensures that the \emph{reverse Cholesky factorisation} -- Cholesky starting from the bottom right instead of the top left -- remains sparse. Moreover, computing and inverting the factorisation almost entirely trivially parallelises across the different elements. By incorporating this approach into an Alternating Direction Implicit (ADI) method \`a la Fortunato and Townsend (2020) we can solve, within a prescribed tolerance, an $hp$-FEM discretisation of the (screened) Poisson equation on a rectangle, in parallel, with quasi-optimal complexity: $O(N^2 \log N)$ operations where $N$ is the maximal total degrees of freedom in each dimension. When combined with fast Legendre transforms we can also solve nonlinear time-evolution partial differential equations in a quasi-optimal complexity of $O(N^2 \log^2 N)$ operations, which we demonstrate on the (viscid) Burgers' equation.

Exceptionally elegant formulae exist for the fractional Laplacian operator applied to weighted classical orthogonal polynomials. We utilize these results to construct a solver, based on frame properties, for equations involving the fractional Laplacian of any power, $s \in (0,1)$, on an unbounded domain in one or two dimensions. The numerical method represents solutions in an expansion of weighted classical orthogonal polynomials as well as their unweighted counterparts with a specific extension to $\mathbb{R}^d$, $d \in \{1,2\}$. We examine the frame properties of this family of functions for the solution expansion and, under standard frame conditions, derive an a priori estimate for the stationary equation. Moreover, we prove one achieves the expected order of convergence when considering an implicit Euler discretization in time for the fractional heat equation. We apply our solver to numerous examples including the fractional heat equation (utilizing up to a $6^\text{th}$-order Runge--Kutta time discretization), a fractional heat equation with a time-dependent exponent $s(t)$, and a two-dimensional problem, observing spectral convergence in the spatial dimension for sufficiently smooth data.

We prove and collect numerous explicit and computable results for the fractional Laplacian $(-\Delta)^s f(x)$ with $s>0$ as well as its whole space inverse, the Riesz potential, $(-\Delta)^{-s}f(x)$ with $s\in\left(0,\frac{1}{2}\right)$. Choices of $f(x)$ include weighted classical orthogonal polynomials such as the Legendre, Chebyshev, Jacobi, Laguerre and Hermite polynomials, or first and second kind Bessel functions with or without sinusoid weights. Some higher dimensional fractional Laplacians and Riesz potentials of generalized Zernike polynomials on the unit ball and its complement as well as whole space generalized Laguerre polynomials are also discussed. The aim of this paper is to aid in the continued development of numerical methods for problems involving the fractional Laplacian or the Riesz potential in bounded and unbounded domains -- both directly by providing useful basis or frame functions for spectral method approaches and indirectly by providing accessible ways to construct computable toy problems on which to test new numerical methods., Comment: 37 pages, 7 tables, 2 figures

We discuss computing with hierarchies of families of (potentially weighted) semiclassical Jacobi polynomials which arise in the construction of multivariate orthogonal polynomials. In particular, we outline how to build connection and differentiation matrices with optimal complexity and compute analysis and synthesis operations in quasi-optimal complexity. We investigate a particular application of these results to constructing orthogonal polynomials in annuli, called the generalised Zernike annular polynomials, which lead to sparse discretisations of partial differential equations. We compare against a scaled-and-shifted Chebyshev--Fourier series showing that in general the annular polynomials converge faster when approximating smooth functions and have better conditioning. We also construct a sparse spectral element method by combining disk and annulus cells, which is highly effective for solving PDEs with radially discontinuous variable coefficients and data.

We study the finite element approximation of the solid isotropic material with penalization method (SIMP) for the topology optimization problem of minimizing the compliance of a linearly elastic structure. To ensure the existence of a local minimizer to the infinite-dimensional problem, we consider two popular regularization methods: $W^{1,p}$-type penalty methods and density filtering. Previous results prove weak(-*) convergence in the space of the material distribution to a local minimizer of the infinite-dimensional problem. Notably, convergence was not guaranteed to \emph{all} the isolated local minimizers. In this work, we show that, for every isolated local or global minimizer, there exists a sequence of finite element local minimizers that strongly converges to the minimizer in the appropriate space. As a by-product, this ensures that there exists a sequence of unfiltered discretized material distributions that does not exhibit checkerboarding.

We develop a sparse spectral method for a class of fractional differential equations, posed on $\mathbb{R}$, in one dimension. These equations can include sqrt-Laplacian, Hilbert, derivative and identity terms. The numerical method utilizes a basis consisting of weighted Chebyshev polynomials of the second kind in conjunction with their Hilbert transforms. The former functions are supported on $[-1,1]$ whereas the latter have global support. The global approximation space can contain different affine transformations of the basis, mapping $[-1,1]$ to other intervals. Remarkably, not only are the induced linear systems sparse, but the operator decouples across the different affine transformations. Hence, the solve reduces to solving $K$ independent sparse linear systems of size $\mathcal{O}(n)\times \mathcal{O}(n)$, with $\mathcal{O}(n)$ nonzero entries, where $K$ is the number of different intervals and $n$ is the highest polynomial degree contained in the sum space. This results in an $\mathcal{O}(n)$ complexity solve. Applications to fractional heat and wave equations are considered.

Topology optimization problems generally support multiple local minima, and real-world applications are typically three-dimensional. In previous work [I. P. A. Papadopoulos, P. E. Farrell, and T. M. Surowiec, Computing multiple solutions of topology optimization problems, SIAM Journal on Scientific Computing, (2021)], the authors developed the deflated barrier method, an algorithm that can systematically compute multiple solutions of topology optimization problems. In this work we develop preconditioners for the linear systems arising in the application of this method to Stokes flow, making it practical for use in three dimensions. In particular, we develop a nested block preconditioning approach which reduces the linear systems to solving two symmetric positive-definite matrices and an augmented momentum block. An augmented Lagrangian term is used to control the innermost Schur complement and we apply a geometric multigrid method with a kernel-capturing relaxation method for the augmented momentum block. We present multiple solutions in three-dimensional examples computed using the proposed iterative solver.

Divergence-free discontinuous Galerkin (DG) finite element methods offer a suitable discretization for the pointwise divergence-free numerical solution of Borrvall and Petersson's model for the topology optimization of fluids in Stokes flow [Topology optimization of fluids in Stokes flow, International Journal for Numerical Methods in Fluids 41 (1) (2003) 77--107]. The convergence results currently found in literature only consider H^1-conforming discretizations for the velocity. In this work, we extend the numerical analysis of Papadopoulos and Suli to divergence-free DG methods with an interior penalty [I. P. A. Papadopoulos and E. Suli, Numerical analysis of a topology optimization problem for Stokes flow, arXiv preprint arXiv:2102.10408, (2021)]. We show that, given an isolated minimizer of the infinite-dimensional problem, there exists a sequence of DG finite element solutions, satisfying necessary first-order optimality conditions, that strongly converges to the minimizer., Comment: arXiv admin note: text overlap with arXiv:2102.10408

The main objective of this work is to utilize state-of-the-art deep learning approaches for the identification of pulmonary embolism in CTPA-Scans for COVID-19 patients, provide an initial assessment of their performance and, ultimately, provide a fast-track prototype solution (system). We adopted and assessed some of the most popular convolutional neural network architectures through transfer learning approaches, to strive to combine good model accuracy with fast training. Additionally, we exploited one of the most popular one-stage object detection models for the localization (through object detection) of the pulmonary embolism regions-of-interests. The models of both approaches are trained on an original CTPA-Scan dataset, where we annotated of 673 CTPA-Scan images with 1,465 bounding boxes in total, highlighting pulmonary embolism regions-of-interests. We provide a brief assessment of some state-of-the-art image classification models by achieving validation accuracies of 91% in pulmonary embolism classification. Additionally, we achieved a precision of about 68% on average in the object detection model for the pulmonary embolism localization under 50% IoU threshold. For both approaches, we provide the entire training pipelines for future studies (step by step processes through source code). In this study, we present some of the most accurate and fast deep learning models for pulmonary embolism identification in CTPA-Scans images, through classification and localization (object detection) approaches for patients infected by COVID-19. We provide a fast-track solution (system) for the research community of the area, which combines both classification and object detection models for improving the precision of identifying pulmonary embolisms., Comment: 16 pages, 6 figures, 1 table, Submitted to the European Radiology journal of Springer

T. Borrvall and J. Petersson [Topology optimization of fluids in Stokes flow, International Journal for Numerical Methods in Fluids 41 (1) (2003) 77--107] developed the first model for topology optimization of fluids in Stokes flow. They proved the existence of minimizers in the infinite-dimensional setting and showed that a suitably chosen finite element method will converge in a weak(-*) sense to an unspecified solution. In this work, we prove novel regularity results and extend their numerical analysis. In particular, given an isolated local minimizer to the infinite-dimensional problem, we show that there exists a sequence of finite element solutions, satisfying necessary first-order optimality conditions, that strongly converges to it. We also provide the first numerical investigation into convergence rates.

Topology optimization problems often support multiple local minima due to a lack of convexity. Typically, gradient-based techniques combined with continuation in model parameters are used to promote convergence to more optimal solutions; however, these methods can fail even in the simplest cases. In this paper, we present an algorithm to perform a systematic exploratory search for the solutions of the optimization problem via second-order methods without a good initial guess. The algorithm combines the techniques of deflation, barrier methods and primal-dual active set solvers in a novel way. We demonstrate this approach on several numerical examples, observe mesh-independence in certain cases and show that multiple distinct local minima can be recovered.

Purpose: Recent work highlights the breadth of reported spin-lattice relaxation rates ($R_1$) for individual tissues. One potential source of variation is the protocol used to determine $R_1$. The methodological dependence of R1 and relaxivity $r_1$ are investigated. Methods: $R_1$ is quantified in gel phantoms with varying concentration of MnCl2, and a small cohort of three healthy volunteers using different acquisition methods. Siemens inversion recovery (IR) and saturation recovery (SR) protocols are applied to phantoms and volunteers. Variable flip angle (VFA) protocols are additionally applied to phantoms. $R_1$ is quantified using single voxel fits, and distributions examined for regions in the thalamus, and cerebellum as well as grey and white matter. Phantoms exclude boundary fits and relaxivity is quantified across the full concentration range. Normality of $R_1$ distributions is assessed by Kolmogorov-Smirnov score, and inter-sequence agreement by two-sample t-test. Results: Phantom relaxivity is found to be 7.16 Hz/mM, 9.22 Hz/mM and 10.65 Hz/mM to 11.91 Hz/mM for IR, SR and VFA methods, respectively. In vivo $R_1$ exhibit low intra-participant variation for IR. SR $R_1$ are lower than IR values with inter- and intra-participant variation on the same order. Brain regions and phantoms mapped with different protocols varied significantly with t-test p-values between 0 and 5E-10. Conclusion: Results suggest a significant protocol dependence of $R_1$, and corresponding relaxivity, suggesting inter-method comparisons should be attempted tentatively, if at all., Comment: 15 pages (5 figures, 5 tables), 3 pages supplementary material (6 figures)

Data from a multi-parametric MRI study of patients with possible early-stage prostate cancer was assessed with a view to creating an efficient clinical protocol. Based on a correlation analysis suggesting that diffusion-weighted imaging (DWI) scores are more strongly correlated with overall PIRADS scores than other modalities such as dynamic contrast enhanced imaging or spectroscopy, we investigate the combination of T2-weighted imaging (T2w) and DWI as a potential diagnostic tool for prostate cancer detection, staging and guided biopsies. Quantification of the noise floor in the DWI images and careful fitting of the data suggests that the mono-exponential model provides a very good fit to the data and there is no evidence of non-Gaussian diffusion for $b$-values up to 1000 s/mm$^2$. This precludes the use of kurtosis or other non-Gaussian measures as a biomarker for prostate cancer in our case. However, the ADC scores for healthy and probably malignant regions are significantly lower for the latter in all 20 but one patient. The results suggest that a simplified mp-MRI protocol combining T2w and DWI may be a good compromise for a cost and time efficient, early-stage prostate cancer diagnostic programme, combining robust MR biomarkers for prostate cancer that can be reliably quantified and appear well-suited for general clinical practice., Comment: to appear in MRI

The controlled modulation of an optical wavefront is required for aberration correction, digital phase conjugation or patterned photostimulation. For most of these applications it is desirable to control the wavefront modulation at the highest rates possible. The digital micromirror device (DMD) presents a cost-effective solution to achieve high-speed modulation and often exceeds the speed of the more conventional liquid crystal spatial light modulator, but is inherently an amplitude modulator. Furthermore, spatial dispersion caused by DMD diffraction complicates its use with pulsed laser sources, such as those used in nonlinear microscopy. Here we introduce a DMD-based optical design that overcomes these limitations and achieves dispersion-free high-speed binary phase modulation. We show that this phase modulation can be used to switch through binary phase patterns at the rate of 20 kHz in two-photon excitation fluorescence applications.

Manipulating the propagation of light through scattering media remains a major challenge for many applications, including astronomy, biomedical imaging and colloidal optics. Wavefront shaping is one of the most promising ways to mitigate scattering and focus through inhomogeneous samples. However, wavefront correction remains accurate over only a limited spatial extent within the scattering medium - a correlation range referred to as the optical memory effect. Here, by selecting only the weakly scattered light for wavefront shaping, we show that the addition of temporal degrees of freedom enhances this correlation range. We investigate spatial scattering correlations by digitally time-gating the early arriving light in the spectral domain. We demonstrate that the range of the translational memory effect for the early arriving light is increased almost fourfold, paving the way for a range of scattering media imaging applications.

The optical memory effect is a well-known type of wave correlation that is observed in coherent fields that scatter through thin and diffusive materials, like biological tissue. It is a fundamental physical property of scattering media that can be harnessed for deep-tissue microscopy or 'through-the-wall' imaging applications. Here we show that the optical memory effect is a special case of a far more general class of wave correlation. Our new theoretical framework explains how waves remain correlated over both space and angle when they are jointly shifted and tilted inside scattering media of arbitrary geometry. We experimentally demonstrate the existence of such coupled correlations and describe how they can be used to optimize the scanning range in adaptive optics microscopes., Comment: 11 pages, 4 figures, supplementary materials included

In this paper, we propose a new technique for two-dimensional phase unwrapping. The unwrapped phase is found as the solution of an inverse problem that consists in the minimization of an energy functional. The latter includes a weighted data-fidelity term that favors sparsity in the error between the true and wrapped phase differences, as well as a regularizer based on higher-order total-variation. One desirable feature of our method is its rotation invariance, which allows it to unwrap a much larger class of images compared to the state of the art. We demonstrate the effectiveness of our method through several experiments on simulated and real data obtained through the tomographic phase microscope. The proposed method can enhance the applicability and outreach of techniques that rely on quantitative phase evaluation.

Acquiring high-contrast optical images deep inside biological tissues is still a challenging problem. Confocal microscopy is an important tool for biomedical imaging since it improves image quality by rejecting background signals. However, it suffers from low sensitivity in deep tissues due to light scattering. Recently, multimode fibers have provided a new paradigm for minimally invasive endoscopic imaging by controlling light propagation through them. Here we introduce a combined imaging technique where confocal images are acquired through a multimode fiber. We achieve this by digitally engineering the excitation wavefront and then applying a virtual digital pinhole on the collected signal. In this way, we are able to acquire images through the fiber with significantly increased contrast. With a fiber of numerical aperture 0.22, we achieve a lateral resolution of 1.5um, and an axial resolution of 12.7um. The point-scanning rate is currently limited by our spatial light modulator (20Hz).

We describe a method for imaging 3D objects in a tomographic configuration implemented by training an artificial neural network to reproduce the complex amplitude of the experimentally measured scattered light. The network is designed such that the voxel values of the refractive index of the 3D object are the variables that are adapted during the training process. We demonstrate the method experimentally by forming images of the 3D refractive index distribution of cells.

Controlling light propagation across scattering media by wavefront shaping holds great promise for a wide range of communications and imaging applications. However, finding the right wavefront to shape is a challenge when the mapping between input and output scattered wavefronts (i.e. the transmission matrix) is not known. Correlations in transmission matrices, especially the so-called memory-effect, have been exploited to address this limitation. However, the traditional memory-effect applies to thin scattering layers at a distance from the target, which precludes its use within thick scattering media, such as fog and biological tissue. Here, we theoretically predict and experimentally verify new transmission matrix correlations within thick anisotropically scattering media, with important implications for biomedical imaging and adaptive optics., Comment: main article (18 pages) and appendices (6 pages)

Light propagation through multimode fibers suffers from spatial distortions that lead to a scrambled intensity profile. In previous work, the correction of such distortions using various wavefront control methods has been demonstrated in the continuous wave case. However, in the ultra-fast pulse regime, modal dispersion temporally broadens a pulse after propagation. Here, we present a method that compensates for spatial distortions and mitigates temporal broadening due to modal dispersion by a selective phase conjugation process in which only modes of similar group velocities are excited. The selectively excited modes are forced to follow certain paths through the multimode fiber and interfere constructively at the distal tip to form a focused spot with minimal temporal broadening. We demonstrate the delivery of focused 500 fs pulses through a 30 cm long step-index multimode fiber. The achieved pulse duration corresponds to approximately 1/30th of the duration obtained if modal dispersion was not controlled. Moreover, we measured a detailed two-dimensional map of the pulse duration at the output of the fiber and confirmed that the focused spot produces a two-photon absorption effect. This work opens new possibilities for ultra-thin multiphoton imaging through multimode fibers.

We present a novel design that exploits the focusing properties of scattering media to increase the resolution and the working distance of multimode fiber based imaging devices. Placing a highly scattering medium in front of the distal tip of the multimode fiber enables the formation of smaller sized foci at increased working distances away from the fiber tip. We perform a parametric study of the effect of the working distance and the separation between the fiber and the scattering medium on the focus size. We experimentally demonstrate submicron focused spots as far away as 800{\mu}m with 532nm light., Comment: 4 pages, 3 figures

We demonstrate Optical-Resolution Photoacoustic Microscopy (OR-PAM), where the optical field is focused and scanned using Digital Phase Conjugation (DPC) through a multimode fiber. The focus is scanned across the field of view using digital means, and the acoustic signal induced is collected by a transducer. Optical-resolution photoacoustic images of a knot made by two absorptive wires are obtained and we report on resolution smaller than 1.5{\mu}m across a 201{\mu}m by 201{\mu}m field of view. The use of a multimode optical fiber for the optical excitation part can pave the way for miniature endoscopes that can provide optical-resolution photoacoustic images at large optical depth., Comment: 10 pages, 3 figures