Randomized numerical algorithms have undergone rapid development in recent years, exhibiting significant potential for practical applications and offering novel methodologies for solving large-scale linear systems. This paper reviewed the state of the art of three primary randomized algorithms for sparse linear systems, examined their features, computational complexity, and key challenges. Based on this analysis, we evaluates the current research progress, identifies gaps between current capabilities and practical application requirements, and proposes future research directions for randomized algorithms targeting large-scale practical problems.

数值线性代数解法器是科学工程计算与工业软件的核心组件，也是影响这些软件计算效率的主要瓶颈，其高效算法设计与性能优化面临复杂应用特征与超级计算机体系结构特征的双重挑战，一直以来都是学术界和工业界广受关注的问题。自2018年开始，国内相关专家发起并组织了解法器快速算法及应用研讨会（Solver会议），该系列会议至今已成功举办了8届，已成为国内该领域研究人员发布研究成果的交流平台，吸引了工业界和实际应用部门的广泛关注。本系列专辑拟邀请Solver会议组织者担任客座编委，不定期组织活跃于该会议的专家和团队，围绕数值代数解法器的快速算法设计、性能优化、自主软件研发和实际应用撰写文章，展示该领域的最新进展，促进该领域在我国的发展。

本专辑由北京应用物理与计算数学研究所徐小文研究员等人（编委名单见后）负责组织，经过严格同行评审，最终录用了8篇文章，涵盖了自主解法器软件与算法库、应用驱动快速算法、面向国产处理器的性能优化、新型算法等主题，一定程度上反映了我国近年来该领域的研究进展。

Large-scale algebraic eigenvalue problems are widely used in fields such as materials science and engineering structure analysis. In multi-scale complex systems, with the problem size increasing and the accumulation of errors induced by scale disparities, traditional eigenvalue algorithms face significant challenges in accuracy, stability and computational efficiency, and in some cases even fail to achieve effective solutions. This paper presents recent progress in the development of the GCGE eigenvalue solver, focusing on the design and implementation of efficient parallel algorithms aimed at enhancing its stability and adaptability in complex systems. To address the performance and stability bottlenecks of conventional solvers when dealing with complex small-scale problems and large-scale ill-conditioned problems, we introduce several improvements to the GCG algorithm. These include a diagonal normalization preconditioning strategy to improve convergence and numerical robustness, and orthogonalization strategies based on L² inner products and stiffness matrix inner products to enhance numerical performance for cases involving indefinite mass matrices. Furthermore, GCGE has been integrated into the SLEPc software framework to support the solution of Hermitian eigenvalue problems and dominant eigenvalue problems, and has been modularly deployed on the Baltam Platform, demonstrating strong potential for engineering applications.

The effect of an oscillating external electric field on ion transport in nanoscale channels is investigated within the framework of the classical Poisson-Nernst-Planck (PNP) model. Three cases are analyzed by multiscale method: (1) an externally applied time-oscillating electric field along the channel direction; (2) a time-oscillating electric field coupled with a spatially periodic electric field along the channel; (3) an oscillating electric field applied along the channel direction in the presence of periodic surface charge distributions on the channel walls. An effective model for high oscillation frequencies is derived by considering the leading order approximation, which shows that the ion distribution and average transport properties within the channel depend only on the time and space averaged properties of the external oscillating electric field. Numerical simulations on a two-dimensional single-channel nanoscale model confirm the validity of these analytical results.

In this paper, firstly, inspired by the ideas of two algorithms: the greedy randomized Kaczmarz algorithm and the geometric probability randomized Kaczmarz algorithm for solving linear equations, two novel greedy randomized algorithms for solving matrix equations are proposed. Then, the convergence of these two methods in matrix equations is proved based on important inequalities. Finally, the numerical experiments are implemented to verify the effectiveness of the proposed methods. The numerical results show that the geometric probability randomized Kaczmarz algorithm outperforms the greedy randomized Kaczmarz algorithm for large-scale systems.

The block algebraic multigrid (BAMG) method is an efficient preconditioning technique for solving discretized systems of partial differential equations and can be applied to multiphase flow models in porous media. In this work, we investigate the impact of different norm selections in block-matrix coarsening and a diagonal-block-based classical interpolation strategy on the convergence behavior of the BAMG method for multiphase, multicomponent reservoir simulations. To address the parallel bottleneck during the setup phase, we introduce the “delayed update” and “structure-preserving update” strategies, and develop an adaptive-setup-based BAMG algorithm guided by a tolerable iteration growth threshold criterion. Furthermore, we perform in-depth optimization and vectorization of performancecritical components in the solve phase, including the smoother and residual computation. Numerical experiments demonstrate the proposed method’s advantages in terms of convergence, efficiency, and parallel scalability. For example, in tests involving grids with over one hundred million elements, the proposed method reduces computation time by 57.6% and improves parallel efficiency by 38.1% when using 8192 CPU cores compared to the traditional method.

The accurate simulation of Fluid-Structure Interaction (FSI) is important for the safety design and operation of nuclear reactors. This paper investigates the simulation of FSI by using a two-way coupled FSI solver based on the high-order discontinuous Galerkin method, which is implemented by the open-source software ExaDG. The solver is verified through the examination of a cantilever beam vibrating in the quiescent water, where the simulated vibration frequencies and damping ratios are compared with theoretical predictions. Furthermore, a numerical example involving turbulent flow is presented. This example employs the Large Eddy Simulation (LES) method to resolve the large scale eddies and utilizes the Synthetic Eddy Method (SEM) to generate inlet boundary conditions, showcasing the potential for applications in engineering and scientific research.

Algebraic multigrid (AMG) is an efficient method for solving linear equation systems as preconditioners. Semi-structured AMG utilizes structured information for efficient computation and supports the presence of unstructured information, thus achieving both high performance and high flexibility, making it widely used in various scenarios of scientific and engineering computing. However, the current mainstream semi-structured AMG solvers still have significant deficiencies in absolute speed and scalability. Therefore, we developed SemiStructMG. On the one hand, it utilizes multidimensional coarsening to reduce complexity, improving single step running speed and scalability; On the other hand, it considers interblock connections in the smoother and interpolation operators, improving convergence in various complex problems. We tested Semi-StructMG in benchmark tests and multiple realworld applications, and achieved speedup of 5.97x, 15.2x, and 3.85x compared to SSAMG, Split and BoomerAMG in hypre.

The X-Solver library aims to implement and optimize Krylov subspace solution methods and preconditioners to solve large and sparse linear systems of equations on distributed memory clusters equipped with many-core GPUs as accelerators. Motivated by the requirements from real-world applications and the trend of hardware integrations, we manage three achievements, i.e., more effective support for distributed calculations, better exploitation of characteristics of modern many-core architectures, and portability across different heterogeneous platforms. Numerical experiments demonstrate advantages over other state-of-the-art libraries in a broad range of applications on the targeting hardware.

Radiation diffusion problems are widespread in multi-physics coupling fields such as astrophysics and inertial confinement fusion. Algebraic Multigrid (AMG) methods based on the physical and algebraic features of a problem have become a hot topic in the field of multigrid research. This paper proposes T2T2-ILU(0)-V-FGMRES solver for efficiently solving linearized discrete systems derived from three-temperature radiation diffusion equations. This solver is a PGMRES solver employing a common unsmoothed aggregation AMG (UA-AMG) preconditioner. Furthermore, several physical and algebraic features are extracted from different discrete systems to enhance the computational performance. Based on these features, we develop an adaptive UA-AMG preconditioned FGMRES solver, termed Adapt-UA-AMG-FGMRES. Numerical experiments show that the new solver exhibits better robustness and computational efficiency. Compared to the T2T2-ILU(0)-V-FGMRES and HMIS-V-FGMRES solvers (best common non-aggregation AMG preconditioners), the CPU time of the new solver is reduced by approximately 49.1% and 25.3%, respectively. It should be noted that the algorithmic design principles can be easily extended to more general model problems, such as multi-group radiative diffusion equations.

The paper integrates low-rank quaternion matrix recovery with the Deep Image Prior method to propose a novel color image inpainting model, LRQMD(Low-Rank Quaternion Matrix Completion with Deep Image Prior). A corresponding numerical algorithm is developed based on the Alternating Direction Method of Multipliers. Furthermore, preliminary analysis of the model's reconstruction capability is conducted using the theories of exact low-rank matrix completion and neural tangent kernels, and the convergence and optimality conditions of the proposed algorithm are established. In terms of network architecture, Inverse Evolution Layers, constructed based on the heat diffusion equation, are introduced. Experiments demonstrate that combining IELs with the LRQMD model improves denoising performance. Finally, numerical experiments compare LRQMD with traditional methods such as LRQMC and QDIP, confirming its superiority in color image restoration tasks.

In large-scale linear system solving, traditional sparse direct solvers typically rely on a single precision computation, making it difficult to flexibly balance computational efficiency and numerical accuracy. To address this issue, a mixed-precision optimization algorithm based on the distributed sparse direct solver PanguLU is proposed, tailored for the heterogeneous multi-zone processor MT-3000. The algorithm dynamically selects storage precision for matrix blocks based on their spatial location and numerical sensitivity, enabling mixedprecision computation during numerical factorization. Meanwhile, a pipelining mechanism that decouples computation precision from storage precision for GEMM subtasks is designed. Experimental results show that the proposed method achieves a speedup between 1.04 and 1.19 times in the numerical factorization phase, while reducing the relative residual by a factor between 1.97 and 4.15 compared to baselines of a single precision, thereby effectively controlling accuracy loss while improving computational performance.

The three-temperature radiation diffusion equations provide an accurate mathematical framework for describing the physical phenomenon where radiation energy propagates through a medium, undergoing scattering, absorption, and emission processes. Based on an overlapping decomposition strategy on physical quantities, we proposed restricted additive and multiplicative Schwarz preconditioning and iterative algorithms in [Yue X, He J, Xu X, Shu S, Wang L. Commun. Comput. Phys., 2022, 32: 829-849], but did not provide any convergence analysis or numerical verification. In this work, we propose reasonable approximation assumptions for the sub-matrix inverses within the involved Schur complements. Accordingly, we construct the inexact restricted additive and multiplicative Schwarz preconditioning algorithms. Motivated by the concepts of norm- and field-of-values-equivalences, we conduct the optimal convergence analyses for their preconditioned generalized minimal residual iterative methods under the stability conditions necessary for ensuring the symmetrypreserving finite volume element discretization (or its associated coefficient matrix) to yield a unique solution. Furthermore, we validate the effectiveness of the theoretical results through experimental test cases derived from realistic simulations of hydrodynamic instability during the deceleration phase of a laser-driven spherical implosion.

This paper presents a modulus-based matrix splitting iteration method for solving a class of vertical linear complementarity problems. By reformulating the vertical linear complementarity problem as an equivalent nonlinear system of equations, a new class of modulus-based matrix splitting iteration methods is created, and the convergence of the algorithm is proved under certain conditions. Finally, two numerical examples are provided to demonstrate the effectiveness of the proposed algorithm.

In this paper, a new hybridization strategy for the sixth order unequal-sized weighted essentially non-oscillatory (HUS-WENO) scheme is designed for solving the nonlinear degenerate parabolic equations on structured meshes. This HUS-WENO scheme only uses the information defined on three unequal-sized stencils to obtain the optimal sixth order accuracy in smooth regions while maintaining stable, non-oscillatory and sharp discontinuity transitions. In addition, the linear weights in the proposed scheme are artificially set to be any random positive numbers with the only requirement that their sum equals one. To reduce CPU time, a hybridization strategy is designed based on the reconstruction polynomial of the six-point stencil in US-WENO scheme, which can accurately, efficiently, and automatically identify the troubled cells, and dose not contain any artificial parameters related to the problems. Finally, some numerical examples are used to verify the performance of the presented WENO scheme in various aspects, such as computational efficiency, low dissipation characteristics, shock capture ability, etc.

Based on the classical spectral collocation method for solving boundary value problems of partial differential equations, the integral and differential operators of fractional order Chebyshev polynomials with integer and fractional order are deduced in vector form, and the solution function that satisfy the consolidation equation and its boundary conditions is established. After substituting the function into the consolidation equation, a system of algebraic equations is deduced with the help of the spectral collocation method. Finally, the numerical algorithm for solving the one-dimensional fractional viscoelastic consolidation problem is presented. It is a direct method for solving the consolidation equation in the time domain, and the excess pore water pressure and the effective stress, as well as settlement are the explicit functions of time, those are composed of finite term series, and it is a useful mathematical tool for solving problems of fractional viscoelastic consolidation with variable loading or multi stage loading. Some numerical examples are shown to illustrate the accuracy and effective of the proposed numerical algorithm.

This paper investigates the problem of recovering the geometric parameters of fractal rough surfaces from measured scattered-field data. By integrating a self-attention mechanism with physical parameter properties, we propose a new neural network architecture, termed the Surface Reconstruction Neural Network (SRNN), for surface parameter estimation. By the structural characteristics of SRNN, we prove its numerical convergence. Numerical experiments further demonstrate the effectiveness of the proposed method.

Acoustic point source inversion is a critical yet ill-posed inverse problem in wave physics, proving particularly challenging when observation data is sparse, limited-aperture, and noisy. To overcome the limitations of traditional methods and existing deep learning models in addressing complex multi-source and strong-interference scenarios, this paper proposes a unified end-to-end inversion network based on Multi-Frequency data fusion and the Transformer (Multi-Frequency Field and Count Transformer, MFFC-Former). The model innovatively aggregates the Multi-Frequency complex response from each measurement point into a feature token. It then leverages the Transformer's powerful self-attention mechanism to capture the global dependencies among all measurement points, thereby synchronously performing source count classification and location indicator field regression within a single network. This end-to-end, multi-task learning paradigm discards complex post-processing steps (such as MCMC) and avoids reliance on information such as noise priors, enhancing solution efficiency and system integration. Numerical experiments under complex conditions, involving up to 6 point sources and high noise levels (up to 20\%), demonstrate that MFFC-Former outperforms a fully connected network (MLP) baseline with a comparable number of parameters in both average localization error and count prediction accuracy. Particularly in the challenging scenario with 6 sources and 20\% noise, where the MLP baseline fails to resolve all sources, MFFC-Former successfully resolves and locates all source points, demonstrating its resolution and robustness in multi-source, strong-interference environments. The results of this study demonstrate that leveraging the Transformer architecture to effectively fuse the intrinsic correlations within Multi-Frequency observation data is a viable pathway for solving such highly ill-posed inverse problems.

Inverse problems for the wave equation aims to reconstruct unknown model coefficients such as sources and media through external observation data. It has wide applications in many engineering fields. This paper is concerned with the inverse moving point source problems, which plays an important role in accurately recovering the motion characteristics and flight parameters of a moving object. First, a mathematical model for the moving source problem in the frequency domain is established based on the the Fourier transform. Then, the trajectory of the moving source is reconstructed using least squares and nonlinear optimization algorithms. For smooth motion trajectory functions, the effects of parameters such as frequency-band, observation directions, and the number of series expansion terms on the inversion results are discussed. In the three-dimensional case, inversion algorithm is performed using multi-frequency far-field and near-field observation data contaminated by white noise. Through extensive numerical experiments, it is found that the algorithm exhibits good accuracy and robustness under a certain level of noise interference. Furthermore, numerical experiments are also conducted for cases where the motion trajectory is a non-smooth function and a multi-scale function, yielding satisfactory results.

This paper combines the Tikhonov regularization method with a sequential quadratic programming(SQP) scheme to propose a new iterative regularization approach for nonlinear inverse problems. Incorporating a line-search strategy, we establish the global convergence and its regularization property in the presence of noise. Several numerical experiments are given to demonstrate its effectiveness.

This article deals with asymptotic solution to a fractional-order dynamic system describing interactions of cell-chalone in micro-environment and inversion of the fractional order. The asymptotic solution is derived by the Laplace-ADM method, and convergence of the asymptotic solution to the exact solution is proved. Based on the computable asymptotic solution with an observation of the chalone at a given time, the inverse order problem is transformed to a nonlinear algebraic equation, and its uniqueness is obtained by monotonicity of the nonlinear function of the order. Furthermore, inverse problem of determining the multi-parameters of the model is discussed also based on the asymptotic solution. Numerical inversions with noisy data are presented showing that the average inversion solution approximates to the exact solution as the noise level goes to small.

This paper concerns the inverse scattering problem for a locally rough surface. We propose a nonlinear integral equation method to reconstruct both the shape and location of the locally rough surface by using scattered field or phaseless total field data for single point source. Based on potential theory, we establish the associated field equation and data equation. The density is solved from the field equation, and then the update of boundary is derived by the linearized data equation. This process continues iteratively until the relative error of the scattered field meets a prescribed tolerance. Furthermore, by employing the reflection principle, we prove the injectivity and dense range property of the Fréchet derivative operator, which provides a theoretical guarantee for the solvability of the linearized data equation. Numerical experiments are presented to demonstrate the effectiveness and robustness of the proposed iterative algorithm.

The Random Feature Method (RFM) is an emerging framework for solving partial differential equations (PDEs). It constructs function spaces for approximating solutions through shallow neural networks, combining the rigor of traditional numerical methods with the flexibility of modern machine learning. In this work, we design and implement a high-performance RFM solver, pyRFM, based on the Python programming language. In terms of framework design, pyRFM provides a complete workflow that covers key components including geometry representation and sampling, feature matrix assembly and equation solving, as well as result visualization, thereby enabling a systematic process from modeling to numerical experiments. Meanwhile, the software is developed with careful consideration of the intrinsic characteristics of RFM, leveraging widely used and well-supported Python scientific computing libraries to ensure ease of use and extensibility. Moreover, pyRFM is fully compatible with both CPU and GPU computing environments, and its entire workflow is free from traditional mesh-based operations, which grants it higher flexibility and efficiency when handling complex geometries.

This paper reviews some of the recent developments in energy-conserving or energystable schemes. It mainly focuses on the relaxation-type schemes. The key is to introduce the relaxation idea and governing equation into some schemes, taking into account the characteristics of the model. The relaxation parameter is determined by the governing equation based on energy conservation or dissipation. On one hand, the governning equation ensures the local existence and uniqueness of the relaxation parameter; on the other hand, it guarantees the structure-preserving properties of the algorithm. Moreover, based on the estimation of the relaxation parameter, the consistency of the relaxation algorithm is derived. Finally, some applications of the algorithm to nonlinear stiff ordinary differential equations and some partial differential equations are discussed.

Machine learning-based methods have advanced electronic structure calculations in groundstate, excited-state, and time-dependent multi-electron systems. For ground states, neural network wave functions with Slater-Jastrow-Backflow forms, trained via variational Monte Carlo, accurately capture electron correlation, achieving precision comparable to or exceeding coupled-cluster approaches. In excited-state calculations, techniques such as state-averaged penalties and natural excited-state variational principles enforce orthogonality and enable accurate prediction of excitation energies and oscillator strengths for atoms and molecules. For time-dependent systems, the time-dependent variational Monte Carlo method, which evolves parameterized wave functions, precisely simulates electron dynamics under strong fields and captures non-equilibrium effects. Integrating pseudopotential with neural networks improves computational efficiency while maintaining accuracy in complex systems, including those with transition metals. These developments highlight the strong representational capacity of neural quantum states and their applicability across diverse quantum chemistry problems, offering effective tools for high-accuracy simulations in physical and chemical sciences.

This study improves the source iteration method to solve the steady-state neutron transport equation and analyze neutron behavior. The method discretizes energy using the group method, angle via discrete ordinates, and space with a discontinuous finite element scheme incorporating an upwind flux. The key improvement involves assembling neutron fluxes for all angles and energy groups into a coupled matrix system, retaining interaction terms between different groups and angles. Validation using KUCA experiments and the NEACRP-L-330 benchmark demonstrates that the enhanced method achieves higher accuracy in calculating the effective multiplication factor and average neutron flux while reducing computational time. This provides an efficient approach for solving neutron transport problems.

This paper proposes the L_{2_∞} norm for generating adversarial examples, which achieves an effective balance between the L₂ and L_∞ norms. Leveraging the geometric interpretation of the L_{2_∞} norm, we analyze its advantages in adversarial attack optimization. The optimization method based on this new norm can stably generate adversarial examples. Numerical experiments demonstrate that the proposed norm significantly improves the transferability and visual imperceptibility of adversarial perturbations. A series of experiments validate the superior performance of perturbations obtained under the L_{2_∞} norm.

The decomposition results are unstable when applying nonnegative matrix factorization to the spectral unmixing of hyperspectral remote sensing images. This paper proposes a perturbation-based minimum volume-constrained nonnegative matrix factorization model to address the problem of weak linear correlation in the column vectors of the endmember matrix (base matrix) in unmixing. Simultaneously introducing a multi-layer structure, a deep minimum volume-constrained nonnegative matrix factorization algorithm was designed. The numerical results indicate that introducing disturbance terms effectively avoids the algorithm working with singular accuracy, and the multi-layer decomposition algorithm obtains more accurate endmember spectral information than the single-layer decomposition algorithm. This study provides an effective solution strategy for spectral unmixing problems with similar endmember spectra. Also, it provides a reference for interpretable machine learning methods by designing regularization terms based on specific problems.

We mainly consider the adaptive coupling of the finite element method and natural boundary element method for the Signorini transmission problem. First, the variational formulation of the original problem and its discrete formulation are presented. Then, a priori and a posteriori error estimates for the coupling method are provided in the new framework. Finally, numerical experiments validate the theoretical analysis results.

Principal component analysis (PCA) has been the most classic and widely used dimensionality reduction method, but it may not be able to accurately capture the essential structure of the data which containing noise or outliers. This paper proposes an adaptive weighted robust sparse principal component analysis method (AWRSPCA), which skillfully integrates the properties of nonconvex penalty function and generalized mean, aiming at inducing sparsity of principal components, i.e., making the principal components more focused on a few key variables, thus enhancing their interpretability. Meanwhile, the method effectively reduces the sensitivity of PCA to noise and outliers by generalized mean, ensuring the robustness of the principal component extraction process. An alternating direction multiplier method is developed for computing the solution of AWRSPCA and the convergence of the algorithm is proved. To verify the effectiveness of the proposed method, we compare the method with several mainstream methods on both synthetic and real datasets. The experimental results show that the new algorithm excels in both sparsity and robustness, not only successfully extracting the key information in the data, but also effectively resisting the interference of noise and outliers.

This paper presents a numerical method based on the charge simulation method for calculating conformal mappings from the bounded multiply connected regions with a rectilinear slit onto the three category canonical slit domains. Firstly, the rectilinear slit is expanded by using the pre-map function, and the bounded multiply connected regions with a rectilinear slit is mapped onto the regular slit domains of straight slit, spiral slit and a disc with spiral domains. On this basis, the corresponding set of constraint equations is established by using the Dirichlet boundary conditions, and for the pathological matrix in the set of constraint equations, it is proposed to solve the charges by using the iterative method of stabilized bi-conjugate gradient based on LU decomposition (LU-BiCGSTAB), which improves the accuracy of the conformal mapping function. Secondly, the conformal mapping computational method proposed in this paper is applied to the simulation of the bypassing of a spiral point vortex around a bounded region with linear obstacles. Finally, corresponding numerical experiments are given to demonstrate the effectiveness of the method.

Neural operators have garnered significant attention in recent years as a powerful framework for learning the mapping between parameter spaces and solution spaces of partial differential equations (PDEs). This paper proposes a novel U-Net variant neural operator based on dilated convolutions, termed U-DNet. The network fully leverages the multi-grid structure of U-Net for multi-scale feature extraction and fusion, while innovatively incorporating dilated convolutions to effectively enhance feature extraction capabilities at each scale. To validate the model’s performance, we conducted systematic experiments on multiple benchmark datasets, including multi-scale elliptic equations, Navier-Stokes equations, and Helmholtz equations. Experimental results demonstrate that U-DNet achieves significant advantages in both prediction accuracy and computational efficiency, highlighting its great potential as a multi-scale operator learning method.