AI探索引力波奥秘

Exploring the Secrets of Gravitational Waves with AI

王赫

2025年1月02日 @超算前沿直播间

hewang@ucas.ac.cn

中国科学院大学 · 国际理论物理中心（亚太地区）

中国科学院大学 · 引力波宇宙太极实验室（北京/杭州）

引力波天文学

什么是引力波

研究引力波的科学意义

蓬勃发展的引力波数据分析

引力波天文学

爱因斯坦于1916年提出广义相对论，并预言了引力波的存在
引力波是广义相对论中的一种强场效应
- 2015年：首次实验探测到双黑洞并合引力波
- 2017年：首个双中子星多信使探测，开启多信使天文学时代
- 2017年：引力波探测成果被授予诺贝尔物理学奖
- 至今：发现了超过 90 个引力波事件
未来：
- 2023-2025年：第四次引力波探测阶段，有希望探测更多不同类型的引力波事件
- 空间引力波探测计划 (LISA/Taiji/Tianqin) / XG (CE/ET)

双星并合系统产生的引力波波源

引力波振幅的测量

地面引力波探测器网络

2017 年诺贝尔物理学奖

多信使天文学

引力波探测打开了探索宇宙的新窗口
不同波源，频率跨越 20 个数量级，不同探测器

引力波研究的科学意义

基础物理学
- 引力子是否有质量, 引力波的传播速度 ...
天体物理学
- 大质量恒星演化模型, 恒星级双黑洞的形成机制 ...
宇宙学
- 哈勃常数的测量, 暗能量 ...

人类成功观测到引力波的五条关键要素：

良好的探测器技术

良好的波形模板

良好的数据分析方法和技术

多个独立探测器间的一致性观测

引力波天文学和电磁波天文学的一致性观测

DOI:10.1063/1.1629411

– 伯纳德·舒尔茨

引力波暂现源星表 (GWTC-3)

首次探测双黑洞并合引力波事件 GW150914

引力波数据处理

匹配滤波技术

困难与挑战

匹配滤波技术

Frequentist hypothesis testing and likelihood princple:
- make some assumptions about signal and noise hypothesis
- write down the likelihood function for a signal in noise
- find the parameters that maximise it
- define a corresponding detection statistic
  \(\rightarrow\) recover the MF
Bayesian hypothesis testing:
- start from the same likelihood
- define a prior over signal parameters
- marginalise over them to arrive at a Bayes factor
- Often the dirty secret: just treat this as a Frequentist detection statistic
  \(\rightarrow\) recover the MF (for certain prior choices)

脉冲响应函数：

h_w[t]

d_w[t]

\rho[t]

线性滤波器

输入序列

输出序列

h_w[t]

高斯且稳态噪声环境下，提取弱信号的最优线性算法
假设：引力波观测数据 = 高斯稳态噪声 + 某引力波信号
\(d(t) = n(t) + h(t)\)
匹配滤波信噪比 (SNR) \(\rho(t)\):

其中

噪声的(单边)功率谱密度 (PSD): \(S_n(f)\) 。
信号处理角度的理解：
- 若某一段时域数据流作为输入，探测统计量 (即匹配滤波信噪比) 是另一段输出的时序数据流，问怎样的线性滤波器可以使得输出结果最大？
Frequentist / Bayesian 角度的理解: 略

\rho^2(t)\equiv \frac{|\langle d|h \rangle(t)|^2} {\langle h|h \rangle}

\langle d|h \rangle (t) = 4\int^\infty_0\frac{\tilde{d}(f)\tilde{h}^*(f)}{S_n(f)}e^{2\pi ift}df

\langle h|h \rangle = 4\int^\infty_0\frac{\tilde{h}(f)\tilde{h}^*(f)}{S_n(f)}df , \, \tilde{h}_w(f) \equiv \frac{\tilde{h}(f)}{ \sqrt{S_n(f)}}

匹配滤波技术的困难与挑战

地面引力波探测科学数据的特点
- 噪声特点：非高斯 + 非稳态
- 信号特点：信噪比低 (约噪声幅度的1/100)
波形模板库的局限性
- 需要大量的精确波形模板以确保无遗漏，至少百万数量级
- 受限于已知引力理论预言的波形模板，难以搜寻超越经典理论的引力波信号
多信使天文学的兴起 + 引力波探测技术的进步
- 低(负)延迟的引力波信号搜寻
- 海量的累积数据和成批的引力波事件，有待高效的仔细分析

真实引力波数据的非高斯性

O1 观测运行时用的波形模板库

在 GW170817 事件后 1.74\(\pm\)0.05s 的伽玛暴 GRB 170817A

机器学习技术

计算智能

机器学习 vs 深度学习

人工智能技术应用于引力波数据处理

科学智能：AI for Science

2016年，AlphaGo 第一版发表在了 Nature 杂志上
2021年，AIphaFold 预测蛋白质结构登上 Science、Nature 年度技术突破
2022年，DeepMind团队通过游戏训练AI发现矩阵乘法算法问题
《达摩院2022十大科技趋势》将 AI for Science 列为重要趋势
- “人工智能成为科学家的新生产工具，催生科研新范式”
2023年，DeepMind发布AI工具GNoME (Nature)，成功预测220万种晶体结构
2023年3月，为贯彻落实国家《新一代人工智能发展规划》，科技部会同自然科学基金委启动“人工智能驱动的科学研究”（AI for Science）专项部署工作，布局“人工智能驱动的科学研究”前沿科技研发体系。
2024.4：美国总统科学技术顾问委员会（PCAST）发布《赋能研究：利用人工智能应对全球挑战》报告
2024.5: 《Science in the age of AI: How artificial intelligence is changing the nature and method of scientific research》 (Royal Soc.)

AlphaGo
围棋机器人

AlphaTensor
发现矩阵算法

AlphaFold
蛋白质结构预测

验证数学猜想

机器学习 vs 深度学习

机器学习：
- 线性回归模型、决策树模型、支撑向量机、马尔科夫链-蒙特卡洛方法 (MCMC) ...
深度学习：
- 用神经网络实现自动特征提取的模型
- 深度神经网络是一个万能的函数拟合器，可以表征任意复杂度的非线性函数映射
- 特点：端到端、数据驱动、过参数化 ...
传统引力波数据分析方法 ~ 传统机器学习方法

人工智能 > 机器学习 > 深度学习

人工智能

机器学习

深度学习

人工智能的一个分支。机器学习是用数据或以往的经验，以此优化计算机程序的性能标准

机器学习的一个分支。基于神经网络结构实现端到端的一种建模方法

任何能实现以人类智能相似的方式做出反应的技术

传统机器学习

深度学习

输入

特征提取

输入

特征

传统机器学习算法

输出

输入

自动特征提取 + 分类

输出

人工智能技术应用于引力波数据处理

深度卷积神经网络模型可以用来搜寻双黑洞并合系统所产生的引力波
- 输入：白化后的时序数据
- 噪声：高斯模拟噪声
- 灵敏度：与匹配滤波方法可比拟
- 执行速度：远胜过匹配滤波方法 (有GPU加持)

PRD, 2018, 97(4): 044039.
PRL, 2018, 120(14): 141103.

人工智能技术应用于引力波数据处理

The most common and direct approach, from Computer Vision (CV) to GW signal processing: pixel point \(\Rightarrow\) sampling point.
Artificial Intelligence (AI) has great potential to revolutionize gravitational wave astronomy by improving data analysis, modeling, and detector development.
Representation and supervised learning crucially extract features from GW signals, autonomously identifying informative features and leveraging labeled data for accuracy.

Exported: Oct, 2023

引力波信号探测

卷积神经网络

匹配滤波-卷积神经网络

卷积神经网络

许多研究者从不同的角度和问题出发，验证在高斯稳态噪声和模拟的理论波形下，卷积神经网络的优异性能。

GW150914

GW151012

GW151226

Krizhevsky, A., Sutskever, I., & Hinton, G. E. (2012)

Convolutional Neural Network (ConvNet or CNN)

将 CNN 模型在真实噪声和真实引力波事件上测试，效果很糟糕 😰

GPS 时间轴

匹配滤波-卷积神经网络

改进并开发神经网络模型，以适应真实的引力波观测数据的任务
匹配滤波算法当中的波形模板 \(\rightarrow\) 卷积层中的卷积核权重参数
匹配滤波感知层 (matched-filtering layer)
可以准确探测到 GWTC-1 中的 11 个真实引力波事件，甚至包括 GW170817
引力波信号处理 \(\rightarrow\) 智能引力波信号处理

GW151012

GW151226

GPS 时间轴

Real-time GW searches for GW150914

Matched-filtering Convolutional Neural Network (MFCNN)

GW150914

He Wang, et al. PRD 101, 10 (2020): 104003

匹配滤波 vs 机器学习

匹配滤波算法 vs 各类机器学习模型，到底孰优孰劣？
引力波信号搜寻算法基准测试 (MLGWSC-1)
PyCBC 和 cWB 是 LIGO-Virgo 主流的两套数据分析流水线
Dataset-4: 采样自 O3a 真实引力波观测数据

在模拟噪声数据上，机器学习算法和 LIGO 最灵敏的信号搜寻流水线相比，是非常有竞争力的
测试的大多数机器学习算法对非高斯的真实噪声背景过于敏感，误报率高
传统信号搜寻算法可以在低误报率水平下识别引力波信号，置信度有保障
测试的机器学习算法对长时信号的识别能力非常有限

http://arxiv.org/abs/2209.11146

引力波数据降噪

引力波数据质量

WaveFormer

引力波数据质量

数据质量的提升是一个非常复杂的问题，超过 20 万个传感器通道的数据会决定引力波科学数据通道的质量
降低引力波数据中非高斯的短时脉冲波干扰 (Glitch)，会有助于减少引力波信号误报率
引力波探测数据中去除 Glitch，是一个多分类问题
- 传统机器学习算法 Powell J, et al. CQG, 2015
- 深度学习算法 Zevin, M, et al. CQG, 2017; Razzano M, Cuoco E. CQG, 2018; Ormiston R, et al. PRR, 2020
与其消除数据的非高斯性，何不直接把信号重构出来？这有助于发现理论预言之外的引力波信号！

真实引力波数据的非高斯性

1400Ripples Air Compressor Blip

Extremely Loud Helix Koi Fish

部分 Glitch 的样例

LIGO-Virgo 数据处理流程

Waveformer

Billion-scale 大模型对全频段的引力波探测数据实现噪声去除和引力波信号的波形重构
Highlights:
- 抑制观测噪声幅度 (包括Glitch) 1-2 个数量级以上
- 对引力波信号幅度和相位信息精确重构

He Wang, et al. MLST. 5, 1 (2024): 015046.

["This", "is", "a", "sample"]

Waveformer

Billion-scale 大模型对全频段的引力波探测数据实现噪声去除和引力波信号的波形重构
Highlights:
- 抑制观测噪声幅度 (包括Glitch) 1-2 个数量级以上
- 对引力波信号幅度和相位信息精确重构

He Wang, et al. MLST. 5, 1 (2024): 015046.

Waveformer

Challenges in Model Interpretability
- The black-box nature of AI models complicates interpretability, challenging the comparison of AI-generated detection statistics with traditional matched filtering chi-square distributions.
- Convincing the scientific community of the pipeline's validity and the statistical significance of new discoveries remains difficult despite the model's ability to identify potential gravitational wave signals.

He Wang, et al. MLST. 5, 1 (2024): 015046.

OURs

LVK. PRD (2016). arXiv:1602.03839

GW151226

GW151012

引力波波源参数反演

贝叶斯统计推断

条件变分自编码器

归一化流模型

Credit: LIGO Magazine.

贝叶斯统计推断

Traditional parameter estimation (PE) techniques rely on Bayesian analysis methods (posteriors + evidence)
Computing the full 15-dimensional posterior distribution estimate is very time-consuming:
- Calculating likelihood function
- Template generation time-consuming
Machine learning algorithms are expected to speed up!

1400Ripples Air Compressor Blip

LIGO-Virgo 数据处理流程

Bayesian statistics

条件变分自编码器

深度生成模型：条件变分自编码器 (CVAE)
基于设计灵敏度的噪声功率谱，高斯模拟噪声
输出一个完整 15 维后验概率分布，约耗时 1s

1400Ripples Air Compressor Blip

耗时对比

一个样例：15 维参数的后验概率分布

归一化流模型

深度生成模型：归一化流模型 (Nflow)
基于 GW150914 附近噪声估计的噪声功率谱
首次实现对真实引力波事件 GW150914 的完整后验参数估计
50,000 后验样本耗时约 8 秒

Wang H, Cao Z, et al. Big Data Mining and Analytics, 2021

归一化流模型

深度生成模型：归一化流模型 (Nflow)
DINGO
- 测试 GWTC-1 的 BBH 事件
- 耗时 < 1 min (≈ 20 s, IMRPhenomPv2)
- 开始为 O4 部署，有望成为新的引力波信号搜寻流水线
DINGO-IS
- 测试 GWTC-3 中 42 BBH 事件
- 耗时 ≲ 1 h (IMRPhenomXPHM)，≈ 10 h (SEOBNRv4PHM, 64 CPU cores)
- 能够计算 evidence

http://arxiv.org/abs/2210.05686

归一化流模型（basic）

1400Ripples Air Compressor Blip

【【机器学习】白板推导系列(三十三) ～流模型(Flow based Model)】

p_{\mathrm{y}}(\mathbf{y})

p_{\mathrm{z}}(\mathbf{z})

\mathbf{z}

\mathbf{y}

T^{-1}

base density

target density

The main idea of flow-based modeling is to express \(\mathbf{y}\in\mathbb{R}^D\) as a transformation \(T\) of a real vector \(\mathbf{z}\in\mathbb{R}^D\) sampled from \(p_{\mathrm{z}}(\mathbf{z})\):

\mathbf{y}=T(\mathbf{z}) \quad \text { where } \quad \mathbf{z} \sim p_{\mathrm{y}}(\mathbf{z})

Note: The invertible and differentiable transformation \(T\) and the base distribution \(p_{\mathrm{z}}(\mathbf{z})\) can have parameters \(\{\boldsymbol{\phi}, \boldsymbol{\psi}\}\) of their own, i.e. \( T_{\phi} \) and \(p_{\mathrm{z},\boldsymbol{\psi}}(\mathbf{z})\).

Change of Variables:

p_{\mathrm{y}}(\mathbf{y})=p_{\mathrm{z}}(\mathbf{z})\left|\operatorname{det} J_{T}(\mathbf{z})\right|^{-1} \quad \text { where } \quad \mathbf{u}=T^{-1}(\mathbf{x}) .

J_{T}(\mathbf{z})=\left[\begin{array}{ccc} \frac{\partial T_{1}}{\partial \mathrm{z}_{1}} & \cdots & \frac{\partial T_{1}}{\partial \mathrm{z}_{D}} \\ \vdots & \ddots & \vdots \\ \frac{\partial T_{D}}{\partial \mathrm{z}_{1}} & \cdots & \frac{\partial T_{D}}{\partial \mathrm{z}_{D}} \end{array}\right]

Equivalently,

The Jacobia \(J_{T}(\mathbf{u})\) is the \(D \times D\) matrix of all partial derivatives of \(T\) given by:

p_{\mathrm{y}}(\mathbf{y})=p_{\mathrm{z}}\left(T^{-1}(\mathbf{y})\right)\left|\operatorname{det} J_{T^{-1}}(\mathbf{y})\right|

归一化流模型（basic）

p_{\mathrm{y}}(\mathbf{y})

p_{\mathrm{z}}(\mathbf{z})

\mathbf{z}

\mathbf{y}

T^{-1}

base density

target density

Rational Quadratic Neural Spline Flows
(RQ-NSF)

Data: target data \(\mathbf{y}\in\mathbb{R}^{15}\) (with condition data \(\mathbf{x}\)).
Task:
- Fitting a flow-based model \(p_{\mathrm{y}}(\mathbf{y} ; \boldsymbol{\theta})\) to a target distribution \(p_{\mathrm{y}}^{*}(\mathbf{y})\)
- by minimizing KL divergence with respect to the model’s parameters \(\boldsymbol{\theta}=\{\boldsymbol{\phi}, \boldsymbol{\psi}\}\),
- where \(\boldsymbol{\phi}\) are the parameters of \(T\) and \(\boldsymbol{\psi}\) are the parameters of \(p_{\mathrm{z}}(\mathbf{z})=\mathcal{N}(0,\mathbb{I})\).
Loss function:
Assuming we have a set of samples \(\left\{\mathbf{y}_{n}\right\}_{n=1}^{N}\sim p_{\mathrm{y}}^{*}(\mathbf{y})\),

Minimizing the above Monte Carlo approximation of the KL divergence is equivalent to fitting the flow-based model to the samples \(\left\{\mathbf{y}_{n}\right\}_{n=1}^{N}\) by maximum likelihood estimation.

\begin{aligned} \mathcal{L}(\boldsymbol{\theta}) &=D_{\mathrm{KL}}\left[p_{\mathrm{y}}^{*}(\mathbf{y}) \| p_{\mathrm{y}}(\mathbf{y} ; \boldsymbol{\theta})\right] \\ &=-\mathbb{E}_{p_{\mathbf{y}}^{*}(\mathbf{y})}\left[\log p_{\mathbf{y}}(\mathbf{y} ; \boldsymbol{\theta})\right]+\text { const. } \\ &=-\mathbb{E}_{p_{\mathbf{y}}^{*}(\mathbf{y})}\left[\log p_{\mathrm{z}}\left(T^{-1}(\mathbf{y} ; \boldsymbol{\phi}) ; \boldsymbol{\psi}\right)+\log \left|\operatorname{det} J_{T^{-1}}(\mathbf{y} ; \boldsymbol{\phi})\right|\right]+\mathrm{const} . \end{aligned}

\mathbb{E}_{p_{\mathbf{y}}^{*}(\mathbf{y})}\left[\log p_{\mathbf{y}}^{*}(\mathbf{y} ; \boldsymbol{\theta})\right]

\mathcal{L}(\boldsymbol{\theta}) \approx-\frac{1}{N} \sum_{n=1}^{N} \log p_{\mathrm{z}}\left(T^{-1}\left(\mathbf{y}_{n} ; \boldsymbol{\phi}\right) ; \boldsymbol{\psi}\right)+\log \left|\operatorname{det} J_{T^{-1}}\left(\mathbf{y}_{n} ; \boldsymbol{\phi}\right)\right|+\mathrm{const.}