Papers

Quantum incommensurate systems have attracted widespread interest due to their unique physical properties. Related studies have made notable progress in recent years. To gain deeper insight, it is both significant and challenging to study a broader range of physical observables for these systems, which requires a comprehensive understanding of their spectral properties and wavefunction behavior. Based on the regularized model recently proposed, this work introduces a regularization framework, rendering physical observables for incommensurate systems mathematically well-defined and computationally accessible with theoretical guarantees. Based on this framework, taking the density of states and the electron density as representative examples, we establish their mathematical characterization, derive their planewave approximations, and provide a rigorous convergence analysis. Numerical experiments validate the effectiveness of our method, demonstrating its practical applicability for generic 1D and 2D incommensurate systems.

We propose $p$ and $rp$-spectral element methods for elliptic boundary layer problems on two dimensional rectangular domains. To resolve boundary layers, we use $p$-version with a fixed mesh and an $rp$-version with a boundary layer mesh consisting of thin needle-like elements near the boundary layer and coarse elements away from the layer. Stability estimates are derived using non-conforming spectral element functions. The numerical scheme for both the methods is based on minimising a least-squares functional in appropriate Sobolev norms. We construct robust preconditioners to manage the condition number of the normal equations derived from the least-squares formulation. The method is able to approximate boundary layers at a rate $O\left(\frac{\log W}{W^2}\right)$ for the $p$-version and at the rate $O({\epsilon}\alpha^{2W})$, uniformly in $\epsilon$, for the $rp$-version, where $0<\epsilon\leq1$ is the boundary layer parameter, $\alpha<1$ is a constant and $W$ denotes the degree of the approximating polynomial. Numerical results are provided for model elliptic boundary layer problems for a range of boundary layer parameters. Simulation results demonstrate the efficiency and robustness of the method in capturing boundary layers on rectangular domains.

For every integer $n \ge 2$, we construct explicit tetrahedral counterexamples showing that repeated longest-edge (LE) n-section refinement does not, in general, preserve shape regularity, contrary to long-standing conjectural expectations. For $n = 2$, the construction disproves the finite-family conjecture put forward by Adler in 1983 and the non-degeneracy conjecture for tetrahedral longest-edge bisection formulated by Rivara and Levin in 1992. It also shows that the sharper higher-dimensional diameter decay suggested by Stynes in 1983,following his 1980 planar finite-similarity result, cannot hold for arbitrary tetrahedra. For $n = 3$, it disproves the published non-degeneracy conjecture for tetrahedral longest-edge trisection formulated in 2011 by Suárez, Abellón, Abad and Plaza on the basis of numerical experiments. The resulting infinite descendant sequences violate both the minimum and maximum angle conditions. For $n \ge 3$, the selected edge is always the unique longest edge at every refinement step. The two-step recurrence used in the construction exhibits two distinct degeneration behaviors: a flat-tetrahedron regime for $2 \le n \le 5$ and a skinny-tetrahedron regime for $n \ge 6$. The degenerating sequences can also be realized within conforming partitions generated by the conforming LE n-section algorithm. For both the classical and the conforming LE n-section algorithms, the maximal element diameter tends to zero as the number of refinement steps tends to infinity, regardless of which longest edge is chosen. Thus diameter convergence, and even conformity, do not prevent shape degeneration. These algorithms should therefore be used with care in applications requiring uniform shape regularity.

The constant-Q transform (CQT) represents audio on a logarithmic frequency axis. Its nonstationary Gabor formulation is exactly invertible, but the unequal numbers of time coefficients in its bands complicate GPU computation. An exact factorisation combines spectral selection, conjugation, windowing, and reordering into a fixed map between one packed Fourier transform and the shorter band inverse transforms. The factors give waveform reconstruction, real adjoints for backpropagation, and bounds on arithmetic depth and block width; overlapping slices permit streaming with bounded memory. Tests on two GPU models show that Flash-CQT reduces analysis-synthesis round-trip time by factors of two to eight relative to a baseline computing the same CQT. The proposed implementation also uses over 30% less peak temporary workspace and reaches a negligible reconstruction error, with a signal-to-noise ratio of about 130 dB, in single-precision floating-point arithmetic. These advances make Flash-CQT a practical, computationally efficient front end for spectral analysis and modern audio machine-learning systems.

This paper develops a compact fourth-order positivity-preserving active flux (AF) method for the one- and two-dimensional compressible Navier–Stokes equations on Cartesian meshes. The method retains the cell averages and shared point values of the standard third-order AF method as its degrees of freedom. To avoid the order reduction that can arise when diffusion is discretized using operators from the standard third-order AF method, while maintaining compactness, the divergence of the viscous flux is discretized directly using compact fourth-order operators. For the inviscid part, incorporating a downwind point value into the biased stencil yields fourth-order accuracy. A monolithic flux limiting blends high-order total numerical fluxes with low-order positivity-preserving counterparts, treating the inviscid and viscous fluxes jointly while maintaining local conservation. Together with a scaling limiter for point values, this procedure preserves density and pressure positivity for both cell averages and point values. Numerical experiments demonstrate fourth-order convergence, positivity preservation, and accurate resolution of shocks and viscous flow structures. For the two-dimensional viscous shock tube, a comparison with a discontinuous Galerkin method shows improved computational efficiency in resolving complex interaction of shock waves and boundary layers.

We develop and analyze HDG methods for two central problems in finite element exterior calculus, the Hodge-Dirac problem and the Hodge-Laplace problem, in arbitrary dimension $n$. Our analysis allows for equal-order polynomial spaces with $\mathcal{O}(1)$ penalty parameters, by contrast with previous work on the Hodge-Laplace problem requiring an underlying conforming complex and $\mathcal{O}(h)$ or $\mathcal{O}(h^{-1})$ penalties. For the Hodge-Dirac problem with degree-$r$ polynomials, our error estimates give optimal order-$(r+1)$ convergence in all form degrees under suitable regularity hypotheses. For the Hodge-Laplace problem, we prove optimal order-$(r+1)$ convergence for the $k$-form solution and order-$(r+\frac{1}{2})$ convergence for the $(k \pm 1)$-form auxiliary variables. We obtain improved order-$(r+1)$ auxiliary-variable estimates when $k=1$ and/or $k=n-1$, which in particular sharpens some recent HDG error estimates for the two-dimensional vector Poisson equation. Our analysis also encompasses cases of lower regularity, including solutions with reentrant-corner singularities on non-convex domains, and requires only mild mesh regularity conditions. The results are illustrated by numerical experiments in dimensions two and three.

This paper is devoted to the non-diffusive computation of solutions to $N$-dimensional hyperbolic systems of conservation laws (HCL) and more specifically the accurate approximation of simple waves. We develop the LeafNet algorithm, which relies on: i) a reformulation of HCL as a coupled system involving smooth solutions, and ii) an accurate physics informed algorithms approximating smooth functions. An error estimate analysis and two-dimensional numerical experiments are proposed to illustrate the presented strategy.

Quasi-brittle fracture is prevalent in concrete, rock, ceramics, composites, and masonry, and its simulation faces a trade-off among accuracy, efficiency, and simplicity. The classical Lattice Spring Model (LSM) captures cracking via bond breakage without remeshing, but its elements are empirical and limited to a few tessellable shapes with fixed Poisson's ratios. We propose a tessellated Isotropic Elastic Lattice Spring Model (IELSM) that discretizes continua into polygonal elements with axial springs and a volumetric constraint, achieving isotropic elasticity on arbitrary polygonal tessellations. Macroscopic isotropy reduces to governing equations whose solvability gives a theoretical criterion for element admissibility, proving conventional tessellable elements and extending to arbitrary regular N-gons and concave elements. Exploiting boundary interpolation compatibility with finite elements, IELSM is assembled by direct node sharing, without interface elements or kinematic constraints. Coupled with an isotropic damage model, a pure bending test and four fracture benchmarks show that the coupling preserves displacement accuracy, yields crack paths and load-displacement curves agreeing with experiments and outperforming standard FEM, and is insensitive to mesh refinement. IELSM can also be restricted to damage-prone regions, with the remainder modeled by finite elements via node sharing. For the benchmarks, this reduces nodes by 41.2-80.9% and CPU time by 34.2-85.2% versus full-domain IELSM. The framework advances IELSM element construction from empirical trial and error to theoretical determination and simplifies coupling to node sharing, offering a balanced route for complex quasi-brittle fracture analysis.

A novel second-order consistent splitting scheme is proposed for the natural convection equations. The scheme is constructed using Taylor expansions about the time level $t^{n+k}$, where $k$ is a parameter to be determined. It is proved to be stable for all $k>4$, with the present study focusing primarily on the case $k=5$. Compared with the conventional splitting scheme based on Taylor expansions about $t^{n+1}$, the proposed scheme exhibits improved stability. By employing the Sobolev inequality and Gronwall's lemma, we rigorously establish stability of the proposed scheme and derive error estimates in both two and three spatial dimensions. Finally, several numerical experiments are presented to verify the stability and accuracy of the proposed scheme and to demonstrate its effectiveness.

A finite difference numerical scheme is proposed and analyzed for the Cahn-Hilliard-Navier-Stokes system, combined with a dynamical boundary condition. Such a physical system has potential applications in the moving contact line problem. The boundary profile is governed by a lower-dimensional energy potential, coupled with a non-homogeneous boundary condition for the phase variable. In the numerical design, a convex-splitting approach is applied to the chemical potential in both bulk and surface levels, which leads to a highly coupled nonlinear system. A semi-implicit discretization is taken in the nonlinear fluid convection, as well as the coupled terms between the fluid motion and phase variable evolution. A careful finite difference approximation and convexity analysis reveals that such a numerical system could be represented as a non-symmetric and monotone mapping associated with the fluid convection. In turn, the unique solvability is valid based on the monotonicity argument. The total energy stability analysis is obtained through a careful summation-by-part calculation. In particular, an optimal rate convergence analysis is theoretically established in this work. The discrete mass conservation of the exact solution is required to preserve the mean-zero property of the error function, so that the associated discrete H_h^{-1} norm is well-defined. A combination of the Fourier projection and an auxiliary function is applied to overcome this difficulty. Furthermore, an approach of rough and refined error estimates concludes the desired convergence result. Some numerical results are presented in this article, which demonstrate the robustness of the proposed numerical scheme. In our knowledge, this work provides a theoretical proof of convergence analysis and error estimate for a numerical scheme to the moving contact line problem, for the first time in the literature.

Anderson acceleration (AA) is an effective technique for accelerating fixed point iterations, but it generally lacks global convergence guarantees. We propose a nonmonotone globalized Anderson acceleration framework that retains the local acceleration of AA while ensuring global convergence. For contractive mappings, without requiring prior knowledge of the contraction factor, we prove that the proposed method reduces exactly to pure AA after finitely many iterations and enjoys global $r$-linear convergence for memory size $m\geq1$ and global residual $q$-linear convergence for $m=1$, with convergence factors no greater than the contraction factor of the underlying fixed point mapping. For nonexpansive mappings, we prove that the fixed point residuals converge globally to zero. A class-agnostic parameter selection rule is further developed to ensure these convergence properties in both settings. Numerical experiments demonstrate the robustness and efficiency of the proposed method.

We study rotational invariance, hyperbolicity, and well-balanced discretization of shallow water linearized moment equations in two horizontal dimensions. Although the direct extension has real characteristic speeds, its directional matrices can lack a complete eigenbasis. We construct a rotationally invariant modification and prove its global hyperbolicity for all states with positive water depth. For quasi-two-dimensional flows over non-flat topography with Navier-slip friction, we develop first- and second-order path-conservative finite-volume schemes that exactly preserve a discrete moving equilibrium constructed by midpoint collocation. Numerical tests confirm equilibrium preservation to roundoff, second-order approximation of smooth stationary branches, and near-second-order accuracy of the second-order scheme for the tested perturbations. Comparisons with three-dimensional two-phase OpenFOAM simulations illustrate the invariant shallow-water limit and show that retaining moments improves the reconstruction of curved rotational velocity profiles, while differences remain in radial motion.

We initiate the study of approximating the top eigenvalue and eigenvector of a random symmetric matrix $ A \in \mathbb{R}^{n\times n} $ using $ q(A)b $ where $q$ is a degree-$d$ polynomial and $b$ is a standard Gaussian vector independent of $A$. For spiked GOE $ Y = \lambda vv^\top + X $, we identify $ d_\star = \frac{\log(n)}{2\log(\lambda)} $ to be the critical degree threshold above which accurate approximation of the top eigenvalue and eigenvector is possible. This sharpens the common belief that spectral methods can be implemented by $ O(\log(n)) $-step power iterations and offers a precise connection between spectral methods and low-degree polynomial algorithms, a popular proxy for all polynomial-time algorithms. For GOE $X$, we identify $ d_\star = n^{1/3+o(1)} $ to be the critical degree threshold for top eigenvector approximation, whereas constant degree suffices for top eigenvalue approximation. Moreover, in the limit where $ d/n^{1/3} $ converges to a positive finite constant, we compute the exact asymptotic eigenvector approximation accuracy in terms of the expected squared overlap. These results significantly improve upon predictions made in randomized numerical linear algebra for deterministic data matrices that the iteration count of power methods with random initialization is governed by the inverse spectral gap. Technically, our analyses leverage extremal properties of Chebyshev polynomials and draw upon the rich literature of random matrix theory.

In recent years a number of finite element methods have been developed for the surface Stokes equations. Constructing standard conforming surface finite element methods for the velocity-pressure formulation of this problem would require simultaneously enforcing $H^1$ conformity and tangentiality of the velocity field, which is not possible on the generalized polyhedral approximating surfaces on which surface FEM are typically posed. Thus methods for this problem weakly enforce either tangentiality or $H^1$ conformity of the velocity field. In this paper we provide an error analysis for an $H({\rm div})$-conforming method based on classical Brezzi-Douglas-Marini elements. In this method continuity of the velocity field is weakly enforced using an interior penalty-type formulation. Advantages of this formulation include exact enforcement of the incompressibility constraint and its construction based on a standard mixed finite element pair. Error analysis of this method in previous work was only valid for lowest-order elements and did not account for “geometric errors”, i.e., variational crimes due to approximation of the continuous surface on which the PDE is posed by an approximating surface on which the finite element method is posed. Here we provide a full error analysis in a natural discontinuous Galerkin energy norm, with geometric errors accounted for. Optimal estimates for the $L_2$ velocity error are also presented. In the case of lowest-order (polyhedral) surface approximations, our proof of optimal $L_2$ estimates requires a stricter than usual condition on the placement of mesh nodes relative to the surface. Numerical experiments are used to illustrate our results.

We study global convergence guarantees of third-order Langevin dynamics for non-convex optimization via simulated annealing with fixed friction and decreasing noise. An explicit three-block distorted entropy transfers dissipation from the noisy auxiliary variable to the full state. Under dissipativity, regularity, and low-temperature functional-inequality assumptions, logarithmic cooling drives the objective values to the global minimum in probability at the barrier-controlled kinetic rate. For the exact-force-integral and midpoint three-stage discretizations, polynomially decreasing steps preserve this rate on the physical time scale. The cubic local endpoint estimate gives a less restrictive sufficient step-size condition than the available frozen-force kinetic result. A comparison with the one-gradient UBU integrator shows how its centered stochastic local error leads, under the same strong-coupling analysis, to a smaller sufficient iteration exponent. Numerical experiments are conducted to illustrate our theory. For a double well objective, third-order Langevin terminal-success point estimates are higher than UBU at both a common horizon and an equal gradient budget. For a high-dimensional nonconvex neural-network objective using synthetic data, independently tuned UBU and third-order Langevin schemes both outperform overdamped Langevin dynamics; the third-order Langevin point estimate is higher. For the same neural-network objective on real data, we show the same point-estimate ordering for best-basin probability and post-quench test accuracy. Numerical code and associated experiment results are publicly available at https://github.com/gagawjbytw/simulated-annealing-third-order-langevin.

The Girard–Hutchinson estimator provides an extremely simple randomized estimate of the Frobenius norm of a matrix $A$ that can only be accessed implicitly via matrix-vector products. In particular, if $\Omega$ is a random Gaussian matrix with $r = O(1/\varepsilon^2)$ columns, than $\frac{1}{r}\|A\Omega\|_F^2$ provides a $(1\pm \varepsilon)$ multiplicative approximation to $\|A\|_F^2$ with high probability. In this work, we introduce a closely related estimator, given by \begin{align*} {\frac{1}{r}\|A\Omega\|_F^2 + \frac{1}{r}\|\Psi^T A\|_F^2 - \frac{1}{r^2}\|\Psi^T A\Omega\|_F^2}, \end{align*} where $\Psi$ is a second, independent random Gaussian matrix with $r$ columns. We prove that this estimator yields a $(1\pm\varepsilon)$ multiplicative approximation to $\|A\|_F^2$ when $r = O(1/\varepsilon)$, a quadratic improvement over Girard–Hutchinson. This dependence on $\varepsilon$ is optimal. Our method, which we call Hutch# (pronounced “Hutch sharp”), matches the complexity of the Hutch++ algorithm [Meyer, Musco, Musco, Woodruff, 2021]. However, unlike Hutch++, Hutch# uses only non-adaptive matrix-vector products with $A$ and $A^T$ and requires no orthogonalization or other advanced linear algebra steps. Thus, Hutch# combines the simplicity of the Girard–Hutchinson estimator and the optimal query complexity of Hutch++.

We propose and analyze a fully discrete finite element scheme for the Beris-Edwards system of nematic liquid crystal dynamics, in which the incompressible Navier-Stokes equations are coupled to a gradient flow for the Landau-de Gennes Q-tensor. The scheme combines a linearly implicit formally second-order accurate backward differentiation formula in time with an incremental Chorin projection step for the incompressibility constraint, conforming finite elements in space, and the invariant energy quadratization approach with mass lumping for the nonlinear bulk potential. Each time step requires the solution of one linear system and one Poisson problem. We show that the scheme is uniquely solvable, that it preserves the symmetry and trace-free structure of the discrete Q-tensor and molecular field, and that it satisfies a discrete energy law without any restriction on the time step. Our main result is that, as the mesh size $h$ and the time step $\Delta t$ tend to zero subject to $h^{2} = o(\Delta t)$, the approximations converge along a subsequence to a weak solution of the Beris-Edwards system. The convergence proof addresses two difficulties: the projection method produces two velocity approximations, only one of which is uniformly bounded in $L^2(0,T;H^1_0(\Omega))$, and the coupling term $\mathcal{H}\nabla Q$ requires strong convergence of $\nabla Q$, which we obtain from the structure of the equation for $\mathcal{H}$ rather than from any discrete $H^2$-bound. Numerical experiments in two dimensions exhibit approximately second-order convergence in space and time, and reproduce the splitting of a $+1$ point defect into two $+1/2$ defects and the transport and deformation of a skyrmion induced by a constant pressure gradient.

In this paper, we introduce a novel second-order, positivity-preserving finite element method for the Vlasov equations using a convex limiting algorithm. The method employs strong-stability-preserving (SSP) Runge-Kutta time integration and a tensor-product construction of the phase-space mesh for efficient high-dimensional implementations. The convex limiting algorithm combines the robust first-order positivity-preserving graph viscosity approach with high-order residual-based viscosity stabilization to obtain a high-order positivity-preserving scheme. Both novel first-order and high-order methods applicable to high-dimensional problems such as the Vlasov system are presented. In addition, we propose a divergence-cleaning technique for Maxwell's equations to ensure that the divergence constraints of the electromagnetic fields are satisfied. Numerical experiments are provided to demonstrate the accuracy and robustness of the proposed methods.

It is well known that the choice of basis functions in hp-FEM heavily influences the stability and the computational cost in order to obtain an approximate solution. For simplicial elements in two and three space dimensions, tensor-product-like basis functions built from Jacobi polynomials with different weights yield optimal properties with respect to condition number and sparsity. In this paper we construct such a high order basis for the Nédélec spaces and for the Raviart-Thomas and Brezzi-Douglas-Marini spaces, modifying existing definitions of Zaglmayr so that the functions belong to the Nédélec spaces of first kind and the Raviart-Thomas spaces. The bases are designed so that the Nédélec space of second kind and the BDM space are extensions of the respective other spaces. In the second part of the paper we introduce biorthogonal basis functions for $H(\text{div}, {\Omega})$ continuing previous research for $H^1({\Omega})$ and $H(\text{curl}, {\Omega})$. These functions can be expressed in closed form as sums of tensor products of Jacobi polynomials, which allows for fast computation of the $L^2$-projection.

We extend the phi-FD finite-difference scheme, which handles domains described by a level-set function on a Cartesian grid, to Neumann boundary conditions. The condition is relaxed near the boundary through the level-set and an auxiliary field, which is eliminated by writing the relaxed condition at a boundary node and at its nearest interior node. The scheme uses only interior nodes and the level-set: no ghost values, flux reconstruction or cut-cell volumes. We prove optimal consistency and show second-order convergence in the L2, H1 and Linf norms on a non-convex domain. Two elementary precautions, an individual normalization of the boundary equations and a rule for the choice of the coupled interior node, make the condition number O(h^-2) and independent of the smallest cut cell, without any stabilization term.

The precise and efficient evaluation of large-scale lattice sums involving power-law kernels is a fundamental computational problem in the simulation of classical and quantum systems with long-range interactions. While methods for spatially isotropic kernels, some based on Epstein zeta functions, have advanced considerably in recent years, the anisotropic case has lagged behind, despite its broad relevance to both fundamental and effective interactions such as the dipole interaction in magnetic materials. In this work, we solve this issue by defining and analyzing anisotropic Epstein zeta functions for which we derive stably computable representations obtained from wave vector derivatives of lattice sums over isotropic interaction kernels. These functions find direct application in the analytical and numerical study of anisotropically interacting lattice systems. Going further, they provide the correction term in an exact equivalence between discrete lattices and their continuous analogs in a recent generalization of the classical Euler-Maclaurin summation formula to lattices and summands involving power-law kernels. Their connection to high-order derivatives of zeta functions can be used to improve convergence rates of numerical algorithms, for instance in micromagnetics, or to provide rapidly convergent expansions suitable for precomputations of generalized zeta functions. We derive a stably computable representation of anisotropic Epstein zeta functions, including the possibility for analytically removing Rayleigh–Wood singularities, and we develop a numerical algorithm for their stable evaluation for any lattice, power-law decay exponent and anisotropy order. We benchmark the algorithm against closed-form identities, direct summation and multi-precision results, obtaining machine precision across various lattices, power-law exponents, and anisotropy orders.

We study the nonparametric statistical inverse problem of recovering the space-dependent permittivity in a coupled Fokker-Planck-Darcy system from discrete, noisy observations of the Fokker-Planck solution. We consider the coupled parabolic-elliptic system on a bounded domain with the physically natural no-flux boundary condition for the Fokker-Planck equation and homogeneous Dirichlet boundary conditions for the Darcy equation. The inverse problem is indirect, since the unknown coefficient enters only through the elliptic equation and is observed only through its effect on the density in the parabolic equation. We first develop the analytical theory of the forward problem required for the statistical analysis, establishing well-posedness and a priori bounds uniform over the admissible set of permittivities. We then derive two stability estimates for the inverse problem, namely a Lipschitz-type forward estimate and a generalized backward estimate. Placing a rescaled Gaussian process prior on the log-shifted permittivity, we show that the posterior contracts around the truth at an explicit polynomial rate in the number of observations, with the posterior mean converging at the same rate. Numerical experiments in one and two dimensions, using preconditioned Crank-Nicolson and ensemble Kalman filter algorithms, complement our theoretical results.

We revisit the derivation of discrete Sobolev and related inequalities on bounded domains for piecewise constant functions on polyhedral meshes. Our results improve the dependence of the Sobolev constant on the Lebesgue exponents and allow us to recover the asymptotics of the continuous case. As a consequence we derive a discrete version of the Trudinger inequality and variants. We extend our results to broken Sobolev spaces and discuss implications of these new or refined discrete functional inequalities in the numerical analysis of non-conforming numerical methods.

We study the stability and consistency of Lloyd's algorithm used in combination with the Lagrangian transport of a density on a Voronoï tessellation. We show that relaxation rates as strong as $O(h^{-1/2})$, where $h$ is the mesh size, give convergence to the continuity equation in Wasserstein distance at the rate $O(h^{1/4})$. Thus the mesh can be kept regular by the correction alone, without the remeshing that Lagrangian methods usually require. An application is proposed for the compressible Euler equations.

Embedded software aims to provide the main characteristics of a cell in real time to manage and optimize battery usage. A complete multidimensional model is not tractable with the low computational resources available in-situ. A class of 0D-in-space models is derived from the P2D Doyle-Fuller-Newman system, with additional assumptions about the functions' shape in space, and strongly reduces the complexity while preserving the main characteristics over time. We propose a conservative, positivity-preserving model built from the P2D system based on the statement that the Butler-Volmer source term is constant in each subdomain (anode, separator, cathode). We develop a mathematically rigorous construction that provides the simplified model and design a scheme to solve the system. We highlight that the model guaranties the conservation property together with the positivity of the concentrations. In particular, the diffusion-limited current density that represents the maximum current that can be used while still preserving the positivity of the cation concentration is obtained explicitly.

The key-value (KV) cache of autoregressive transformers can be viewed as a fourth-order tensor spanning attention heads, tokens, features, and grouped layers. We measure the singular-value spectra of all four mode unfoldings on Mistral-7B-v0.3 and LLaMA-2-13B and compare four standard tensor decompositions: Tucker, CP, tensor train, and t-SVD, at matched storage. The spectra partition the four axes into two classes. The token and feature modes carry low-rank structure, particularly for keys. The head and layer modes are nearly full-rank and resist compression at any practical error level. Among the four decompositions, Tucker achieves the lowest reconstruction error at every compression ratio from $2\times$ to $5\times$, because it can leave the full-rank modes untouched. Comparisons with two-dimensional unfolding baselines show that the preferred representation differs between keys and values: 2D methods achieve lower key error, while four-way Tucker achieves lower value error at matched storage. A mode-pinning theorem certifies the full-rank preservation from the measured spectra alone. Two further spectral properties affect the compressible modes without touching the full-rank ones: values reach a higher error floor than keys at every ratio, and post-RoPE keys lose $41\%$ - $64\%$ of their pre-RoPE compressibility on both models.

In this paper, we are concerned with matrices formed by block-diagonal factors interleaved with fixed permutations – a flexible family of structured matrices. This class has recently drawn interest in deep learning architectures for its balanced expressivity-efficiency trade-off, yet efficient computational strategies for working with it remain to be found. We approach this problem through Riemannian geometry and examine under what conditions this class admits a smooth manifold structure. For the practically important case of orthogonal two-factor matrices, we derive the essential Riemannian tools and propose efficient algorithms for their implementation. The algorithms leverage automatic differentiation, support parameter sharing within each factor, and avoid explicit dense matrix construction. We test them within the Riemannian optimization framework on the best matrix approximation problem and for parameter-efficient fine-tuning of large language models. Beyond the two-factor setting, we study the geometric and matrix-theoretic properties of factorizations with a larger number of block-diagonal factors.

We study a Kirchhoff-type nonlinear integro-differential equation in two spatial dimensions whose coefficients are allowed to depend on time, and we construct conservative discretizations in time for the associated initial–boundary value problem. We consider two symmetric three-layer schemes of Crank–Nicolson type, a locally linear one and a genuinely nonlinear one, and we show that each of them preserves a discrete analogue of the total mechanical energy of the homogeneous problem with constant coefficients. For the nonlinear scheme we establish uniform apriori bounds on the discrete solution and on its discrete velocity by working directly with the discrete energies, without invoking a nonlinear discrete Grönwall inequality, the constants still grow exponentially in the final time, and we prove local second-order convergence in time, both for the solution and for the central-difference approximation of its first time derivative. The nonlinear system arising at each time level is solved by a fixed-point iteration: given iterates at the two preceding levels that satisfy the apriori bounds, that equation has exactly one solution and the iteration converges to it at a geometric rate once the time step is small enough. Since the apriori bounds are index-local, alternating them with that one-step solver constructs the trajectory stepwise, for time-dependent coefficients as well, on the local interval on which those bounds hold. Numerical experiments, carried out in a setting in which the spatial discretization contributes no error, exhibit the conservation of the discrete invariants, confirm the second order in time and verify the geometric convergence of the fixed-point iteration.

The spectral limited-memory preconditioner (sLMP) clusters leading eigenvalues of symmetric positive definite matrices to accelerate conjugate gradient (CG) convergence. In practice, the cluster point is often chosen to be unity. In some cases, however, this choice can fail to accelerate convergence relative to unpreconditioned CG, even when highly accurate spectral information is available. Alternative cluster points have been proposed based on exact-arithmetic convergence analysis, but such analysis does not explain this finite-precision behaviour. We study how the cluster point influences two sources of numerical error in sLMP-preconditioned CG. First, we analyse the propagation of floating-point rounding errors during application of the preconditioner and derive computable relative-error bounds. For the dominant subspace (spanned by the eigenvectors associated with the leading eigenvalues of the unpreconditioned system) and its orthogonal complement (spanned by the remaining eigenvectors), these bounds are minimized by a weighted median and a weighted arithmetic mean of the leading eigenvalues, respectively. Our analysis explains why small cluster points can strongly amplify errors in the dominant subspace. Second, we investigate sensitivity to perturbations in the dominant spectral information when constructing the preconditioner. The resulting perturbation bound is minimized by a weighted median of the perturbed dominant eigenvalues, with weights determined by the eigenvector perturbation magnitudes. Numerical experiments on synthetic problems illustrate the predicted rounding-error and sensitivity behaviour. Together, these results show that cluster-point selection in finite precision should account for exact-arithmetic convergence, rounding errors, and inaccuracies in the available spectral information.

The design of profile extrusion dies remains a challenging task due to the complex rheological behavior of polymer melts and the geometric intricacies of flow channels. Traditional manual optimization approaches, which rely heavily on human experience, are inefficient and often employ unvalidated heuristics. To address these challenges, we present a deterministic and explainable framework for automatic die design based on adjoint-based shape optimization. This approach enables the computation of sensitivities that directly indicate beneficial modifications to the flow channel geometry. A major difficulty in such optimization processes lies in generating boundary-conforming meshes that evolve consistently with changing geometries. To overcome this issue, we employ non-boundary conforming geometry representation methods that eliminate the need for an explicit surface representation along physical boundaries. A dedicated reconstruction technique is developed to recover accurate sensitivity information at the virtual interface between fluid and solid regions. The proposed algorithm is demonstrated on 3D geometries with varying complexity, including realistic extrusion die flow channels. Several objective functionals relevant to industrial applications, such as flow balance at the outflow, are considered. The results highlight significant improvements in performance metrics while maintaining numerical robustness. This work showcases the potential of adjoint-based techniques for automated die design in a domain still largely governed by manual trial-and-error procedures, establishing a foundation for data-efficient, sustainable manufacturing workflows using computational rheology.