【问题标题】:How to efficiently compare two sets of unordered vectors (`np.ndarray`)?如何有效地比较两组无序向量(`np.ndarray`)?
【发布时间】:2023-04-08 15:20:01
【问题描述】:

我试图让我的问题更清楚。在这最后还能找到老题外话。

问题:我有一个n x 3 矩阵A (np.ndarray) 的n 点在三个维度上。如果这些点作为无序集是静止的,则我们说这些点相对于变换R(3 x 3 矩阵)是对称的。

这意味着 AA @ R.T 相差 1. 最多一个排列和 2. 校正排列后,两个矩阵可能会因数值容差参数而有所不同:np.allclose(A, permuted(A @ R.T)) == True(我不知道@987654329事先@,这肯定取决于R)。

问题:如何创建以下函数:

def is_symmetric(A, R, atol=1e-5) -> bool:
    # checks symmetry as defined above, considering both numerical noise
    # and permutation of vectors.

(关于可能的方法的一些讨论以及我的怀疑和尝试在下面找到。)


老题外话

我想检查点集合中的对称性,表示为向量。这意味着检查这些点在空间中矩阵变换的应用(例如旋转或平面镜像)是否是不变的。

import numpy as np
R = np.array([[0, -1, 0],  # 90° rotation around z as an example
              [1,  0, 0],
              [0,  0, 1]])

重点是我接受向量的排列:只要将某个初始位置转换为其他某个预先存在的位置,我就可以了。这意味着检查从一个到另一个转换的向量对。

朴素的解决方案将遍历A @ R.T 的行(其中A 是一个矩阵,其行是点位置)并尝试为@ 的每个初始行匹配一个变换向量987654335@,它似乎随着列数呈二次方增长。

另一种可能性是对向量进行预排序(例如,通过它们的坐标值):

A = np.array([[1,   0,   0],  # some points
              [0,   1,   0],
              [0,   0,   1],
              [0.5, 0.5, 0]])
A = np.array(sorted(A, key=lambda v: (v[2], v[1], v[0])))  # sort by z, y, x values
# array([[1. , 0. , 0. ],
#        [0.5, 0.5, 0. ],
#        [0. , 1. , 0. ],
#        [0. , 0. , 1. ]])

A_rotated = np.array(sorted(A @ R.T, key=lambda v: (v[2], v[1], v[0])))
# array([[-1. ,  0. ,  0. ],   # no match
#        [-0.5,  0.5,  0. ],   # no match
#        [ 0. ,  1. ,  0. ],   # match
#        [ 0. ,  0. ,  1. ]])  # match

(这种方法有两种方法,所以 O(n ln(n))?) 第三个想法是比较从原始向量和旋转向量创建的集合。我有一种直觉,这和比较排序的向量一样好。

但还有一件事:如何处理近似比较?我想接受两个向量 vw 相等,如果例如np.allclose(v, w) == True 或同等名称(即abs(v - w) < eps 或类似名称):

np.allclose([1, 0, 0], [1, 0, 0])
# True
np.allclose([1, 0, 0], [1 + 1e-5, 0, 0], atol=1e-5)
# True
np.allclose([1, 0, 0], [1 + 1e-4, 0, 0], atol=1e-5)
# False

那么问题来了:我如何(有效地)比较两组(无序)向量的相等性,同时考虑数值近似(例如np.allclose)?

【问题讨论】:

  • 请准确一点,解释一下isclose()是什么意思,并举例说明。
  • @CatalinaChircu 我编辑以更好地表达我的意思。我对isclose 产生了一些混淆,我的意思是allclose(我将两个不同的ndarrays 作为单个实体进行比较)。

标签: python numpy vector comparison numpy-ndarray


【解决方案1】:

这很简单,您将向量的所有值(<atol)转换为 0。您可以使用 numpy.vectorize 执行此操作,查找文档 here

import numpy as np

def transform(x, atol):
    if x>0 and x-1<atol:
        return 1
    else:
        return x

myarray = np.array([[1, 0, 0], [1 + 1e-4, 0, 0]])
vf = np.vectorize(transform)
new_array = vf(myarray, atol=1e-5)

现在您可以进行对称验证了。

如果我回答了你的问题,请告诉我。

更新

我看到你在numpy 中已经有了函数isclose(),你应该使用它,而不是从头开始使用我的函数。

尽管如此,我不会删除我的答案,以保持与它相关的 cmets。

【讨论】:

  • 我不知道要比较的配对也是问题的一部分。此外,您的函数在您的示例中返回相同的给定数组。我将编辑我的问题以使其更清楚。
  • 对不起,我更正了。事实上,必须改变的是接近1的值。
  • 我在问题中添加了一个(希望)更好的解释。找到向量对也是问题的一部分。
  • 如果我理解得很好,你已经有了你的置换算子(假设切换 90 度),你必须对初始矩阵 A 进行置换,从而获得第二个矩阵 A1。然后你比较它们以检查它们是否对称。您需要的是对称运算符。我做对了吗?
  • 寻找排列并不是一个好的解决方案,因为它通常会随着向量的数量而阶乘增长。请参阅我对问题的建议答案。您能帮我找到更好的解决方案吗?
【解决方案2】:

这是一个使用np.lexsort 工作的函数:

def is_symmetric(A, R, *args, **kwargs):
    A = np.asanyarray(A)
    A = A[np.lexsort(A.T)]

    A_t = A @ np.asanyarray(R).T
    A_t = A_t[np.lexsort(A_t.T)]
    return np.allclose(A, A_t, *args, **kwargs)

一些结果:

R = np.array([[0, -1, 0],  # 90° rotation as an example
              [1,  0, 0],
              [0,  0, 1]])

is_symmetric([[0, 0, 0]], R)
# True
is_symmetric([[1, 0, 0],
              [0, 1, 0],
              [0, 0, 0]], R)
# False
is_symmetric([[1, 0, 0],
              [0, 1, 0],
              [0, 0, 0],
              [-1, 0, 0]], R)
# False
is_symmetric([[1, 0, 0],
              [0, 1, 0],
              [0, 0, 0],
              [-1, 0, 0],
              [0, -1, 0]], R)
# True

100000 个随机向量的性能似乎还不错:

A = np.random.rand(100000, 3)
%timeit is_symmetric(A, R)
# 82.2 ms ± 75.4 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

【讨论】:

  • 排列后是否可以丢失或添加行?你的初始数组是 (3,3),我看到你得到大小 (1,3) 和 (5,3)。哪个是你的初始数组,哪个是排列后的数组? 90 度排列后的 A 可能是大小 A 或大小 A.T.
  • 对我来说,R 旋转 90 度给出:[ [0, 1, 0] [0, 0, -1] [1, 0, 0]]。你的意思是矩阵的旋转?
  • @CatalinaChircu,您的矩阵围绕不同的轴旋转。我的回答显示了可变数量向量的用例。
  • 确实你提到了 3 个维度,但我只能在你的矩阵中看到两个维度。那你为什么不提供 3D 数组呢?我知道你的意思是当你有 2D 时看不到的 z 轴。好的。
  • 我所有的矩阵都有三列,因为A 矩阵中的每一行都有一个点。这就是为什么A 需要右乘以R.T
猜你喜欢
  • 2018-11-17
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2022-01-10
  • 2011-05-10
  • 1970-01-01
相关资源
最近更新 更多