python中的list和numpy中的矩阵分析

python中的list和numpy中的矩阵分析

• Author : Jasper Yang
• School : Bupt

preface

由于之前在做GIbbsLDA++的源码学习，并且将其c++的源码翻译成了python的版本。后来有朋友用我的实现在大数据量的情况下内存跑崩溃了，仔细去网上一查，才发现了python中的list的实现方式是一种很泛化的面对各种类型的数据结构，这个结构用来做二位数组比numpy中的narrays需要占用更多更过的内存，后面我会详细分析。（但是我发现好像并不是如此）

正文

他两在使用上的区别

	l=[1,2L,3.0,'a'] # 这是一个放了不同数据类型的list
	  
	a=np.array([1,2,3]) # list -> array   
	b=np.array([[1,2,3],[4,5,6]])   
	c=list(a)   # array到list的转换  
	print l
	print a,np.shape(a)
	print b,np.shape(b)  
	print c,np.shape(c)  

程序输出如下：

	[1, 2L, 3.0, 'a']
	[1 2 3] (3,)
	[[1 2 3]
	 [4 5 6]] (2, 3)
	[1, 2, 3] (3,)

可以看到list都有逗号来隔开，因为list中为每一个数据分配了一个指针，正因为这种实现方式所以可以面对不同的数据类型。但因如此，l用了4个指针和4个数据的内存空间。

	# numpy 定义矩阵
	array=([[1,2,3],
        [4,5,6],
        [7,8,9]])
        
	print array[1]
	print array[1][1]

	5
	[4, 5, 6]

下面我想看看tuple，list，array占用内存的情况，用到了sys里的函数

	import sys
	a = (1,2,3,4,5)
	b = [1,2,3,4,5]
	c = np.array([1,2,3,4,5])
	d = np.array((1,2,3,4,5))
	print(sys.getsizeof(a))
	print(sys.getsizeof(b))
	print(sys.getsizeof(c))
	print(sys.getsizeof(d))

output:

	# bytes 为单位
	96
	112
	136
	136

这么一看，有点太奇怪了，怎么会是用了array之后比list占用的内存还要大~？

再试试

	s1 = ()
	s2 = []
	s3 = zero(0) # 创建0矩阵，用array的方法
	print(sys.getsizeof(s1))
	print(sys.getsizeof(s2))
	print(sys.getsizeof(s3))

output:

	56
	72
	96

这里发现tuple初始化费用最小，list次之，array的最大，那么说明了array并没有比占用更小的内存空间，但是，没错，凡事都有但是。array在矩阵操作的性能上必然是大大大于list的，见网上总结如下。

• Python 的 list 是动态类型，可以包含不同类型的元素，所以没有支持诸如点乘等数学函数，因为要为 list 实现这些操作会牺牲性能。
• Numpy 数组是 静态类型 并且 齐次。 元素类型在数组创建的时候就已经确定了。
• Numpy 数组节约内存。(这点我不能赞同，难道是getsizeof统计不到指针？)
• 由于是静态类型，对其数学操作函数（如矩阵乘法，矩阵加法）的实现可以使用 C 或者 Fortran 完成。

而且通过array的itemsize属性发现每次增加数据，array的开销是更多的，也就是说数组越大，array占用的内存比list多越多。

我从Python中优化NumPy包使用性能的教程这里看到了一个很有意思的解释。

为什么NumPy数组如此高效？

一个NumPy数组基本上是由元数据（维数、形状、数据类型等）和实际数据构成。数据存储在一个均匀连续的内存块中，该内存在系统内存（随机存取存储器，或RAM）的一个特定地址处，被称为数据缓冲区。这是和list等纯Python结构的主要区别，list的元素在系统内存中是分散存储的。这是使NumPy数组如此高效的决定性因素。

为什么这会如此重要？主要原因是：

1. 低级语言比如C，可以很高效的实现数组计算（NumPy的很大一部分实际上是用C编写）。例如，知道了内存块地址和数据类型，数组计算只是简单遍历其中所有的元素。但在Python中使用list实现，会有很大的开销。(这里也是提到了很大开销，我很想知道怎么才能看到那很大的开销在哪？)
2. 内存访问模式中的空间位置访问会产生显著地性能提高，尤其要感谢CPU缓存。事实上，缓存将字节块从RAM加载到CPU寄存器。然后相邻元素就能高效地被加载了（顺序位置，或引用位置）。
3. 数据元素连续地存储在内存中，所以NumPy可以利用现代CPU的矢量化指令，像英特尔的SSE和AVX，AMD的XOP等。例如，为了作为CPU指令实现的矢量化算术计算，可以加载在128,256或512位寄存器中的多个连续的浮点数。

此外，说一下这样一个事实：NumPy可以通过Intel Math Kernel Library (MKL)与高度优化的线性代数库相连，比如BLAS和LAPACK。NumPy中一些特定的矩阵计算也可能是多线程，充分利用了现代多核处理器的优势。

总之，将数据存储在一个连续的内存块中，根据内存访问模式，CPU缓存和矢量化指令，可以确保以最佳方式使用现代CPU的体系结构。

在逛知乎里，我又发现了很多关于为什么numpy这么快的讨论，很有意思。

numpy的许多函数不仅是用C实现了，还使用了BLAS（一般Windows下link到MKL的，Linux下link到OpenBLAS）。基本上那些BLAS实现在每种操作上都进行了高度优化，例如使用AVX向量指令集，甚至能比你自己用C实现快上许多，更不要说和用Python实现的比。。

作者：Zhipeng
链接：https://www.zhihu.com/question/30823702/answer/49696394
来源：知乎
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。


numpy底层使用BLAS做向量，矩阵运算。像求平均值这种vector operation，很容易使用multi-threading或者vectorization来加速。比如MKL就有很多优化。

p.s. 一个小小的experiment可以看看（摘抄--大大的橙子）

	from sys import getsizeof
	   
	class A(object): pass
	class B: pass
	
	for x in (None, 1, 1L, 1.2, 'c', [], (), {}, set(), B, B(), A, A()):
	  print "{0:20s}\t{1:d}".format(type(x).__name__, sys.getsizeof(x))
	
	
	
	NoneType                16
	int                     24
	long                    28
	float                   24
	str                     34
	list                    64
	tuple                   48
	dict                    272
	set                     224
	classobj                96
	instance                64
	type                    896
	A                       56

一份numpy常用函数总结

• arange

    # create a range
    x = arange(0, 10, 1) # arguments: start, stop, step    
    x
    
    => array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
    
    
    x = arange(-1, 1, 0.1)
    x
  

  output:

      array([ -1.00000000e+00,  -9.00000000e-01,  -8.00000000e-01,
              -7.00000000e-01,  -6.00000000e-01,  -5.00000000e-01,
              -4.00000000e-01,  -3.00000000e-01,  -2.00000000e-01,
              -1.00000000e-01,  -2.22044605e-16,   1.00000000e-01,
               2.00000000e-01,   3.00000000e-01,   4.00000000e-01,
               5.00000000e-01,   6.00000000e-01,   7.00000000e-01,
               8.00000000e-01,   9.00000000e-01])
  

• mgrid

    x, y = mgrid[0:5, 0:5] 
    print(x)
    print(y)
  

output:

	array([[0, 0, 0, 0, 0],
          [1, 1, 1, 1, 1],
          [2, 2, 2, 2, 2],
          [3, 3, 3, 3, 3],
          [4, 4, 4, 4, 4]])

	array([[0, 1, 2, 3, 4],
          [0, 1, 2, 3, 4],
          [0, 1, 2, 3, 4],
          [0, 1, 2, 3, 4],
          [0, 1, 2, 3, 4]])

• linspace 和 logspace

    linspace(0, 10, 20)
    logspace(0, 10, 10, base=e)
  

output:

	array([  0.        ,   0.52631579,   1.05263158,   1.57894737,
	         2.10526316,   2.63157895,   3.15789474,   3.68421053,
	         4.21052632,   4.73684211,   5.26315789,   5.78947368,
	         6.31578947,   6.84210526,   7.36842105,   7.89473684,
	         8.42105263,   8.94736842,   9.47368421,  10.        ])


	array([  1.00000000e+00,   3.03773178e+00,   9.22781435e+00,
	         2.80316249e+01,   8.51525577e+01,   2.58670631e+02,
	         7.85771994e+02,   2.38696456e+03,   7.25095809e+03,
	         2.20264658e+04])

• random data

    random.rand(4,4)
    random.randn(4,4)
  

output:

	array([[ 0.2659394 ,  0.59798318,  0.16067613,  0.21704906],
	       [ 0.98081628,  0.10926225,  0.86342814,  0.83130784],
	       [ 0.83563301,  0.18313372,  0.73531414,  0.67861801],
	       [ 0.38322499,  0.02943103,  0.43471297,  0.87099786]])

	array([[ 0.18834898,  0.28862928,  0.58291415,  0.57712703],
	       [ 0.75071525, -0.39247518, -0.35748584, -0.42839121],
	       [-1.12789581, -0.67267265, -0.05525579, -0.89712592],
	       [-0.14731484, -0.72237449, -0.16594984,  0.62914291]])

• diag

    diag([1,2,3])
    diag([1,2,3], k=1) 
  

output

	array([[1, 0, 0],
	       [0, 2, 0],
	       [0, 0, 3]])
	       
	array([[0, 1, 0, 0],
          [0, 0, 2, 0],
          [0, 0, 0, 3],
          [0, 0, 0, 0]])

• zeros 和ones

    zeros((2,2))
    ones((2,2))
  

  output

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

做了上面这份总结之后使用各种数据类型的创建就不用总是查找了，虽然都很简单。

numpy的索引很有意思，基本分成了以下三类

• 花式索引 （arange）
• 索引列表 （take） -----> 最快
• 布尔掩码 （compress）

一些优化，在 stackoverflow 上看的，为了更好的记住这些trick我决定记在这篇博客里。

	import numpy
	x = numpy.array([0] * 1000000)
	for i in range(1,len(x)):
	  x[i] = x[i-1] + i
	  
	a[0] += 1232234234234324353453453
	OverflowError: Python int too large to convert to C long

可以变成

import numpy as np
x = np.arange(1000000).cumsum()
同时上面出错数字超出了范围可以如下解决

a = np.array([0], dtype=object)
a[0] += 1232234234234324353453453

使用 cumsum() 这样做叫做 vectorized operations 。

全文完，希望能帮到你学到姿势~
另外，给大家推荐一篇我的关于python中list的实现的博客 ---> python中list的实现

补充 -----> 浅拷贝和深拷贝以及广播

numpy中的矩阵属于浅拷贝

a = np.zeros((2,2))
b = a
b[1][1] = 1
print(a)

a的值会被改变，这是因为这是浅拷贝，b获得的是a的指针，数据区域还是一样的

list则是深拷贝

a = [1,2,3,4]
b = a
b[2] = 6
print(a)

a的值不会被改变，这是因为深拷贝，将数据区域包括指针都复制了一遍过来

还有一个numpy矩阵的广播，这个是一个很神奇的使用方式

a = np.arange(0, 50, 10).reshape(-1, 1)  
b = np.arange(0, 4)  
print a  
print b  

这时的输出

[[ 0]
 [10]
 [20]
 [30]
 [40]]
[0 1 2 3]

当我们使用广播

np.add(a,b)

输出变成了

array([[ 0,  1,  2,  3,  4],
   [10, 11, 12, 13, 14],
   [20, 21, 22, 23, 24],
   [30, 31, 32, 33, 34],
   [40, 41, 42, 43, 44],
   [50, 51, 52, 53, 54]])

是不是很神奇呢，用法一目了然。 see also --> numpy.broadcast

广播用以描述numpy中对两个形状不同的阵列进行数学计算的处理机制。较小的阵列“广播”到较大阵列相同的形状尺度上，使它们对等以可以进行数学计算。广播提供了一种向量化阵列的操作方式，因此Python不需要像C一样循环。广播操作不需要数据复制，通常执行效率非常高。然而，有时广播是个坏主意，可能会导致内存浪费以致计算减慢。

paper done 2017/4/20