SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2709-A2736, August 2024. Abstract. The present study is an extension of the work done by Peddle, Haut, and Wingate [SIAM J. Sci. Comput., 41 (2019), pp. A3476–A3497] and Haut and Wingate [SIAM J. Sci. Comput., 36 (2014), pp. A693–A713], where a two-level Parareal method with mapping and averaging is examined. The method proposed in this paper is a multilevel Parareal method with arbitrarily many levels, which is not restricted to the two-level case. We give an asymptotic error estimate which reduces to the two-level estimate for the case when only two levels are considered. Introducing more than two levels has important consequences for the averaging procedure, as we choose separate averaging windows for each of the different levels, which is an additional new feature of the present study. The different averaging windows make the proposed method especially appropriate for nonlinear multiscale problems, because we can introduce a level for each intrinsic scale of the problem and adapt the averaging procedure such that we reproduce the behavior of the model on the particular scale resolved by the level. The method is applied to nonlinear differential equations. The nonlinearities can generate a range of frequencies in the problem. The computational cost of the new method is investigated and studied on several examples.
SIAM 科学计算期刊》，第 46 卷第 4 期，第 A2709-A2736 页，2024 年 8 月。 摘要本研究是对 Peddle、Haut 和 Wingate [SIAM J. Sci. Comput., 41 (2019), pp. A3476-A3497] 以及 Haut 和 Wingate [SIAM J. Sci. Comput., 36 (2014), pp.本文提出的方法是一种具有任意多层次的多层次 Parareal 方法，它并不局限于两层情况。我们给出了一个渐近误差估计值，在只考虑两级的情况下，该估计值与两级估计值相减。引入两个以上的水平对平均过程有重要影响，因为我们为每个不同的水平选择了不同的平均窗口，这是本研究的另一个新特点。不同的平均窗口使所提出的方法特别适用于非线性多尺度问题，因为我们可以为问题的每个内在尺度引入一个层次，并调整平均程序，从而重现模型在该层次所解析的特定尺度上的行为。该方法适用于非线性微分方程。非线性会在问题中产生一系列频率。新方法的计算成本在几个例子中进行了调查和研究。
{"title":"Multilevel Parareal Algorithm with Averaging for Oscillatory Problems","authors":"Juliane Rosemeier, Terry Haut, Beth Wingate","doi":"10.1137/23m1547123","DOIUrl":"https://doi.org/10.1137/23m1547123","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2709-A2736, August 2024. <br/> Abstract. The present study is an extension of the work done by Peddle, Haut, and Wingate [SIAM J. Sci. Comput., 41 (2019), pp. A3476–A3497] and Haut and Wingate [SIAM J. Sci. Comput., 36 (2014), pp. A693–A713], where a two-level Parareal method with mapping and averaging is examined. The method proposed in this paper is a multilevel Parareal method with arbitrarily many levels, which is not restricted to the two-level case. We give an asymptotic error estimate which reduces to the two-level estimate for the case when only two levels are considered. Introducing more than two levels has important consequences for the averaging procedure, as we choose separate averaging windows for each of the different levels, which is an additional new feature of the present study. The different averaging windows make the proposed method especially appropriate for nonlinear multiscale problems, because we can introduce a level for each intrinsic scale of the problem and adapt the averaging procedure such that we reproduce the behavior of the model on the particular scale resolved by the level. The method is applied to nonlinear differential equations. The nonlinearities can generate a range of frequencies in the problem. The computational cost of the new method is investigated and studied on several examples.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"22 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-08-20","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"142217605","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
Gabriel Barrenechea, Ernesto Castillo, Douglas Pacheco
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2660-A2682, August 2024. Abstract. This article investigates different implicit-explicit (IMEX) methods for incompressible flows with variable viscosity. The viscosity field may depend on space and time alone or, for example, on velocity gradients. Unlike most previous works on IMEX schemes, which focus on the convective term, we propose also treating parts of the diffusive term explicitly, which can reduce the coupling between the velocity components. We present different IMEX alternatives for the variable-viscosity Navier–Stokes system, analyzing their theoretical and algorithmic properties. Temporal stability is proven for all the methods presented, including monolithic and fractional-step variants. These results are unconditional except for one of the fractional-step discretizations, whose stability is shown for time-step sizes under an upper bound that depends solely on the problem data. The key finding of this work is a class of IMEX schemes whose steps decouple the velocity components and are fully linearized (even if the viscosity depends nonlinearly on the velocity) without requiring any CFL condition for stability. Moreover, in the presence of Neumann boundaries, some of our formulations lead naturally to conditions involving normal pseudotractions. This generalizes to the variable-viscosity case what happens for the standard Laplacian form with constant viscosity. Our analysis is supported by a series of numerical experiments.
{"title":"Implicit-explicit Schemes for Incompressible Flow Problems with Variable Viscosity","authors":"Gabriel Barrenechea, Ernesto Castillo, Douglas Pacheco","doi":"10.1137/23m1606526","DOIUrl":"https://doi.org/10.1137/23m1606526","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2660-A2682, August 2024. <br/> Abstract. This article investigates different implicit-explicit (IMEX) methods for incompressible flows with variable viscosity. The viscosity field may depend on space and time alone or, for example, on velocity gradients. Unlike most previous works on IMEX schemes, which focus on the convective term, we propose also treating parts of the diffusive term explicitly, which can reduce the coupling between the velocity components. We present different IMEX alternatives for the variable-viscosity Navier–Stokes system, analyzing their theoretical and algorithmic properties. Temporal stability is proven for all the methods presented, including monolithic and fractional-step variants. These results are unconditional except for one of the fractional-step discretizations, whose stability is shown for time-step sizes under an upper bound that depends solely on the problem data. The key finding of this work is a class of IMEX schemes whose steps decouple the velocity components and are fully linearized (even if the viscosity depends nonlinearly on the velocity) without requiring any CFL condition for stability. Moreover, in the presence of Neumann boundaries, some of our formulations lead naturally to conditions involving normal pseudotractions. This generalizes to the variable-viscosity case what happens for the standard Laplacian form with constant viscosity. Our analysis is supported by a series of numerical experiments.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"81 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-08-19","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"142217626","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2683-A2708, August 2024. Abstract. Building on the successes of local kernel methods for approximating the solutions to partial differential equations (PDEs) and the evaluation of definite integrals (quadrature/cubature), a local estimate of the error in such approximations is developed. This estimate is useful for determining locations in the solution domain where increased node density (equivalently, reduction in the spacing between nodes) can decrease the error in the solution. An adaptive procedure for adding nodes to the domain for both the approximation of derivatives and the approximate evaluation of definite integrals is described. This method efficiently computes the error estimate at a set of prescribed points and adds new nodes for approximation where the error is too large. Computational experiments demonstrate close agreement between the error estimate and actual absolute error in the approximation. Such methods are necessary or desirable when approximating solutions to PDEs (or in the case of quadrature/cubature), where the initial data and subsequent solution (or integrand) exhibit localized features that require significant refinement to resolve and where uniform increases in the density of nodes across the entire computational domain is not possible or too burdensome.
{"title":"Adaptivity in Local Kernel Based Methods for Approximating the Action of Linear Operators","authors":"Jonah A. Reeger","doi":"10.1137/23m1598052","DOIUrl":"https://doi.org/10.1137/23m1598052","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2683-A2708, August 2024. <br/> Abstract. Building on the successes of local kernel methods for approximating the solutions to partial differential equations (PDEs) and the evaluation of definite integrals (quadrature/cubature), a local estimate of the error in such approximations is developed. This estimate is useful for determining locations in the solution domain where increased node density (equivalently, reduction in the spacing between nodes) can decrease the error in the solution. An adaptive procedure for adding nodes to the domain for both the approximation of derivatives and the approximate evaluation of definite integrals is described. This method efficiently computes the error estimate at a set of prescribed points and adds new nodes for approximation where the error is too large. Computational experiments demonstrate close agreement between the error estimate and actual absolute error in the approximation. Such methods are necessary or desirable when approximating solutions to PDEs (or in the case of quadrature/cubature), where the initial data and subsequent solution (or integrand) exhibit localized features that require significant refinement to resolve and where uniform increases in the density of nodes across the entire computational domain is not possible or too burdensome.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"16 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-08-19","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"142217607","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page C479-C507, August 2024. Abstract. Investigating blood flow in the cardiovascular system is crucial for assessing cardiovascular health. Computational approaches offer some noninvasive alternatives to measure blood flow dynamics. Numerical simulations based on traditional methods such as finite-element and other numerical discretizations have been extensively studied and have yielded excellent results. However, adapting these methods to real-life simulations remains a complex task. In this paper, we propose a method that offers flexibility and can efficiently handle real-life simulations. We suggest utilizing the physics-informed neural network to solve the Navier–Stokes equation in a deformable domain, specifically addressing the simulation of blood flow in elastic vessels. Our approach models blood flow using an incompressible, viscous Navier–Stokes equation in an arbitrary Lagrangian–Eulerian form. The mechanical model for the vessel wall structure is formulated by an equation of Newton’s second law of momentum and linear elasticity to the force exerted by the fluid flow. Our method is a mesh-free approach that eliminates the need for discretization and meshing of the computational domain. This makes it highly efficient in solving simulations involving complex geometries. Additionally, with the availability of well-developed open-source machine learning framework packages and parallel modules, our method can easily be accelerated through GPU computing and parallel computing. To evaluate our approach, we conducted experiments on regular cylinder vessels as well as vessels with plaque on their walls. We compared our results to a solution calculated by finite element methods using a dense grid and small time steps, which we considered as the ground truth solution. We report the relative error and the time consumed to solve the problem, highlighting the advantages of our method.
{"title":"A Meshless Solver for Blood Flow Simulations in Elastic Vessels Using a Physics-Informed Neural Network","authors":"Han Zhang, Raymond H. Chan, Xue-Cheng Tai","doi":"10.1137/23m1622696","DOIUrl":"https://doi.org/10.1137/23m1622696","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page C479-C507, August 2024. <br/> Abstract. Investigating blood flow in the cardiovascular system is crucial for assessing cardiovascular health. Computational approaches offer some noninvasive alternatives to measure blood flow dynamics. Numerical simulations based on traditional methods such as finite-element and other numerical discretizations have been extensively studied and have yielded excellent results. However, adapting these methods to real-life simulations remains a complex task. In this paper, we propose a method that offers flexibility and can efficiently handle real-life simulations. We suggest utilizing the physics-informed neural network to solve the Navier–Stokes equation in a deformable domain, specifically addressing the simulation of blood flow in elastic vessels. Our approach models blood flow using an incompressible, viscous Navier–Stokes equation in an arbitrary Lagrangian–Eulerian form. The mechanical model for the vessel wall structure is formulated by an equation of Newton’s second law of momentum and linear elasticity to the force exerted by the fluid flow. Our method is a mesh-free approach that eliminates the need for discretization and meshing of the computational domain. This makes it highly efficient in solving simulations involving complex geometries. Additionally, with the availability of well-developed open-source machine learning framework packages and parallel modules, our method can easily be accelerated through GPU computing and parallel computing. To evaluate our approach, we conducted experiments on regular cylinder vessels as well as vessels with plaque on their walls. We compared our results to a solution calculated by finite element methods using a dense grid and small time steps, which we considered as the ground truth solution. We report the relative error and the time consumed to solve the problem, highlighting the advantages of our method.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"6 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-08-19","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"142217606","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2634-A2659, August 2024. Abstract. The multigrid V-cycle method is a popular method for solving systems of linear equations. It computes an approximate solution by using smoothing on fine levels and solving a system of linear equations on the coarsest level. Solving on the coarsest level depends on the size and difficulty of the problem. If the size permits, it is typical to use a direct method based on LU or Cholesky decomposition. In settings with large coarsest-level problems, approximate solvers such as iterative Krylov subspace methods, or direct methods based on low-rank approximation, are often used. The accuracy of the coarsest-level solver is typically determined based on the experience of the users with the concrete problems and methods. In this paper, we present an approach to analyzing the effects of approximate coarsest-level solves on the convergence of the V-cycle method for symmetric positive definite problems. Using these results, we derive coarsest-level stopping criterion through which we may control the difference between the approximation computed by a V-cycle method with approximate coarsest-level solver and the approximation which would be computed if the coarsest-level problems were solved exactly. The coarsest-level stopping criterion may thus be set up such that the V-cycle method converges to a chosen finest-level accuracy in (nearly) the same number of V-cycle iterations as the V-cycle method with exact coarsest-level solver. We also utilize the theoretical results to discuss how the convergence of the V-cycle method may be affected by the choice of a tolerance in a coarsest-level stopping criterion based on the relative residual norm. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://doi.org/10.5281/zenodo.11178544.
SIAM 科学计算期刊》，第 46 卷第 4 期，第 A2634-A2659 页，2024 年 8 月。 摘要多网格 V 循环法是求解线性方程组的一种常用方法。它通过在精细级上进行平滑处理并在最粗级上求解线性方程组来计算近似解。在最粗级上求解取决于问题的大小和难度。如果规模允许，通常采用基于 LU 或 Cholesky 分解的直接方法。在处理大型最粗级问题时，通常会使用近似求解器，如迭代克雷洛夫子空间方法或基于低阶近似的直接方法。最粗级求解器的精度通常根据用户对具体问题和方法的经验来确定。在本文中，我们提出了一种分析近似最粗级求解对对称正定问题 V 循环方法收敛性影响的方法。利用这些结果，我们得出了最粗级停止准则，通过该准则，我们可以控制使用近似最粗级求解器的 V 循环方法计算出的近似值与精确求解最粗级问题时计算出的近似值之间的差异。因此，可以设定最粗级停止准则，使 V 循环方法在（几乎）与使用精确最粗级求解器的 V 循环方法相同的 V 循环迭代次数内收敛到所选的最细级精度。我们还利用理论结果讨论了 V 循环方法的收敛性如何受到基于相对残差规范的最粗级停止准则中容限选择的影响。计算结果的可重复性。本文被授予 "SIAM 可重复性徽章"：代码和数据可用"，以表彰作者遵循了 SISC 和科学计算界重视的可重现性原则。读者可通过 https://doi.org/10.5281/zenodo.11178544 获取代码和数据，以重现本文中的结果。
{"title":"The Effect of Approximate Coarsest-Level Solves on the Convergence of Multigrid V-Cycle Methods","authors":"Petr Vacek, Erin Carson, Kirk M. Soodhalter","doi":"10.1137/23m1578255","DOIUrl":"https://doi.org/10.1137/23m1578255","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2634-A2659, August 2024. <br/> Abstract. The multigrid V-cycle method is a popular method for solving systems of linear equations. It computes an approximate solution by using smoothing on fine levels and solving a system of linear equations on the coarsest level. Solving on the coarsest level depends on the size and difficulty of the problem. If the size permits, it is typical to use a direct method based on LU or Cholesky decomposition. In settings with large coarsest-level problems, approximate solvers such as iterative Krylov subspace methods, or direct methods based on low-rank approximation, are often used. The accuracy of the coarsest-level solver is typically determined based on the experience of the users with the concrete problems and methods. In this paper, we present an approach to analyzing the effects of approximate coarsest-level solves on the convergence of the V-cycle method for symmetric positive definite problems. Using these results, we derive coarsest-level stopping criterion through which we may control the difference between the approximation computed by a V-cycle method with approximate coarsest-level solver and the approximation which would be computed if the coarsest-level problems were solved exactly. The coarsest-level stopping criterion may thus be set up such that the V-cycle method converges to a chosen finest-level accuracy in (nearly) the same number of V-cycle iterations as the V-cycle method with exact coarsest-level solver. We also utilize the theoretical results to discuss how the convergence of the V-cycle method may be affected by the choice of a tolerance in a coarsest-level stopping criterion based on the relative residual norm. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://doi.org/10.5281/zenodo.11178544.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"23 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-08-14","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"142217608","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page B554-B578, August 2024. Abstract. We present an algorithm to solve the dispersive depth-averaged Serre–Green–Naghdi equations using patch-based adaptive mesh refinement. These equations require adding additional higher derivative terms to the nonlinear shallow water equations. This has been implemented as a new component of the open source GeoClaw software that is widely used for modeling tsunamis, storm surge, and related hazards, improving its accuracy on shorter wavelength phenomena. We use a formulation that requires solving an elliptic system of equations at each time step, making the method implicit. The adaptive algorithm allows different time steps on different refinement levels and solves the implicit equations level by level. Computational examples are presented to illustrate the stability and accuracy on a radially symmetric test case and two realistic tsunami modeling problems, including a hypothetical asteroid impact creating a short wavelength tsunami for which dispersive terms are necessary. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/rjleveque/ImplicitAMR-paper and in the supplementary materials (ImplicitAMR-paper.zip [174KB]).
{"title":"Implicit Adaptive Mesh Refinement for Dispersive Tsunami Propagation","authors":"Marsha J. Berger, Randall J. LeVeque","doi":"10.1137/23m1585210","DOIUrl":"https://doi.org/10.1137/23m1585210","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page B554-B578, August 2024. <br/> Abstract. We present an algorithm to solve the dispersive depth-averaged Serre–Green–Naghdi equations using patch-based adaptive mesh refinement. These equations require adding additional higher derivative terms to the nonlinear shallow water equations. This has been implemented as a new component of the open source GeoClaw software that is widely used for modeling tsunamis, storm surge, and related hazards, improving its accuracy on shorter wavelength phenomena. We use a formulation that requires solving an elliptic system of equations at each time step, making the method implicit. The adaptive algorithm allows different time steps on different refinement levels and solves the implicit equations level by level. Computational examples are presented to illustrate the stability and accuracy on a radially symmetric test case and two realistic tsunami modeling problems, including a hypothetical asteroid impact creating a short wavelength tsunami for which dispersive terms are necessary. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/rjleveque/ImplicitAMR-paper and in the supplementary materials (ImplicitAMR-paper.zip [174KB]).","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"26 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-08-14","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"142217612","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2607-A2633, August 2024. Abstract. Tikhonov regularization is a widely used technique in solving inverse problems that can enforce prior properties on the desired solution. In this paper, we propose a Krylov subspace based iterative method for solving linear inverse problems with general-form Tikhonov regularization term [math], where [math] is a positive semidefinite matrix. An iterative process called the preconditioned Golub–Kahan bidiagonalization (pGKB) is designed, which implicitly utilizes a proper preconditioner to generate a series of solution subspaces with desirable properties encoded by the regularizer [math]. Based on the pGKB process, we propose an iterative regularization algorithm via projecting the original problem onto small dimensional solution subspaces. We analyze the regularization properties of this algorithm, including the incorporation of prior properties of the desired solution into the solution subspace and the semiconvergence behavior of the regularized solution. To overcome instabilities caused by semiconvergence, we further propose two pGKB based hybrid regularization algorithms. All the proposed algorithms are tested on both small-scale and large-scale linear inverse problems. Numerical results demonstrate that these iterative algorithms exhibit excellent performance, outperforming other state-of-the-art algorithms in some cases.
{"title":"A Preconditioned Krylov Subspace Method for Linear Inverse Problems with General-Form Tikhonov Regularization","authors":"Haibo Li","doi":"10.1137/23m1593802","DOIUrl":"https://doi.org/10.1137/23m1593802","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2607-A2633, August 2024. <br/> Abstract. Tikhonov regularization is a widely used technique in solving inverse problems that can enforce prior properties on the desired solution. In this paper, we propose a Krylov subspace based iterative method for solving linear inverse problems with general-form Tikhonov regularization term [math], where [math] is a positive semidefinite matrix. An iterative process called the preconditioned Golub–Kahan bidiagonalization (pGKB) is designed, which implicitly utilizes a proper preconditioner to generate a series of solution subspaces with desirable properties encoded by the regularizer [math]. Based on the pGKB process, we propose an iterative regularization algorithm via projecting the original problem onto small dimensional solution subspaces. We analyze the regularization properties of this algorithm, including the incorporation of prior properties of the desired solution into the solution subspace and the semiconvergence behavior of the regularized solution. To overcome instabilities caused by semiconvergence, we further propose two pGKB based hybrid regularization algorithms. All the proposed algorithms are tested on both small-scale and large-scale linear inverse problems. Numerical results demonstrate that these iterative algorithms exhibit excellent performance, outperforming other state-of-the-art algorithms in some cases.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"3 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-08-14","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"142227248","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2582-A2606, August 2024. Abstract. We introduce a fast direct solver for variable-coefficient elliptic PDEs on surfaces based on the hierarchical Poincaré–Steklov method. The method takes as input an unstructured, high-order quadrilateral mesh of a surface and discretizes surface differential operators on each element using a high-order spectral collocation scheme. Elemental solution operators and Dirichlet-to-Neumann maps tangent to the surface are precomputed and merged in a pairwise fashion to yield a hierarchy of solution operators that may be applied in [math] operations for a mesh with [math] degrees of freedom. The resulting fast direct solver may be used to accelerate high-order implicit time-stepping schemes, as the precomputed operators can be reused for fast elliptic solves on surfaces. On a standard laptop, precomputation for a 12th-order surface mesh with over 1 million degrees of freedom takes 10 seconds, while subsequent solves take only 0.25 seconds. We apply the method to a range of problems on both smooth surfaces and surfaces with sharp corners and edges, including the static Laplace–Beltrami problem, the Hodge decomposition of a tangential vector field, and some time-dependent nonlinear reaction-diffusion systems. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: code and data available”, as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/danfortunato/surface-hps-sisc.
{"title":"A High-Order Fast Direct Solver for Surface PDEs","authors":"Daniel Fortunato","doi":"10.1137/22m1525259","DOIUrl":"https://doi.org/10.1137/22m1525259","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2582-A2606, August 2024. <br/> Abstract. We introduce a fast direct solver for variable-coefficient elliptic PDEs on surfaces based on the hierarchical Poincaré–Steklov method. The method takes as input an unstructured, high-order quadrilateral mesh of a surface and discretizes surface differential operators on each element using a high-order spectral collocation scheme. Elemental solution operators and Dirichlet-to-Neumann maps tangent to the surface are precomputed and merged in a pairwise fashion to yield a hierarchy of solution operators that may be applied in [math] operations for a mesh with [math] degrees of freedom. The resulting fast direct solver may be used to accelerate high-order implicit time-stepping schemes, as the precomputed operators can be reused for fast elliptic solves on surfaces. On a standard laptop, precomputation for a 12th-order surface mesh with over 1 million degrees of freedom takes 10 seconds, while subsequent solves take only 0.25 seconds. We apply the method to a range of problems on both smooth surfaces and surfaces with sharp corners and edges, including the static Laplace–Beltrami problem, the Hodge decomposition of a tangential vector field, and some time-dependent nonlinear reaction-diffusion systems. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: code and data available”, as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/danfortunato/surface-hps-sisc.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"23 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-08-13","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"142217609","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page C448-C478, August 2024. Abstract. We studied the least-squares ReLU neural network (LSNN) method for solving a linear advection-reaction equation with discontinuous solution in [Z. Cai et al., J. Comput. Phys., 443 (2021), 110514]. The method is based on a least-squares formulation and uses a new class of approximating functions: ReLU neural network (NN) functions. A critical and additional component of the LSNN method, differing from other NN-based methods, is the introduction of a properly designed and physics preserved discrete differential operator. In this paper, we study the LSNN method for problems with discontinuity interfaces. First, we show that ReLU NN functions with depth [math] can approximate any [math]-dimensional step function on a discontinuity interface generated by a vector field as streamlines with any prescribed accuracy. By decomposing the solution into continuous and discontinuous parts, we prove theoretically that the discretization error of the LSNN method using ReLU NN functions with depth [math] is mainly determined by the continuous part of the solution provided that the solution jump is constant. Numerical results for both two- and three-dimensional test problems with various discontinuity interfaces show that the LSNN method with enough layers is accurate and does not exhibit the common Gibbs phenomena along discontinuity interfaces.
SIAM 科学计算期刊》，第 46 卷第 4 期，第 C448-C478 页，2024 年 8 月。 摘要我们在[Z. Cai et al., J. Comput. Phys., 443 (2021), 110514]一文中研究了求解具有不连续解的线性平流反应方程的最小二乘 ReLU 神经网络（LSNN）方法。该方法基于最小二乘法，并使用了一类新的近似函数：ReLU 神经网络 (NN) 函数。与其他基于 NN 的方法不同，LSNN 方法的一个关键和额外的组成部分是引入了一个经过适当设计并保留了物理特性的离散微分算子。在本文中，我们研究了针对不连续界面问题的 LSNN 方法。首先，我们证明了深度为[math]的 ReLU NN 函数可以以任意规定的精度将矢量场产生的不连续界面上的任意[math]维阶跃函数近似为流线。通过将解分解为连续部分和不连续部分，我们从理论上证明了使用深度[数学]ReLU NN 函数的 LSNN 方法的离散化误差主要由解的连续部分决定，前提是解的跳跃是恒定的。对具有各种不连续界面的二维和三维测试问题的数值结果表明，具有足够多层次的 LSNN 方法是精确的，不会在不连续界面上出现常见的吉布斯现象。
{"title":"Least-Squares Neural Network (LSNN) Method for Linear Advection-Reaction Equation: Discontinuity Interface","authors":"Zhiqiang Cai, Junpyo Choi, Min Liu","doi":"10.1137/23m1568107","DOIUrl":"https://doi.org/10.1137/23m1568107","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page C448-C478, August 2024. <br/> Abstract. We studied the least-squares ReLU neural network (LSNN) method for solving a linear advection-reaction equation with discontinuous solution in [Z. Cai et al., J. Comput. Phys., 443 (2021), 110514]. The method is based on a least-squares formulation and uses a new class of approximating functions: ReLU neural network (NN) functions. A critical and additional component of the LSNN method, differing from other NN-based methods, is the introduction of a properly designed and physics preserved discrete differential operator. In this paper, we study the LSNN method for problems with discontinuity interfaces. First, we show that ReLU NN functions with depth [math] can approximate any [math]-dimensional step function on a discontinuity interface generated by a vector field as streamlines with any prescribed accuracy. By decomposing the solution into continuous and discontinuous parts, we prove theoretically that the discretization error of the LSNN method using ReLU NN functions with depth [math] is mainly determined by the continuous part of the solution provided that the solution jump is constant. Numerical results for both two- and three-dimensional test problems with various discontinuity interfaces show that the LSNN method with enough layers is accurate and does not exhibit the common Gibbs phenomena along discontinuity interfaces.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"95 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-08-12","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"142217611","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
Aleksandr Y. Aravkin, Robert Baraldi, Dominique Orban
SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2557-A2581, August 2024. Abstract. We develop a Levenberg–Marquardt method for minimizing the sum of a smooth nonlinear least-squares term [math] and a nonsmooth term [math]. Both [math] and [math] may be nonconvex. Steps are computed by minimizing the sum of a regularized linear least-squares model and a model of [math] using a first-order method such as the proximal gradient method. We establish global convergence to a first-order stationary point under the assumptions that [math] and its Jacobian are Lipschitz continuous and [math] is proper and lower semicontinuous. In the worst case, our method performs [math] iterations to bring a measure of stationarity below [math]. We also derive a trust-region variant that enjoys similar asymptotic worst-case iteration complexity as a special case of the trust-region algorithm of Aravkin, Baraldi, and Orban [SIAM J. Optim., 32 (2022), pp. 900–929]. We report numerical results on three examples: a group-lasso basis-pursuit denoise example, a nonlinear support vector machine, and parameter estimation in a neuroscience application. To implement those examples, we describe in detail how to evaluate proximal operators for separable [math] and for the group lasso with trust-region constraint. In all cases, the Levenberg–Marquardt methods perform fewer outer iterations than either a proximal gradient method with adaptive step length or a quasi-Newton trust-region method, neither of which exploit the least-squares structure of the problem. Our results also highlight the need for more sophisticated subproblem solvers than simple first-order methods.
{"title":"A Levenberg–Marquardt Method for Nonsmooth Regularized Least Squares","authors":"Aleksandr Y. Aravkin, Robert Baraldi, Dominique Orban","doi":"10.1137/22m1538971","DOIUrl":"https://doi.org/10.1137/22m1538971","url":null,"abstract":"SIAM Journal on Scientific Computing, Volume 46, Issue 4, Page A2557-A2581, August 2024. <br/> Abstract. We develop a Levenberg–Marquardt method for minimizing the sum of a smooth nonlinear least-squares term [math] and a nonsmooth term [math]. Both [math] and [math] may be nonconvex. Steps are computed by minimizing the sum of a regularized linear least-squares model and a model of [math] using a first-order method such as the proximal gradient method. We establish global convergence to a first-order stationary point under the assumptions that [math] and its Jacobian are Lipschitz continuous and [math] is proper and lower semicontinuous. In the worst case, our method performs [math] iterations to bring a measure of stationarity below [math]. We also derive a trust-region variant that enjoys similar asymptotic worst-case iteration complexity as a special case of the trust-region algorithm of Aravkin, Baraldi, and Orban [SIAM J. Optim., 32 (2022), pp. 900–929]. We report numerical results on three examples: a group-lasso basis-pursuit denoise example, a nonlinear support vector machine, and parameter estimation in a neuroscience application. To implement those examples, we describe in detail how to evaluate proximal operators for separable [math] and for the group lasso with trust-region constraint. In all cases, the Levenberg–Marquardt methods perform fewer outer iterations than either a proximal gradient method with adaptive step length or a quasi-Newton trust-region method, neither of which exploit the least-squares structure of the problem. Our results also highlight the need for more sophisticated subproblem solvers than simple first-order methods.","PeriodicalId":49526,"journal":{"name":"SIAM Journal on Scientific Computing","volume":"34 1","pages":""},"PeriodicalIF":3.1,"publicationDate":"2024-08-12","publicationTypes":"Journal Article","fieldsOfStudy":null,"isOpenAccess":false,"openAccessPdf":"","citationCount":null,"resultStr":null,"platform":"Semanticscholar","paperid":"142217628","PeriodicalName":null,"FirstCategoryId":null,"ListUrlMain":null,"RegionNum":2,"RegionCategory":"数学","ArticlePicture":[],"TitleCN":null,"AbstractTextCN":null,"PMCID":"","EPubDate":null,"PubModel":null,"JCR":null,"JCRName":null,"Score":null,"Total":0}
扫码关注我们
求助内容：
应助结果提醒方式：
京公网安备 11010802042870号