Galmoss： GPU 加速的 Galaxy Profile Fitting 软件包

摘要

我们介绍 galmoss，这是一个基于 python、torch 驱动的工具，用于二维拟合星系轮廓。 通过无缝启用 GPU 并行化，galmoss 满足了大规模星系调查的高计算需求，将星系轮廓拟合置于 LSST 时代。 它包含广泛使用的配置文件，例如 Sérsic、指数盘、费雷尔、King、高斯和莫法特配置文件，并允许轻松集成更复杂的模型。 使用单个 NVIDIA A100 GPU 对斯隆数字巡天 (SDSS) g 波段的 8,289 个星系进行了测试，galmoss 在大约 10 分钟内完成了经典的 Sérsic 剖面拟合。 基准测试表明，galmoss 的计算速度比默认实现快 6 $\times$。

1简介

星系是宇宙的组成部分，由不同的恒星组成部分组成，例如核球、盘、棒、旋臂和核星团。 河外天文学的主要驱动力之一是根据星系的光度图像研究星系的结构和形态特性，这已被证明与星系的形成和演化路径相关（van der Wel，2008；Conselice， 2014；迪毛罗等人，2022）。

已经开发了几种用于星系结构和形态分析的方法。 经典的眼球形态分类，例如哈勃序列（Hubble，1926），通常是描述性的并且依赖于视觉检查。 非参数形态分析（Ferrari 等人，2015；Rodriguez-Gomez 等人，2019）采用定量指标来表示浓度、不对称性和平滑度等属性（Conselice，2003）. 这些指标被认为是稳健的，因为它们不依赖于任何底层模型。 相比之下，星系表面亮度剖面的分析涉及到将光分布与参数模型进行拟合，以利用特定的模型参数来解释和量化形态（Nantais等人，2013；庄和Ho，2022）。 这种轮廓拟合方法已经从一维（de Vaucoulurs，1958；Sersic，1968）发展到二维分析，通过结合点扩散函数和非轴对称分量等因素来提高精度（Andredakis 等人，1995；Schade 等人，1995；Byun 和 Freeman，1995）。

星系中二维光剖面的拟合变得越来越复杂，需要更快、更可扩展的算法。 两种主要方法占主导地位：马尔可夫链蒙特卡罗 (MCMC) 和基于梯度的方法。 基于 MCMC 的项目，例如 profit (Robotham 等人, 2017) 和 autogalaxy (Nightingale 等人, 2023)，迭代采样以近似目标分布，拟合值和不确定性源自马尔可夫链的统计数据。 相反，基于梯度的方法，由 galfit (Peng 等人，2010) 和 imfit (Erwin, 2015) 等工具使用，通过直接导航到损失函数的局部最小值来提供更快的解决方案。 然而，这种速度会损害参数空间的探索和不确定性估计的稳健性，而基于 MCMC 的方法中不存在这种权衡（Chen 等人，2016）。

随着斯隆数字巡天（SDSS，约克等人，2000）、暗能量巡天（DES，协作：等人，2016）等天文观测数据量的不断扩大、中国巡天太空望远镜(CSST, Zhan, 2011)、遗留时空巡天(LSST, Abell 等人, 2009)、以及 EUCLID (Laureijs等人，2011)，星系轮廓拟合任务出现了重大挑战。 基于梯度的方法和 MCMC 方法在有效处理大量星系方面都面临着困难。

为了应对这些挑战，越来越多地采用深度学习策略。 例如，deeplegato (Tuccillo 等人, 2018) 和 galnets (Li 等人, 2022; Qiu 等人, 2023) 训练卷积神经网络从模拟星系图像中预测轮廓拟合参数。 这些深度学习方法在星系图像与其轮廓参数之间建立了隐式联系，大大加速了参数估计的过程（另请参阅 gamornet，Ghosh，2019）。 然而，由于众所周知的“黑匣子”问题（Castelvecchi，2016），这些方法不能轻易适应新的数据域，这意味着模型生成的隐式映射通常缺乏清晰的解释。

解决速度挑战的另一种方法涉及基于并行梯度的方法，该方法完全依赖于 GPU 加速的矩阵并行计算。 pytorch (Paszke 等人, 2019)、tensorflow (Abadi 等人, 2016)< 等流行框架都支持这一点/t3> 和 jax (Bradbury 等人, 2018)，支持自动微分的结合。 这种方法的应用范围从星系运动学建模（Bekiaris，2017）到宇宙学N体模拟（Modi等人，2021）。

为了加速星系轮廓拟合，将其带入 LSST 时代，并为星系形成和演化的文献做出贡献，我们引入了 galmoss。 我们的目标是通过提供一个框架来利用传统和现代方法的优势，该框架能够实现无缝并行化，同时保留更传统方法的可解释性。 galmoss 将批量星系光度图像处理为多维数组，从而促进 GPU 上的高效拟合。 这种基于张量的方法允许通过调整批量大小和优化 GPU 利用率来动态管理计算期间的数据大小。 对于拟合过程，galmoss 由于其相对较低的计算和内存要求而使用梯度下降。 最后，galmoss提供了两种不确定性估计方法。 在星系轮廓拟合领域，astropot （Stone 等人，2023） 在其串行优化过程中利用 GPU 加速和自动微分，同时保留物理可解释性。 虽然我们的软件包专注于处理单个星系的大量图像，但 astrophoto 已针对在单个拥挤的图像中拟合多个星系进行了优化。

本文的结构如下：2 部分概述了 galmoss 工作流程。 3 部分描述了我们用来展示软件性能的数据。 4部分以Sérsic配置文件为例详细介绍了模型图像的生成，并介绍了一组内置配置文件以及高级应用程序。 5部分讨论了图像生成后的拟合过程。 最后，6 部分通过案例研究展示了galmoss 的表现，并在7 部分中提出了结论性意见。

2 galmoss的一般工作流程

图1： galmoss工作流程：从参数初始化和图像处理到拟合过程和最终数据产品生成，包括拟合模型和残差，全部以FITS格式存储。

图1展示了galmoss的工作流程。 它首先读取由用户定义的初始参数值组成的向量$p=\{p_G,p_L\}$，以及数据图像、西格玛图像（数据不确定性）、掩模图像和点扩散函数（以下简称PSF）图像。 第一个子集$p_G$描述了几何属性，它定义了椭圆的等光线的中心位置、长轴/短轴以及位置角。 第二个子集$p_L$读取星系径向亮度分布$\zeta (𝐫)$的初始值。

建立这些初始值后，galmoss 会根据选定的配置文件生成一系列图像模型。 然后，这些模型与图像的 PSF 进行卷积，并通过向所有像素添加平均天空值来进行增强。 此过程产生图像模型 $M_{i,j}$，然后将其与图像数据 $D_{i,j}$ 和相应的 sigma 图像 ${\sigma }_{i,j}$ 进行比较。 拟合过程采用 ${\chi }^2$ 损失函数以及掩模 ${\mathrm{\Omega }}_{i,j}$ 来确定优化中包含哪些像素。 在内部，galmoss 采用嵌入在 pytorch 中的梯度下降方法。

完成拟合过程后，该包将数据产品（包括图像块和拟合参数值）保存到指定的输出路径。 以FITS格式保存的图像块文件包括原始图像、最佳拟合模型图像和残差。 每个组件都存储在给定的 FITS 头部数据单元 (HDU) 中，在安装多个组件（例如凸块 + 圆盘）的情况下，每个组件也可以与其他数据产品一起单独访问。

3数据

为了展示 galmoss 功能，我们从 SDSS 目录中收集了星系图像($128\times 128$ 像素）和 PSF 图像($41\times 41$ 像素）。 我们的数据选择以独立星系参数的可用性为指导，由 MaNGA PyMorph 光度增值目录 （MPP-VAC-DR17，Domínguez Sánchez 等人，2022） 提供，涵盖来自ManNGA 最终版本（Bundy 等人，2014）。 所选样本主要由平均红移 z $\sim$ 0.03 的星系组成，并且具有相对均匀的恒星质量分布（协作，2022）。 该目录提供了通过二维表面轮廓拟合获得的光度参数。 这些参数包括来自单个 Sérsic 和组合 Sérsic + Exponential 模型的结果，使用 pymorph 管道（Vikram 等人，2010） 计算。 pymorph 集成了 sextractor 和 galfit，前者简化了初始参数的估计和图章的生成，后者旨在拟合银河系轮廓。

在我们的选择过程中，我们排除了位于板块边缘的星系，因为在这些星系中不可能获得完整的 $128\times 128$ 像素图像。 我们还遗漏了 MPP-VAC-DR17 目录中单次 Sérsic 拟合过程中失败的星系。 最终样本包含 8,289 个星系。 我们仅使用 g 波段图像来确保整个样本具有良好的信噪比。 在 pymorph 之后，我们使用 sextractor （Bertin 和 Arnouts，1996） 生成蒙版图像和初始值。

4 图像生成和内置配置文件

对于给定的一组星系或轮廓参数，galmoss 根据所选轮廓生成表示每个像素的强度或表面亮度的模型图像。 galmoss 提供六种内置径向剖面：Sérsic、Exponential disk、Ferrer、King、Gaussian 和 Moffat，以及 Flat 天空模型（可以选择包含更复杂的天空模型）。 通过假设圆形对称，这些径向轮廓被转换成二维图像。 此过程有效地将一维轮廓投影到二维圆形区域上以创建最终图像。 为了更具体地说明这个投影过程，考虑到 Sérsic 剖面在建模星系剖面中的广泛应用，我们以 Sérsic 剖面为例。

此外，galmoss支持不同配置文件的组合和集成。 有关代码实际使用的详细说明可以在B、C和D中找到。

4.1服务

图2： J162123.19+322056.4的g波段galmoss拟合结果，显示出良好的拟合质量。 左边四张图显示了模型、模型残差（数据模型）、数据和残差分布。 右侧显示一维投影，其中上图用黑线表示模型分布，带有误差条的橙色填充标记表示星系图像的通量分布及其不确定性。 下图表示通量残留，带有蓝色填充标记和误差条。

Sérsic剖面（Sersic，1968）通常用于拟合椭圆星系的表面亮度分布，以及其他星系类型的盘和核球部分。 强度 $I(r)$ 作为半径 $r$ 的函数由下式给出

$$I(r)=I_{\mathrm{e}}\exp \left\{-{\nu }_n\left[{\left(\frac{r}{r_{\mathrm{e}}}\right)}^{\frac{1}{n}}-1\right]\right\},$$ (1)

其中$I_{\mathrm{e}}$表示有效半径$r_{\mathrm{e}}$处的表面亮度，它包含总轮廓通量的一半。 Sérsic 指数 $n$ 指示轮廓的曲率，参数 ${\nu }_n$ 是根据伽玛分布的逆累积分布函数进行数值计算的。 Sérsic 剖面的多功能性在于它能够通过调整 Sérsic 指数 $n$ 来模拟各种星系形态。例如，设置 $n=4$ 会将其转换为 de Vaucoulurs 配置文件，称为 $r^{1/4}$ 配置文件。 使用 $n=1$ 获得指数磁盘分布，$n=0.5$ 产生高斯分布，使用户能够通过简单地改变 $n$ 来探索这些不同的分布。

在一维环境中，半径$r$是指距星系中心的距离。 在二维中，给定星系轮廓中心$(x_{\mathrm{c}},y_{\mathrm{c}})$、椭圆轮廓的位置角$\theta$、轴比$q$、圆半径$r$ 由 定义（参见例如 Robotham 等人，2017）：

$r$ $={\left[r_{\mathrm{maj}}^{(B+2)}+{\left({\displaystyle \frac{r_{\min }}{q}}\right)}^{(B+2)}\right]}^{\frac{1}{B+2}},$ (2)

在哪里

$r_{\mathrm{maj}}$ $=\left|\cos (\theta )(x-x_{\mathrm{c}})\right|+\left|\sin (\theta )(y-y_{\mathrm{c}})\right|,$ (3)

$r_{\min }$ $=\left|-\sin (\theta )(x-x_{\mathrm{c}})\right|+\left|\cos (\theta )(y-y_{\mathrm{c}})\right|.$ (4)

箱形参数$B$引入了椭圆形状的灵活性（Athanassoula等人，1990），允许更通用的椭圆轮廓。 $B=0$ 对应于标准椭圆，正值会导致更多的盒状形状，负值会导致更多的盘状形状。 这种灵活性在星系建模中特别有用，其中一系列椭圆形状可以代表在星系中观察到的不同形态。

不是直接使用 Sérsic 中的表面亮度 $I_{\mathrm{e}}$，而是使用幅度 $m$ 及其相应的零点幅度 $m_0$ 来指定强度级别。 这种方法与广泛使用的星系拟合工具所采用的方法一致，例如 galfit (Peng 等人, 2010) 和 profit (Robotham 等人, 2017)，便于与观测数据进行比较。 两个量的关系如下：

$$I_{\mathrm{e}}=\frac{r_{\mathrm{box}}{\nu }_n^{2n}{10}^{-0.4(m-m_0)}}{2\pi q{r}_{\mathrm{e}}^{2}n\mathrm{\Gamma }(2n)\exp ({\nu }_n)},$$ (5)

在哪里

$$r_{\mathrm{box}}=\frac{\pi (2+B)}{2\beta (\frac{1}{2+B},\frac{1}{2+B})}.$$ (6)

在此公式中，$r_{\mathrm{box}}$ 用作 $I_{\mathrm{e}}$ 中的几何校正因子，以考虑受圆度或箱形水平影响的完美椭圆的偏差。 具体来说，当$B$ = 0时，表示完美椭圆，$r_{\text{box}}$ = 1时，表示没有几何校正。 beta 函数 $\beta (a,b)$ 使用关系 $\beta (a,b)=\frac{\mathrm{\Gamma }(a)\mathrm{\Gamma }(b)}{\mathrm{\Gamma }(a+b)}$ 计算，其中 $\mathrm{\Gamma }$ 表示 Gamma 函数。

图2展示了使用星系J162123.19+322056.4（我们目录中的星系之一）的Sérsic剖面的拟合结果。 该模型的质量由残差证明，突出了拟合的精度。 对于那些对实现细节感兴趣的人，A 中提供了用于此分析的代码。

4.2 其他可用的配置文件

在本节中，我们将简要讨论 Sérsic 之外的内置配置文件。 除平坦天空轮廓外，所有这些轮廓都具有遵循方程 2 至方程 4 的圆形半径 $r$。

4.2.1 指数磁盘配置文件

在指数盘剖面中，强度定义如下：

$$I(r)=I_0\exp \left(-\frac{r}{r_{{}_{\mathrm{s}}}}\right),$$ (7)

其中 $I_0$ 是中心轮廓表面的亮度（在径向距离 $r=0$ 处），$r_{{}_{\mathrm{s}}}$ 是圆盘刻度长度。

幅度($m,m_0)$、$r_{{}_{\mathrm{s}}}$ 和 $I_0$ 之间的关系由下式给出

$$I_0=\frac{r_{\mathrm{box}}{10}^{-0.4(m-m_0)}}{2\pi q{r}_{\mathrm{s}}^2},$$ (8)

其中$q$是轴比，$r_{\mathrm{box}}$如方程6中所定义。

4.2.2 修改的 Ferrer 配置文件

改进的费雷尔轮廓，其特征是核心接近平坦，外围形状迅速截断，最初由Ferrers（1877）提出，后来由Laurikainen等人（2005）修改。 强度定义如下：

$$I(r)=I_0{\left[1-{\left(\frac{r}{r_{\mathrm{out}}}\right)}^{2-b}\right]}^a,$$ (9)

其中 $I_0$ 也是中心表面亮度参数，$r_{\mathrm{out}}$ 是外截断半径，$b$ 和 $a$ 是控制斜率的参数分别是截断和核心。

幅度 ($m$, $m_0$) 和剖面模型参数之间的关系由下式给出

$$I_0=\frac{r_{\mathrm{box}}{10}^{-0.4(m-m_0)}}{\pi q{r}_{\mathrm{out}}^{2}a\beta (a,1+\frac{2}{2-b})},$$ (10)

其中 $q$ 也是轴比，$\beta$ 是方程 6 中定义的 beta 函数。 修改后的费雷尔剖面特别适合模拟星系中的杆结构（Blázquez-Calero等人，2020；Dalla Bontà等人，2018）。

4.2.3 经验（修改）King 配置文件

自从 King (1962) 首次给出说明以来，经验 King 轮廓已被广泛用于拟合银河系和河外球状星团 （BAOlab222https://github.com/soerenslarsen/baolab 提出Larsen，2014；Chies-Santos 等人，2007；Tripathi 等人，2023）。 该轮廓的强度定义如下：

(11)

这里，$I_0$也是中心表面亮度参数。 核心半径 $r_{\mathrm{c}}$ 表示密度开始偏离均匀性的尺度，而截断半径 $r_{\mathrm{t}}$ 则标记簇的边界。 全局幂律因子 $\alpha$ 决定了密度随着距中心距离的增加而下降的速率。 球状星团中恒星的浓度可以使用参数$r_{\mathrm{t}}$和$r_{\mathrm{c}}$来定义，用$c=\log \left(\frac{r_{\mathrm{t}}}{r_{\mathrm{c}}}\right)$表示。

图3： 面板比较了所选 $\sim$ 8,000 个 SDSS 星系的测量 Sérsic 剖面参数以及 $R^2$ 估计。 每个面板显示星系星等($m$)、位置角($\theta$)、轴比($q$)、有效半径($r_e$) ) 和 Sérsic 索引 ($n$)。 点的透明度说明了它们的浓度。 蓝线显示 galfit 结果等于 galmoss 结果。 顶部和右侧的直方图分别说明了 galmoss 和 galfit 测量值的分布。 比较显示所有面板具有适当的一致性。

4.2.4 高斯分布

在高斯分布中，强度定义如下：

$$I(r)=I_0\exp \left(-\frac{r^2}{2{\sigma }^2}\right),$$ (12)

其中$I_0$也是中心表面亮度参数，$\sigma$是径向色散。 在 galmoss 实现中，使用半高全宽 (FWHM = 2.354$\sigma$)。 高斯分布的经典应用包括对简单的 PSF 和点源进行建模。

4.2.5 莫法特简介

莫法特轮廓（莫法特，1969）通常用于模拟真实的望远镜PSF，定义其强度如下：

$$I(r)=I_0{\left[1+{\left(\frac{r}{r_{\mathrm{d}}}\right)}^2\right]}^{-n},$$ (13)

在哪里

$$r_{\mathrm{d}}=\frac{\mathrm{FWHM}}{2\sqrt{2^{\frac{1}{n}}-1}}.$$ (14)

$I_0$也是中心密度。 集中指数 $n$ 指示分布是更像洛伦兹分布 ($n$ =1) 还是更像高斯分布 ($n\to \mathrm{\infty }$)。

当考虑轴比$q$和箱形参数（方程6中的$r_{\mathrm{box}}$)的椭圆时，观测星等$m$之间的关系剖面模型参数为

$$I_0=\frac{r_{\mathrm{box}}(n-1){10}^{-0.4(m-m_0)}}{\pi q{r}_{\mathrm{d}}^2}.$$ (15)

与高斯分布类似，莫法特分布也可用于对点源建模。

4.2.6 平坦的天空

与其他径向光轮廓不同，在 galmoss 中，天空轮廓由所有没有径向矩阵的像素的天空平均值 $I_{\mathrm{sky}}$ 控制：

$$I=I_{\mathrm{sky}}.$$ (16)

在galmoss中，sky Profile中唯一需要的参数是$I_{\mathrm{sky}}$。 如果天空轮廓随位置（半径）变化，则可以将其包含为用户定义的轮廓。

5图像拟合和评估

生成模型图像后，galmoss 利用 python 和 pytorch 功能实现扩展的 $\chi 2$ 似然和梯度下降优化。 此外，还可以进行不确定性估计。

5.1 最大似然估计

galmoss 包采用 ${\chi }^2$ 可能性：

$\ln ℒ=-{\displaystyle \frac{k}{2}}\ln 2-\ln \mathrm{\Gamma }(k/2)+\left({\displaystyle \frac{k}{2}}-1\right)\ln {\chi }^2-{\displaystyle \frac{{\chi }^2}{2}},$ (17)

${\chi }^2={\displaystyle \sum _{i=1}^{N_i}}{\displaystyle \sum _{j=1}^{N_j}}{\mathrm{\Omega }}_{i.j}{\displaystyle \frac{{(D_{i,j}-M_{i,j})}^2}{{\sigma }_{i,j}^2}},$ (18)

其中$M_{i,j}$表示像素$i,j$处的模型，${\sigma }_{i,j}$表示$D_{i,j}$的不确定性。 在 ${\chi }^2$ 分布中，自由度 $k$ 是未包含在掩模中的具有星系通量的像素数与模型拟合参数数之间的差。

galmoss内部采用梯度下降进行并行拟合。 尽管收敛速度比其他传统方法（例如 Levenberg-Marquardt）慢，但它的计算和内存需求较低，可在 GPU 上实现高效的并行拟合。

5.2 置信区间

Galmoss采用两种策略来估计星系图像中的不确定性，主要使用基于高斯分布假设的协方差矩阵（Hogg等人，2010）。 参数不确定性是根据协方差矩阵的对角线计算的，为了提高计算效率，使用雅可比矩阵（例如，Gavin，2022）导出。 该方法将不确定性近似为 ${\sigma }_{\mathrm{p}}=\sqrt{\mathrm{diag}{\left[J^{T}WJ\right]}^{-1}}$，其中 $J$ 是雅可比矩阵，$W$ 是与 $W_{i,j}=1/{\sigma }_{i,j}^2$ 的对角线，提供 1-$\sigma$ 参数的置信水平。 不确定性估计的第二种方法是引导法，它涉及对每个星系图像进行多次替换重采样（通常是 100 次迭代）并重新拟合它们。 不确定性由 ${\sigma }^2=\frac{1}{M}{\sum }_{i=1}^M{\left[m_i-\overline{m}\right]}^2$ 给出，其中 $\overline{m}$ 是估计参数的平均值。 图4说明了星系随机子集的不确定性估计方法的比较。 我们显示通过引导分析确定的不确定性（红色误差条）和从协方差矩阵导出的不确定性（蓝色误差条）。 一般来说，自举不确定性涵盖的范围更广，表明估计更为保守。

图4： 使用两种方法对拟合星系参数的不确定性估计进行说明性比较：协方差矩阵（用蓝色误差条表示）和自举法（用红色误差条表示），用于随机选择样本中的星系。 通常，观察到从自举得出的不确定性大于从协方差矩阵获得的不确定性。

图5： galmoss 和 galfit 之间的加速性能评估。 每条线代表由 galfit 的运行时间归一化的加速时间。 结果表明，当批量大小略超过 10 时，galmoss 的速度超过 galfit，当批量大小为 1,000 时，速度快约 6 $\times$。 该基准测试是在具有以下规格的机器中执行的：CPU – 2.2GHz Intel Xeon Silver 4210； GPU – NVIDIA A100 (80GB)；操作系统 – Ubuntu Linux 18.04 64 位；内存 – 3.9 TB。

6结果

为了根据已建立的方法评估galmoss的性能，我们通过拟合来自同一数据集的单个星系剖面并将结构参数（例如Sérsic指数和有效半径）与使用以下方法获得的参数进行比较来进行验证测试galfit，如 MPP-VAC-DR17 中的目录所示。 这种比较充当基本的“健全性检查”，以验证 galmoss 拟合结果的可靠性。 图 3 显示关键星系参数的总体一致性 - 星等 ($m$)、位置角 ($\theta$)、轴比 ($q$ )、有效半径 ($r_{\mathrm{e}}$) 和 Sérsic 指数 ($n$) - 通过 $R^2$ 决定系数进行量化。

幅度、位置角和轴比的面板显示出很强的线性相关性，表明与 galfit 结果高度一致。 具体来说，在位置角面板中，显着偏离识别线的对象是那些具有大轴比($q\sim 1$)的对象。 这表明，在这些情况下，光分布对位置角 ($\theta$) 的变化相对不敏感。

众所周知，参数 $r_{\mathrm{e}}$ 和 $n$ 很难精确拟合（Trujillo 等人，2001），它们表现出线性关系，尽管与到其他参数。 值得注意的是，在 Sérsic 指数面板中，galmoss 值通常低于在 $n$ 值较高时从 galfit 导出的值。这种差异可能归因于图像质量的变化以及 pytorch 中使用的优化算法的特征。 为了深入分析这种观察到的偏差，我们引导读者阅读E部分。

除了上述测试之外，我们还使用相同的数据集和初始值使用 galfit 和 galmoss 进行了独立的速度比较。 使用 galmoss 拟合 8,289 个星系大约需要 10 分钟，比串行 galfit 快六倍。 图5展示了不同批量大小之间的效率增益，表明随着批量大小的增加，galmoss的速度优势变得更加明显，直至达到GPU容量的极限。

7结论

在这项研究中，我们引入了 galmoss，这是一个开源软件包，专门为在大型数据集上拟合星系光轮廓而设计，非常适合 LSST 时代。 galmoss 基于 torch 框架构建，为批量星系图像提供高效的 2D 表面亮度轮廓拟合，并具有 GPU 加速的额外优势。 它具有用户友好的界面，可以轻松定义参数及其范围。 此外，galmoss能够通过协方差矩阵和自举方法有效地量化不确定性。 该软件包还支持新配置文件的集成，使其成为星系结构研究中统计分析的适应性强的多功能工具。

基准测试表明，与 galfit 相比，galmoss 可以实现高达 6 倍的加速增益，具体取决于样本大小。 galmoss 的一个新颖之处是它使用本机并行化来处理单个星系的多个观测场作为单个多维张量的元素，从而通过高效的 PyTorch GPU 提高其速度 -加速矩阵计算和内存使用（Paszke等人，2019）。 这种方法与astropot形成鲜明对比，后者专注于具有多个星系或联合多波段拟合的较大区域。 galmoss 当前的限制是需要手动选择剖面模型和初始参数，以及其对于小星系批次的 GPU 利用率次优——未来版本需要改进的领域。

galmoss 可在 GitHub 上免费获取 ，泽诺多 并在 Python 包索引 readthedocs 包含示例用法以及包的概述。

致谢

我们感谢陈功宇为论文的逻辑和语言流程提供的帮助。 我们感谢沉孟亭帮助准备 SDSS 数据集。 我们感谢 Zhu Chen 和 Ziqi Ma 在学习 sextractor 的使用过程中提供的帮助。 我们感谢 Qiqi Wu 和 Shuairu Zhu 在 galmoss 应用中提供的帮助。

ACS 感谢巴西机构 Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) 和南里奥格兰德研究基金会 (FAPERGS) 通过赠款 CNPq-11153/2018-6 和 FAPERGS/CAPES 19/ 提供的资金2551-0000696-9。 SS感谢中国载人航天工程的研究资助 CMS-CSST-2021-A07，国家重点研发计划（No. 2022YFF0503402, 2019YFA0405501,），国家自然科学基金（No. 12073059 & 12141302）和上海市学术/技术研究带头人（22XD1404200）。

斯隆数字巡天 IV 的资金由阿尔弗雷德·P·斯隆基金会、美国能源部科学办公室和参与机构提供。 SDSS 感谢犹他大学高性能计算中心的支持和资源。 SDSS 网站是 www.sdss4.org。

附录A如何使用galmoss：单个Sérsic案例

在这里，我们演示如何使用 galmoss 包将单个 Sérsic 配置文件拟合到 SDSS 图像数据。

首先，我们需要加载必要的包。

⬇

1import Galmoss as gm

接下来，我们需要定义参数对象并将它们与配置文件实例关联。 星系参数的初始估计由sextractor提供。 值得注意的是，我们在这个简单的示例中没有包含 boxiness 参数，尽管它在 galmoss 框架中可用。

⬇

1# define parameter objects and profile

2

3sersic = gm.lp.Sersic(

4 cen_x=gm.p(65.43),

5 cen_y=gm.p(64.95),

6 pa=gm.p(-81.06, angle=True),

7 axis_r=gm.p(0.64),

8 eff_r=gm.p(7.58, pix_scale=0.396),

9 ser_n=gm.p(1.53, log=True),

10 mag=gm.p(17.68, M0=22.5)

11)

可以利用图像集（星系图像、西格玛图像、PSF图像、掩模图像）以及所选的轮廓来制定综合数据集对象。

⬇

1dataset = gm.Dataset(

2 galaxy_index="J162123.19+322056.4",

3 image_path="./J162123.19+322056.4_image.fits",

4 sigma_path="./J162123.19+322056.4_sigma.fits",

5 psf_path="./J162123.19+322056.4_psf.fits",

6 mask_path="./J162123.19+322056.4_mask.fits"

7 mask_index=2,

8 img_block_path="./test_repo",

9 result_path="./test_repo"

10)

11

12dataset.define_profiles(sersic=sersic)

拟合过程中初始化超参数后，就可以开始拟合了。 随后，我们运行不确定性估计过程。

⬇

1fitting = gm.Fitting(dataset=dataset,

2 batch_size=1,

3 iteration=1000)

4fitting.fit()

5fitting.uncertainty(method="covar_mat")

拟合过程完成后，拟合结果和img_blocks将保存在相应的路径中。

附录B如何使用galmoss：凸块+圆盘盒

图6： J100247.00+042559.8的g带galmoss分解结果，显示出良好的拟合质量。 左边六张图显示了总模型、圆盘模型、凸出模型、数据、模型残差（数据模型）和残差分布。 右侧显示一维投影。 上图用黑线表示模型分布，带有误差线的橙色点表示星系图像的通量分布及其不确定性。 红色虚线和蓝色虚线分别表示凸出部分和圆盘部分的分布。 下图表示带有误差条的蓝点的通量残留。

在这里，我们演示如何使用两个 Sérsic 配置文件的组合，使用 galmoss 包对 SDSS 图像数据进行圆盘和凸块分解。

⬇

1import Galmoss as gm

导入包后，后续步骤需要定义参数对象。 为了确保两个配置文件中的中心参数保持相同，只需指定一次中心参数，然后将其合并到各种配置文件中即可。

⬇

1xcen = gm.p(65.97)

2ycen = gm.p(65.30)

为了快速启动，我们让圆盘和凸出轮廓共享来自 sextractor 的初始值，凸出部分的初始 Sérsic 索引为 1，圆盘部分的初始 Sérsic 索引为 4。

⬇

1bulge = gm.lp.Sersic(

2 cen_x=xcen,

3 cen_y=ycen,

4 pa=gm.p(58.70, angle=True),

5 axis_r=gm.p(0.75),

6 eff_r=gm.p(4.09, pix_scale=0.396),

7 ser_n=gm.p(4),

8 mag=gm.p(17.97, M0=22.5)

9)

10

11disk = gm.lp.Sersic(

12 cen_x=xcen,

13 cen_y=ycen,

14 pa=gm.p(ini_value=58.70, angle=True),

15 axis_r=gm.p(0.75),

16 eff_r=gm.p(ini_value=4.09, pix_scale=0.396),

17 ser_n=gm.p(ini_value=1),

18 mag=gm.p(ini_value=17.97, M0=22.5)

19)

与单个配置文件的情况相比，我们只需要更改配置文件定义的代码片段。 这里我们选择使用bootstrap来计算不确定性。

⬇

1dataset = gm.Data_Box(

2 galaxy_index="J100247.00+042559.8"],

3 image_path="./J100247.00+042559.8_image.fits",

4 sigma_path="./J100247.00+042559.8_sigma.fits",

5 psf_path="./J100247.00+042559.8_psf.fits",

6 mask_path="./J100247.00+042559.8_mask.fits"

7 img_block_path="./test_repo",

8 result_path="./test_repo"

9)

10

11dataset.define_profiles(bulge=bulge, disk=disk)

12fitting = gm.Fitting(dataset=dataset,

13 batch_size=1,

14 iteration=1000)

15fitting.fit()

16fitting.uncertainty(method="bstrap")

图6显示了分解结果，残差显示了高质量的拟合。

附录C定义一组星系剖面的示例代码

在这里，我们提供一个示例，展示如何定义配置文件组合并从图像函数中提取模型图像以进行初步审查。 图7显示了一个模型图像，其特征是中央星系具有盘状、核球和条状，每个星系分别由两个Sérsic剖面和一个费雷尔剖面定义。 此外，还有一个由 Sérsic 剖面定义的侧星系、由 Moffat 剖面定义的两个点源以及由 King 剖面定义的疏散星团。 所有这些都与高斯分布生成的 PSF 图像进行卷积。

⬇

1import numpy as np

2import Galmoss as gm

3import torch_optimizer as optim

导入包后，第一个子模型代表中央星系，包括圆盘、凸起和条形。

⬇

1# The galaxy center

2xcen = gm.p(65.)

3ycen = gm.p(64.)

4

5# Define parameters and profiles

6Bisk = gm.lp.Sersic(

7 cen_x=xcen,

8 cen_y=ycen,

9 pa=gm.p(55., angle=True),

10 axis_r=gm.p(0.44),

11 eff_r=gm.p(13, pix_scale=0.396),

12 ser_n=gm.p(1),

13 mag=gm.p(17, M0=22.5),

14 box=gm.p(0.1)

15)

16

17Bulge = gm.lp.Sersic(

18 cen_x=xcen,

19 cen_y=ycen,

20 pa=gm.p(0., angle=True),

21 axis_r=gm.p(1),

22 eff_r=gm.p(6, pix_scale=0.396),

23 ser_n=gm.p(6),

24 mag=gm.p(21.5, M0=22.5)

25)

26

27Bar = gm.lp.King(

28 cen_x=xcen,

29 cen_y=ycen,

30 pa=gm.p(-40., angle=True),

31 axis_r=gm.p(0.2),

32 mag=gm.p(21.5, M0=22.5),

33 trunc_r=gm.p(5, pixScale=0.396),

34 trunc_a=gm.p(0.5),

35 trunc_b=gm.p(1.9),

36 box=gm.p(0.1)

37)

第二个子模型是侧星系。

⬇

1Sersic_s = gm.lp.Sersic(

2 cen_x=gm.p(45.43),

3 cen_y=gm.p(100.95),

4 pa=gm.p(-60., angle=True),

5 axis_r=gm.p(0.7),

6 eff_r=gm.p(13, pix_scale=0.396),

7 ser_n=gm.p(1.6),

8 mag=gm.p(17.5, M0=22.5),

9 box=gm.p(-0.6)

10)

第三个子模型是两个点源。 它们与理想化 PSF 图像共享色散 $\sigma$ = $1.515\mathrm{″}$，两者均由高斯分布产生。

⬇

1# Two point-souce

2P1 = gm.lp.Gaussian(

3 cen_x= gm.p(75),

4 cen_y= gm.p(80),

5 pa=gm.p(0., angle=True),

6 axis_r=gm.p(1),

7 inten=gm.p(0.1),

8 fwhm=gm.p(1.2,pix_scale=0.396),

9)

10

11P2 = gm.lp.Gaussian(

12 cen_x=gm.p(88.),

13 cen_y=gm.p(40.),

14 pa=gm.p(0., angle=True),

15 axis_r=gm.p(1),

16 inten=gm.p(0.1),

17 fwhm=gm.p(1.2, pix_scale=0.396),

18)

19

20# The idealized PSF

21PSF = gm.lp.Gaussian(

22 cen_x=gm.p(20),

23 cen_y=gm.p(20),

24 pa=gm.p(0., angle=True),

25 axis_r=gm.p(1),

26 inten=gm.p(1/(2*np.pi*((0.6/0.396)**2))),

27 fwhm=gm.p(1.2, pix_scale=0.396),

28)

没有卷积的总模型图像是三个子分量的总和（见图C.7(a))。 要提取模型图像，需要首先定义网格，因为网格生成仅在拟合对象内的拟合过程中自动发生。

⬇

1x = torch.linspace(0.5, 128-0.5, 128)

2y = torch.linspace(128-0.5, 0.5, 128)

3grid = torch.meshgrid(y, x, indexing=’ij’)

⬇

1bulge = Bulge.image_via_grid_from(grid, mode="initial_model")

2disk = Disk.image_via_grid_from(grid, mode="initial_model")

3bar = Bar.image_via_grid_from(grid, mode="initial_model")

4sersic_s = Sersic_s.image_via_grid_from(grid, mode="initial_model")

5p1 = P1.image_via_grid_from(grid, mode="initial_model")

6p2 = P2.image_via_grid_from(grid, mode="initial_model")

7result = bulge + disk + bar + sersic_s + p1 + p2

为了模拟视觉效果，将从相对较小的网格生成理想化的 PSF 图像，然后与模型图像（点源除外）进行卷积。 结果如图C.7(b)所示。

⬇

1#Grid generation

2x_psf = torch.linspace(0.5, 40-0.5, 40)

3y_psf = torch.linspace(40-0.5, 0.5, 40)

4grid_psf = torch.meshgrid(x_psf, y_psf, indexing=’ij’)

5

6# PSF image

7psf = PSF.image_via_grid_from(grid_psf, mode="initial_model")

8

9#Convolving

10result_conv = torch.squeeze(torch.nn.functional.conv2d((bulge+disk+bar+sersic_s).reshape(1,1, 128,128), psf.reshape(1, 1, 40, 40), padding="same", groups=1)) + p1 + p2

(a) The Model image without convolving a PSF image.

(b) The Model image with a PSF image convolved.

图7： 提供模型图像以供初步审查。 该模型由一个中心星系组成，该星系具有圆盘、核球和条，每个星系分别由两个 Sérsic 剖面和一个 Ferrer 剖面定义。 此外，还有一个由 Sérsic 剖面定义的侧星系、由高斯剖面定义的两个点源。

附录D用户定义的配置文件

在4.2.6节中，我们介绍了天空剖面。 在这里，我们展示如何集成更真实的天空轮廓作为用户定义轮廓的示例。 在这个新的天空剖面中，天空强度定义如下：

$$I=I_0+k_{{}_{\mathrm{x}}}(x-x_0)+k_{{}_{\mathrm{y}}}(y-y_0),$$ (19)

其中$(x_0,y_0)$是图像的几何中心，$I_0$是$(x_0,y_0)$处的天空值，$(k_{{}_{\mathrm{x}}},k_{{}_{\mathrm{y}}})$描述天空的变化沿 x 和 y 轴的值。

为了集成这个新的配置文件，我们需要定义一个新类，并将其命名为$NewSky$。 此类应包含一个指导根据给定参数生成图像模型的函数。 这里有一些注意事项。 首先，每个轮廓类都应该继承自基础类LightProfile，以访问通用的灯光轮廓功能。 接下来，应将参数加载到 __init__ 函数中的配置文件类中，并将其设置为属性（例如 self.sky_0 = sky_0）。

⬇

1import galmoss as gm

2

3class NewSky(gm.LightProfile):

4 def __init__(self, sky_0, grad_x, grad_y):

5 super().__init__()

6 self.psf = False

7 self.sky_0 = sky_0

8 self.grad_x = grad_x

9 self.grad_y = grad_y

10

11 def image_via_grid_from(self,

12 grid,

13 mode="updating_model"):

14 return (self.sky_0.value(mode)

15 + self.grad_x.value(mode)

16 * (grid[1]-(grid[0].shape[1] + 1)/2)

17 + self.grad_y.value(mode)

18 * (grid[0]-(grid[0].shape[0] + 1)/2))

轮廓方程在 $image\mathrm{\_}via\mathrm{\_}grid\mathrm{\_}from$ 函数中定义。 参数值在模式值之后提取，默认值为updating_model。 这种模式需要在整个拟合过程中不断更新的值，并且已经被广播为适合多维矩阵计算的形状。

新定义的配置文件可以按如下方式使用：

⬇

1sky = NewSky(sky_0=gm.p(0.3),

2 grad_x=gm.p(2),

3 grad_y=gm.p(3))

附录E高Sérsic指数端的拟合偏差

如图3所示，我们在比较结果中观察到高Sérsic指数($n$)端存在普遍偏差。 在galfit结果中$n$值较大（例如$n>6$）的星系，在galmoss结果中大多显示出较低的值（例如）。

为了进一步研究高$n$星系的类型，我们使用MDLM-VAC-DR17目录提供的形态分类。 我们将8,289个星系识别为2,849个ETG和5,440个LTG，并进一步将ETG识别为1,891个椭圆星系、793个S0星系和165个未定义星系。 我们发现 ETG 占整个星系样本的 34.37 $\%$，在偏向区域代表 88.78 $\%$。 为了进一步研究，我们绘制了 E 和 S0 中 $n$ 值的直方图（图 8)。 在图中，S0 星系中 galmoss 拟合结果和 galfit 拟合结果的直方图几乎一致。 然而，椭圆星系中galmoss拟合结果与galfit拟合结果的直方图分布有很大差异。 该结果表明，椭圆星系对拟合 $n$ 值的不一致贡献最大。

观察到的结果有几个原因。 其中之一是观察的局限性。 大多数椭圆星系都有一个陡峭的中央核心和一个平坦的外翼。 由于平翼的信噪比较差，${\chi }^2$对$n$的变化不敏感，在优化过程中起到了一个小梯度的作用。 另一个可能的原因与优化器的属性有关。 pytorch 中最流行的优化器使用自适应学习率，在低梯度条件下降低学习率以加快收敛速度​​。 因此，与galfit相比，$n$的灵敏度在高值端进一步降低，这增加了本例中两个软件之间结果的差异。

图8： 直方图显示了E和S0星系中$n$值的分布。 实线代表galmoss结果，虚线代表galfit结果。 紫色和粉色分别对应S0和E。 从galmoss和galfit拟合值的不一致可以看出，椭圆星系对$n$拟合结果高端不一致的影响最大。