The study of inverse problems has demonstrated significant value and broad application prospects across multiple domains. However, inverse problems in geophysics remain highly challenging due to their inherent complexity and uncertainty. In recent years, the rapid development and application of quantum computing technologies have opened new avenues for addressing this difficult problem. Based on the two currently prominent paradigms of quantum computing, this paper first introduces an iterative inversion algorithm based on quantum annealing. This algorithm integrates the global search capabilities of quantum annealing with the local optimization characteristics of iterative inversion. Through both linear and nonlinear tomography experiments, the proposed algorithm achieves inversion accuracy comparable to that of classical methods. In addition, this paper explores the potential of variational quantum linear system algorithm in solving inverse problems. We introduce a preliminary attempt using this algorithm in solving inverse problems, demonstrating promising application prospects by constructing and optimizing quantum circuits tailored for these challenges. The research presented in this paper provides new ideas and directions for the application of quantum computing in solving inverse problems in geophysics. With the continuous development and application of quantum computing technology, it is expected that quantum computing methods will play a significant role in solving large-scale geophysical inverse problems in the future.

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

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

Projection algorithms play an important role in tackling large-scale constrained optimization problems, especially when the orthogonal projection onto the constraint set is easy to implement. In projection algorithms, the choice of step size greatly affects the convergence rates. However, theoretically “optimal” step sizes often have poor effects in practical scenarios. Existing literature has demonstrated that relaxing the optimal step size for unconstrained optimization can achieve acceleration for steepest descent algorithm. Building on this insight, our work extends these relaxation strategies to projection algorithms, including both the fundamental scheme and the “prediction-correction” type, specifically for constrained convex quadratic optimization problems. Our numerical experiments indicate that these relaxation strategies significantly enhance the performance of both types of projection algorithms with optimal step sizes, demonstrating advantages across various problem sizes.

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.

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.

As the SPD stiffness equation’s scale expands, the memory requirements for direct solver constantly increase. Utilizing abundant hard disk resources to assist in reducing memory demand has become a solution. Three out-of-core(OOC) solving strategies using block supernode Cholesky method are proposed. 1) Fastest OOC strategy ensures that each supernode is read and saved once; 2) Minimum memory OOC strategy ensures that except for supernodes that still need to participate in factorization, all other supernodes are freed; 3) Restricted memory OOC strategy decides whether to adopt strategy 1) or 2) based on the given memory size. Each strategy uses OpenMP to implement asynchronous parallelism based on task pool in a shared memory environment. MUMPS failed to solve a 4984362 dimensional stiffness equation with OOC mode, while the solver using strategy 1) successfully solved with 2328.07s and only 18.2GB memory, saving 87% of 142.85GB memory required for IC solving; PARDISO OOC took 4643.18s and used 21.3GB memory.

An algorithm for ill-conditioned linear equations is proposed based on the support vector machine (SVM). The ill-conditioned linear equations are transformed to the regression problems for the hyper-plane in high dimensional space. And it can be solved by classic SVM and Least squares SVM (LS-SVM). Compared to the other methods of neural networks, SVM can obtain the numerical solutions with global and unique solution by small errors, since that it is convex quadratic programming algorithm. Numerical examples show that the presented method can achieve numerical solutions with relative high accuracy and stability.

With the development of computer science and information technology, computerized symbolic computation provides a more powerful approach in complex calculations and numerical simulation. Describing the nonlinear mechanism of shallow water wave, one type of (2 + 1)-dimensional generalized variable-coefficient shallow wave equation is investigated in this paper. With the Hirota bilinear method, bilinear form, soliton solutions, bilinear Bäacklund transformation and lump solution are derived and solved. The nonlinear mechanism of shallow water waves described by the soliton solution and lump solution are both illustrated with numerical simulations. One different and novel characteristics of the special solution to the (2 + 1)-dimensional generalized variable-coefficient shallow wave equation is that the lumps move along the y-axis with zero velocity in the x-axis direction.

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.

Geological parameters in petroleum reservoirs are highly heterogeneous, but the exploration approaches are very limited. In consequnce, geological models are usually highly uncertain. No matter conventional history matching or data space inversion, it is necessary to predict future reservoir dynamics under the uncertainty of geological models. This is conventionally done by reservoir simulation. Yet, for the large ensemble of realisattions to reflect the uncertainty of geological models, reservoir simulation is computaionally expensive. For this problem, we extend dynamic mode decomposition (DMD) to the parameter space, and conduct prediction based on a series of geological realisations reflecting the uncertainty of formation parameters. The change of reservoir flow variables tend to become gradual over time. Therefore, we build training and testing data sets, and determine the time after which the dynamic data approximately satisfies local linearity. Then DMD can be used to conduct prediction as a surrogate for reservoir simulation, in order to enhance the efficiency of prediction. The new method is validated using single and two-phase transient Darcy flow test cases.

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.

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.

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.

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.

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.

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 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.

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.

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.

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.