有限度量嵌入：好的算法吗？

Question

我有一个有限的度量空间，它是一个对称的 k x k 距离矩阵。我想找一个算法，把这个空间大致地嵌入到欧几里得空间 R^(k-1) 中。虽然通过解决距离给出的方程组并不总是能做到这一点，但我希望找到一个可以在某个（非常小的）可控误差范围内进行嵌入的解决方案。

目前我使用多维尺度法（MDS），并将输出维度设置为 (k-1)。我意识到，通常情况下，MDS 可能更适合于将嵌入维度减少到小于 (k-1) 的情况（通常是 2 或 3），而对于我这种特定情况，可能有更好的算法。

问题是：有什么好的、快速的算法可以在 R^{k-1} 中实现大小为 k 的度量空间，使用欧几里得距离？

一些参数和提示：

(1) 我的 k 值相对较小。比如 3 < k < 25。

(2) 我其实不在乎是否嵌入到 R^{k-1}。如果能简化或加快速度，任何 R^N 都可以，只要是等距的。如果我增加到 R^k 或 R^(2k+1)，能有更快的算法或更小的误差，我也很满意。

(3) 如果你能指向一个 Python 实现，我会更高兴。

(4) 任何比 MDS 更好的方法都可以。

Answer 1

好的，按照承诺，这里有一个简单的解决方案：

说明一下：用 d_{i,j}=d_{j,i} 表示点 i 和 j 之间的平方距离。N 是点的数量。p_i 表示点 i，p_{i,k} 是点的第 k 个坐标。

希望我这次能正确推导出算法。之后会有一些解释，方便你检查推导过程（我讨厌出现太多索引的时候）。

这个算法使用动态规划来找到正确的解决方案。

Initialization:
p_{i,k} := 0 for all i and k

Calculation:
for i = 2 to N do
   sum = 0
   for j = 2 to i - 1 do
     accu = d_{i,j} - d_{i,1} - p_{j,j-1}^2
     for k = 1 to j-1 do
       accu = accu + 2 p_{i,k}*p_{j,k} - p_{j,k}^2
     done
     p_{i,j} = accu / ( 2 p_{j,j-1} )
     sum = sum + p_{i,j}^2
   done
   p_{i,i-1} = sqrt( d_{i,0} - sum )
done

如果我没有犯什么严重的索引错误（我通常会犯），这应该就能解决问题。

这个想法是：

我们把第一个点随便放在原点，这样会让我们的工作更简单。注意，对于点 p_i，我们从来不会在 k > i-1 的情况下设置坐标 k。也就是说，对于第二个点，我们只设置第一个坐标；对于第三个点，我们只设置第一个和第二个坐标，依此类推。

现在假设我们已经有了所有点 p_{k'} 的坐标，其中 k'<i，我们想计算 p_{i} 的坐标，以满足所有 d_{i,k'} 的条件（我们暂时不关心 k>k' 的点的任何约束）。如果我们写下这些方程，我们会得到：

d_{i,j} = \sum_{k=1}^N (p_{i,k} - p_{j,k} )^2

因为对于 k>k'，p_{i,k} 和 p_{j,k} 都等于零，所以我们可以简化为：

d_{i,j} = \sum_{k=1}^{k'} (p_{i,k} - p_{j,k} )^2

还要注意，由于循环不变性，当 k>j-1 时，所有 p_{j,k} 都会是零。所以我们可以把这个方程拆分为：

d_{i,j} = \sum_{k=1}^{j-1} (p_{i,k} - p_{j,k} )^2 + \sum_{k=j}^{k'} p_{i,j}^2

对于第一个方程，我们得到：

d_{i,1} = \sum_{k=1}^{j-1} (p_{i,k} - p_{j,k} )^2 + \sum_{k=j}^{k'} p_{i,i,1}^2

这个方程稍后需要特别处理。

现在如果我们解决正常方程中的所有二项式，我们得到：

d_{i,j} = \sum_{k=1}^{j-1} p_{i,k}^2 - 2 p_{i,k} p_{j,k} + p_{j,k}^2 + \sum_{k=j}^{k'} p_{i,j}^2

从这个方程中减去第一个方程，我们得到：

d_{i,j} - d_{i,1} = \sum_{k=1}^{j-1} p_{j,k}^2 - 2 p_{i,k} p_{j,k}

对于所有 j > 1。

如果你看看这个，你会发现 p_i 的所有坐标平方都消失了，而我们需要的平方已经知道了。这是一个线性方程组，可以很容易地用线性代数的方法解决。实际上，这组方程还有一个特别的地方：这些方程已经是三角形形式，所以你只需要最后一步来传播解。对于最后一步，我们只剩下一个简单的二次方程，可以通过开平方来解决。

希望你能跟上我的思路。现在有点晚了，我的脑袋有点晕，索引太多了。

编辑：是的，确实有索引错误。已经修正。我会在有时间的时候尝试用 Python 实现这个算法，以便测试一下。

Answer 2

为什么不试试 LLE 呢？你也可以在那找到代码。