在Python中实现3D向量:numpy vs x,y,z字段

Question

我正在用Python实现3D Vector类.我的向量有坐标x,y和z(所有浮点数),我需要决定如何存储这些信息.我在这里至少可以看到三个选项:

1)制作三个独立的浮点字段:self.x,self.y,self.z

class Vector:

  def __init__(self, x, y, z):
    self.x = x
    self.y = y
    self.z = z

2)制作一个列表,比如说self.data,有三个元素.如果对象可以是常量,我也可以使用元组.

class Vector:

  def __init__(self, x, y, z):
    self.data = [x,y,z]

3)制作一个numpy数组,比如self.data,有三个元素.

import numpy as np    

class Vector:

  def __init__(self, x, y, z):
    self.data = np.array([x,y,z])

对于选项(2)和(3),我可以实现属性和设置器来访问单个坐标

@property
def x(self):
  return self.data[0]

4)为什么不进行冗余？我可以同时拥有一个列表(或元组或numpy数组)和单独的字段x,y和z.

该类用于执行常见操作,例如向量添加,内积,叉积,旋转等.需要考虑这些操作的性能.

是否有我更喜欢的解决方案,为什么？

Answer 1

这个问题有不同的方面,我可以给你一些关于如何解决这些问题的提示.请注意,这些都是建议,你肯定需要看看你最喜欢哪一个.

支持线性代数

您提到要支持线性代数,例如向量加法(逐元素加法),交叉积和内积.这些是可用的,numpy.ndarray所以你可以选择不同的方法来支持它们:

1. 只需使用numpy.ndarray并且不要为自己的课程烦恼:

   import numpy as np
   vector1, vector2 = np.array([1, 2, 3]), np.array([3, 2, 1])
   np.add(vector1, vector2)      # vector addition
   np.cross(vector1, vector2)    # cross product
   np.inner(vector1, vector2)    # inner product
   

   没有定义内置矢量旋转,numpy但有几个可用的源,例如"3D矢量的旋转".所以你需要自己实现它.

2. 您可以创建一个类,独立于_how存储属性并提供__array__方法.这样你可以支持(所有)numpy函数,就好像你的实例是numpy.ndarray自己的:

   class VectorArrayInterface(object):
       def __init__(self, x, y, z):
           self.x, self.y, self.z = x, y, z
   
       def __array__(self, dtype=None):
           if dtype:
               return np.array([self.x, self.y, self.z], dtype=dtype)
           else:
               return np.array([self.x, self.y, self.z])
   
   vector1, vector2 = VectorArrayInterface(1, 2, 3), VectorArrayInterface(3, 2, 1)
   np.add(vector1, vector2)      # vector addition
   np.cross(vector1, vector2)    # cross product
   np.inner(vector1, vector2)    # inner product
   

   这将返回与第一种情况相同的结果,因此您可以为numpy函数提供一个接口,而无需使用numpy-array.如果您有存储在您的类numpy的阵列的__array__方法可以简单地返回它,这可能是用于存储你的论点x,y并z为numpy.ndarray内部(因为这基本上是"免费").

3. 你可以继承np.ndarray.我不会在这里详细介绍,因为这是一个高级主题,可以很容易地证明整个答案本身.如果您真的考虑到这一点,那么您应该查看"Subclassing ndarray"的官方文档.我不推荐它,我参与了几个做子类的类,np.ndarray并且在这条路上有几个"粗糙的egdes".

4. 您可以自己实施所需的操作.这是重新发明的轮子,但它具有教育性和趣味性 - 如果只有少数几个.我不推荐这个用于严肃的制作,因为这里也有几个已经在numpy函数中得到解决的"粗糙边缘".例如溢出或下溢问题,功能的正确性,......

   可能的实现(不包括旋转)可能如下所示(这次是内部存储的列表):

   class VectorList(object):
       def __init__(self, x, y, z):
           self.vec = [x, y, z]
   
       def __repr__(self):
           return '{self.__class__.__name__}(x={self.vec[0]}, y={self.vec[1]}, z={self.vec[2]})'.format(self=self)
   
       def __add__(self, other):
           x1, y1, z1 = self.vec
           x2, y2, z2 = other.vec
           return VectorList(x1+x2, y1+y2, z1+z2)
   
       def crossproduct(self, other):
           x1, y1, z1 = self.vec
           x2, y2, z2 = other.vec
           return VectorList(y1*z2 - z1*y2,
                             z1*x2 - x1*z2,
                             x1*y2 - y1*x1)
   
       def scalarproduct(self, other):
           x1, y1, z1 = self.vec
           x2, y2, z2 = other.vec
           return x1*x2 + y1*y2 + z1*z2
   

   注意:您可以实现这些可编码的方法并实现__array__我之前提到的方法.这样你就可以支持任何期望的功能numpy.ndarray,也可以使用你自己开发的方法.这些方法并不是唯一的,但是你会得到不同的结果,上面的方法会返回一个标量或者一个Vector但是如果你经历过,__array__你会得到numpy.ndarray回报.

5. 使用包含3D矢量的库.从某种意义上说,这是其他方面最简单的方法,它可能非常复杂.从好的方面来说,现有的类可能是开箱即用的,它可能在性能方面进行了优化.另一方面,您需要找到一个支持您的用例的实现,您需要阅读文档(或通过其他方式弄清楚它是如何工作的),并且您可能会遇到对您的项目来说非常糟糕的错误或限制.啊,你得到一个额外的依赖,你需要检查许可证是否与您的项目兼容.另外,如果您复制实现(检查许可证是否允许!),您需要维护(即使它只是同步)外部代码.

性能

在这种情况下,性能很棘手,所提到的用例非常简单,每个任务应该是微秒级 - 所以你应该能够每秒执行几千到几百万次操作.假设你没有引入不必要的瓶颈!但是,您可以微观优化操作.

让我从一些一般的tipps开始:

• 避免numpy.ndarray< - > list/ float操作.这些都很贵!如果大多数操作使用numpy.ndarrays,则不希望将值存储在列表中或作为单独的属性存储.同样,如果您想要访问Vector这些值的各个值或迭代这些值或对它们执行操作,list然后将它们存储为列表或单独的属性.

• 使用numpy三个值进行操作相对低效.numpy.ndarray对于大数组来说非常棒,因为它可以更有效地存储值(空间)并且比纯Python操作更好地扩展.然而,这些优点对于小型阵列来说有一些重要的开销(例如length << 100,这是一个有根据的猜测,而不是一个固定的数字!).python解决方案(我使用上面已经介绍过的解决方案)可以比这种小型数组的numpy解决方案快得多:

  class VectorArray:
      def __init__(self, x, y, z):
          self.data = np.array([x,y,z])
  
  # addition: python solution 3 times faster
  %timeit VectorList(1, 2, 3) + VectorList(3, 2, 1)  
  # 100000 loops, best of 3: 9.48 µs per loop
  %timeit VectorArray(1, 2, 3).data + VectorArray(3, 2, 1).data  
  # 10000 loops, best of 3: 35.6 µs per loop
  
  # cross product: python solution 16 times faster
  v = Vector(1, 2, 3)
  a = np.array([1,2,3])  # using a plain array to avoid the class-overhead
  %timeit v.crossproduct(v)  
  # 100000 loops, best of 3: 5.27 µs per loop
  %timeit np.cross(a, a)     
  # 10000 loops, best of 3: 84.9 µs per loop
  
  # inner product: python solution 4 times faster
  %timeit v.scalarproduct(v)  
  # 1000000 loops, best of 3: 1.3 µs per loop
  %timeit np.inner(a, a)      
  # 100000 loops, best of 3: 5.11 µs per loop
  

  但是就像我说的那样,这些时间是微秒级,所以这就是微观优化.但是,如果您专注于课程的最佳表现,那么使用纯python和自我实现的功能可以更快.

  一旦尝试进行大量线性代数运算,就应该利用numpys向量化运算.其中大多数与您描述的类不兼容,并且完全不同的方法可能是合适的:例如,以与numpys函数正确接口的方式存储数组向量数组(多维数组)的类!但我认为这个答案超出了范围,并且不会真正回答你的问题,这个问题仅限于只存储3个值的类.

• 我用不同的方法使用相同的方法做了一些基准测试,但这有点作弊.通常,您不应该为一个函数调用计时,您应该测量程序的执行时间.在程序中,被称为数百万次的函数中的微小速度差异可以比仅仅被称为几次的方法中的大速度差异产生更大的整体差异....或者不是!我只能为函数提供时间,因为您没有共享程序或用例,因此您需要找出哪种方法最适合您(正确性和性能).

结论

还有其他几个因素需要考虑哪种方法最好,但这些因素更多是"元"因素,与您的计划没有直接关系.

• 重新发明轮子(自己实现功能)是一个学习的机会.你需要确保它正常工作,你可以计时,如果它太慢,你可以尝试不同的方法来优化它.你开始考虑算法复杂性,常数因素,正确性......而不是考虑"哪个函数将解决我的问题"或"我如何使numpy函数正确解决我的问题".

• 使用NumPy进行长度为3的阵列可能就像"用苍蝇拍摄大炮",但这是一个很熟悉numpy功能的好机会,将来你会更多地了解NumPy的工作原理(矢量化,索引,广播, ...)即使NumPy不适合这个问题和答案.

• 尝试不同的方法,看看你有多远.我在回答这个问题时学到了很多东西,尝试这些方法很有趣 - 比较差异的结果,调整方法调用的时间并评估它们的局限性!