优化∣浅析黑洞成像背后的稀疏重构算法原理

优化∣浅析黑洞成像背后的稀疏重构算法原理

编者按:本文讨论了最近在网络上得到热议的黑洞成像背后的稀疏重构算法原理,涉及到优化、机器学习和统计理论相关的知识。内容整理自Katie Bouman博士(现为加州理工学院教授)读博期间的一系列工作。

文章作者: @覃含章
责任编辑:@覃含章
文章发表于微信公众号【运筹OR帷幄】:优化∣浅析黑洞成像背后的稀疏重构算法原理
欢迎原链接转发,转载请私信@运筹OR帷幄获取信息,盗版必究。
敬请关注和扩散本专栏及同名公众号,会邀请全球知名学者发布运筹学、人工智能中优化理论等相关干货、知乎Live及行业动态:『运筹OR帷幄』大数据人工智能时代的运筹学

今天才关注到黑洞照片出炉这个事件,作为一个(伪)科幻迷心情激动,不过我注意到网络上绝大多数内容还是在讨论这张照片在物理学上的意义(如用来验证爱因斯坦的相对论,以及黑洞观测的历史等等)。

我个人呢,也是在了解了一些细节之后发现黑洞照片成像其实除了来自全世界实验物理学家的努力,其实也涉及到了大量(稀疏)图像处理、统计理论、非凸优化算法、混合高斯模型、隐马尔科夫模型、贝叶斯推断的问题。通过查阅资料,我发现,给黑洞拍照其实对统计学、计算机科学和运筹学/优化理论这些领域实际上也带来了很多有意思的、具有挑战性的问题。

以及,21世纪的科学边界的发展应当是跨学科的。

这里主要想给大家介绍一位叫做Katie Bouman的小姐姐关于黑洞照片的图像重构(reconstruction)方面的一系列工作。Katie小姐姐17年获得了MIT CS PhD学位,博士期间的主要工作便是“基于物理模型反问题的极端成像:从角落里一窥黑洞”(也是其博士论文的题目,Extreme Imaging via Physical Model Inversion: Seeing Around Corners and Imaging Black Holes),现在在Caltech的CMS系当助理教授。

本文基于Katie小姐姐个人主页上的一系列公开资料,包括她关于黑洞成像的Ted talk,博士答辩视频,她制作的精美的PPT和她的博士论文。

Katie是在一开始对天文物理学零了解的情况下开始做这个项目的,很传奇的经历

黑洞成像中的稀疏重构问题

本节我们先讨论为什么给黑洞拍照这件事情要牵扯到很多数学问题呢?如果你对天文物理学和人类目前的天文物理技术设备现状完全不了解,可能这是你的第一个问题:“拍照”不就是对着目标“咔擦”一声就完事了?而这次的黑洞照片,从2017年4月就开始拍摄,过了约摸两年才能公布出照片,科学家们是不是太拖沓了一点?

你可能会觉得,给黑洞拍照不就是手机“咔擦”一声的事情...

这个想法主要的问题是,黑洞离我们地球都很远很远。根据报道,这次的黑洞照片的目标位置在M87星系中心离地球大约5500万光年处。因此,观测难度相当于人在地球,却要在月球表面寻找一个橘子。然而,依据人类现有的望远镜的分辨率是万万达不到的,如见下图,这是一张有1w多个像素点的月球表面局部照片,而每个像素点的大小可以放得下100多万个橘子,在这种精度下想直接获得黑洞的照片(假设它看的到)无疑于痴人说梦。。。而且,以目前的航天实力我们也没有办法派航天飞船飞到黑洞附近去(靠的太近了也会直接被黑洞引力扯碎)。。。

这是目前分辨率最高的月球表面照片之一(通过不断放大局部得到)

实际上,我们目前的射电望远镜的分辨率和它的体积/大小差不多是成比例增长的,要达到上述的分辨率,容易经过计算以我们现有的技术需要至少跟地球那么大的望远镜才行。。

嗯,差不多这么大的望远镜就可以直接观测黑洞了...

那么假如我们有了地球这么大的望远镜,我们是不是能看见黑洞呢?这次终于有个好消息,我们可以看到!因为虽然黑洞本体连光都不放过,但是我们可以来观察黑洞的“事件视界”(Event Horizon)。也就是一个在黑洞附近的物质盘,其中的气体剧烈运动而发出强烈的辐射,这也正是我们常常在科幻片中所能看到的那围绕着黑洞轮廓的碟形光环。

好了,接下来的问题就是,我们现在并没有那么大的望远镜,但能否通过将一系列的小望远镜组合起来,“近似”地看到黑洞呢?答案就是,这便是本次黑洞“摄猎”团队所做的事情。观测事件视界的望远镜就叫做Event Horizon Telescope(EHT),因此团队也就叫做EHT团队。

全世界各地的一些EHT望远镜


因为虽然我们并没有那么大的一台望远镜可以给我们准确的图像(稠密的信号),现实当中我们只能利用全世界团队们手中的那几台望远镜得到一些稀疏的信号。因此,其实最终我们是希望通过这些稀疏的信号去反推(重构)那个完整的图像(见下图)。这是什么意思呢?我们现在看到的这张黑洞照片是实际上通过统计模型推断(基于“无线电之光”稀疏成像)出来的,以极高概率和原始黑洞照片非常非常像的照片,但并不是什么实际直接拍到的照片!

左边:如果我们有一台覆盖全地球的射电望远镜(得到稠密图像,因为我们可以得到各个频率分布下需要的光) 右边:现实...(只能得到“稀疏”图像)

由此,在观测完毕之后(射电望远镜们得到的事件视界的波长数据),我们就知道我们的拍照问题变成了一个算法和数学问题。海量的数据从观测团队那里被汇入了位于MIT和马普研究所的实验室内,开始进行统计/优化模型和算法的训练和重构。这之中,科学家们需要对海量数据进行处理(有一个细节,据说5天的观测就累积了差不多3500TB的数据,最后都是装到硬盘里物理运输到数据中心的...),并不断地修正提升自己的统计推断模型和优化算法(数据处理和算法实现都在超算集群上完成),并且两个实验室之间还进行了严格的结果比对,这才是为什么从最初拍摄到我们能看到花了那么久的时间(所谓的“冲洗”两年)...

当然,实际重构问题并没有前面那张图那么简单,因为地球无时不刻在运动(自转+公转),所以我们对观测数据还需要一些额外处理...但这不见得是坏事,因为地球的转动让我们得到了更多的信号(从不同的自转/公转角度)


基于VLBI的静态源稀疏图像重构

本节我们先不考虑黑洞图像可能会随着时间迁移改变,我们指出上节的图像重构问题可以简单看成下图的形式,即我们在数学上可以简单看成我们所需要的黑洞照片的像素点 x 和视界事件传来的无线电光(信号/测量值)y 呈函数关系 f(x)=y,那么如果是我们有稠密的信号源(地球般大的射电望远镜),实际上我们只需要取反函数,即 f (y),就可以直接得到我们想要的图像(像素点阵) 。

图像重构是一个反问题(inverse problem)


当然,主要的问题也如前所述,我们并没有稠密的信号源,所以我们只能基于稀疏信号对 x 进行重构,得到一个近似的照片(像素点阵)x 这里我们先提一句传统方法CLEAN,一种基于逆傅里叶变换的贪心重构算法。这其实就是一个所谓的反卷积(deconvolution)算法。可惜的是,这些传统方法在黑洞成像里没那么管用,最主要的问题就是由于无线电光到达地球上的不同EHT望远镜的时间根据地球公转/自转的不同阶段和望远镜自身的方位会有不同程度的延时,具体的延时长度可以认为是随机的...

两个EHT望远镜收到信号的时间会有延时,具体的延时长度可以看成是随机(random)的...

因此,Katie等人的研究小组提出使用贝叶斯模型来解决这个问题。为了阐述清楚他们的模型,我们先考虑最简单的高斯统计模型如下。

注意到这里 x 所代表的是黑洞图像的先验(prior)分布,我们可以用一些其它黑洞/星系/甚至日常图片的仿真/近似/拼图图片集 {Zn:n=1: ...,N} 来训练先验分布的参数:(样本均值和经验协方差矩阵)

这样我们就完成了对x先验分布的训练,注意到和经典模型不同的是,y 我们现在认为是一个均值为 f(x) 的随机变量。那么根据贝叶斯公式,我们知道图像 的后验分布满足

于是,我们的重构图像 就可以通过最大化log后验似然函数(最小化-log后验似然函数)来得到:

就是这么简单!任何有过工科统计学本科背景的同学都应该很熟悉这个模型,但谁会想到这就会是第一张黑洞图像重构算法的基础呢?

注意,上式也可以写成常见的一般形式:

当然,实际当中直接用这个高斯模型是不行的,因为它对图像的拟合能力有限。在黑洞图像重构中,Katie小组使用的是混合高斯模型(Gaussian Mixture Model, GMM)。

具体来说,我们现在认为 Zn 来自于 C 个聚类组成的一个混合高斯分布:

同样,只要有本科统计/机器学习课程的基础,我们就知道可以用著名的EM(Expectation-Maximization)算法对这类GMM模型进行训练。而这也就是Katie小组所使用的方法!

贝叶斯框架下的反问题


三基于VLBI的动态源稀疏图像重构

如前所述,之前的静态模型并没有考虑黑洞的观测图像其实应该是随着时间不断变化的,也就是说之前的高斯模型其实从准确性的角度来说更应该建模成动态模型

我们可以利用下图来简单理解这个动态成像模型。在一个离散的时间区间 t 中,我们观察到基于图像 Xt 的带随机误差的信号yt 。我们认为 Xt 和相邻的图片 Xt-1,Xt+1 是相似的(同一分钟内黑洞图像的变化不会很大),然后图像在整体时间轴上的动态变化就由矩阵A 来刻画。

动态成像模型图示

这边Katie小组选择的 A 的形式如下:

这里我们不仔细阐述这个矩阵的物理含义和里面各种符号的定义,只是说明这个矩阵可以通过调整参数来刻画如旋转、断裂、线性变形、等比例缩放等各种基本的运动。

同样的,学过本科统计/机器学习课程的同学马上又知道了,如果我们先认为 事先已知,这类随时间变化的混合高斯模型其实可以看成所谓的高斯隐马尔科夫模型(Gaussian Hidden Markov Model)。这里具体细节我们略去不表。

backward更新步骤

然后我们就知道可以用EM算法来训练我们的模型了!对具体推导和细节感兴趣的同学,建议去看Katie的文章。而这个模型,就是黑洞图像背后的算法原理。


四最后的一点碎碎念

毫无疑问,黑洞照片的问世是激动人心的。但是对我个人而言,不仅是看到照片本身让我激动,更让我感到激动的是这背后的数学原理和算法实现。实际上,黑洞成像的成功不仅仅让物理学家们激动,对于做图像处理、统计、计算机科学、优化算法等各个方向的研究者们来说可能是更加激动的一件事情。

因为它其实引出了更多的问题。

从我个人的角度,最让我兴奋的是Katie没有去使用一个黑箱子(blackbox)模型,比如现在大火的各种深度学习里面的卷积神经网络等等(当然这也是目前深度学习研究者们趋之若鹜的一个方向)。相反,Katie采用的是一个更加传统的统计模型对图像进行建模。这个高斯隐马尔科夫模型虽然传统,但其实对这个模型的训练/算法设计还远没有做到很好。这也是Katie本人在博士答辩的时候指出的一个未来的方向(见下图)。

截自Katie博士答辩的PPT

原因当然有很多,最显而易见的一点就在于优化问题本身高度非凸(EM算法只是一个针对该优化问题的一个heuristics,只能保证找到一个局部最优,但对解离全局最优的质量基本没有保证)。因此,这类GMM模型直到今天在算法理论方面仍遗留有大量还未解决的问题。

与此同时,这种和时间序列相关的动态稀疏成像问题本身也值得更多的关于数据本身的研究。本质来说,我们到底如何通过有随机噪声的稀疏信号去进行图像重构,Katie等人的高斯隐马尔科夫模型当然是一种成功的思路,但是否针对不同的问题我们还可以设计出更好的统计模型呢?还有,针对这种超大规模的图像重构问题是否可以有更好的训练算法(基于超算集群)?

另外,即使是这个高斯隐马尔科夫模型,如何选择先验分布?如何选择EM算法的初始点?如何选择正则项?等等等等,都需要一个更好的理论指导呢。

黑洞的照片已经拍出来了,但是这些问题可能才是作为非物理学家/算法工程师的我们更应该去思考的。


参考文献:

Bouman, Katherine Louise. Extreme imaging via physical model inversion: seeing around corners and imaging black holes. Diss. Massachusetts Institute of Technology, 2017.


相关文章推荐

文为知乎问题的精华回答选编,回答者均是知乎上该问题涉及到的领域里非常活跃和具有一定影响力的用户。回答从各个角度优化理论能给深度学习带来的方方面面的影响。

《优化理论能给深度学习带来怎样的革命?


温馨提示

可以在 本公众号后台 回复关键词:“黑洞成像”获取Katie小姐姐的Ted演讲视频(中文字幕),如果觉得有用, 请勿吝啬你的留言和赞哦!~


板块招聘信息:

【优化】版块目前招募副主编,要求:

1. 国内外运筹学、优化理论等相关专业博士在读或以上

2. 有时间有热情有能力去分享自己的想法,做一些运筹学或者人工智能的科普工作

3. 对优化理论其中任何一个领域比较熟悉即可

4. 有工作热情,每周能有时间付出

5. 具有创造力,能够有独立创作专题系列文章的能力


同时也欢迎相关的背景的同学加入【优化】板块担任责任编辑、原创文章/电子书系列志愿者~

关于我们:『运筹OR帷幄』团队掠影

请将简历发送至:operations_r@163.com

欢迎加入我们这个大家庭!


扫二维码关注『运筹OR帷幄』公众号:

点击查看『运筹OR帷幄』志愿者招募介绍及加入方式

编辑于 2019-04-11

文章被以下专栏收录