Python 机器学习——线性代数和矩阵运算：从matlab迁移到python

诚然，没有一门语言能够撼动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)

)。

关于矩阵运算从python到matlab的迁移，更详细的信息请见：http://wiki.scipy.org/NumPy_for_Matlab_Users