Python 機器學習——線性代數和矩陣運算：從matlab遷移到python

阿新 • • 發佈：2019-01-26

誠然，沒有一門語言能夠撼動matlab的矩陣或科學計算在學術圈的地位，因其簡潔的語法（matrix是其基本資料型別），因其矩陣運算的便捷，因其術業有專攻（matlab：為科學計算而生），因其名字matlab：matrix laboratory，所在的公司名mathworks：math works。我在寫過一些matlab和python的程式碼之後，油然發過一句感慨“沒有一門語言能比matlab還更具數學感”。然而，因為numpy、pandas和matplotlib 等一眾優秀的開發者（不排除從matlab陣營溜出來的）貢獻的一眾優秀的開源的庫，讓python具備了和matlab一樣的功能，為工程而生的python從此因為有了數學家的參與就相當初的matlab一樣，也學術起來，工程學術通吃。

本文試圖回答的問題包括：

為什麼矩陣運算要從matlab遷移到python？
如何進行遷移，其中會涉及哪些基本程式設計理念的差異？
遷移的過程中需要注意哪些細節？

python矩陣運算，更準確地說，是numpy矩陣運算，為了更為方便地使用numpy庫，如下文使用的那樣，我們需要匯入numpy庫並重命名為np

import numpy as np

零、程式設計理念的對比

0.1 程式設計正規化

matlab是面向過程的程式設計方式，而python既支援面向過程又支援面向物件，是一種多正規化（multi paradigms）的程式語言。因此，不難理解python程式語言中廣泛存在的以下的兩種等價實現方式：

np.dot(X, w)    # 呼叫全域性函式，面向過程的程式設計方式
X.dot(w)        # 呼叫物件的成員函式，面向物件的程式設計方式

0.2 matlab從1開始計數，python從0開始

0.2.1 對矩陣而言：

r = size(A, 1);         % 表示的是行數
c = size(A, 2);         % 表示的是列數
i = 0:size(A, 1) - 1;
A(i, :)             % error： 矩陣的下標索引必須是正整數型別(>=1)或邏輯型別(true/false)

0.2.2 對python而言：

A.shape[0]      # 行數， axis = 0
A.shape[1]      # 列數,   axis = 1
# 等價的（或者叫面向過程的）表達方式
np.shape(A)[0]
np.shape(A)[1]
# 可從0開始索引矩陣的行
r = A.shape[0]
A[0:r, :]

0.3 matlab矩陣索引用的是小括號，python是中括號

小括號對於一種支援面向物件、支援運算子過載的語言來說，具有特別的意義，過載括號運算子是仿函式實現的命門。如果瞭解c++的運算子過載機制，一個類如果過載了括號運算子，便可稱作仿函式，把一個類當做函式來使用。恰好，python也支援運算子過載。

class Prob(object):
    def __init__(self, lhs):
        self.lhs = lhs
    def __call__(self, rhs):
        return self.lhs * rhs
    # def __call__(self, lhs, rhs)
    #   return lhs * rhs
if __name__ == '__main__':
    p = Prob(2)     # 呼叫的是類建構函式，也即__init__
    print(p(5))     # 呼叫的是例項的括號運算子，也即__call__

0.4 切片的端點值

matlab中的切片（也即：表示式），是包含兩個端點值的，也即是一個閉區間

1:5     % 1, 2, 3, 4, 5

而python中的切片是一個左閉右開的區間：

A[0:2, :]   # 索引的是第零行和第一行，而不包括第二行

0.5 維度的順序

% matlab
>> zeros(2, 3, 4)
ans(:,:,1) =

   0   0   0
   0   0   0

ans(:,:,2) =

   0   0   0
   0   0   0

ans(:,:,3) =

   0   0   0
   0   0   0

ans(:,:,4) =

   0   0   0
   0   0   0

# python
>>> np.zeros((2, 3, 4))
array([[[ 0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.]],

       [[ 0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.]]])

所以如果一幅256×256的彩色影象在matlab中的shape為256×256×3，在python中的shape為3×256×256。

一、矩陣的基本操作

1.1 建立矩陣

1.1.1 matlab/octave

>> A = [1, 2, 3; 4, 5, 6; 7, 8, 9]
> A =
     1     2     3
     4     5     6
     7     8     9

分號表示一行的結束

1.1.2 python

A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
A = 
    array([[1, 2, 3],
           [4, 5, 6], 
           [7, 8, 9]])

以Python的基本資料型別list作為np.array的引數，生成numpy.ndarray，也即多維陣列

1.2 獲取矩陣的維度

1.2.1 matlab/octave

[r, c] = size(A)    % r = 3, c = 3
size(A, 1)          % 3 
size(A, 2)          % 3
size(A, 3)          % 1
numel(A)            % 返回的是A的元素個數，即A的全部軸的乘積

1.2.2 python

shapes = A.shape    # 是tuple型別的返回值

而：

A.size          # 返回的是A的元素數，即行*列
A.size == np.prod(A.shape)

1.3 索引矩陣的行和列

1.3.1 matlab/octave

A(1, :)     % A的第一行
A(1:2, :)   % A的前兩行
A(:, 1)     % A的第一列
A(:, 1:2)   % A的前兩列
A(end, :)   % A的最後一行

1.3.2 python

A[0, :]     # 第一行
A[0:2, :]   # 第一行，第二行
A[:, 0]     # 第一列
A[:, 0:2]   # 第一列，第二列
A[-1, :]        # 最後一行

1.4 通過斷言（predicate）提取矩陣的行和列

1.4.1 matlab/octave

mod(A, 4) == 0
    0   0   0
    1   0   0
    0   1   0
A(mod(A， 4) == 0)   % 返回矩陣中是4的倍數的元素，

因取模運算元%在矩陣中有著特殊的意義（註釋），故取模運算用了一個函式來替代。

1.4.2 python

A%4 == 0
    array([[False, False, False],
           [ True, False, False],
           [False,  True, False]], dtype=bool)

A[A%4==0]           % 返回矩陣中是4的倍數的元素

1.5 獲取特定位置上的元素

以第一個元素為例：

1.5.1 matlab/octave

A(1, 1)

1.5.2 python

A[0, 0]

1.6 向量的操作

1.6.1 matlab

a = [1; 2; 3]   % 建立列向量
b = [1, 2, 3]   % 建立行向量, b = [1 2 3]
b = [1 2 3]'    % 行轉換為列

1.6.2 python

a = np.array([[1], [2], [3]])
b = np.array([1, 2, 3])
b = b[:, np.newaxis]    # 或者 b = b[np.newaxis].T

1.7 矩陣變維

將一個3*3矩陣轉換為一個行向量

1.7.1 matlab/octave

B = reshape(A, 1, numel(A))

1.7.2 python

B = A.reshape(1, A.size)    # array([[1, 2, 3, 4, 5, 6, 7, 8, 9]])
B = A.ravel()               # array([1, 2, 3, 4, 5, 6, 7, 8, 9])

1.8 矩陣拼接

1.8.1 matlab/octave

A(3, :) = []    % 移除矩陣的某一行（第三行），可見matlab中矩陣操作之強大
A = 
    1   2   3
    4   5   6

B = [7, 8, 9; 10, 11, 12]
C = [A; B]      % 列方向上的拼接
    1   2   3
    4   5   6
    7   8   9
    10  11  12

C = [A B]       % 行方向上的拼接
    1   2   3   7   8   9
    4   5   6   10  11  12

1.8.2 python

A = A[0:2, :]       # 移除矩陣的第三行
B = np.array([[7, 8, 9], [10, 11, 12]])
C = np.concatenate((A, B), axis=0)  # 列方向的拼接
C = np.concatenate((A, B), axis=1)  # 行方向的拼接

機器學習的演算法實踐的過程中，如單層神經網路，支援向量機或者神經網路，常常會遇到對輸入X從d維增廣到1+d維的情況（每個輸入的第一維值為1）：

N = 100, d = 5
X = np.random.randn(N, d)
x0 = np.ones((X.shape[0], 1))
X = np.concatenate((x0, X), axis=1)

1.9 向量（矩陣）層疊

1.9.1 matlab/octave

a = [1, 2, 3]
b = [4, 5, 6]
c = [a', b']
c = [a; b]

1.9.2 python

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
np.c_[a, b]                 # c_: column
    array([[1, 4],
       [2, 5],
       [3, 6]])

np.r_[a, b]                 # r_: row
    array([1, 2, 3, 4, 5, 6])

二、幾種特殊矩陣

2.1 隨機矩陣（均勻分佈）

2.1.1 matlab/octave

rand(3, 2)

2.1.2 python

np.random.rand(3, 2)

無論是matlab中的rand()函式，還是numpy.random中的rand()函式，生成的隨機數都是服從[0-1]的均勻分佈（uniformed distributed），如何生成任意區間的隨機數呢？例[2, 5]區間上的隨機數 2*rand()+3

2.2 0矩陣, 全一矩陣，單位矩陣

符合大小的0矩陣，常常用以申請空間，初始化矩陣, 預先分配記憶體，提高執行的速度

2.2.1 matlab/octave

zeros(3, 2)
ones(3, 2)
eye(3)

2.2.2 python

np.zeros((3, 2))
np.ones((3, 2))
np.eye(3)

2.3 對角矩陣

2.3.1 matlab/octave

a = [1, 2, 3]
diag(a)
diag(diag(a)    % [1, 2, 3]

2.3.2 python

a = np.array([1, 2, 3])
np.diag(a)
np.diag(np.diag(a))     % array([1, 2, 3])

三、矩陣運算

3.1 矩陣與標量運算

3.1.1 matlab/octave

A = [1, 2, 3; 4, 5, 6; 7, 8, 9]
A*2
A+2
A-2
A/2         % 最後統統轉換為A中的每一個元素（element-wise）同標量的運算

3.1.2 python

A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
A*2
A/2
A-2
A+2

3.2 矩陣與矩陣運算

無論是matlab中的乘法運算（*），還是numpy中的np.dot()運算，本質上執行的都是矩陣的乘法運算，都需要滿足矩陣A的列數等於矩陣B的行數。

3.2.1 matlab/octave

A = [1, 2, 3; 4, 5, 6]  % 2*3
B = [1, 2; 3, 4; 5, 6]  % 3*2
A*B

3.2.2 python

A = np.array([[1, 2, 3], [4, 5, 6]])
B = np.array([[1, 2], [3, 4], [5, 6]])
A.dot(B)    # np.dot(A, B)

3.3 矩陣與向量運算

矩陣與向量的乘積，是矩陣與矩陣乘積的特例

3.4 按位矩陣與矩陣運算

3.4.1 matlab/octave

A.*A
A.+A
A.-A
A./A
A.^A

3.4.2 python

A*A
A+A
A-A
A/A

3.5 矩陣元素的冪乘

3.5.1 matlab/octave

A.^2

3.5.2 python

np.power(A, 2)

3.6 矩陣的冪乘

方陣才有冪乘運算

3.6.1 matlab/octave

A^2

3.6.2 python

np.linalg.matrix_power(A, 2)

3.7 矩陣轉置

3.7.1 matlab/octave

A'

3.7.2 python

A.T

3.8 矩陣行列式

3.8.1 matlab/octave

det(A)

3.8.2 python

np.linalg.det(A)

3.9 矩陣求逆

3.9.1 matlab/octave

inv(A)

3.9.10 python

A_inv = np.linalg.inv(A)
assert(np.dot(A, A_inv).all() == (np.eye(2)).all())

3.10 計算矩陣協方差矩陣

協方差矩陣刻畫的是屬性間的關係，標準協方差矩陣的求法：

協方差(i,j)=（第i列所有元素-第i列均值）*（第j列所有元素-第j列均值）/（樣本數-1）(分母減去1是為了實現無偏估計)

3.10.1 matlab/octave

>> x1 = [4.0000 4.2000 3.9000 4.3000 4.1000]’
>> x2 = [2.0000 2.1000 2.0000 2.1000 2.2000]'
>> x3 = [0.60000 0.59000 0.58000 0.62000 0.63000]’
>> cov( [x1,x2,x3] )
ans =
   2.5000e-02   7.5000e-03   1.7500e-03
   7.5000e-03   7.0000e-03   1.3500e-03
   1.7500e-03   1.3500e-03   4.3000e-04

3.10.1 python

>> x1 = np.array([ 4, 4.2, 3.9, 4.3, 4.1])
>> x2 = np.array([ 2, 2.1, 2, 2.1, 2.2])
>> x3 = np.array([ 0.6, 0.59, 0.58, 0.62, 0.63])
>> np.cov([x1, x2, x3])
array([[ 0.025  ,  0.0075 ,  0.00175],
       [ 0.0075 ,  0.007  ,  0.00135],
       [ 0.00175,  0.00135,  0.00043]])

3.11 計算特徵值和特徵向量

3.11.1 matlab/octave

[eig_vec, eig_val] = eig(A)

3.11.2 python

eig_val, eig_vec = np.linalg.eig(A)

3.12 生成高斯分佈的資料集

3.12.1 matlab/octave

mu = [0, 0]
cov = [2, 0; 0, 2]
X = mvnrnd(mu, cov, 1000)

3.12.2 python

mu = np.array([0, 0])
cov = np.array([[2, 0], [0, 2]])
X = np.random.multivariate_normal(mu, cov, 100)

四、 numpy中的matrix與array對比

雖然numpy也存在matrix型別，且matrix語法更類似與matlab中的矩陣運算，而numpy中的陣列，進行矩陣運算時，和matlab中的用法存在較大的差異，比如matlab中的*表示矩陣相乘，numpy中的陣列卻表示按位相乘。但絕大多數人仍然推薦numpy中的陣列，因為numpy中的絕大多數函式的返回型別都是numpy.ndarray()，如果堅持使用numpy中的matrix型別的話，需要進行繁瑣的型別轉換(np.mat(A))。