对齐¶

align 是一个可以用于对齐两个不同数据集的模块。特别是，目前存在的三个类，即 SignFlips、常规的 OrthogonalProcrustes 和 SeedlessProcrustes，都旨在校正数据的正交变换，其确切形式未知。这样做的动机是正交不可识别性，这在处理各种嵌入方法时很常见，无论是在统计图还是其他领域。需要注意的是，如果两个图是使用 omnibus 嵌入的——它们不需要对齐。

回想一下，正交矩阵是任何满足 \(Q^T Q = Q Q^T = I\) 的矩阵 \(Q\)。对于所有对齐方法，主要思想是从输入数据集 \(X\) 和 \(Y\) 中推断出一个正交矩阵 \(Q\)，然后使用 \(Q\) 将 \(X\) 对齐到 \(Y\)，通过将 \(X\) 转换为 \(\hat{X} = X Q\)。align 中类的主要区别在于它们适用的设置，以及对矩阵 \(Q\) 所做的额外假设。

[1]:

import numpy as np
import matplotlib.pyplot as plt
from graspologic.plot import heatmap
from scipy.stats import special_ortho_group, ortho_group

%matplotlib inline

简单旋转示例¶

首先，让我们构建两个简单的数据集 \(X\) 和 \(Y\)，其中 \(X\) 的所有条目来自 \(Uniform(0, 1)\) 分布，而 \(Y\) 是 \(X\) 的副本，随后使用随机生成的矩阵 \(Q\) 围绕原点旋转，如下所示：

[2]:

np.random.seed(314)
X = np.random.uniform(0, 1, (15, 2))
Q = special_ortho_group.rvs(2)
Y = X @ Q

[3]:

heatmap(Q, figsize=(4,4), vmin=-1, vmax=1)
Q

[3]:

array([[-0.94924938, -0.31452442],
       [ 0.31452442, -0.94924938]])

../../_images/tutorials_aligning_aligning_5_1.png

[4]:

fig, ax = plt.subplots(1, 1, figsize=(6, 6))
ax.set_title("Simple rotation example")
ax.scatter(X[:,0], X[:,1], label=r"$X$", marker='x')
ax.scatter(Y[:,0], Y[:,1], label=r"$Y$", s=100, alpha=0.5)
ax.set(xlim=(-1.20, 1.20), ylim=(-1.20, 1.20))
ax.set_xticks([0])
ax.set_yticks([0])
ax.legend()
ax.grid();

../../_images/tutorials_aligning_aligning_6_0.png

符号翻转¶

SignFlips 假设 \(Q\) 不仅是正交的，而且是对角矩阵，这限制了矩阵的对角线上只有 \(\pm 1\)，其他地方都是 \(0\)。本质上，它沿着某些轴翻转数据集，使得两个数据集最终位于同一个正交象限（广义象限）。数据集“位于”哪个象限由某些标准定义，默认情况下是维度方向的中位数。如果数据集需要沿某个维度翻转，则 \(Q\) 中的相应条目将为 \(-1\)，否则为 \(1\)。

[5]:

from graspologic.align import SignFlips
aligner_SF = SignFlips()
X_prime_SF = aligner_SF.fit_transform(X, Y)

由于 \(X\) 完全位于第一象限，其在两个维度上的中位数都大于 \(0\)。同时，\(Y\) 主要位于第三象限，因此其维度中位数都小于 \(0\)。因此，为了使两个数据集到达同一象限，\(X\) 需要沿 x 轴和 y 轴翻转。所以，算法找到的 \(Q_{SF}\) 是

[6]:

heatmap(aligner_SF.Q_, figsize= (4,4), vmin=-1, vmax=1)
aligner_SF.Q_

[6]:

array([[-1,  0],
       [ 0, -1]])

../../_images/tutorials_aligning_aligning_10_1.png

[7]:

fig, ax = plt.subplots(1, 1, figsize=(6, 6))
ax.set_title("Simple rotation example, Sign Flips")
ax.scatter(X[:,0], X[:,1],
           label=r"$X$", marker='x')
ax.scatter(Y[:,0], Y[:,1],
           label=r"$Y$", s=100, alpha=0.5)
ax.scatter(X_prime_SF[:,0],
           X_prime_SF[:,1],
           label=r"$X_{SF}$",
           marker='$F$')
ax.set(xlim=(-1.20, 1.20), ylim=(-1.20, 1.20))
ax.set_xticks([0])
ax.set_yticks([0])
ax.legend()
ax.grid();

../../_images/tutorials_aligning_aligning_11_0.png

SignFlips 正确地将 \(X\) 带入第三象限，但由于它是基于启发式的方法，对齐并不完美。请注意，它并不假设 \(X\) 和 \(Y\) 的条目之间存在任何匹配关系。实际上，这两个数据集甚至不需要具有相同的大小。

如果提供单位矩阵作为第二个数据集，此类也可用于将数据集带到第一象限。

正交Procrustes¶

OrthogonalProcrustes 通过解决经典的正交Procrustes问题获得 \(Q_{OP}\)。

正交Procrustes问题是线性代数中的一个矩阵近似问题，要求找到一个正交矩阵\(Q\)，它能够最接近地将给定的矩阵\(X\)映射到另一个给定的矩阵\(Y\)。形式上，它寻求找到一个矩阵\(Q_{OP}\)，使得\(|| X Q_{OP} - Y ||_F\)最小化。

由于这个公式强烈依赖于给定矩阵行之间的对应关系，因此只有当输入数据集 \(X\) 和 \(Y\) 具有相同的形状 (\(n\), \(d\)) 时，才能应用这个类，并且假设 \(X\) 中的第 \(i^{th}\) 个条目与 \(Y\) 中的第 \(i^{th}\) 个条目的顶点值相同，除了某些噪声。在图设置中，这意味着两个图的潜在位置是等价的。

OrthogonalProcrustes 有一个广为人知的闭式解。实际上，这个类只是简单地封装了 scipy.linal.orthogonal_procrustes()。

[8]:

from graspologic.align import OrthogonalProcrustes
aligner_OP = OrthogonalProcrustes()
X_prime_OP = aligner_OP.fit_transform(X, Y)

请注意，\(Q_{OP}\) 是一个通用的正交矩阵，而不是对角矩阵。

[9]:

heatmap(aligner_OP.Q_, figsize= (4, 4), vmin=-1, vmax=1)
aligner_OP.Q_

[9]:

array([[-0.94924938, -0.31452442],
       [ 0.31452442, -0.94924938]])

../../_images/tutorials_aligning_aligning_20_1.png

[10]:

fig, ax = plt.subplots(1, 1, figsize=(6, 6))
ax.set_title("Simple rotation example, OrthogonalProcrustes")
ax.scatter(X[:,0], X[:,1], label=r"$X$", marker='x')
ax.scatter(Y[:,0], Y[:,1], label=r"$Y$", s=100, alpha=0.5)
ax.scatter(X_prime_OP[:,0],
           X_prime_OP[:,1],
           label=r"$X_{OP}$",
           marker='$O$')
ax.set(xlim=(-1.20, 1.20), ylim=(-1.20, 1.20))
ax.set_xticks([0])
ax.set_yticks([0])
ax.legend()
ax.grid();

../../_images/tutorials_aligning_aligning_21_0.png

正如预期的那样，OrthogonalProcrustes 在 \(X\) 是 \(Y\) 的旋转副本时获得精确解。

无籽Procrustes¶

SeedlessProcrustes 是一种算法，它将常规的 Procrustes 扩展到条目之间没有已知对应关系的情况。事实上，即使两个数据集中的条目数量不同，它也可以使用。

形式上，如果 \(X\) 有 \(n\) 个条目，且 \(Y\) 有 \(m\) 个条目，它试图求解一个正交矩阵 \(Q_{SP}\)，以及一个 \(n \times m\) 的矩阵 \(P\)，并满足一个特殊条件，即行之和为 \(1/m\)，列之和为 \(1/n\)，使得量 \(|| X Q_{SP} - P Y ||_F\) 最小化。

与常规的procrustes不同，这个问题没有已知的闭式解。SeedlessProcrustes尝试通过一种迭代算法来解决这个优化问题，该算法交替使用最优传输和常规procrustes。有关此算法的更多详细信息，请参阅SeedlessProcrustes的参考资料中的注释部分。

[11]:

from graspologic.align import SeedlessProcrustes
aligner_SP = SeedlessProcrustes()
X_prime_SP = aligner_SP.fit_transform(X, Y)

[12]:

heatmap(aligner_SP.Q_, figsize=(4,4), vmin=-1, vmax=1)
aligner_SP.Q_

[12]:

array([[-0.94472887, -0.32785266],
       [ 0.32785266, -0.94472887]])

../../_images/tutorials_aligning_aligning_27_1.png

[13]:

fig, ax = plt.subplots(1, 1, figsize=(6, 6))
ax.set_title("SeedlessProcrustes, Simple rotation example")
ax.scatter(X[:,0], X[:,1], label=r"$X$", marker='x')
ax.scatter(Y[:,0], Y[:,1], label=r"$Y$", s=100, alpha=0.5)
ax.scatter(X_prime_SP[:,0],
           X_prime_SP[:,1],
           label=r"$X_{SP}$",
           marker='$S$')
ax.set(xlim=(-1.20, 1.20), ylim=(-1.20, 1.20))
ax.set_xticks([0])
ax.set_yticks([0])
ax.legend()
ax.grid();

../../_images/tutorials_aligning_aligning_28_0.png

SignFlips、OrthogonalProcrustes、SeedlessProcrustes之间的比较¶

为了便于比较，简单旋转示例中三个类的结果被绘制在同一张图上。

[14]:

fig, ax = plt.subplots(1, 1, figsize=(6, 6))
ax.set_title("Aligning, Simple rotation example")
ax.scatter(X[:,0], X[:,1], label=r"$X$", marker='x')
ax.scatter(Y[:,0], Y[:,1], label=r"$Y$", s=100, alpha=0.5)
ax.scatter(X_prime_SF[:,0],
           X_prime_SF[:,1],
           label=r"$X_{SF}$",
           marker='$F$')
ax.scatter(X_prime_OP[:,0],
           X_prime_OP[:,1],
           label=r"$X_{OP}$", marker='$o$')
ax.scatter(X_prime_SP[:,0],
           X_prime_SP[:,1],
           label=r"$X_{SP}$", marker='$S$')
ax.set(xlim=(-1.20, 1.20), ylim=(-1.20, 1.20))
ax.set_xticks([0])
ax.set_yticks([0])
ax.legend()
ax.grid();

../../_images/tutorials_aligning_aligning_31_0.png

[15]:

fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, figsize = (16, 4))
heatmap(Q, title=r'$Q$', vmin=-1, vmax=1, ax=ax1)
heatmap(aligner_SF.Q_, title=r'$Q_{SF}$', vmin=-1, vmax=1, ax=ax2)
heatmap(aligner_OP.Q_, title=r'$Q_{OP}$', vmin=-1, vmax=1, ax=ax3)
heatmap(aligner_SP.Q_, title=r'$Q_{SP}$', vmin=-1, vmax=1, ax=ax4);

../../_images/tutorials_aligning_aligning_32_0.png

在这个简单的例子中，OrthogonalProcrustes 和 SeedlessProcrustes 在 Frobenius 范数下将 \(X\) 与 \(Y\) 对齐得非常好，而 SignFlips 只能将它们带到同一个象限。事实上，当一个数据集是另一个数据集的旋转副本时，OrthogonalProcrustes 可以将其对齐到数值 epsilon 精度。

[16]:

norm_SignFlips = np.linalg.norm(X_prime_SF - Y)
norm_Orthogonal = np.linalg.norm(X_prime_OP - Y)
norm_Seedless = np.linalg.norm(X_prime_SP - Y)
print('Simple example')
print('Difference in Frobenious norm after Sign Flips ', norm_SignFlips)
print('Difference in Frobenious norm after Orthogonal Procrustes ', norm_Orthogonal)
print('Difference in Frobenious norm after Seedless Procrustes ', norm_Seedless)

Simple example
Difference in Frobenious norm after Sign Flips  1.047909310874878
Difference in Frobenious norm after Orthogonal Procrustes  5.800542824167306e-16
Difference in Frobenious norm after Seedless Procrustes  0.046291889999320164

随机示例¶

回想一下，与其他两种方法不同，OrthogonalProcrustes 严重依赖于两个数据集的顶点匹配的事实。如果顶点被打乱，即使在这个简单的例子中，其性能也会下降，如下图所示。

[17]:

from sklearn.utils import shuffle

[18]:

X_shuffled = shuffle(X)

[19]:

aligner_SF_shuffled = SignFlips()
X_prime_SF_shuffled = aligner_SF_shuffled.fit_transform(X_shuffled, Y)
aligner_OP_shuffled = OrthogonalProcrustes()
X_prime_OP_shuffled = aligner_OP_shuffled.fit_transform(X_shuffled, Y)
aligner_SP_shuffled = SeedlessProcrustes()
X_prime_SP_shuffled = aligner_SP_shuffled.fit_transform(X_shuffled, Y)

[20]:

fig, ax = plt.subplots(1, 1, figsize=(6, 6))
ax.set_title("Aligning, Shuffled example")
ax.scatter(X_shuffled[:,0], X_shuffled[:,1], label=r"$X$", marker='x')
ax.scatter(Y[:,0], Y[:,1], label=r"$Y$", s=100, alpha=0.5)
ax.scatter(X_prime_SF_shuffled[:,0],
           X_prime_SF_shuffled[:,1],
            label=r"$X_{SF}$", marker='$F$')
ax.scatter(X_prime_OP_shuffled[:,0],
           X_prime_OP_shuffled[:,1],
            label=r"$X_{OP}$", marker='$o$')
ax.scatter(X_prime_SP_shuffled[:,0],
           X_prime_SP_shuffled[:,1],
            label=r"$X_{SP}$", marker='$S$')
ax.set(xlim=(-1.20, 1.20), ylim=(-1.20, 1.20))
ax.set_xticks([0])
ax.set_yticks([0])
ax.legend()
ax.grid();

../../_images/tutorials_aligning_aligning_40_0.png

[21]:

fig, (ax1, ax2, ax3, ax4) = plt.subplots(1, 4, figsize = (16, 4))
heatmap(Q, title=r'$Q$', vmin=-1, vmax=1, ax=ax1)
heatmap(aligner_SF_shuffled.Q_, title=r'$Q_{SF}$', vmin=-1, vmax=1, ax=ax2)
heatmap(aligner_OP_shuffled.Q_, title=r'$Q_{OP}$', vmin=-1, vmax=1, ax=ax3)
heatmap(aligner_SP_shuffled.Q_, title=r'$Q_{SP}$', vmin=-1, vmax=1, ax=ax4);

../../_images/tutorials_aligning_aligning_41_0.png

在这个打乱的例子中，尽管数据集 \(X_{shuffled}\) 看起来与图中的数据集 \(X\) 相同，但 OthogonalProcrustes 的性能比上面简单旋转示例中的要差。这是因为在从 \(X\) 打乱后，\(X_{shuffled}\) 中的第 \(i^{th}\) 个条目的值不再与 \(Y\) 的第 \(i^{th}\) 个条目相对应。

未匹配的示例¶

如上所述，SeedlessProcrustes 可以应用于具有不同条目数量的两个数据集。考虑以下示例：\(X\) 与之前相同，但新的 \(Y\) 是一个不同的、更大的数据集，从相同的分布生成，然后进行旋转。

[22]:

Y_unmatched = np.random.uniform(0, 1, (30, 2)) @ Q

[23]:

aligner_SF_unmatched = SignFlips()
X_prime_SF_unmatched = aligner_SF_unmatched.fit_transform(X, Y_unmatched)
aligner_SP_unmatched = SeedlessProcrustes()
X_prime_SP_unmatched = aligner_SP_unmatched.fit_transform(X, Y_unmatched)

[24]:

fig, ax = plt.subplots(1, 1, figsize=(6, 6))
ax.set_title("Unmatched input example, Seedless Procrustes")
ax.scatter(X[:,0], X[:,1], label=r"$X$", marker='x')
ax.scatter(Y_unmatched[:,0], Y_unmatched[:,1], label=r"$Y$", s=100, alpha=0.5)
ax.scatter(X_prime_SF_unmatched[:,0],
           X_prime_SF_unmatched[:,1],
           label=r"$X_{SF}$",
           marker='$F$')
ax.scatter(X_prime_SP_unmatched[:,0],
           X_prime_SP_unmatched[:,1],
           label=r"$X_{SP}$",
           marker='$S$')
ax.set(xlim=(-1.20, 1.20), ylim=(-1.20, 1.20))
ax.set_xticks([0])
ax.set_yticks([0])
ax.legend()
ax.grid();

../../_images/tutorials_aligning_aligning_47_0.png

[25]:

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize = (12, 4))
heatmap(Q, title=r'$Q$', vmin=-1, vmax=1, ax=ax1)
heatmap(aligner_SF_unmatched.Q_, title=r'$Q_{SF}$', vmin=-1, vmax=1, ax=ax2)
heatmap(aligner_SP_unmatched.Q_, title=r'$Q_{SP}$', vmin=-1, vmax=1, ax=ax3);

../../_images/tutorials_aligning_aligning_48_0.png

[26]:

print("Mean of the dataset X:      ", np.mean(X, axis=0))
print("Mean of the dataset Y:      ", np.mean(Y, axis=0))
print("Mean of the dataset X_{SF}: ", np.mean(X_prime_SF_unmatched, axis=0))
print("Mean of the dataset X_{SP}: ", np.mean(X_prime_SP_unmatched, axis=0))

Mean of the dataset X:       [0.58510607 0.48290266]
Mean of the dataset Y:       [-0.4035269 -0.6424252]
Mean of the dataset X_{SF}:  [-0.58510607 -0.48290266]
Mean of the dataset X_{SP}:  [-0.38535747 -0.65348582]

SeedlessProcrustes 仍然能够对齐分布，而 SignFlips 方法只是一个好的启发式方法。请注意，OthogonalProcrustes 在这种情况下根本不适用。