python學(xué)習(xí)筆記 --- python中的list和numpy中的矩陣分析

DobbyKim 發(fā)布于2019-07-30 14:18 / 3276人閱讀

摘要：中的和中的矩陣分析由于之前在做的源碼學(xué)習(xí)，并且將其的源碼翻譯成了的版本。在逛知乎里，我又發(fā)現(xiàn)了很多關(guān)于為什么這么快的討論，很有意思。作者鏈接來源知乎著作權(quán)歸作者所有。

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

Author : Jasper Yang

School : Bupt

preface

由于之前在做GIbbsLDA++的源碼學(xué)習(xí)，并且將其c++的源碼翻譯成了python的版本。后來有朋友用我的實現(xiàn)在大數(shù)據(jù)量的情況下內(nèi)存跑崩潰了，仔細去網(wǎng)上一查，才發(fā)現(xiàn)了python中的list的實現(xiàn)方式是一種很泛化的面對各種類型的數(shù)據(jù)結(jié)構(gòu)，這個結(jié)構(gòu)用來做二位數(shù)組比numpy中的narrays需要占用更多更過的內(nèi)存，后面我會詳細分析。（但是我發(fā)現(xiàn)好像并不是如此）

正文 他兩在使用上的區(qū)別

    l=[1,2L,3.0,"a"] # 這是一個放了不同數(shù)據(jù)類型的list
      
    a=np.array([1,2,3]) # list -> array   
    b=np.array([[1,2,3],[4,5,6]])   
    c=list(a)   # array到list的轉(zhuǎn)換  
    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中為每一個數(shù)據(jù)分配了一個指針，正因為這種實現(xiàn)方式所以可以面對不同的數(shù)據(jù)類型。但因如此，l用了4個指針和4個數(shù)據(jù)的內(nèi)存空間。

    # 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占用內(nèi)存的情況，用到了sys里的函數(shù)

    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占用的內(nèi)存還要大~？

再試試

    s1 = ()
    s2 = []
    s3 = zero(0) # 創(chuàng)建0矩陣，用array的方法
    print(sys.getsizeof(s1))
    print(sys.getsizeof(s2))
    print(sys.getsizeof(s3))

output:

    56
    72
    96

這里發(fā)現(xiàn)tuple初始化費用最小，list次之，array的最大，那么說明了array并沒有比占用更小的內(nèi)存空間，但是，沒錯，凡事都有但是。array在矩陣操作的性能上必然是大大大于list的，見網(wǎng)上總結(jié)如下。

Python 的 list 是動態(tài)類型，可以包含不同類型的元素，所以沒有支持諸如點乘等數(shù)學(xué)函數(shù)，因為要為 list 實現(xiàn)這些操作會犧牲性能。

Numpy 數(shù)組是靜態(tài)類型并且齊次。元素類型在數(shù)組創(chuàng)建的時候就已經(jīng)確定了。

Numpy 數(shù)組節(jié)約內(nèi)存。(這點我不能贊同，難道是getsizeof統(tǒng)計不到指針？)

由于是靜態(tài)類型，對其數(shù)學(xué)操作函數(shù)（如矩陣乘法，矩陣加法）的實現(xiàn)可以使用 C 或者 Fortran 完成。

而且通過array的itemsize屬性發(fā)現(xiàn)每次增加數(shù)據(jù)，array的開銷是更多的，也就是說數(shù)組越大，array占用的內(nèi)存比list多越多。

我從Python中優(yōu)化NumPy包使用性能的教程這里看到了一個很有意思的解釋。

為什么NumPy數(shù)組如此高效？

一個NumPy數(shù)組基本上是由元數(shù)據(jù)（維數(shù)、形狀、數(shù)據(jù)類型等）和實際數(shù)據(jù)構(gòu)成。數(shù)據(jù)存儲在一個均勻連續(xù)的內(nèi)存塊中，該內(nèi)存在系統(tǒng)內(nèi)存（隨機存取存儲器，或RAM）的一個特定地址處，被稱為數(shù)據(jù)緩沖區(qū)。這是和list等純Python結(jié)構(gòu)的主要區(qū)別，list的元素在系統(tǒng)內(nèi)存中是分散存儲的。這是使NumPy數(shù)組如此高效的決定性因素。

為什么這會如此重要？主要原因是：

低級語言比如C，可以很高效的實現(xiàn)數(shù)組計算（NumPy的很大一部分實際上是用C編寫）。例如，知道了內(nèi)存塊地址和數(shù)據(jù)類型，數(shù)組計算只是簡單遍歷其中所有的元素。但在Python中使用list實現(xiàn)，會有很大的開銷。(這里也是提到了很大開銷，我很想知道怎么才能看到那很大的開銷在哪？)

內(nèi)存訪問模式中的空間位置訪問會產(chǎn)生顯著地性能提高，尤其要感謝CPU緩存。事實上，緩存將字節(jié)塊從RAM加載到CPU寄存器。然后相鄰元素就能高效地被加載了（順序位置，或引用位置）。

數(shù)據(jù)元素連續(xù)地存儲在內(nèi)存中，所以NumPy可以利用現(xiàn)代CPU的矢量化指令，像英特爾的SSE和AVX，AMD的XOP等。例如，為了作為CPU指令實現(xiàn)的矢量化算術(shù)計算，可以加載在128,256或512位寄存器中的多個連續(xù)的浮點數(shù)。

此外，說一下這樣一個事實：NumPy可以通過Intel Math Kernel Library (MKL)與高度優(yōu)化的線性代數(shù)庫相連，比如BLAS和LAPACK。NumPy中一些特定的矩陣計算也可能是多線程，充分利用了現(xiàn)代多核處理器的優(yōu)勢。

總之，將數(shù)據(jù)存儲在一個連續(xù)的內(nèi)存塊中，根據(jù)內(nèi)存訪問模式，CPU緩存和矢量化指令，可以確保以最佳方式使用現(xiàn)代CPU的體系結(jié)構(gòu)。

在逛知乎里，我又發(fā)現(xiàn)了很多關(guān)于為什么numpy這么快的討論，很有意思。

numpy的許多函數(shù)不僅是用C實現(xiàn)了，還使用了BLAS（一般Windows下link到MKL的，Linux下link到OpenBLAS）。基本上那些BLAS實現(xiàn)在每種操作上都進行了高度優(yōu)化，例如使用AVX向量指令集，甚至能比你自己用C實現(xiàn)快上許多，更不要說和用Python實現(xiàn)的比。。

作者：Zhipeng
鏈接：https://www.zhihu.com/question/30823702/answer/49696394
來源：知乎
著作權(quán)歸作者所有。商業(yè)轉(zhuǎn)載請聯(lián)系作者獲得授權(quán)，非商業(yè)轉(zhuǎn)載請注明出處。

numpy底層使用BLAS做向量，矩陣運算。像求平均值這種vector operation，很容易使用multi-threading或者vectorization來加速。比如MKL就有很多優(yōu)化。

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}	{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常用函數(shù)總結(jié)

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.]])

做了上面這份總結(jié)之后使用各種數(shù)據(jù)類型的創(chuàng)建就不用總是查找了，雖然都很簡單。

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

花式索引（arange）

索引列表（take） -----> 最快

布爾掩碼（compress）

一些優(yōu)化，在 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()
同時上面出錯數(shù)字超出了范圍可以如下解決

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

使用 cumsum() 這樣做叫做 vectorized operations 。

全文完，希望能幫到你學(xué)到姿勢~
另外，給大家推薦一篇我的關(guān)于python中l(wèi)ist的實現(xiàn)的博客 ---> python中l(wèi)ist的實現(xiàn)

補充 -----> 淺拷貝和深拷貝以及廣播

numpy中的矩陣屬于淺拷貝

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

a的值會被改變，這是因為這是淺拷貝，b獲得的是a的指針，數(shù)據(jù)區(qū)域還是一樣的

list則是深拷貝

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

a的值不會被改變，這是因為深拷貝，將數(shù)據(jù)區(qū)域包括指針都復(fù)制了一遍過來

還有一個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中對兩個形狀不同的陣列進行數(shù)學(xué)計算的處理機制。較小的陣列“廣播”到較大陣列相同的形狀尺度上，使它們對等以可以進行數(shù)學(xué)計算。廣播提供了一種向量化陣列的操作方式，因此Python不需要像C一樣循環(huán)。廣播操作不需要數(shù)據(jù)復(fù)制，通常執(zhí)行效率非常高。然而，有時廣播是個壞主意，可能會導(dǎo)致內(nèi)存浪費以致計算減慢。