一种利用自动机和小波分析解决Monge–Kantorovich问题的方法 Armando Sánchez-Nungaray Marcelo Pérez-Medel Carlos González-Flores Raquiel R. López-Martínez Martín Solís-Pérez

《Mathematical and Computational Applications》：A Method for Solving the Monge–Kantorovich Problem Using an Automaton and Wavelet Analysis Armando Sánchez-Nungaray, Marcelo Pérez-Medel, Carlos González-Flores, Raquiel R. López-Martínez and Martín Solís-Pérez

　　摘要：本文介绍了一种自动机，旨在改进Monge–Kantorovich（MK）问题的可行解，特别是在成本函数连续的情况下效果更为显著。为了提升其性能，使用离散小波变换获得了一个良好的初始解。具体而言，解决了一个运输问题，其中成本矩阵由变换的近似系数组成，将变量数量减少到原始离散问题的四分之一。通过使用细节系数，扩展了这个简化问题的解，从而得到了原始问题的一个可行解。这个解作为调整自动机的初始状态，自动机的最终状态提供了运输问题最优解的近似值。

1. 引言
Monge–Kantorovich（MK）问题被认为是最优运输理论的基础公式之一，为包括图像配准、质量重新分配问题、测度的数值逼近以及现代机器学习模型在内的广泛应用提供了理论框架[1,2,3]。高效解决这个问题需要近似在高维空间上定义的测度和成本函数，这带来了重大的计算挑战，尤其是在使用精细离散化或寻求逐步精细化的解决方案时[4,5,6,7]。
Monge–Kantorovich问题的表述如下。设X和Y是两个紧子集。我们用μ表示定义在X和Y上的有限测度族。对于每个-可测集合A和B，其在X和Y上的边际测度分别由(1)和(2)定义。设c是一个定义在X和Y上的实值函数，μ?和μ?分别是X和Y上的有限测度。那么Monge–Kantorovich质量传输问题由(3)给出，这构成了一个在具有预定边际的正测度空间上的线性优化问题[5,6]。

基于小波的多分辨率分析（MRA），特别是Haar MRA，为构建MK问题的分层近似提供了一个自然的框架。在这种方法中，涉及的测度和成本函数被投影到嵌套的近似空间上，产生了一系列分辨率逐渐增大的离散问题。这种多尺度结构显著降低了问题的计算复杂性，同时保留了其基本特征[8,9,10,11]。然而，即使在这种框架内，在精细分辨率下获得最优解仍然可能需要解决大规模的运输问题。
基于小波分解的多分辨率方案已经成功应用于Monge–Kantorovich问题。特别是，基于Haar的近似被用来构建MK问题的离散版本，并为其解决方案推导出高效的数值方案[12,13,14]。这些方法利用小波变换的分层特性在粗略层次获得良好的初始近似，然后随着分辨率的提高进行精细化。

为了进一步提高计算效率，我们提出用一种专门设计用于操作排列的自动机来补充基于小波的近似。该自动机利用局部成本信息来指导可行解之间的转换，确保每次转换都保持可行性，同时改进或维持目标函数的价值。自动机的底层结构与经典的分配和运输问题密切相关[15,16,17]，其行为可以被解释为一种元启发式精细化策略[4]。
在我们的方法中，多分辨率分析提供了一个在每个分辨率层次上接近最优的初始可行解。具体来说，我们将离散小波变换应用于成本函数，并解决相应的运输问题，其成本矩阵由近似系数组成。这个简化问题只涉及原始离散问题变量的四分之一，从而大幅降低了计算成本。类似策略在之前的研究中被证明是有效的[12,13,14]。
在近似层次上获得的运输问题的解随后通过结合小波分解的细节系数扩展为原始问题的可行解。这个扩展解作为调整自动机的初始状态。自动机达到的最终状态被视为Monge–Kantorovich问题最优解的近似值。这种近似程序的收敛性和稳定性得到了关于近似方案和测度收敛性的经典结果的支持[5,6,18]。

本文的内容如下：第2节，我们介绍了一系列关于多分辨率分析（MRA）和小波变换的初步结果，这些结果将用于获得Monge–Kantorovich问题的近似值。第3节，我们介绍了一个用于排列群的调整自动机。这一工具在本文的发展中将至关重要，因为它将使我们能够改善运输问题的解决方案。第4节，我们介绍了Monge–Kantorovich问题的离散化，指出可行解是排列，并定义了一个在自动机每个状态（排列和位置）上的函数，这将确保下一个状态优于或至少等于当前状态。然后我们通过从随机排列开始的一些示例来说明自动机的工作原理。第5节，我们给出了几个经典的Monge–Kantorovich问题示例，并基于使用小波变换和排列调整自动子在多分辨率分析的每个层次上解决运输问题的方法来给出解决方案的近似值。

2. 基础知识
本节介绍了所提出方法论背后的理论基础。特别是，它介绍了基于Haar小波的多分辨率分析（MRA）及其在二维情况下的扩展，以及通过多尺度近似进行测度投影的属性。
在实数线上定义的多分辨率分析（MRA）被定义为函数空间中的一组子空间序列。这些子空间满足以下属性：
(1) ... 对于每个i，|f_i(x)| ≤ |f_j(x)| ≤ L。
(2) ... 对于所有i ≠ j，f_i(x) · f_j(x) = 0。
(3) ... 存在一个函数φ，称为缩放函数，使得{φ(x - a)/λ_i, φ(x + a)/λ_i}构成一组正交的平移函数。
(4) ... 对于所有i，φ(x)是连续的。
(5) ... 存在一个函数g，使得{φ(x - a)/λ_i, φ(x + a)/λ_i}构成一个关于[0,1]的正交系。
存在一个函数σ，使得{φ(x - λ_i·x, φ(x + λ_i·x)}构成一组正交的平移和缩放函数。
设f_i和f_j是定义在[0,1]上的函数。令A是[0,1]的某个子集。对于每对i, j，我们定义函数g_ij：g_ij(x) = φ(x - λ_i·x - σ·(x - a), φ(x + λ_i·x + σ·a)。因此，基于Haar缩放函数生成的子空间序列构成了Haar多分辨率分析（Haar-MRA）。此外，集合{φ(x - a)/λ_i, φ(x + a)/λ_i}构成了[0,1]的一个基，对于每个i, j，空间X_i可以表示为X_i = φ(x - a/λ_i, x + λ_i/λ_j)的直和。

基于小波的多分辨率分析（MRA）在实数线上可以扩展到二维情况。在这种扩展中，小波函数被定义为定义在[0,1]上的Haar小波函数，这些函数可以用于构造MK问题的分层近似。这样，问题的计算复杂性大大降低，同时保留了其基本特征。

3. 排列上的调整自动机
本节介绍了一种在排列群上定义的有限自动机，旨在细化Monge–Kantorovich（MK）问题的离散可行解。该自动机在后续章节中作为一个核心工具，它能够进行系统的局部调整，从而减少相关成本并提高解决方案的质量。
自动机是一种计算的理论模型，用于描述和分析根据预定状态集和转换规则处理输入并产生输出的系统。正式地，自动机定义如下：
定义1. 自动机被定义为一个5元组：(Q, Σ, Δ, φ, ω)，其中：
Q是有限的状态集；
Σ是称为字母表的有限输入符号集；
Δ是转换函数；
φ是状态转换规则；
ω是初始状态；
F是接受（或最终）状态集。
组件的解释：
1. 状态（Q）：在任何给定时间，自动机占据有限的状态之一。这些状态编码了系统的内部配置，可以解释为它的“内存”。
2. 字母表（Σ）：输入字母表是自动机能够读取的有限符号集。例如，在二元输入的情况下，字母表可以是{0, 1}。
3. 转换函数（Δ）：转换函数指定了自动机如何响应输入符号从一个状态转移到另一个状态。正式地，它是一个映射f: Σ → Q。
4. 初始状态（ω）：计算从特定的初始状态开始，自动机从该状态开始处理输入。
5. 接受状态（F）：子集F被指定为接受（或最终）状态集。如果整个输入被处理后，自动机终止在属于F的状态中，则称输入字符串被自动机接受。

有关自动机的更多信息，请参阅[19,20]。现在我们回顾以下概念：
定义2. 设σ_n是在n上的一个排列。n上的所有排列以及函数组合构成一个称为n次对称群的群，用S_n表示。
一个排列可以表示为一个长度为n的数组（或列表），其中第i个元素表示排列在σ_n下的像。排列的逆排列σ(Collections)n是这样一个排列：σ.reverse(n) = σ^(n-1)。
可以通过交换数组中的索引和相应的值来计算排列的逆排列。关于排列的详细讨论可以参见[21]。

示例1。设σ = (1, 2, 3)。那么它的逆排列σ^[?1] = (3, 2, 1)。
示例2。设σ = (0, 1, 2)。它的矩阵表示为[0, 1, 2]。
注意，对于每个固定的排列σ和固定的索引k，与位置k相关的元素构成数组中的相应位置。第二行表示排列在数组中的位置，k表示自动机作用的位置。

示例3。为了示例上述段落，我们考虑以下关于位置的排列。逆排列σ^[?1] = (3, 2, 1)。现在，计算相对于σ的位置。最后，矩阵表示为[1, 2, 3]。在第5节中呈现的图形中，表示排列（可行解）的图表是从下到上、从左到右编号的。此外，1的位置用有色点表示。相比之下，在之前的矩阵表示中，编号是从上到下、从右到左进行的。例如，考虑到第5节中使用的行重排序，之前的矩阵对应于以下排列。最后，使用有色点表示排列的方式如下：邻居的解释如下。在矩阵表示中，位置的邻居是指在行和列中包含值1的位置，因为是固定的，我们将位置简单地称为位置k。具体来说，这些位置是：为了简单起见，我们将用来表示位置k的邻居，如图2所示。图2. 排列矩阵中位置k的邻居。定义3. 我们按照以下方式从元素构造了一个调整后的自动机：我们考虑5元组(11)，其中：是一个有限的状态集，是一个称为字母表的有限输入符号集，是转移函数，由定义，其中是初始状态，是的任意元素。是接受（或最终）状态的集合。4. 通过排列群上的调整自动机改进Monge–Kantorovich问题的解在本节中，我们提出了一种混合方法，将基于Haar小波的多分辨率分析（MRA）与作用于排列的有限自动机结合起来，旨在获得Monge–Kantorovich（MK）问题的最优解。这种方法利用MRA的层次结构构造出一个初始的近似解，然后通过前一节中描述的自动机进行局部调整来 refined 这个解。如果度量满足(3)并且配对是有限的，则称该度量为问题的一个可行解。如果存在一个达到目标函数最优值的可行度量，则称问题是可解的。在这种情况下，被称为(3)的一个最优解。此外，假设相对于Lebesgue度量在上是绝对连续的，并且和相对于Lebesgue度量在上是绝对连续的，则问题可以通过在上的Haar多分辨率分析进行自然离散化。回想一下，根据方程(3)，在j层上的Monge–Kantorovich问题由给出：对于每个-可测集和，其中，和是将度量、和分别投影到j层的Haar MRA的结果。特别是，我们考虑具有成本函数、基础集和边际度量的问题。由于实际应用通常涉及离散化公式，因此在上的多分辨率分析应用于我们得到以下目标问题：在这里，表示位于x轴上区间内的初始质量部分，该部分被分配给y轴上的区间。我们将j离散的单位正方形称为由正方形形成的网格（见(9)），它将集合划分为个块，每个块自然地与点相对应。上述公式是一个分配问题，是运输问题的一个特例。问题实例有一定数量的代理和一定数量的任务。任何代理都可以被分配去执行任何任务，可能会根据代理-任务分配产生一些成本。匈牙利方法是一种组合优化算法，它可以在多项式时间内解决分配问题，并且先于后来的原始-对偶方法被提出。它由Harold Kuhn在1955年开发并发布，他将其命名为“匈牙利方法”，因为该算法在很大程度上基于两位匈牙利数学家Dénes K?nig和Jen? Egerváry的早期工作。更多细节见[15,16,17]。我们假设存在一个简单的解决方案(12)。也就是说，是一个可行解，对于任何，它必然满足对于所有，以及对于所有。所有(12)的可行解都对应于排列。因此，我们使用第3节中介绍的自动机，在这个设置中。我们的目标是改进给定的(12)解决方案，其中表示相关的成本矩阵。基于这个成本矩阵和自动机的当前状态，我们引入以下量。定义4. 我们定义与矩阵相关联的状态的评估函数如下：评估函数的解释如下：如果，则任何改变都会改善可行解。如果，则在中的位置和k的改变是改善可行解的更好局部改变。如果，则在位置和k的改变是改善可行解的更好局部改变。如果，则在位置和k的改变是改善可行解的更好局部改变。如果，则在位置和k的改变是改善可行解的更好局部改变。根据构造，如果我们考虑一个状态并应用带有字母表符号的自动机，我们得到以下结果。显然，和都是(12)的可行解。此外，根据在(14)中定义的函数f的构造，我们有等价于其中和分别表示与和相关的矩阵的条目。总之，我们得到以下结果。定理1. 考虑由(12)给出的运输问题，并考虑一个状态，其中是一个排列和运输问题的可行解。设f是一个根据成本矩阵(14)评估状态的函数。如果，则等价于给定自动机的初始状态，我们应用自动机来改进可行解如下。因此，我们发现自动机的最终状态由给出。因此，自动机的最终状态提供了比初始状态更好的解决方案。示例4. 图3展示了将调整自动机应用于成本函数在3层离散化的Monge–Kantorovich问题，产生了大小为的 cost矩阵。图3. 优化过程中候选解的演变。背景颜色代表成本矩阵的值。绿色点表示当前正在评估的交换中涉及的分配，而黄色点对应于在该步骤中保持不变的分配。图像序列展示了自动机如何逐步改进候选解，从初始的随机可行分配开始，直到使用匈牙利算法计算出最优解。第一幅图像对应于从随机排列获得的随机生成的可行解。其余图像显示了自动机在评估相邻元素之间可能的交换时进行的连续调整。在每一步中，彩色标记表示正在评估的候选解。绿色点表示当前正在测试的交换中涉及的元素，而黄色点对应于在该评估步骤中保持不变的分配。评估过程按行顺序进行。从第0行开始，算法评估与相邻元素的可能交换；然后继续处理第1行，依此类推，直到检查所有行。最后一幅图像对应于使用匈牙利算法获得的最优解。在这个例子中，所提出的自动机产生的解决方案与最优解一致。5. 基于小波和自动机的提出方法的示例自动机构成了改进可行解的一个可行替代方案。然而，希望其初始状态是一个已经接近最优解的Monge–Kantorovich问题的可行解。为了构建这样的初始解，我们采用了[13]中提出的方案的第一阶段，然后将其作为自动机的输入。我们的近似方案包括以下步骤：1. 成本函数的离散化：成本函数在j层上离散化，得到的离散化成本表示为。2. 小波变换的应用：对应用小波变换，我们得到：一个低通分量，表示为；三个高通分量，表示为，和。3. 问题的解决：使用成本函数和[12]中提出的方法，我们计算问题的最优解。4. 问题的可行解构建：使用与高通分量一起，我们构建在j层上的可行解，具有以下属性：(a) (b) (c) (d)这里，表示符号函数。换句话说，的支持包含在的支持中。组分和不参与构建，因为它们影响问题的边界条件。最后，的符号决定了解决方案缩放的方向。5. 自动机的实现：我们使用可行解作为自动机(11)的输入。自动机的输出得到最终近似。5.1. 示例1在这个例子中，我们考虑成本函数。在图4中，左侧表示第5层的离散化，其中成本函数在成本矩阵上转换为，右侧表示成本函数的三维图形。图4. 成本函数及其离散化。根据[12]，这一层的最优解在图5中说明，其中红色点表示运输问题的最优分配。图5. 使用匈牙利算法得到的最优解。下面，我们实现了本节开头描述的算法。将小波变换应用于之前的函数，我们得到图6中的分解。图6. 小波分量。小波分解分量：近似和细节。这里，左上角图示了较低层次的近似，右上角图示了水平细节，左下角图示了垂直细节，右下角图示了对角细节。通过解决由较低层次近似给出的成本函数运输问题，我们得到度量。目标是从和组分构建一个可行解（第5层），如图7所示。图7. 分配比较。原始分配与近似分配之间的比较。图7的描述。左侧的图示了运输问题的最优解，其中成本矩阵是从图6左上角面板中的基于小波的近似得到的。中间的图示了这个解决方案到更高级别的自然扩展；然而，它不构成运输问题的一个可行解。右侧的图示了在更高层次上的运输问题的一个可行解。这个解决方案是通过选择中间图中每个块的一个对角线得到的，选择是由对角线细节矩阵的符号指导的（见图6的右下角面板）。下面，我们实现了在(11)中定义的自动机，如图8所示：图8. 自动机的实现，其初始状态由图7中的右侧图示定义。图8的描述：上面的图示了当初始状态由图7的右侧面板给出时，调整自动机执行的操作。绿色点表示自动机进行局部优化的位置。为了清晰起见，省略了没有发生变化的步骤。左下角的图示了我们的程序获得的最终解决方案，它包括应用小波变换、解决带有成本矩阵（第4层）的运输问题、扩展解决方案，然后应用自动机（第5层）。右下角的图示了直接使用成本矩阵（第5层）解决运输问题得到的解决方案，如[12]中所述。我们算法得到的解决方案和直接应用匈牙利算法得到的解决方案并不完全相同；然而，它们都产生了相同的最小值，因此可以被视为等效解决方案。在本工作中提出的数值实验中，所提出的自动器始终收敛到一个最优分配。然而，对该方法的收敛性质进行完整的理论描述，包括确定最优解唯一性的精确条件，仍然是一个未解决的问题。确定这些条件并建立正式的收敛保证是未来研究的一个有趣方向。

5.2. 示例2
在这个示例中，我们考虑成本函数。在某一层次上的离散化函数如下图9所示：
图9. 成本函数及其离散化。
根据[12]，该层次上的可行解如下图10所示：
图10. 最优分配。使用匈牙利算法在原始成本矩阵上的最优分配。

另一方面，通过对某一层次上的离散化应用小波变换，我们得到了图11中的分解结果，这与图6中显示的分解类似。
图11. 小波分量。小波分解组件的合成图：近似值和细节。
通过解决由较低层次近似值给出的成本函数的运输问题，我们得到了测量值。根据这个解，我们得到了一个可行的解（图12），其形式与图7中呈现的类似。
图12. 分配比较。原始分配与近似分配之间的比较。
在这种情况下，从该测量值构建的自动机没有进行任何修改，如图13所示，这与上一个示例中进行的比较类似。

5.3. 示例3
在这个示例中，我们考虑成本函数。在某一层次上的离散化函数如下图14所示：
图14. 成本函数及其离散化。
根据[12]，该层次上的可行测量值在图15中示意。
图15. 最优分配。使用匈牙利算法在原始成本矩阵上的最优分配。
通过对某一层次上的离散化应用小波变换，我们得到了图16中的分解结果，这与图6中显示的分解类似。
图16. 小波分量。小波分解组件：近似值和细节。
由较低层次近似值给出的成本函数对应的运输问题求解产生了测量值。根据这个结果，我们得到了一个可行的解（图17），其形式与图7中描绘的类似。
图17. 分配比较。原始分配与近似分配之间的比较。
图18展示了自动机对可行测量值进行的更改，从而得到了可行测量值，以及自动机在应用缩放程序后生成的候选解序列。该过程从缩放后获得的近似值开始，并进行初步调整，每个后续图像代表在评估可能的交换时生成的新近似值。倒数第二图像表示自动机获得的最终解，而最后一幅图像显示了用于比较的匈牙利算法计算出的最优解。

6. 结论
这项工作的一个重要部分是在排列群上引入了一个自动机，这被证明是一个有用的工具，可以改进与MK问题相关的运输问题的可行解。这代表了将自动机理论应用于优化问题的一个重要贡献。
在所有与经典Monge–Kantorovich问题相关的运输问题（分配问题）示例中，这些问题都具有代理的数量。因此，使用匈牙利算法解决这类问题的计算成本为。
我们解决由此产生的运输问题的方法包括以下步骤：
1. 对大小为的 cost matrix 应用小波变换，其计算成本为。
2. 在小波变换的近似系数上解决运输问题（分配问题）。在这个较粗的分辨率层次上，问题涉及到 agents，这比原始的分配问题更简单。这个简化后的分配问题使用匈牙利算法来解决，以便为较小的 cost matrix 获得最优分配。
3. 将层次上的解扩展到层次 j 的计算成本为，因为对于大小为的数组中的每个元素，只涉及固定数量的步骤。
4. 线性自动机在数组大小上的计算成本为。
因此，这种方法解决与MK问题相关的运输问题的计算成本显著低于直接应用匈牙利算法的成本。由于我们算法的第1、3和4步都是的阶，该方法的主要计算成本集中在第2步。在这项工作中，这一步是通过将匈牙利算法应用于只涉及一半代理的简化问题来解决的。
尽管匈牙利算法对于这个简化问题的阶为，但在最坏的情况下具有相同的渐近阶，但实际上所需的操作更少。此外，为了解决具有 agents 的问题，可以再次应用我们的算法，将问题简化为具有 agents 的问题，依此类推。研究这些技术的极限和范围是未来研究的一个方向。
此外，这种方法产生的结果与使用匈牙利算法解决与MK问题相关的运输问题得到的结果相当。为了说明这一点，我们在层次上展示了三个示例，清楚地证明了获得了等效的结果。
在本文考虑的运输问题中，成本矩阵是从连续函数的离散化中得到的，这为矩阵系数提供了一定的局部稳定性。未来的工作将旨在探索这种技术在成本矩阵熵较高的情况下的极限和范围。
在图3中，过程从一个随机初始解开始，直接应用自动机来解决所考虑的问题。在这种情况下，自动机在层次3成功找到了最优解。然而，没有研究自动子在更高分辨率层次上的应用。作为未来研究的一个方向，分析将自动机应用于更一般分配问题的范围和局限性将是有趣的，特别是在成本矩阵不是从连续函数的离散化中得到的情况下。
此外，通过将更多的多分辨率层次应用于成本函数来评估该技术的性能和鲁棒性也将是研究的内容。
本文提出的基于自动机的精细化策略也可以应用于最优部分运输问题。在这种设置中，只运输了一部分总质量，这在分配结构中引入了额外的约束。自动机的局部调整机制通过顺序交换分配来操作，可以适应于处理部分匹配，允许未分配的元素或松弛变量代表未使用的质量。尽管这样的扩展需要对自动机的状态表示和评估规则进行修改，但迭代局部细化的基本思想仍然适用。研究这种扩展并分析其数值行为是未来研究的一个有趣方向。