图片来自 Unsplash，由作者修改
·
关注 发表在 Towards Data Science ·14 分钟阅读·2023 年 8 月 4 日
--
魔方是一个具有巨大状态空间的规划问题原型，且只有一个解决方案。这正是“针在干草堆中”这一概念的定义。如果没有指导（即使你每秒可以旋转面 100 次），你可能会在整个宇宙的时间里也无法成功。关于它的一切似乎都涉及到巨大的数字。
我们在这里计算的量是一个例外。通过它，你将获得对一个困难问题（以及任何类似的规划问题）的简单视角。我们需要两个要素，一个随机过程和一个最佳求解器。后者是一个设备（真实的或理想的），可以在任何初始状态下使用最少的移动次数来解决魔方（或类似问题）。我们将完全回答以下问题：
如果一个已解魔方经历了 N 次随机转动， p(d|N) 的概率是多少，即一个最佳求解器需要 d 次移动才能将其恢复到原始状态？
在正常情况下，如果有人要求你解魔方，你只会得到一个打乱的魔方，没有任何参考或标签。在这里，我们有关于打乱状态的一条信息：它是在从已解状态开始的N次随机移动后获得的。这条信息很有用！
在计算上，你可以尝试以不同的方式解密魔方。一个魔方项目的目标可能在于以次优方式解决任何或某些状态，或者以最佳方式解决每个可能的状态（例如，这将需要著名的35 CPU 年）。一个魔方求解器通常涉及两个方面，一个搜索算法和一个启发式函数。通过选择这两个方面，我们参数化了我们的方法的难度、效率或计算要求。
在启发式函数领域，即搜索引导，总是存在创新的空间。历史上，魔方的启发式方法是将对打乱的魔方面片相对于其已解状态位置的曼哈顿距离估算相结合。直到最近，神经网络才被用作启发式方法。
神经网络 = 魔方的新启发式方法
神经网络的工作很简单（一个分类器）：你输入一个魔方的状态x，然后预测该状态的深度d。状态d的深度定义为从该状态开始解决魔方所需的最少移动次数。请注意以下几点。如果我们有一个知道任何状态深度的设备，我们实际上就有了一个最佳求解器，因为每次我们可以选择一个使状态深度更低的移动，直到达到深度 = 0（已解状态）。
这里的问题是如何训练该网络。或者，具体来说，如何获得准确的训练数据。除非你已经拥有一个最优求解器，否则没有简单的方法知道一个打乱状态x的真实深度d。我们没有最优求解器，或者，我们不想使用计算代价高昂的求解器。我们想从头开始构建一个近似且高效的最优求解器，并尽量减少人工输入，同时也需要准确的训练数据：
training_data = (x , d).
正如我们所说，d的准确性很难获得，但将某个特定打乱状态与数字N关联却很容易：即通过对已解决状态进行N次随机移动生成的状态。然后
p(d|N) 估计 d, 给定 N.
***p(d|N)***将用于提高该训练数据的准确性。前述论文的作者建立了第一个魔方深度分类神经网络。他们的训练数据形式为：
training_data = (x , N).
他们将d视为N. 这个选择通过在训练过程中使用类似 Bellman 的循环动态提高标签的准确性来进行补偿。这里计算的概率p(d|N)为训练数据的准确性提供了一个更好的起点（仅通过随机旋转已解决状态N次即可获得大量数据）。
计算p(d|N)相当于问一个随机游走者在N步之后会离d多远。不是在方格网格上行走，而是在一个拥有 10 的 19 次方节点（立方体状态）和相似数量连接（合法移动）的巨大魔方图上行走。如果我们选择一个布局，将节点按深度组织：将已解决状态置于中心，深度为d的状态位于距离中心d的半径上，则图将看起来非常对称。径向（深度）方向提供了一个非常简单的视角。
在这里，我们采用所谓的 3x3x3 魔方的四分之一转度量，其中一次移动涉及 90 度的面旋转，无论是顺时针还是逆时针。在这种度量下，有十二种可能的移动。如果我们选择了不同的度量，例如半转度量（也将 180 度的面旋转作为一次移动），那么***p(d|N)***的表达式将会有所不同。
要获得p(d|N)，我们需要使用某种领域知识，但我们不想处理图、模式数据库或群体理论。我们将使用一些更“基础”的东西：
包含深度为d的魔方状态数量的列表
这个列表（由 2012 年“上帝的数字”论文的作者提供）没有指定哪些状态在某个特定深度，只提供了它们的总数，并且没有提及任何N.
# Depth population list
# number of cubes' states at a depth d in the quarter-turn metric
D_d={
# depth number of states
0: 1,
1: 12,
2: 114,
3: 1068,
4: 10011,
5: 93840,
6: 878880,
7: 8221632,
8: 76843595,
9: 717789576,
10: 6701836858,
11: 62549615248,
12: 583570100997,
13: 5442351625028,
14: 50729620202582,
15: 472495678811004,
16: 4393570406220123,
17: 40648181519827392,
18: 368071526203620348,
19: 3000000000000000000, # approximate
20: 14000000000000000000, # approximate
21: 19000000000000000000, # approximate
22: 7000000000000000000, # approximate
23: 24000000000000000, # approximate
24: 150000, # approximate
25: 36,
26: 3,
27: 0
}
深度对状态数量的对数刻度图
首先，深度大于 26 时没有状态（在四分之一转动度量中，上帝的数字是 26）。其次，列表中报告了19到24之间d的状态的近似数量。我们稍后需要对此保持谨慎。第三，如果我们绘制对数尺度图，我们可以看到大多数深度（除了接近末端的那些）呈线性增长。这意味着状态数量D(d)以指数方式增长。将对数图的线性部分拟合成一条直线，我们发现d = 3和d = 18之间，状态数量增长为
在3 < d < 18的深度上，平均来说，9.34个12次移动会使你远离已解决状态，而2.66次会使你更接近已解决状态。
要找出p(d|N)，我们可以把深度类看作马尔可夫过程的站点。让我解释一下：
随机旋转立方体面被描述为深度类之间的马尔可夫过程（一维随机游走）。作者提供的图像。
一个深度类d是指所有在深度d下的立方体状态（到达已解决状态的最小步数）。如果我们在深度类d中随机选择一个状态，并用随机的动作旋转一个随机的面，这将以概率p_d给我们一个深度为d + 1的状态，或者以概率q_d给我们一个深度为d - 1的状态。在四分之一转动度量中，没有自类转换。
这定义了一个马尔可夫过程，其中特定的站点是一个完整的深度类。在我们的例子中，只有相邻的d类是一步跳转连接的。准确地说，这是一个离散时间出生-死亡马尔可夫链。由于站点数量是有限的，因此该链也是不可约且遍历的，并且存在唯一的平稳分布。
我们假设在每次选择随机动作时概率是均匀分布的。这会产生一些深度类之间的转移概率p_d, q_d（待计算）。随机动作的数量N是马尔可夫过程的离散时间。这也是一个一维随机游走：在每个站点（深度类编号d）中，前进的概率是p_d，后退的概率是q_d。这个一维链，粗略来说，是鲁比克图中的“径向”方向（按深度-径向布局组织）。
任何马尔可夫过程都可以用转移矩阵M编码。M的**(i,j)项是从站点i跳到站点j**的概率。在我们的例子中，只有以下项不同于零：
这里 p_0 = 1: 从深度等级0**（仅包含已解决状态）我们只能跳到深度等级1（不存在等级***-1***）。同样，q_26 = 1: 从深度等级26我们只能跳到深度等级25**（不存在等级27）。出于同样的原因，**p_26 或 **q_0 不存在。
我们将立方体的随机移动作用映射为一个一维深度等级随机游走者，以概率 q_d 和 p_d 来回跳动。长时间的行走会发生什么？或者，游走者在长时间的行走后访问特定位置的次数是多少？在现实生活中：当立方体经历随机旋转时，深度等级的访问频率是多少？
从长远来看，无论起点是什么，游走者在深度等级 d 上花费的时间与该深度等级的人口 D(d) 成正比。这是这里的重点：
（归一化的）深度人口列表 D(d) 应被解释为表示我们深度等级马尔可夫过程的平稳分布的向量。
从数学上讲，D(d) 是 M 的左特征向量
这个矩阵方程将给出26个线性方程，我们将从中得到 p_i’ 和 q_i*’*。
考虑到 **p_0 = q_26 = 1, 我们可以将这些重新写作
详细平衡方程。图像由作者提供。
这些被称为 详细平衡方程：流量，定义为站态位置人口与跳跃概率的乘积，在两个方向上是相同的。解为：
并且 p_i 是通过 p_i + q_i = 1. 获得的。
这些解有趣的地方在于，因为 q_i 是一个概率，我们应该有
这转化为分布 D_k 的以下条件：
这是深度人口 D_k 应该满足的一个不等式塔。明确地，它们可以组织为：
特别是，最后两个不等式是
因为 D_27 = 0, 我们得到下限和上限相等，所以
或者：
偶数位置的总人口应该等于奇数位置的总人口！
我们可以将其视为偶数和奇数站点之间的详细平衡：每一步总是到达不同且相邻的深度类。任何跳跃都会将你从奇数深度类（所有奇数深度类的类）带到偶数深度类（所有偶数深度类的类）。因此，奇数到偶数类的跳跃发生的概率为 1（反之亦然）。由于两个方向的概率都是 1，它们的数量应该通过详细平衡来相等。
出于同样的原因，马尔可夫过程将达到一个周期为二的“平稳分布”，在每次移动后在偶数和奇数站点之间切换（离散时间N）。
我们计划使用的数据的深度人口D_d在source中报告的是大致的，对于d = 19,20,21,22,23,24\。因此，没有保证它会满足所有这些条件（不等式）。如果我们得到一些概率q_i*超出了[0,1]范围（如情况所示！）。特别是，如果我们尝试检查最后一个条件（偶数-奇数人口平衡），它差距很大！（更新：见末尾注释）
奇数类似乎人口不足（这是作者选择报告数据的近似的结果）。为了使结果有效（使概率在[0,1]范围内），我们决定将之前的大数字添加到深度类别 21 的人口中（具有最大人口的奇数类，或者，最不容易注意到这个大增量的类）。通过这个修正，所有得到的概率似乎都是正确的（这意味着不等式也得到了满足）。
跳跃概率为：
p_i =
{1., 0.916667, 0.903509, 0.903558, 0.903606, 0.903602, 0.90352, 0.903415,
0.903342, 0.903292, 0.903254, 0.903221, 0.903189, 0.903153, 0.903108,
0.903038, 0.902885, 0.902409, 0.900342, 0.889537, 0.818371, 0.367158,
0.00342857, 6.24863*1e-12, 0.00022, 0.0833333}
# i from 0 to 25
q_i =
{0.0833333, 0.0964912, 0.0964419, 0.096394, 0.0963981, 0.0964796,
0.096585, 0.096658, 0.0967081, 0.0967456, 0.0967786, 0.0968113,
0.0968467, 0.0968917, 0.0969625, 0.0971149, 0.0975908, 0.0996581,
0.110463, 0.181629, 0.632842, 0.996571, 1., 0.99978, 0.916667, 1.}
# i from 1 to 26注意几乎所有前p_i（直到i = 21）都接近1。这些是远离已解决状态的概率。接近已解决状态的概率（q_i）对于i大于21几乎为1。这使我们认识到为何解决魔方困难：随机行走者（或魔方的随机移动者）将“永远困在”深度类21的邻域中。
将p_i, q_i数值代入转移矩阵M中，马尔可夫过程得到完整描述。特别是，如果我们以概率一从站点0开始，经过N步后，随机行走者将以概率到达站点d：
这是我们寻找的概率：
从数值上看，我们了解到p(d|N)只有在N和d具有相同奇偶性时才非零（这是某些魔方学者所熟知的）。下面我们为不同的N绘制了一些p(d|N)：
一些概率 p(d|N)。作者图像。
例如：经过 N = 18 次随机移动（绿色曲线），结果立方体的状态在深度 d = 17 的可能性比在深度 d = 19 更高。我们还观察到，在 N = 31 或 32 时，p(d|N) 与平稳分布 D(d) 非常接近（除了它在偶数和奇数位置之间来回切换）。这是另一个回答多少步足以说明我们真的打乱了魔方的问题。
请注意，我们解决了一个逆问题。我们从平稳分布中得到了转移概率。这对于一般的马尔科夫过程是不可能的。特别是，对于半转度量，通过这里描述的方法不能找到 p(d|N)。半转度量是不同的，因为我们可以在移动后停留在同一深度类别（按他们的移动定义）。这些自深度类别跳跃在转移矩阵的对角线上引入了额外的概率 r_i，我们将有更多的变量而不是方程。
即使从计算角度来看，魔方是一个 35 CPU 年问题，它的许多方面仍然可以通过分析或适度的数值计算来描述。我们在这里计算的概率就是一个例子。我们所说的一切都可以很容易地推广到更复杂的魔方-兄弟。一个很酷的推广是进入更多维度：S = 4D, 5D, 6D，……维度的魔方。在这些情况下，状态空间随着 S 的增加而呈指数增长。因此，有些魔方的马尔科夫链可以长到我们想要的长度。换句话说，我们有类似的谜题，其中上帝的数字可以大到我们想要的程度（粗略地说，上帝的数字是状态数量的对数，而状态数量随着维度 S 的增加而增加）。在这些情况下，我们可以采取某些极限来解释我们3D情况的某些方面，就像下一个：
在大上帝数 G 的极限情况下， p(d|N) 的概率应该接近二项分布
很容易看出为什么会这样。如果 G 很大，D_d 的指数增长将在大多数 d 下非常稳定。这是一个大胆但并非过于疯狂的猜测：
对于远离 0 和 G 的 k。如我们所说，在 S = 3D 的情况下 b = 9.34。对于更高的 S，b 应该增加（拥有更多面会增加分支因子 b）。这转化为以下概率值：
当i远离原点（i >> 1）且上帝的数字（i << G）时，q_i会接近一个常数值1/(b+1)。在这个范围内，p_i也将是常数。你可以看到，这里为3D情况计算的p_i和q_i的值对于i = 3, …, 15几乎是常数，并且q_i大约等于1/(b+1)，其中b = 9.34。对于0 << i << G，我们将得到一个具有常数返回概率（q）和前进概率（p）的一维随机行走者。在这种情况下，行走者的位置将由类似二项分布的分布描述。
随机行走者在经过N次试验后，向右走k步（成功率p）和向左走N-k步（成功率q）的概率为Binomial(k,N,p)。在经过N步后行走的距离将是
d = k - (N - k)
从中我们得到
从这里我们可以得到最可能的d的解析估计
最可能的d（在二项分布范围内）随N线性增长，斜率依赖于“有效”分支因子b。随着立方体维度S的增加，分支因子b也会增加，而最可能的深度d实际上接近于N。
更新说明。在这个故事发布后，我联系了 Tomas Rokicki 和 Morley Davidson（2012 年证明上帝的数字在半转度量中为 20的两位作者）关于他们报告的D(d)以及我使用这些数据得到的负概率。他们与我分享了更准确的数据，包括深度 d = 19, …, 24 的上下界。他们的数据与这里得到的不等式完全兼容，并且与偶数深度类别的数量应等于奇数深度类别的数量（在四分之一转度量中）的事实相一致。使用这些新数据计算的概率具有微不足道的修正。