Papers

We revisit the modelling of the material decomposition problem in multispectral computed tomography and propose to leverage the known absorption spectra of different materials. With this we both derive new uniqueness conditions for the material decomposition and also offer a three-step reconstruction model. Our uniqueness criterion only requires certain characteristics of the bin sensitivity functions and the material attenuation spectra and hence, can be computed a-priori before measurements are taken.

Complex-valued neural networks (CVNNs) are increasingly adopted for complex-valued data; however, they are often trained with first-order optimizers inherited from the real-valued case. The efficiency of these methods depends largely on the step size, and their step-size rules ignore the angular information available in the complex plane. We address step-size adaptation in the complex domain by introducing AURA (Angular Update Rate Adaptation), a per-parameter step-size adaptation that can be added on top of any first-order optimizer, and removed from it, without altering its update direction. AURA measures the agreement between consecutive updates of each complex parameter, in length, alignment, and sense of rotation, and enlarges the step when they are consistent and reduces it when they are not. It requires no additional gradient evaluations and only inexpensive vector operations per step. We combine AURA with Adam and Muon and compare the resulting methods with well-known first-order optimizers on four test cases of increasing complexity, ranging from the approximation of scalar complex functions to physics-informed training. Fully connected neural networks are used throughout this work. All hyperparameters other than the step size are held fixed across test cases; for one case, we also tune the hyperparameters of each optimizer under the same budget. Our empirical tests show that AURA improves the convergence of its base optimizer in most cases with a small per-step overhead, and we identify the conditions under which it fails to do so.

In this paper, we investigate the elastic inverse scattering problem in a two-dimensional cavity with an impedance boundary condition. We prove that the scattered fields inside the cavity satisfy a reciprocity relation, and prove that both the shape of cavity and the impedance parameter on the boundary can be uniquely determined from the scattered fields measured on an interior curve. Numerically, we first employ the linear sampling method to reconstruct the shape, and then recover the impedance parameter by combining potential theory with the least squares method. Several numerical experiments are presented to demonstrate the feasibility and effectiveness of our methods.

The multirate partial differential equation (PDE) technique is considered as one of the most powerful methods for the simulation of radio frequency (RF) circuits. The method reformulates the underlying ordinary differential algebraic equation (DAE) by a PDE. The DAE solution is obtained along a characteristic curve of the PDE. Numerical solution techniques comprise (envelope modulated) Harmonic Balance (HB), multistep finite difference (FD) and spline/wavelet techniques. All techniques have disadvantages when mixed lumped/distributed devices occur and the waveforms exhibit pulse shaped behavior. Furthermore, communication systems are modeled and simulated by the Equivalent Complex Baseband (ECB) approach. Signals and impulse responses are described by their baseband counterparts, employing the Hilbert transformation. A coupled system/circuit simulation requires therefore novel numerical techniques for the solution of integro-differential equations.

The tensor product multilevel method (TPML) approximates high-dimensional functions from scattered data on general bounded domains by combining Smolyak's sparse grid construction with the kernel-based multilevel interpolation, coupling the kernel scale to the fill distance at each level so that condition numbers stay uniformly bounded. Its original formulation depended only implicitly on the target function, which precluded numerical evaluation. Using a recent nodal representation of the low-dimensional multilevel operators, we derive an explicit, computable representation of the high-dimensional interpolant. On this basis we prove error estimates in mixed-regularity Sobolev and continuity norms, extending the existing $ L_2 $-theory, and we characterise the full range of anisotropy weights that attain the optimal rate against the number of degrees of freedom. Three explicit weight strategies follow, equilibrating the accuracy, the degrees of freedom, or the cost-benefit ratio. Numerical experiments in up to ten dimensions illustrate the theory.

We study various discretisation schemes for the regular Sturm–Liouville operator on a bounded interval and for the Laplace operator on an arbitrary planar domain, in both cases subject to zero Dirichlet boundary conditions. The objective is to identify a discretisation scheme such that the corresponding discretised operator—a matrix of size $N \times N$—produces $N$ eigenvalues that approximate as closely as possible the first $N$ eigenvalues of the corresponding operator at the continuous level. By means of numerical experiments we examine several conventional discretisation schemes, and we corroborate the known fact that the conventional discretisations fail to achieve the objective, with the failure attributable to the poor approximation of high-index eigenvalues, and, as a result, to the non-uniform spectrum approximation. In contrast, the newly proposed discrete sine transform based discretisation scheme is designed in such a way that it replicates the asymptotic behaviour of high-index eigenvalues, thereby providing the sought uniform spectrum approximation. Of equal significance is the fact that the proposed discrete sine transform based scheme can be, unlike many Fourier transform based methods, applied to non-rectangular domains.

The Scott-Vogelius element is a classical divergence-free mixed finite element for the Stokes problem that has attracted decades of research attention, yet its theoretical framework remains incomplete. In two dimensions, the inf-sup stability on Freudenthal and other regular meshes has been rigorously established. In three dimensions, Zhang established inf-sup stability on Freudenthal meshes for $k\ge6$ in 2011, while numerical evidence indicates that inf-sup stability remains true for $k=4,5$. In this paper, we strengthen Zhang's approach and prove that the Scott–Vogelius element is inf-sup stable on Freudenthal meshes for every velocity degree $k\ge 4$, thereby resolving a conjecture proposed by Farrell, Mitchell, and Scott in 2024. The proof proceeds by explicit constructions on local patches and does not rely on computer verification.

In connection with eigenvalue perturbation bounds, Li [Linear Algebra Appl., 278 (1998), pp. 317-326] formulated three minimization problems involving Hadamard products. The third asks how to relate $$ \min_{W\ \mathrm{nonsingular}}\|W\circ G\|_\mathrm{ui}\,\|W^{-1}\circ G^\mathrm{T}\|_\mathrm{ui} \quad\mbox{and}\quad \min_{P\ \mathrm{permutation}}\|P\circ G\|_\mathrm{ui}^2, $$ where $\|\cdot\|_\mathrm{ui}$ denotes a unitarily invariant norm and $\circ$ denotes the Hadamard product, i.e., entrywise multiplication. The first quantity never exceeds the second, and a reverse inequality with a constant independent of the matrix size was conjectured. Romeo and Tilli [Linear Algebra Appl., 326 (2001), pp. 161-172] showed that if such a bound holds for every $G$, then $\|\cdot\|_\mathrm{ui}$ must be uniformly equivalent to the Frobenius norm $\|\cdot\|_\mathrm{F}$. We prove that the two quantities coincide when $\|\cdot\|_\mathrm{ui}=\|\cdot\|_\mathrm{F}$: for every complex matrix $G$, $$ \min_{W\ \mathrm{nonsingular}} \|W\circ G\|_\mathrm{F}\,\|W^{-1}\circ G^\mathrm{T}\|_\mathrm{F} =\min_{P\ \mathrm{permutation}}\|P\circ G\|_\mathrm{F}^2. $$ The associated conjecture therefore holds with the optimal constant $c=1$.

A controlled comparative study of a mesh-free Kolmogorov–Arnold Physics-Informed Neural Network (KAN-PINN) applied to nonlinear strain-limiting partial differential equations is presented. Three distinct training runs across varying material-parameter pairs $(\alpha, \beta)$ are evaluated within a fixed computational pipeline comprising an exact hard boundary ansatz, a robust pseudo-Huber residual loss, and a multi-stage Adam-to-L-BFGS optimizer schedule. A fundamental metric inversion is revealed by the analysis: while training loss reflects operator-specific residual scaling and constitutive stiffness, the operator-independent relative $L_2$ error demonstrates that intermediate parameter regimes achieve superior solution accuracy and enhanced generalization.

In this paper, we develop a high-order variable-scaled energetic variational neural network (VS-EVNN) method for phase-field gradient flows. Each discretization stage uses an updated-state network whose spatial inputs are centered and divided by the interface scale, while the sampled energy and movement penalties remain unchanged relative to the baseline EVNN formulation. This rescaling is motivated by a coordinate change in the energy that normalizes the gradient coefficient of the transformed functional. We combine the scaled representation with high-order, multi-stage variational extrapolation schemes and a heuristic adaptive time-stepping controller. Each stage is solved as a network minimization problem, in which the earlier stages serve as fixed anchors on the quadrature grid, yielding a neural realization of these multi-stage variational schemes, including a mobility-weighted, mass-conserving version for Cahn–Hilliard flows. The method is tested on one- and two-dimensional Allen–Cahn and Cahn–Hilliard flows, the latter in the mass-conserving metric. Under the reported optimization settings, the scaled networks yield lower solution errors than the unscaled EVNN baseline, and the higher-order schemes improve the accuracy in the test examples.

We prove that the original Barrett–Garcke–Nürnberg (BGN) method for mean curvature flow of curves converges in the manifold distance with order $h^2$ for the time step $\tau=h^2$. The key new technique is backward error analysis, which stems from the numerical analysis of ordinary differential equations. By combining the expansions of the discrete velocity in space and time, we construct an approximation flow whose interpolation satisfies the BGN method with an improved defect, of order $h^4$ with tangential part of order $h^6$. This construction also identifies the tangential velocity of the BGN method and shows that $\tau=O(h^2)$ is the correct relation between $\tau$ and $h$. With the improved defect, we then prove $H^1$ superconvergence of order $h^4$ of the numerical solution to the approximation flow using techniques of evolving surface finite element methods, which yields the convergence in the manifold distance. Numerical experiments confirm the predicted orders of defects.

We consider the problem of reconstructing a kernel $k$ and image $u$ in a blind deconvolution problem from convolutional data $g = k*u$. We particularly focus on the 3D case, i.e., on volumetric images. The problem is approached by imposing a semiparametric CP decomposition for the kernel and a variational TV penalty for the image. We also discuss a motion blur ambiguity effect, i.e., a nonuniqueness issue of the kernel in the blind deconvolution problem that may appear in certain videos of moving objects. The blind deconvolution algorithm uses an alternating minimization approach with TV regularization for $u$ and a CP decomposition step for the kernel. These are implemented with positivity projection, kernel normalization, and\textbackslash or causal support projection steps. Several numerical examples for 3D MRI images, hyperspectral images, and a grayscale test video show the effectiveness of the method.

We study an unsteady fluid-poroelastic structure interaction (FPSI) problem, in which the free fluid region is governed by the incompressible Stokes equations and the poroelastic region by the fully dynamic Biot system. The original physical coupling conditions are equivalently reformulated as Robin-type interface conditions and an auxiliary interface variable is introduced as common Robin data shared by the two subsystems. This reformulation decomposes the coupled FPSI problem into fluid and poroelastic subproblems posed on their respective subdomains. The problem is discretized in space by the finite element method (FEM), while the backward Euler scheme and the BDF2 scheme combined with the second-order interface extrapolation are employed for the temporal discretization. By exploiting the resulting discrete dynamic relations, we derive an update formula for the auxiliary interface variable and obtain fully discrete first- and second-order Robin partitioned schemes in which the two subproblems are solved independently at each time step. The unconditional stability of both schemes is rigorously established. Numerical experiments verify the temporal convergence orders of the proposed schemes and investigate the influence of the Robin parameter. A classical pressure-wave benchmark verifies the accuracy of the proposed methods, while supplementary nonlinear moving-domain simulations illustrate their applicability to more complex FPSI configurations.

Safety evaluation in sampling-based stochastic model predictive control often requires numerical estimation of exit functionals. The approximation of first-exit times and exit indicators is therefore a key numerical bottleneck, and discretization error in these quantities directly affects the resulting controller. This paper studies how existing higher-order methods for strong approximation of exit times can be brought into safe control. Two cases are highlighted. For general noncommutative dynamics, an adaptive order-1 Milstein discretization is used together with Lévy-area simulation via Wiktorsson's method. For commutative dynamics, an adaptive order-1.5 construction achieves a stronger exit-time rate. Under a local anti-concentration condition on the exit-time law, we show that strong exit-time approximation transfers to strong approximation of the failure indicator. The methods are then studied in the context of chance-constrained path integral control, which provides an exact continuous-time representation of safety through exit events. Numerical experiments compare the two cases in terms of strong exit-time error, failure-indicator error, and closed-loop constraint satisfaction, showing improvement over Euler-Maruyama and thereby enabling existing and future techniques whose applicability depends on improved strong approximation.

Variational autoencoders (VAEs) offer an efficient approach to amortized Bayesian inference for inverse problems, but posterior accuracy can depend strongly on the choice of variational regularization, particularly when the inverse problem contains weakly identified parameter directions. This study investigates three objectives: a reverse Kullback–Leibler formulation (VAE-KL), an asymmetric Jensen–Shannon formulation (VAE-JS), and a Jensen–Shannon–Wasserstein formulation (VAE-JSWA), which replaces the reverse Kullback–Leibler regularizer with the squared 2-Wasserstein distance while retaining forward-Kullback–Leibler posterior supervision. A full-covariance Gaussian encoder and a pre-trained physics-based surrogate are used for amortized posterior inference. A local linear–Gaussian analysis in the generalized Fisher basis is developed to characterize the variance-dependent gradients of the three objectives. The formulations are first evaluated using linear–Gaussian benchmarks with known posterior solutions and subsequently tested on nonlinear physics-based inverse problems, including an inverse problem governed by a linear ODE and two PDE-constrained problems. VAE-KL performs slightly better than the other formulations in the well-conditioned benchmark, where all three approaches yield comparable posterior approximations, whereas VAE-JSWA provides substantially lower posterior errors in the strongly ill-conditioned benchmark. The nonlinear physics-based problems exhibit a similar conditioning-dependent trend, with JS-based formulations providing greater benefit as posterior ill-conditioning increases. These results indicate that posterior conditioning is an important factor in selecting variational objectives and motivate geometry-adaptive variational inference for Bayesian inverse problems.

Multiparameter eigenvalue problems arise in boundary value problems, stability analysis, and delay-differential equations. Despite their importance, existing methods either require solving extremely large global problems or rely on local iterative techniques that only recover a few eigenvalues at a time. In this work we develop the first contour method for analytic multiparameter eigenvalue problems. The key theoretical ingredient is a new residue formula for multivariate matrix-valued analytic functions, extending Beyn's Keldysh-based residue theorem for meromorphic operator functions to several complex variables. Using this, we derive a multidimensional analogue of Beyn's contour method that computes all the eigenvalues of a multiparameter eigenvalue problem contained in a prescribed region of $\mathbb{C}^d$. The resulting algorithm targets eigenvalues in a region without constructing preposterously large matrices and is embarrassingly parallelizable. Numerical experiments demonstrate that we can now successfully solve large multiparameter problems arising in applications where existing approaches struggle.

We study the squared Rayleigh-residual functional on complex projective space and its negative gradient flow. For an arbitrary complex matrix, every positive-residual critical point has a strictly descending tangent direction; consequently, the local minima are precisely the eigenlines. For normal matrices, we describe the critical set in terms of circles through spectral points and identify the amplitude dynamics with a replicator equation whose payoff matrix has rank at most four. This gives explicit logarithmic first integrals and a reduced representation of interior trajectories. We compute transverse escape rates near critical components and analyze how polynomial reweighting changes unstable coordinates. Finally, we distinguish local saddle trapping from concentration of Haar-distributed initial states and derive the density of the filtered ensemble. Numerical examples with random spectra and a Jordan block illustrate convergence, nonmonotone Rayleigh motion, and the improvement of finite-time spectral exploration by low-degree random polynomial filters.

High-order differential operators and nonlinear coupling make it challenging to construct conservative and energy-stable schemes for coupled gradient-flow systems. We present GradAgent, a knowledge-guided multi-agent system that coordinates three agents across model analysis, algorithm design and proofs, and numerical implementation and validation. Independent audits strengthen reliability by uncovering mathematical errors and proof gaps, guiding revisions, and maintaining consistency across stages. In Reconstruction Mode, GradAgent reconstructs 20 published studies and organizes audited knowledge in an extensible knowledge graph (KG), GradAgent-KG, linking model structures, discretization strategies and proofs, implementations, and numerical evidence. In Design Mode, the agents assess the applicability of retrieved knowledge and develop new schemes informed by relevant discretization strategies. Applied to the fully coupled multicomponent vesicle phase-field-fluid model, GradAgent yields three first-order and three second-order schemes across three algorithmic families, including four linear, decoupled schemes. Under stated assumptions, all six schemes conserve membrane component mass and vesicle volume and dissipate their respective temporally discrete energies unconditionally. Comparisons with and without GradAgent-KG show that it promotes diversity in structure-preserving scheme design for this target model. Numerical tests confirm second-order spatial accuracy, the expected temporal orders, conservation, and temporally discrete energy dissipation, while three-dimensional shear-flow simulations agree qualitatively with experiments. These results demonstrate GradAgent's ability to combine reusable knowledge, coordinated reasoning, and independent auditing to develop and validate structure-preserving algorithms for complex coupled systems.

A domain derivative for time-dependent electromagnetic scattering from perfect conductors is established. By proceeding through the Laplace domain, frequency-dependent bounds on solutions to Maxwell's equations are derived. These bounds are used both for establishing time regularity properties of the domain derivative and for proving convergence of the proposed Runge–Kutta convolution quadrature semi-discretization in time. A full convergence analysis is also carried out for pointwise evaluations of the domain derivative, when this time discretization is combined with a Galerkin method in space. Numerical experiments validate the convergence orders in time and space. Eventually, the domain derivative is applied in an iterative shape reconstruction algorithm, in which the electric near field is measured at some receiver positions away from the perfect conductor. Numerical examples show the feasibility of this algorithm and in particular highlight its robustness, when additional noise is applied to the data.

A ghost-point immersed boundary method is coupled with auxiliary differential equations to impose a locally-reacting impedance condition on a wall that does not conform to a Cartesian grid, in the framework of the linearized Euler equations. The ghost-point reconstruction follows the normal-stencil approach: an elliptical least squares cloud of fluid nodes is used to build a high-order polynomial along the local wall normal, which is then extrapolated onto the ghost points by one-dimensional La grange interpolation. The impedance, modeled as a massspringdamper oscillator, is advanced in time through two auxiliary state variables per ghost point, coupled to the reconstruction at every RungeKutta stage. The method is validated against exact and semi-analytical references on a flat immersed wall, a cylindrical wall reflecting an acoustic pulse, and a cylinder scattering a harmonic monopole. Convergence orders ranging from about three on the curved wall to about four and a half on the flat wall are measured, and a cross-check against an exact harmonic reference shows that a persistent offset observed against a semi-analytical reference is a property of that ref erence rather than of the coupled solver. A spectral stability analysis of the linearized semi-discrete operator further shows that the normal-stencil reconstruction reduces, without eliminating, a known instability that occurs when the immersed wall becomes tangent to the grid.

We develop an adaptive localized orthogonal decomposition (ALOD) method for Helmholtz problems with parametrized loads. Under standard resolution, oversampling, and discrete-stability conditions, we establish, to our knowledge, the first wavenumber-robust residual a posteriori error estimate for a Helmholtz LOD discretization. The estimate is two-sided relative to the fine-grid solution and includes a computable bound for the corrector-localization error. In addition, we design a novel family-aware adaptive algorithm using additive Schwarz corrections to construct a shared regional enrichment space for Helmholtz problems with multiple loads. Energy-weighted proper orthogonal decomposition compresses this space. We establish stability and quasi-optimality of the enriched formulation, with bounds independent of the enrichment space. Numerical experiments demonstrate the effectiveness of ALOD for multiple loads and large wavenumbers.

The result of Greenbaum, Pták, and Strakoš that for a given set of eigenvalues, any convergence curve is possible [SIMAX 1996] and the subsequent parameterization of such matrix-right-hand side pairs $(A,\mathbf{b})$ of Arioli, Pták, and Strakoš [BIT 1998] demonstrated that the behavior of the \gmres could not be completely characterized by the eigenvalues of $A$ alone. In this paper, we consider how to use this theory to understand the admissible and attainable \gmres behavior for matrices with constrained structure, focussing on non-Hermitian (non-symmetric) tridiagonal Toeplitz matrices. We show that Toeptliz structure necessarily constrains the how the theory from these papers can manifest but that a continuum of admissible behaviors is still attainable.

Discrete gradient methods provide a structure-preserving approach to energy-consistent time integration of dissipative Hamiltonian systems. In their conventional formulation, however, they generally attain at most second-order convergence. Extending the perturbation principle previously developed for Poisson systems to the dissipative setting, we derive fourth-order average vector field discrete gradient schemes. A congruence-based completion preserves the positive semi-definiteness of the dissipation matrix, thereby guaranteeing unconditional energy dissipativity for all positive step sizes. This completion concept provides a systematic route to higher-order energy-consistent schemes. For purely dissipative systems, an alternative exponential completion is derived, together with structure-preserving rational approximations. In particular, we characterize Padé approximants that preserve positive semi-definiteness and can be incorporated into energy-associated splitting schemes. Numerical experiments on a damped physical pendulum and a dissipative Fermi–Pasta–Ulam system demonstrate the importance of the positive semi-definite completion for avoiding nonphysical numerical instabilities and capturing the correct qualitative dynamics. They further show that the proposed methods achieve the desired accuracy and energy consistency at a computational cost competitive with existing Galerkin-type approaches. Moreover, the fourth-order schemes outperform the conventional second-order discrete gradient method in computational efficiency over the considered accuracy range, despite their higher per-step cost.

Katz centrality is a well-established measure to identify and rank the most important nodes in complex networks by means of a linear system solve. Recent works have developed notions of Katz centrality for temporal, i.e., time-evolving networks. Their drawback is that small changes in the network structure may drastically change centralities across large parts of the network. Moreover, an unproportional effort in the network science community has been devoted to the study of node-based quantities while edge-based measures are far less explored. In this manuscript, we introduce a novel failure-aware temporal multiplex network model for computing edge-based Katz centralities for spatial networks. This class of networks admits the use of the line graph as network representation. We use a block-triangular and -banded supra-adjacency matrix representation, modeling inter-layer connections as low-rank matrices, which assigns a spatially constrained emphasis on up- or downdated edges. These could represent, e.g., breakages of water/gas pipes or obstructions in street networks in real-time infrastructure monitoring applications. We analyze the structure of the of the block matrix inverse with respect to entry decay in the blocks above the block-diagonal. This gives rise to two truncation approaches that allow drastic computational runtime reductions at the cost of introducing a controlled truncation error. Numerical experiments on a range of real-world spatio-temporal networks of size up to $2\cdot 10^8$ illustrate accuracy and efficiency with runtime gains of up to a factor of $679$.

We establish arbitrary-order a priori regularity estimates for the reverse heat flow in the Ornstein–Uhlenbeck setting, that is, for the probability-flow ODE of score-based diffusion models. For compactly supported, possibly singular targets satisfying a dimension-uniform doubling condition on supporting-cap masses, we bound the iterated material derivatives of the flow and their first spatial derivatives pointwise, with constants independent of the dimension $d$ and explicit in time and support radius. The proof rests on two ingredients: a spatially uniform bound on normal-direction posterior fluctuations, which is the only point where the geometry of the target enters, and a centered posterior-moment normal form closed under material differentiation. This algebraic structure, together with a weighted regularity calculus, yields quantitative bounds at every order, including the mixed space–time directional derivatives required by Runge–Kutta stages. As an application, for exact-score sampling stopped at forward time $\delta>0$, order-$p$ Taylor and explicit Runge–Kutta schemes reach accuracy $\varepsilon$ in $\widetilde O(d^{1/p}\varepsilon^{-1/p})$ steps in total variation, and the Taylor scheme in $\widetilde O(d^{1/(2p)}\varepsilon^{-1/p})$ steps in $W_2$.

Emulation of high-accuracy floating-point matrix multiplication, as in the Ozaki scheme, splits the inputs into low-precision components and multiplies them pairwise. These products must be error-free, and each is an integer matrix product times a scale factor. FP4 Tensor Cores are the fastest on the NVIDIA B200 and B300 but cannot hold INT8 operands. The FP4 values scaled by 2 form the set $S = \{0, \pm1, \pm2, \pm3, \pm4, \pm6, \pm8, \pm12\}$, which contains every residue modulo 13, so with carries any integer splits into base-13 digits that FP4 can store. Prior work splits each INT8 operand into 3 such digits (limbs) with weights $(1, 13, 169)$ and multiplies them pairwise, 9 FP4 matrix multiplications (GEMMs) for INT8$\times$INT8. The classical ways to reduce products, such as the Karatsuba and Toom–Cook methods, do not apply as they stand: sums of limbs reach $\pm 24$ and leave $S$. This paper asks how many FP4 GEMMs are needed for one integer matrix product. We propose Adaptive Weight Encoding (AWE): the limbs take freely chosen integer weights, the stored planes are linear combinations of limbs, and the exact product is the sum of the FP4 GEMMs scaled by reconstruction coefficients. For each input range, we searched these choices for encodings with fewer products and found INT8$\times$INT8 in 6 products and INT4$\times$INT8 in 4. The formulation also holds modulo $m$, which covers the residue number systems of Ozaki scheme II: for the FP64 significand, the 75 products of prior work are reduced to 59. The boundary in product count between encodings with and without residues lies near input width 15. We release the encodings found.

CUR approximation is a popular approach for constructing an interpretable low-rank approximation of a matrix. There are many approaches for constructing the CUR approximation, many of which are based on the SVD. Results in the literature bounding the quality of the approximation usually rely on the assumption that all or part of the SVD is computed exactly. In this work, we relax this assumption and show how the quality of the approximation is affected by an inexact SVD. Moreover, we analyze the effect of finite precision arithmetic on the DEIM-CUR method. Our analysis suggests guidelines for selecting the working precisions for the SVD and the DEIM index selection algorithm. In particular, it suggests that in many cases one can use very low precisions for these computations without significantly impacting the approximation quality. A set of numerical experiments is presented to demonstrate the theoretical results.

In this paper, we extend numerical schemes based on parameterized Young measures and linear programming to two space dimensions and apply them to the Kelvin-Helmholtz instability governed by the two-dimensional Euler equations of gas dynamics. We construct first-, second-, third-, fifth-, seventh-, and ninth-order Young-measure schemes and compare them with corresponding standard local Lax-Friedrichs (LLF) flux-splitting schemes and LLF schemes equipped with an entropy-based local characteristic decomposition (LLF-ELCD). We examine instantaneous and time-averaged density profiles, cumulative averages across spatial reconstruction orders, regional empirical density distributions, averaged density marginals, and several selection criteria for dissipative weak solutions. The numerical results reveal scheme-dependent flow patterns, particularly in the small-scale structures generated during the roll-up of the shear layers. At every reconstruction order considered, the Young-measure schemes yield the largest time-averaged physical entropy. By contrast, the LLF-ELCD schemes do not systematically yield larger entropy than the standard LLF schemes and therefore do not constitute a consistent maximum-entropy selection mechanism. The averaged density marginals of the Young measures are concentrated on several neighboring density states, and the later cumulative averages are closer for all Young-measure schemes. These results demonstrate the feasibility of the two-dimensional Young-measure formulation and show that the objective function in the linear-programming problem can act as an effective selection mechanism. The observed entropy preference is consistent with the local optimization of the expected physical entropy. It cannot be reproduced by merely incorporating entropy into a conventional numerical construction.

The transposition method known from the theory of second-order elliptic partial differential equations is formulated for saddle-point problems under certain structural assumptions known from the Brezzi splitting plus mapping properties that reflect elliptic regularity in concrete instances of differential operators with boundary conditions. A framework is devised that grants approximation by standard, unmodified Galerkin-type methods, which may be nonconforming with respect to the more general data allowed by the transposed operator. That the data need not be regularized is an advantage particular to the saddle-point approach. Examples covered by the framework include the Laplacian, the Stokes operator, linear elasticity, the bi-Laplacian, and the stationary Maxwell problem with $L^2$ boundary data.

In this paper, we propose and analyze a novel three-field symmetric formulation for optimal control problems governed by linear thermo-poroelasticity. The state variables are the solid displacement $\boldsymbol{u}$, the fluid pressure $p$, and the temperature $\theta$, and both the distributed fluid source $m_p$ and the distributed heat source $m_\theta$ serve as control variables. We prove well-posedness of the symmetric differentiated formulation, establish the existence and uniqueness of an optimal control pair $(\bar{m}_p,\bar{m}_\theta)$, and derive the first-order necessary optimality conditions via a coupled three-field adjoint system. For the fully discrete scheme based on dG(0) time discretization, we employ a variational discretization approach for both controls and derive a priori error estimates of order $\mathcal{O}(h^s + \Delta t)$ for the state, adjoint, and control errors. We also present manufactured-solution experiments that validate the fully discrete optimality system and illustrate the storage-degenerate regime covered by our analysis.