import numpy as np import pandas as pd一、NumPy矩阵运算基础
在机器学习基础阶段，需要掌握的矩阵及线性代数基本理论包括：
矩阵的形变及特殊矩阵的构造方法：包括矩阵的转置、对角矩阵的创建、单位矩阵的创建、上/下三角矩阵的创建等；矩阵的基本运算：包括矩阵乘法、向量内积、矩阵和向量的乘法等；矩阵的线性代数运算：包括矩阵的迹、矩阵的秩、逆矩阵的求解、伴随矩阵和广义逆矩阵等；矩阵分解运算：特征分解、奇异值分解和SVD分解等。补充知识：
1.NumPy中的矩阵表示
在NumPy中，二维数组(array)和matrix类型对象都可以用于表示矩阵，并且也都具备矩阵的代数学方法。
利用数组创建矩阵
A = np.array([[1, 2], [1, 1]]) A
array([[1, 2],
[1, 1]])
type(A)
numpy.ndarray
利用mat创建矩阵
AM = np.mat(A) AM
matrix([[1, 2],
[1, 1]])
type(AM)
numpy.matrix
关于两种对象类型的选取，此处进行简单说明：
NumPy中的matrix类型对象和MATLAB中的matrix类型等价，和NumPy中数组类型对象底层基本结构不同；在NumPy中，针对大规模数据，数组类型对象的计算速度要快于矩阵类型对象；矩阵类型对象可以通过运算符直接进行矩阵乘法，而二维数组要进行矩阵乘法（及其他矩阵运算），则必须要使用包括linalg（线性代数运算）模块在内的相关函数。
AM * AM
matrix([[3, 4],
[2, 3]])
A.dot(A)
array([[3, 4],
[2, 3]])
# 新版NumPy也支持使用符号进行矩阵乘法 A @ A
array([[3, 4],
[2, 3]])
为了执行更高效的计算、以及确保代码整体基本对象类型统一，课程如无说明，将统一使用二维数组表示矩阵。
2.NumPy中特殊矩阵构造方法在实际线性代数运算过程中，经常涉及一些特殊矩阵，如单位矩阵、对角矩阵等，相关创建方法如下：
| 函数 | 描述 |
|---|---|
| a.T | 数组a转置 |
| np.eye(n) | 创建包含n个分量的单位矩阵 |
| np.diag(a1) | 以a1中各元素，创建对角矩阵 |
| np.triu(a) | 取矩阵a中的上三角矩阵 |
| np.tril(a) | 取矩阵a中的下三角矩阵 |
# 创建一个2*3的矩阵 a1 = np.arange(1, 7).reshape(2, 3) a1
array([[1, 2, 3],
[4, 5, 6]])
# 转置 a1.T
array([[1, 4],
[2, 5],
[3, 6]])
矩阵的转置就是每个元素行列位置互换
# 创建单位矩阵 np.eye(3)
array([[1., 0., 0.],
[0., 1., 0.],
[0., 0., 1.]])
单位矩阵之所以被称为“单位”，核心原因在于单位矩阵和任何矩阵相乘，都将返回原矩阵。
a = np.arange(5) a
array([0, 1, 2, 3, 4])
np.diag(a)
array([[0, 0, 0, 0, 0],
[0, 1, 0, 0, 0],
[0, 0, 2, 0, 0],
[0, 0, 0, 3, 0],
[0, 0, 0, 0, 4]])
# 对角线向上偏移一位 np.diag(a, 1)
array([[0, 0, 0, 0, 0, 0],
[0, 0, 1, 0, 0, 0],
[0, 0, 0, 2, 0, 0],
[0, 0, 0, 0, 3, 0],
[0, 0, 0, 0, 0, 4],
[0, 0, 0, 0, 0, 0]])
# 对角线向下偏移一位 np.diag(a, -1)
array([[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0],
[0, 0, 2, 0, 0, 0],
[0, 0, 0, 3, 0, 0],
[0, 0, 0, 0, 4, 0]])
a1 = np.arange(9).reshape(3, 3) a1
array([[0, 1, 2],
[3, 4, 5],
[6, 7, 8]])
# 取上三角矩阵 np.triu(a1)
array([[0, 1, 2],
[0, 4, 5],
[0, 0, 8]])
# 上三角矩阵向左下偏移一位 np.triu(a1, -1)
array([[0, 1, 2],
[3, 4, 5],
[0, 7, 8]])
# 上三角矩阵向右上偏移一位 np.triu(a1, 1)
array([[0, 1, 2],
[0, 0, 5],
[0, 0, 0]])
# 下三角矩阵 np.tril(a1)
array([[0, 0, 0],
[3, 4, 0],
[6, 7, 8]])
3.NumPy中矩阵基本运算
由于NumPy中我们使用二维数组来表述矩阵，因此二维数组也就具备了数组和矩阵的两重属性。其中数组属性决定的基本运算相对简单，基础运算（如加减乘除）就是对应位置元素进行逐元素计算，而矩阵属性决定的运算则稍显复杂，当然矩阵的相关线性代数运算将在下一小节讨论，在基础运算上，矩阵和数组核心的区别在于乘法运算。
当然，从另一个角度考虑，其实对于向量和矩阵这种具备一定结构的对象，有很多种容易混淆的计算规则。对于常用的计算规则，我们通过将其划分成三类以帮助大家理解：
| 描述 | 解释/函数 |
|---|---|
| 逐元素相乘 | 向量、矩阵通用 |
| 每个对应位置元素相乘 | * |
| 逐元素相乘后相加 | 也被称为点积（内积），向量，矩阵通用 |
| 向量点积 | vdot、dot、inner |
| 矩阵点积 | vdot |
| 矩阵乘法 | 代数学意义的矩阵相乘 |
| 矩阵乘法 | dot、matmul、@ |
* :逐元素相乘
- matrix和np.array的区别
a = np.mat(np.array([[-3, 0,1,5], [2, -1,4,7], [1,3,0,6]])) a
matrix([[-3, 0, 1, 5],
[ 2, -1, 4, 7],
[ 1, 3, 0, 6]])
b = np.mat(np.array([[7,-1,0],[-2,4,0],[0,5,3],[1,-3,8]])) b
matrix([[ 7, -1, 0],
[-2, 4, 0],
[ 0, 5, 3],
[ 1, -3, 8]])
a*b
matrix([[-16, -7, 43],
[ 23, -7, 68],
[ 7, -7, 48]])
a = np.array([[-3, 0,1,5], [2, -1,4,7], [1,3,0,6]]) b = np.array([[7,-1,0],[-2,4,0],[0,5,3],[1,-3,8]]) a*b
--------------------------------------------------------------------------- ValueError Traceback (most recent call last)in 1 a = np.array([[-3, 0,1,5], [2, -1,4,7], [1,3,0,6]]) 2 b = np.array([[7,-1,0],[-2,4,0],[0,5,3],[1,-3,8]]) ----> 3 a*b ValueError: operands could not be broadcast together with shapes (3,4) (4,3)
a = np.arange(4) a
array([0, 1, 2, 3])
a * a
array([0, 1, 4, 9])
A = a.reshape(2, 2) A
array([[0, 1],
[2, 3]])
A * A
array([[0, 1],
[4, 9]])
向量点积
所谓点积（也被称为内积），指的是向量或矩阵对应位置元素相乘后相加。向量点积有三种实现方法，分别是dot、vdot和ineer。
a
array([0, 1, 2, 3])
np.dot(a, a)
14
a.dot(a)
14
(a * a).sum()
14
np.vdot(a, a)
14
np.inner(a, a)
14
矩阵点积
值得注意的是，矩阵内积只有vdot一种方式实现。
A
array([[0, 1],
[2, 3]])
np.vdot(A, A)
14
(A * A).sum()
14
注意，高维数组的inner并不是内积，而是一种类似tensordot的沿着尾轴实现和积的计算过程，该方法并不通用，此处暂不做介绍。
矩阵乘法
NumPy中，我们可以使用诸多方法实现矩阵乘法，包括dot、@、matmul等。
a1 = np.arange(1, 7).reshape(2, 3) a1
array([[1, 2, 3],
[4, 5, 6]])
a2 = np.arange(1, 10).reshape(3, 3) a2
array([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
# 矩阵乘法 np.matmul(a1, a2)
array([[30, 36, 42],
[66, 81, 96]])
此处也简单回顾矩阵乘法运算，上述相乘过程如下所示：
4.NumPy中矩阵代数运算值得注意的是，矩阵相乘要求左乘矩阵列数和右乘矩阵行数相同，而内积计算过程则严格要求两个向量/矩阵形状完全一致。
如果说矩阵的基本运算是矩阵基本性质，那么矩阵的线性代数运算，则是我们利用矩阵数据类型在求解实际问题过程中经常涉及到的线性代数方法，具体相关函数如下：
| 函数 | 描述 |
|---|---|
| np.trace(A) | 矩阵的迹 |
| np.linalg.matrix_rank(A) | 矩阵的秩 |
| np.linalg.det(A) | 计算矩阵A的行列式 |
| np.linalg.inv(A) | 矩阵求逆 |
同时，由于线性代数所涉及的数学基础知识较多，从实际应用的角度出发，我们将有所侧重的介绍实际应用过程中需要掌握的相关内容，并通过本节末尾的实际案例，来加深线性代数相关内容的理解。
NumPy中的linalg是linear algebra（线性代数）的简写，也是NumPy中保存线性代数相关计算函数的模块。
矩阵的迹（trace）
矩阵的迹的运算相对简单，就是矩阵对角线元素之和，在NumPy中，可以使用trace函数进行计算。
A = np.array([[1, 2], [4, 5]]) A
array([[1, 2],
[4, 5]])
np.trace(A)
6
当然，对于矩阵的迹来说，计算过程不需要是方正
B = np.arange(1, 7).reshape(2, 3) B
array([[1, 2, 3],
[4, 5, 6]])
无[3,3]的情况下，默认取[1,1]和[2,2]
np.trace(B)
6
矩阵的秩(rank)
矩阵的秩（rank），是指矩阵中行或列的极大线性无关数，且矩阵中行、列极大无关数总是相同的，任何矩阵的秩都是唯一值，满秩指的是方阵（行数和列数相同的矩阵）中行数、列数和秩相同，满秩矩阵有线性唯一解等重要特性，而其他矩阵也能通过求解秩来降维，同时，秩也是奇异值分解等运算中涉及到的重要概念。
所谓线性相关，其实也就是线性表示，如果 y = w x + b y=wx+b y=wx+b，我们则称y可以由x线性表示，二者线性相关，反之则线性无关。类似，如果 y = w 1 x 1 w 2 x 2 + b y=w_1x_1w_2x_2+b y=w1x1w2x2+b，则我们称y可以由 x 1 、 x 2 x_1、x_2 x1、x2线性表示，y与 x 1 、 x 2 x_1、x_2 x1、x2线性相关。
matrix_rank计算矩阵的秩
A = np.array([[1, 3, 4], [2, 1, 3], [1, 1, 2]]) A
array([[1, 3, 4],
[2, 1, 3],
[1, 1, 2]])
np.linalg.matrix_rank(A)
2
对于矩阵A来说，第三列明显可以由第一列和第二列相加得出，因此极大线性无关组只有两列。
B = np.array([[1, 3, 4], [2, 1, 3], [1, 1, 10]]) B
array([[ 1, 3, 4],
[ 2, 1, 3],
[ 1, 1, 10]])
np.linalg.matrix_rank(B)
3
矩阵的行列式(det)
所谓行列式，我们可以简单将其理解为矩阵的一个基本性质或者属性，通过行列式的计算，我们能够知道矩阵是否可逆，从而可以进一步求解矩阵所对应的线性方程。当然，更加专业的解释，行列式的作为一个基本数学工具，实际上就是矩阵进行线性变换的伸缩因子。
对于任何一个n维方正，行列式计算过程如下：
更为简单的情况，如果对于一个2*2的矩阵，行列式的计算就是主对角线元素之积减去另外两个元素之积
A = np.array([[1, 2], [4, 5]]) A
array([[1, 2],
[4, 5]])
np.linalg.det(A)
-2.9999999999999996
A的秩计算过程如下：
对于行列式的计算，要求矩阵必须是方阵，也就是行列数必须一致。
B = np.arange(1, 7).reshape(2, 3) B
array([[1, 2, 3],
[4, 5, 6]])
np.linalg.det(B)
--------------------------------------------------------------------------- LinAlgError Traceback (most recent call last)in ----> 1 np.linalg.det(B) <__array_function__ internals> in det(*args, **kwargs) ~AppDataRoamingPythonPython37site-packagesnumpylinalglinalg.py in det(a) 2153 a = asarray(a) 2154 _assert_stacked_2d(a) -> 2155 _assert_stacked_square(a) 2156 t, result_t = _commonType(a) 2157 signature = 'D->D' if isComplexType(t) else 'd->d' ~AppDataRoamingPythonPython37site-packagesnumpylinalglinalg.py in _assert_stacked_square(*arrays) 201 m, n = a.shape[-2:] 202 if m != n: --> 203 raise LinAlgError('Last 2 dimensions of the array must be square') 204 205 def _assert_finite(*arrays): LinAlgError: Last 2 dimensions of the array must be square
A矩阵不满秩，所以他存在相关性
A = np.array([[1, 3, 4], [2, 1, 3], [1, 1, 2]]) A
array([[1, 3, 4],
[2, 1, 3],
[1, 1, 2]])
np.linalg.det(A)
0.0
矩阵的逆
对于满秩的方正来说，可以求其逆矩阵。从基本定义上来看，如果矩阵B和矩阵A相乘能够得到单位矩阵，即：
B ⋅ A = E B cdot A = E B⋅A=E
则称B为A的逆矩阵，也可将B写作 A − 1 A^{-1} A−1。当然，逆矩阵的性质是相互的，我们也可称A为B的逆矩阵，或者A和B互为逆矩阵。
A = np.array([[1, 1], [3, 1]]) A
array([[1, 1],
[3, 1]])
np.linalg.det(A)
-2.0000000000000004
然后使用inverse函数进行逆矩阵求解
np.linalg.inv(A)
array([[-0.5, 0.5],
[ 1.5, -0.5]])
简单试探逆矩阵的基本特性
A.dot(np.linalg.inv(A))
array([[1.00000000e+00, 1.11022302e-16],
[2.22044605e-16, 1.00000000e+00]])
当然，对于逆矩阵，还有很多其他理解角度。例如，从方程组求解角度来看，逆矩阵的存在就代表着方程组存在唯一解，并且逆矩阵本身也是方程组求解的关键；从矩阵分解角度来看，逆矩阵是一种最为基础的矩阵分解的形式。关于这些相关内容，我们都将在后续课程中逐渐介绍。
二、矩阵方程与向量求导方法另外，在本节内容中，我们还将介绍solve（方程组求解）、lstsq（最小二乘法）相关函数的使用。
在铺垫了基本矩阵和线性代数相关知识后，接下来，我们尝试将Lesson 1中的方程组表示形式转化为矩阵表示形式，并借助矩阵方法进行相关方程的求解。并且，在Lesson 1中，我们已经简单讨论了最小二乘法这一优化算法的基本思路，最小二乘法一个最基础的优化算法，无论是其背后的数学推导还是实际应用，都值得继续深究。因此，本节开始，我们先从矩阵方程入手，先进行矩阵运算的相关方法的回顾、以及进行矩阵求导方法的讲解，再从一个更加严谨的数学角度出发，讨论最小二乘法的基本原理。
1.方程组求解与矩阵方程求解
至此我们就将方程组转化为了矩阵方程，并且，借助矩阵运算，我们可以直接在矩阵方程中对参数向量 X X X进行求解。首先我们利用NumPy基础知识，通过创建二维张量去表示上述矩阵方程中的A和B
A = np.array([[20, 8], [8, 4]]) A
array([[20, 8],
[ 8, 4]])
B = np.array([[28, 12]]).T B
array([[28],
[12]])
注，此时B也是二维张量，可以使用矩阵乘法。
B.ndim
2
然后通过行列式计算结果，简单验证A是否满秩：
np.linalg.matrix_rank(A)
2
当然，也可以通过观察A的行列式计算结果是否为0，来判断A是否满秩
np.linalg.det(A)
15.999999999999991
对于满秩矩阵，我们可以求其逆矩阵
np.linalg.inv(A).dot(A)
array([[1., 0.],
[0., 1.]])
然后在矩阵方程左右两端同时左乘其逆矩阵，即可解出X的取值
A − 1 A X = A − 1 B A^{-1}AX=A^{-1}B A−1AX=A−1B
X = A − 1 B X=A^{-1}B X=A−1B
np.matmul(np.linalg.inv(A), B)
array([[1.],
[1.]])
# 也可以使用dot方法，对于二维数组，dot就是执行矩阵乘法 np.linalg.inv(A).dot(B)
array([[1.],
[1.]])
即 X = [ w b ] = [ 1 1 ] X = left [begin{array}{cccc} w \ b \ end{array}right] =left [begin{array}{cccc} 1 \ 1 \ end{array}right] X=[wb]=[11]
当然，此外，NumPy中还提供了一种求解矩阵方程的函数，类似于上述 A ∗ X T = B A*X^T=B A∗XT=B的矩阵方程，我们还可以通过下式进行求解：
np.linalg.solve(A, B)
array([[1.],
[1.]])
2.向量求导运算
由于在编程实践层面上更倾向于使用矩阵/向量而不是方程组的形式进行计算，因此包括最小二乘法在内的一系列优化方法和算法的理论讲解，我们也将采用矩阵/向量作为基本数据结构进行概念讲解和数学公式推导。在正式讲解最小二乘法的数学原理之前，我们需要先补充一些关于向量求导的相关知识。
2.2 常见向量求导公式在前期学习中，数学理论推导涉及到的求导以向量变元求导居多，因此，除了掌握基本的向量求导方法以外，我们还需要推导几个常用的向量求导公式，在这些公式中，向量求导结果都能通过一些已经定义的向量简洁表示。
关于矩阵函数和矩阵方程二者概念的区别：
矩阵方程：指变量为矩阵的方程；矩阵函数：同函数矩阵，指自变量和因变量都是n阶矩阵的函数，也可以简单理解成由函数构成的矩阵，并且每个函数的变量都是矩阵。 三、最小二乘法的推导过程及使用方法
从数学角度讨论最小二乘法的基本理论，并见尝试简单实现最小二乘法求解损失函数的一般过程。
1.模型及方程组的矩阵形式改写
其中
w ^ hat w w^：方程系数所组成的向量，并且我们将自变量系数和截距放到了一个向量； x ^ hat x x^：方程自变量和1共同组成的向量； X ^ hat X X^：样本数据特征构成的矩阵，并在最后一列添加一个全为1的列； y y y：样本数据标签所构成的列向量; y ^ hat y y^：预测值的列向量。 2.构造损失函数
在方程组的矩阵表示基础上，我们可以以SSE作为损失函数基本计算流程构建关于 w ^ hat w w^的损失函数：
S S E L o s s ( w ^ ) = ∣ ∣ y − X w ^ ∣ ∣ 2 2 = ( y − X w ^ ) T ( y − X w ^ ) SSELoss(hat w) = ||y - Xhat w||_2^2 = (y - Xhat w)^T(y - Xhat w) SSELoss(w^)=∣∣y−Xw^∣∣22=(y−Xw^)T(y−Xw^)
需要补充两点基础知识：
向量的2-范数计算公式
上式中，
∣
∣
y
−
X
w
^
T
∣
∣
2
||y - Xhat w^T||_2
∣∣y−Xw^T∣∣2为向量的2-范数的计算表达式。向量的2-范数计算过程为各分量求平方和再进行开平方。例如
a
=
[
1
,
−
1
,
]
a=[1, -1,]
a=[1,−1,]，则
∣
∣
a
∣
∣
2
=
1
2
+
(
−
1
)
2
=
2
||a||_2= sqrt{1^2+(-1)^2}=sqrt{2}
∣∣a∣∣2=12+(−1)2
=2
。
向量的1-范数为各分量绝对值之和。值得注意的是，矩阵也有范数的概念，不过矩阵的范数计算要比向量复杂得多。
2-范数计算转化为内积运算
向量的2-范数计算结果其实就是向量（将其是做矩阵）的交叉乘积计算结果后开平方。例如，
a
=
[
1
,
−
1
]
a=[1, -1]
a=[1,−1]，则
a
a
a的交叉乘积为
a
⋅
a
T
=
[
1
,
−
1
]
⋅
[
1
−
1
]
=
2
a cdot a^T = [1, -1] cdot left [begin{array}{cccc} 1 \ -1 \ end{array}right]=2
a⋅aT=[1,−1]⋅[1−1]=2，开平方后等于其2-范数计算结果。
3.最小二乘法求解损失函数的一般过程
代码实现过程：
X = np.array([[1, 1], [3, 1]]) X
array([[1, 1],
[3, 1]])
y = np.array([2, 4]).reshape(2, 1) y
array([[2],
[4]])
np.linalg.inv(X.T.dot(X)).dot(X.T.dot(y))
array([[1.],
[1.]])
np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
array([[1.],
[1.]])
即解得 w = 1 , b = 1 w=1,b=1 w=1,b=1，即模型方程为 y = x + 1 y=x+1 y=x+1。至此，我们即完成了最小二乘法的推导以及简单实现。其中，相比于记忆最终最小二乘法求解损失函数的计算公式，更加重要的是理解并掌握矩阵表示方程组的方法以及矩阵求导计算方法。
最小二乘法计算函数
当然，除了借助NumPy中矩阵运算的函数进行最小二乘法计算以外，NumPy中也有独立的用于计算最小二乘法的函数：np.linalg.lstsq。通过该方法，我们能够在直接输入特征矩阵和标签数组的情况下进行最小二乘法的计算，并在得出最小二乘法计算结果的同时，也得到一系列相应的统计指标结果。
np.linalg.lstsq(X, y, rcond=-1)
(array([[1.],
[1.]]),
array([], dtype=float64),
2,
array([3.41421356, 0.58578644]))
其中，返回的元组中第一个元素是最小二乘法计算结果，第二个元素是SSE的计算结果，第三个元素是矩阵X的秩，最后一个结果是矩阵X的奇异值。
