Papers

In this paper, we consider the convergence of some popular Kaczmarz-type methods, including the randomized Kaczmarz (RK) method, the randomized extended Kaczmarz (REK) method, the randomized block Kaczmarz (RBK) method, and the randomized double block Kaczmarz (RDBK) method for doubly noisy linear system. To the best of our knowledge, this is the first work that considers and compares the convergence of the Kaczmarz-type methods for consistent linear systems, as well as that for inconsistent linear systems to solve double noisy linear system. We prove that these methods converge to a neighborhood of the smallest norm least squares solution of the underlying noiseless system. Specifically, RK and RBK that are designed for consistent linear systems, often outperform REK and RDBK that are designed for inconsistent linear systems, respectively, in terms of the convergence rate or the convergence horizon. Moreover, we establish the convergence of the more general sketch-and-project framework for the doubly noisy linear system. The analysis provides more comprehensive results without imposing any additional restrictions. Numerical experiments are performed to validate our theoretical results.

We establish sharp Barron regularity for Coulombic many-electron wave functions after extraction of the universal cut-off Jastrow factors. Following the factorization of Fournais et al.~\cite[Definition~1.4]{FournaisEtAl2005}, for a Coulombic eigenfunction $\psi$ we define the successive quotients by \[ \phi=e^{-F_{2,\mathrm{cut}}}\psi \quad\text{and}\quad \phi_3=e^{-F_{3,\mathrm{cut}}}\phi=e^{-(F_{2,\mathrm{cut}}+F_{3,\mathrm{cut}})}\psi. \] Then \[ \phi,\phi_3\in\mathcal{B}^s(\mathbb{R}^{3N}) \qquad\text{for every }s<2. \] This range is optimal among universal factorizations. No factor depending only on the particle number and the nuclear data, but not on the eigenfunction or its eigenvalue, can make every corresponding quotient belong to $\mathcal{B}^2$. We also determine the exact endpoint growth. Writing $\varepsilon=2-s$, we prove that, for either $u=\phi$ or $u=\phi_3$, there is a computable constant $M$ independent of $\varepsilon$ such that \[ \left\|u\right\|_{\mathcal{B}^{2-\varepsilon}}\leq\frac{M}{\varepsilon^2}\left\|u\right\|_{\mathcal{B}^1}. \] For the unperturbed two-electron atom we prove, with a constant independent of $\varepsilon$, \[ \left|\left\|\phi_3\right\|_{\mathcal{B}^{2-\varepsilon}}-\frac{32\pi Z\lvert\phi_3(0,0)\rvert}{\varepsilon^2}\right|\leq\frac{C}{\varepsilon}. \] Hence the quadratic rate in the upper bound is sharp whenever $\lvert\phi_3(0,0)\rvert\neq0$, as is the case for the ground state.

Quasiperiodic parabolic equations (QPEs) describe various physical processes in systems with long-range order without decay, such as heat conduction, diffusion, and transport, and are naturally posed on the whole space. However, practical computations could be performed on finite domains, where preserving the global quasiperiodic structure through appropriate boundary conditions remains a key challenge. Traditional periodic boundary conditions (PBCs) cannot capture quasiperiodic information beyond the computational domain and suffer from Diophantine errors caused by rational approximations of irrational numbers. To address this issue, we propose a class of quasiperiodic boundary conditions (QBCs) based on a homomorphism between the physical domain and a high-dimensional torus, enabling the long-range quasiperiodic structure to be captured at the boundaries. We further establish rigorous convergence results for numerical discretizations of finite-size QPEs with QBCs. Numerical experiments demonstrate that QBCs avoid the influence of Diophantine errors and enable accurate and efficient computations for both high- and low-regularity cases, with improved convergence over PBCs.

In this paper, we present an $hp$-version time-stepping spectral Monte Carlo method for solving semi-linear parabolic equations. The key innovation lies in constructing an exponentially accurate stochastic algorithm that integrates a residual iteration scheme on Gauss-type nodes in both temporal and spatial directions with a reconstruction strategy rooted in spectral methods. To address the long-time simulations and initial singularities that are often challenging for traditional stochastic algorithms (e.g., walk-on-spheres method), we further develop an $hp$-version time-stepping framework that employs multiple time steps and, respectively, geometric time partitions with linearly increasing polynomial degrees to handle these difficulties. Notably, the proposed algorithm bypasses the need to solve linear systems required by traditional spectral methods and remarkably supports parallel computation at both temporal and spatial grid points. We rigorously establish exponential convergence rates for the multistep method within a finite number of iterations. Extensive numerical experiments are conducted to demonstrate the spectral accuracy and computational efficiency of the proposed method in long-time simulations, problems with initial singularities, and a five-dimensional problem, thereby validating the theoretical results.

In this paper, we derive optimal order a posteriori error estimates for the variable step-size explicit two-step exponential Adams (E-Adams2) method for parabolic problems. We begin by introducing an E-Adams2 approximation, deffned by the piecewise linear approximate solutions, which leads to suboptimal error estimates. To recover optimal order error estimates, we introducean appropriate reconstruction of the approximation with the second order residual for the explicit E-Adams2 method, which plays key roles in deriving optimal order a posteriori error estimates for the proposed explicit method for linear and semilinear parabolic equations. Various numerical experiments are carried out to verify the correct convergence rates of the a posteriori quantities, and the high efffciency of the adaptive algorithm.

Across implicit-neural-representation (INR) architectures and analytic benchmarks we observe that a thoroughly tuned Adam (especially its learning rate (lr), e.g. in a hyperparameter sweep from $lr = 0.05$ to $10^{-8}$) can potentially reach a very low loss even on ill-conditioned loss landscape or converge at a plateau far above the loss attained by second-order methods. This report defines the measured metrics that help determine if Adam can mitigate the ill-conditioning on a given loss landscape. We provide the indicators by which each outcome is determined, that are: the condition number of the Hessian and of the Adam-preconditioned Hessian $D^{-1/2}HD^{-1/2}$ (with the derivation from Adam's update rule), the diagonal mass $\rho$ that distinguishes axis-aligned from cross-coupled ill-conditioning, the negative spectral mass estimated by stochastic Lanczos quadrature, and the gradient energy fractions over curvature bands, including the flat fraction that indicates the Adam stall. A worked out $2\times 2$ example and an illustration show the reasons why a diagonal preconditioning by Adam can remove axis-aligned ill-conditioning by rescaling and why it cannot do the same if the ill-conditioning is cross coupled. In addition, we present a case study of FINER image fitting architecture that goes over the whole loss landscape analysis framework: the fitting architecture description, reasons due to which its landscape stalls Adam at saddles, the measured PSNR values through our tuned baselines to the $120$–$134$\,dB results of the blockwise second order methods, the error maps behind those numbers, and description of the benefits such image fitting accuracy gives in practice.

In this work, a new two-grid method presented for the elliptic partial differential equations is generalized to the time-dependent linear parabolic partial differential equations. The new two-grid waveform relaxation method uses the numerical method of lines, replacing any spatial derivative by a discrete formula, obtained here by the finite element method. A convergence analysis in terms of the spectral radius of the corresponding two-grid waveform relaxation operator is also developed. Moreover, the efficiency of the presented method and its analysis are tested, applying the two-dimensional heat equation.

In this work, a semialgebraic mode analysis (SAMA) is proposed for investigating the convergence of a multigrid waveform relaxation method applied to the Finite Element (FE) discretization of the heat equation in two and three dimensions. This analysis for finite element methods is more involved and more general than that for Finite Difference (FD) discretizations, since mass matrix must be considered. The proposed analysis results in a very useful tool to study the behaviour of the multigrid waveform relaxation method depending on the parameters of the problem.

Overparameterized networks often have ill-conditioned, saddle rich loss landscapes. They are usually trained with first-order methods that struggle on those landscapes, while many second-order methods stay infeasible at scale. Some scalable cubic methods replace the Hessian by a computationally cheap approximation and drop the exact block curvature. This paper introduces ARC-block, a blockwise adaptive regularization with cubics (ARC) optimizer. For each parameter tensor it minimizes a cubic model over the true block Hessian, with a per-block adaptive cubic constant and a monotone guard on the full loss. Large tensors are handled matrix-free in a Lanczos-built Krylov subspace. ARC-block's CubicKrylov subsolver's descent is monotone with a per-block cubic decrease bound. For exact steps, and for Krylov steps that meet the Cartis–Gould–Toint residual test, we provide an $\mathcal{O}(\sqrt{B}\varepsilon^{-3/2})$ sweep bound above a cross-block coupling floor. Another subsolver, ARC-$\varphi_1$, uses a trust-controlled exponential relaxation step. Among the cubic Newton methods evaluated here, these are the only ones whose steps stay exact on every parameter block (within the Krylov subspace), including on a $91.4$M-parameter implicit neural representation (INR). Run to convergence on FINER 2D image-fitting INR network, ARC-$\varphi_1$ reaches $133.5$ dB peak signal-to-noise ratio (PSNR), with no post peak collapse. Tuned Adam peaks at $78.2$ dB at $2{,}629$ s (step ${\sim}625$k, about $44$ minutes), further steps degrade its PSNR. On the same benchmark, SOAP, L-BFGS, and the scalable diagonal cubic method AdaCubic stay at or below Adam's $78.2$ dB peak. Over the benchmark's 20 images, run to a plateau rule, ARC-$\varphi_1$'s per-image PSNR has median $132.4$ and mean $131.5\pm3.3$ dB against tuned Adam's $72.8\pm4.5$ dB.

This paper investigate the conforming virtual element method for the optimal control problem governed by the linear poroelasticity equation. Both distributed control and Neumann boundary controls are considered with the optimize-then-discretize approach. We establish the well-posedness of the discrete state and adjoint problems in three-field formulation by deriving the first-order optimality condition. The optimal a priori error estimates are developed through a novel poroelastic projection for both control problems, and these estimates are uniform with respect to the relevant physical parameters of the model problem. Moreover, the numerical experiments are carried out using the primal dual active set strategy, and demonstrate the support for the derived theoretical results.

Newton's iteration is a fundamental tool for root-finding and numerical solutions of systems of equations. The iteration rapidly refines the initial approximation to the exact root, and in general the convergence is quadratic. Since the method requires finding the function value and its derivative at each iteration, in some cases, it may not converge. This is because the value of the derivative gets close to zero. In this paper, we introduce a modified and stable algorithm of Newton's iteration method that addresses this issue, reduces computational cost, and improves efficiency. In addition. we analyze the convergence properties of the modified method to demonstrate its effectiveness.

We present a data-driven reduced-order modeling (ROM) framework for the zeroth angular moment of the solution to the radiative transfer equation (RTE) rather than the full phase-space solution. Our construction is based on the Peierls integral formulation of the angularly averaged density. For media with isotropic scattering, the density satisfies a closed second-kind Fredholm equation with a globally attenuated, weakly singular kernel. We project this equation directly. For media with anisotropic scattering, we utilize the average-fluctuation decomposition to derive a closed system for a projection-based ROM. Numerical simulations are presented to illustrate the effectiveness of the ROMs we implemented.

Regional climate risk assessment is critical for applications such as infrastructure design, disaster forecasting, and insurance resource allocation. However, estimating regional (i.e., high-spatial-resolution) risk with global climate models (GCMs) remains computationally prohibitive, which has driven the development of downscaling methods for coarse GCM outputs. Downscaling is vital for rare events, since quantifying their extreme properties requires high spatial resolution and very long GCM simulations. These methods non-intrusively increase GCM resolution while correcting statistical biases from unresolved fine-scale processes, thereby improving the accuracy of extreme event statistics with long return periods. A key challenge is preserving dynamical consistency, as freely evolving GCM trajectories are not expected to track the observational dataset used for training the correction operator. This is critical for causal extreme event analyses, where storyline-based risk assessment, i.e., extreme event catalogs, is necessary for effective planning. We address this challenge by introducing Dynamically and Statistically Consistent downscaling (DySCo), a non-intrusive framework yielding high-resolution climate projections consistent with coarse GCM dynamics. DySCo relies on a data-driven reformulation of nudging to create dynamically paired training trajectories without intrusive GCM modifications. Using these paired trajectories, we train a dynamically and statistically consistent, two-stage operator. We evaluate the method by downscaling the Community Earth System Model v2 Large Ensemble (LENS2) in time and space towards historical reanalysis. Results show DySCo achieves superior dynamical consistency with the coarse GCM trajectories, essentially applying a minimal, causal correction to the GCM, preserving top statistical performance comparable to state-of-the-art unsupervised models.

We study the smallest singular value of nonuniform Fourier matrices in two settings: clustered nodes and perturbations of an equispaced grid. By reducing the problem to spectral norm estimates for periodic nonuniform interpolation matrices, we obtain nearly optimal bounds in both cases. For clustered nodes, we derive the first local separation condition in which each required gap depends only on the sizes of the two neighboring clusters. For perturbations with the bound \(1/4\leq L<1/2\), our result confirms the conjecture of Austin and Trefethen on the \(2\)-norm Lebesgue constant up to a logarithmic factor.

In this paper, we propose and analyze a diffusion time-changed susceptible-infected-susceptible (SIS) epidemic model driven by time-changed Brownian motion. We prove that the proposed model admits a unique global positive solution for any initial value in $(0,N)$. The extinction and persistence of the disease are then investigated. To approximate the diffusion time-changed SIS model, we construct a positivity-preserving logarithmic Euler-Maruyama (LEM) method. Assuming that the time-changed is given by the inverse of a standard $\alpha$-stable subordinator with $\alpha\in(0,1)$, we prove that the numerical solution converges strongly to the exact solution with order $\alpha$. Finally, numerical experiments are provided to confirm the predicted convergence rates and illustrate the positivity-preserving property of the proposed method.

We introduce weighted Laplacian flow (WLF), a deterministic particle-flow framework for sampling from a target distribution known only up to normalization. The proposed method evolves the logarithmic density ratio through a transport equation and determines the particle velocity from a target-weighted Poisson problem, resulting in a nonlocal and kernel-free mechanism for redistributing mass. For bounded smooth domains, we establish global well-posedness of the flow for Lipschitz initial data and prove that the flow contracts the discrepancy between the evolving and target densities at an exact exponential rate in terms of the oscillation of the log-density ratio, as well as in the associated projective metric. The relative entropy also satisfies a precise dissipation relation involving the two directional Kullback–Leibler divergences. In addition, we identify a variational interpretation of WLF. The method can be viewed as a gradient flow of the reverse Kullback–Leibler divergence under a target-anchored metric. Consequently, the method achieves an explicit relaxation scale without requiring log-concavity or spectral-gap assumptions on the target distribution. Numerical experiments on multimodal, heavy-tailed, and high-dimensional targets demonstrate the effectiveness of the proposed approach in capturing long-range mass transport.

Data-driven methods for earthquake and tsunami early warning rely on large ensembles of rupture scenarios and their resulting waveforms, but generating such datasets with repeated high-fidelity seismic and tsunami simulations is prohibitively expensive. We exploit the linear time-invariant structure of both dynamics to precompute elastic Green's functions and acoustic-gravity adjoint responses, reducing the source-to-waveform map to two consecutive convolution operators. We evaluate these convolutions with a distributed, FFT-accelerated GPU pipeline that partitions the large seafloor grid across GPUs and directly generates the final observation waveforms. We demonstrate the scalability of this pipeline for the Cascadia Subduction Zone with 963 subfaults, 2,416,530 seafloor grid points, 64 observation locations, and 256 timesteps, requiring 9.45 TiB of aggregate GPU memory. On 64 GB200 GPUs within one NVL72 domain, the pipeline generates waveforms in 24 ms per rupture once the response operators are resident, enabling large rupture ensembles to be evaluated within minutes.

We introduce a continuous data assimilation method for particle-in-cell simulations of the Vlasov-Poisson equation when only hydrodynamic moments are observed. The forecast state is an empirical measure on phase space, whereas the observed fields (density, bulk velocity, and temperature) constrain only a few velocity moments and leave the velocity-space shape of the distribution undetermined. We construct the moment feedback as a Wasserstein gradient flow of a moment-mismatch functional over phase-space measures. The resulting drift acts directly on particle positions and velocities, couples the density, momentum, and energy residuals through a single variational structure, and vanishes on the entire moment-compatible set. Under the standard Wasserstein metric, the energy residual produces a position correction that grows quadratically with the particle speed, and the particle system falls outside standard well-posedness theory. Our primary formulation pairs the quadratic moment mismatch with a velocity-weighted Wasserstein metric that penalizes spatial transport at large peculiar velocity relative to the observed bulk flow, which removes this growth. A direction-split variant retains the plain metric instead. Under the same weighted metric, an alternative moment-relative-entropy functional yields an affine, shape-preserving velocity correction and explicit global moment balances for the space-inhomogeneous system. We prove that the finite-particle scheme with linear Lenard-Bernstein collisions is globally well posed. In 1D1V and 2D2V experiments with several collision models, the nudged formulations reduce bulk-velocity and temperature errors by up to two orders of magnitude relative to an unassimilated run.

The randomized extended Kaczmarz method is an effective iterative framework for solving large-scale inconsistent linear systems. In this paper, we extend this framework to third-order inconsistent tensor linear systems under the t-product and propose the Tensor Greedy Double Block Extended Kaczmarz (TGDBEK) method. At each iteration, TGDBEK dynamically constructs active blocks of row and column slices via a residual-based greedy selection strategy, prioritizing the slices associated with the largest residual norms. Unlike existing tensor block Kaczmarz variants that rely on static, predefined partitions –the tensor randomized extended block Kaczmarz (TREBK) method, its greedy counterpart (TREGBK), and the tensor randomized extended average block Kaczmarz (TREABK) method – TGDBEK adapts the active block sizes dynamically at each step using a single intuitive threshold parameter $\eta$. We establish the theoretical linear convergence of TGDBEK to the unique minimum-norm least-squares solution $\mathcal{A}^\dagger * \mathcal{B}$. Extensive numerical benchmarks on synthetic dense and sparse tensor systems, as well as multidimensional multichannel color and 3D volumetric MRI image deblurring problems, demonstrate that TGDBEK substantially outperforms state-of-the-art tensor Kaczmarz solvers in both iteration count and CPU running time.

Designing a high-order structure-preserving discretization is a constrained mathematical search: conservation, a positive discrete inner product, physical spectral behavior, boundary accuracy, bandwidth, and partial differential equation (PDE) error must hold simultaneously. We test whether large language models (LLMs) can help while remaining non-authoritative. The motivating MOLE implementation of the Corbino-Castillo staggered operators satisfies a general discrete Gauss identity and conserves exactly, yet its order-six and order-eight Dirichlet blocks develop four non-real boundary-localized modes. An endpoint-supported positive-definite identity would instead force a real non-positive spectrum, so the search changes the closure and norm architecture. An LLM proposes only a typed construction program; a deterministic linear-program compiler generates coefficients; an independent verifier tests algebra, positivity, physical modes, conditioning, and manufactured PDEs; and coupled rational reconstruction provides exact certificates. Across 1,200 solver evaluations, an externally fixed verifier accepted 55.0% of full-metric-feedback proposals and 53.3% of illumination-archive proposals, versus 13.3% for uniform random search. Four leading LLM-originated programs were reconstructed exactly. The strongest order-six-interior, order-four-boundary candidate lowers the prior positive-diagonal spectral-radius constant by 25.4% and its PDE error 62.7-fold. Its certified heat-equation energy is contractive, whereas the order-six reference exhibits 6.6% transient growth; both have the same RK4 stability limit. The LLM proposes structural hypotheses; deterministic mathematics determines validity.

We introduce a large language model (LLM)-driven approach for generating dynamic algorithmic dispatch heuristics in high-performance linear algebra. By combining prompt engineering with LLaMA 3 and a curated performance database, the model learns to synthesize selection heuristics that exploit structural patterns to identify fast algorithmic choices. A case study on LU factorization demonstrates the model's ability to replicate expert-designed strategies. This work, developed as part of the DARPA-MIT SmartSolve project, highlights the promise of LLMs for algorithmic discovery and the development of more adaptive, fast linear algebra software.

The classical procedure for solving hydraulic networks relies on the assumption of constitutive equations, which state the relationship between the pressure gradient and fluxes along an edge. In this paper, we propose a data-driven framework that bypasses these constitutive models, formulating the incompressible flow problem directly on the graph topology. By assigning discrete measured data points to network edges, the problem is cast as a mixed-integer quadratic optimization over nodal pressures, edgewise states, and data assignments, accommodating both laminar (convex) and turbulent (non-convex) regimes. To solve this, we evaluate three algorithms: a GPU-accelerated Brute Force method, the Alternating Direction Method (ADM), and Deterministic Annealing (DA). The Brute Force method certifies global optima for small networks, establishing a good baseline for the iterative solvers. We demonstrate that ADM is highly sensitive to its initialization, requiring a faithful surrogate model to avoid local minima. In contrast, DA eliminates this dependence through unsupervised clustering. By annealing the data assignment from the centroid to strict nearest-neighbor projections, DA consistently reaches the global optimum without prior manifold reconstruction. Furthermore, numerical experiments reveal that DA is robust to noisy data, remains thermodynamically admissible on all but the coarsest and noisiest datasets, and sustains its convergence rate on larger networks where Brute Force is intractable and ADM degrades. Finally, the framework is successfully validated on complex configurations, including mixed-component networks and a $958$-edge arteriovenous bed featuring a non-Newtonian Carreau–Yasuda model, demonstrating its scalability and practical applicability.

The primary challenge of conducting PDE-constrained optimization for high-dimensional problems, such as kinetic equations, is the often prohibitive memory cost. Computing gradients using the adjoint state method would require the storage of the entire time history of the forward solution. For such problems, where the memory cost for storing a single instance of the forward solution can already be a limiting factor, this is clearly not feasible. In this paper, we propose a memory-efficient adjoint state method that compresses the forward and adjoint solution with a dynamical low-rank approximation (a model order reduction technique) and bypasses the need to store the entire forward solution by employing a time-reversible low-rank integrator. The dynamical low-rank approach introduces a number of challenges: reversibility can fail in the rank-deficient case and the low-rank trajectories can show chaotic behavior. In particular, the latter has a number of important consequences for the optimization problem. We address those challenges and show that our method can drastically reduce the memory requirement for gradient-based optimization of kinetic equations. In particular, we consider two examples from kinetic plasma physics: optimizing beam profiles to suppress a bump-on-tail instability and shaping a particle beam using external electric fields.

The retrieval dynamics of a modern Hopfield network is the gradient flow of a log-sum-exp energy, while the attention update is its exact difference-of-convex minimization step. We study which time discretizations preserve not only energy decay and equilibria but also basins of attraction. We introduce energy cells, connected components of sublevel sets containing one attractor and no other critical point. Our main theorem shows that every finite energy cell below the escape energy is contained simultaneously in the basin of the continuous flow, every relaxed attention map $\Psi_\theta=(1-\theta)\,\mathrm{id}+\theta\,\mathrm{attention}$ for $0<\theta<2$, and implicit Euler throughout its uniqueness regime. A parameter-uniform unit-curvature majorant yields unconditional dissipation and a monotone interpolation of each discrete step. We also derive explicit local contraction bounds near well-separated patterns, with a certified optimal slight overrelaxation; characterize proximal tunneling and overshoot beyond the preservation regimes; compare first-order error constants; establish an order barrier for scalar reparametrizations of the relaxed family; construct a second-order scalar-auxiliary-variable scheme; and extend cell preservation to damped difference-of-convex iterations in Bregman geometry, including a certified overrelaxed window under bounded asymmetry. Nine numerical campaigns test the bounds and failure mechanisms. In two-dimensional basin experiments, all observed disagreements between continuous and discrete retrieval occur above the attractor-specific numerically inferred escape level.

We give a self-contained finite-precision analysis of the symmetric Lanczos algorithm without reorthogonalization. In particular, we derive the perturbed three-term recurrence, Paige's loss-of-orthogonality identity, containment of all computed Ritz values, and localization of stabilized Ritz values. We then prove a Greenbaum-type backward stability result, exhibiting a nearby problem on which exact Lanczos produces the computed tridiagonal matrix. Our proofs simplify those of Paige and Greenbaum, at the cost of hiding polynomial factors in the iteration count.

We present a comprehensive framework for tangent vector fields and differential operators on the rotation group $\mathrm{SO}(3)$ using harmonic series expansions, which provides a mathematical foundation for their implementation in the crystallographic texture analysis software MTEX. The central idea is to employ the standard left- and right-invariant frames as global orthonormal frames of the tangent bundle, thereby avoiding the numerical instabilities of the classical tangent space basis derived from the Jacobian of the Euler angle parametrization. While the choice of frames is primarily motivated by their geometric and numerical properties, their full potential emerges in the harmonic setting. Representing tangent vector fields through harmonic expansions of their frame components, we derive explicit frequency domain formulas for the gradient, divergence, and curl, allowing these operators to be applied directly to the harmonic coefficients. Moreover, we show that these differential operators preserve harmonic band-limitedness. In addition, the left- and right-invariant representations of tangent vector fields can be transformed into one another directly in the frequency domain, with an increase in harmonic bandwidth of at most one degree.

These lecture notes form the first part of a master's-level course on advanced numerical linear algebra. Their aim is not only to present the classical algorithms, but to show why the subject has become considerably more central than it was a generation ago. Numerical linear algebra grew up alongside the numerical solution of partial differential equations, and for a long time that is where its large sparse systems came from. Ranking the nodes of a network, assimilating observations into a weather forecast, and fitting a model to a large noisy data set now lead to problems of the same kind: too large to factorise, structured, and accessible only through matrix-vector products. Strikingly few ideas are needed for all of them. Each chapter therefore develops a standard topic and then puts it to work outside its original setting. We treat norms, factorisations, conditioning and floating-point arithmetic; sparse matrices arising from finite differences, from graphs and from machine learning; stationary iterations and the smoothing property; the conjugate gradient and Lanczos methods, with spectral clustering and regularisation by early stopping; Arnoldi and GMRES, with PageRank and large least squares; and finally preconditioning, Schwarz domain decomposition and multigrid. We assume a first course in linear algebra. Every section closes with a summary of what should be retained and every chapter with exercises, several drawn from past examinations. Accompanying Python code reproduces the numerical illustrations.

Pairwise reciprocal matrices are fundamental to the Analytic Hierarchy Process (AHP),a decision-making model. While the Direct Least Squares (DLS) method provides an intuitive mechanism for deriving priority vectors without complex transformations, it is susceptible to solution non-uniqueness. Under high levels of inconsistency, such as severe cyclic contradictions, the DLS optimization landscape becomes non-convex, yielding multiple distinct global minima. Consequently, priority rankings become unstable and critically dependent on initial algorithmic guesses. Furthermore, established prioritization operators (POs), including normalization techniques, the Eigenvector method, Singular Value Decomposition, Cosine Maximization, and the Pseudo-Inverse Gram Matrix (the closed-form solution of Weighted Least Squares), frequently generate disparate outcomes. To overcome these structural deficiencies, this paper introduces the Anchored Regularized Direct Least Squares (ARDLS) optimization model as a harmonizing framework. ARDLS integrates uniquely determined established POs as theoretical anchors within a regularization penalty. This integration systematically breaks mathematical symmetries and tilts the optimization landscape to guarantee convergence upon a single, unique global minimum. By minimizing the root mean square variance (RMSV) of the initial baseline vectors, ARDLS effectively unifies these divergent solutions. Comprehensive numerical experiments validate that the framework successfully fine-tune the solution of established POs by reducing RMSV while ensuring strict mathematical uniqueness. The practical utility of the method is further demonstrated through a numerical case study resolving an innovation fund dilemma in FinTech project selection. The proposed ARDLS approach offers a robust alternative to classical AHP across a wide range of decision-making domains.

This paper studies the numerical solution of elastic wave propagation on networks, modeled by elastodynamic equations posed on each edge, coupled at the nodes through suitable transmission conditions. We propose and analyze a hybridizable discontinuous Galerkin method that exploits the network structure to reduce the global problem at each time step to a linear system whose size depends only on the number of network nodes and not on the polynomial degree of the discretization. Combining it with an energy-conservative implicit time discretization, we derive a priori error estimates of optimal order in space and time. The implicit time discretization avoids the severe CFL restriction caused by the large variation in fiber segment lengths. To efficiently solve the resulting, typically ill-conditioned global system, we introduce a two-level overlapping additive Schwarz preconditioner. Under suitable assumptions on the network, we establish uniform convergence of the resulting preconditioned conjugate gradient method. Numerical experiments confirm the theoretical findings.

It has been conjectured that the operator norm of the Lyapunov operator induced by the Frobenius norm is always attained at a symmetric matrix. The conjecture is known to hold for all matrices of order at most five. We give an integer matrix of order seven for which the skew-symmetric restricted norm is strictly larger than the symmetric restricted norm. A rational separator and exact-arithmetic certificates establish the strict inequality without relying on floating-point computations. A direct-sum construction yields counterexamples in every order $n \geq 7$; the case $n = 6$ remains open.