Lesson 2.矩阵运算基础、矩阵求导与最小二乘法

Lesson 2.矩阵运算基础、矩阵求导与最小二乘法

在Lesson 1中，我们介绍了关于机器学习的一般建模流程，并且在基本没有数学公式和代码的情况下，简单介绍了关于线性回归的一般实现形式。不过这只是在初学阶段、为了不增加基础概念理解难度所采取的方法，但所有的技术最终都是为了解决实际问题的，因此，接下来，我们就在之前的基础上更进一步，从一个更加严谨的理论体系出发、来尝试进行一种更加贴合实际应用所采用的一般方法的建模方法的学习。

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

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中特殊矩阵构造方法

在实际线性代数运算过程中，经常涉及一些特殊矩阵，如单位矩阵、对角矩阵等，相关创建方法如下：

# 创建一个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)
# 上三角矩阵向右上偏移一位
np.triu(a1, 1)
# 上三角矩阵向右上偏移一位
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中我们使用二维数组来表述矩阵，因此二维数组也就具备了数组和矩阵的两重属性。其中数组属性决定的基本运算相对简单，基础运算（如加减乘除）就是对应位置元素进行逐元素计算，而矩阵属性决定的运算则稍显复杂，当然矩阵的相关线性代数运算将在下一小节讨论，在基础运算上，矩阵和数组核心的区别在于乘法运算。

当然，从另一个角度考虑，其实对于向量和矩阵这种具备一定结构的对象，有很多种容易混淆的计算规则。对于常用的计算规则，我们通过将其划分成三类以帮助大家理解：

*:逐元素相乘

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。

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中矩阵代数运算

如果说矩阵的基本运算是矩阵基本性质，那么矩阵的线性代数运算，则是我们利用矩阵数据类型在求解实际问题过程中经常涉及到的线性代数方法，具体相关函数如下：

同时，由于线性代数所涉及的数学基础知识较多，从实际应用的角度出发，我们将有所侧重的介绍实际应用过程中需要掌握的相关内容，并通过本节末尾的实际案例，来加深线性代数相关内容的理解。

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]])
np.trace(B)
#6

矩阵的秩(rank)
  矩阵的秩（rank），是指矩阵中行或列的极大线性无关数，且矩阵中行、列极大无关数总是相同的，任何矩阵的秩都是唯一值，满秩指的是方阵（行数和列数相同的矩阵）中行数、列数和秩相同，满秩矩阵有线性唯一解等重要特性，而其他矩阵也能通过求解秩来降维，同时，秩也是奇异值分解等运算中涉及到的重要概念。

所谓线性相关，其实也就是线性表示，如果 푦=푤푥+푏 ，我们则称y可以由x线性表示，二者线性相关，反之则线性无关。类似，如果푦=푤1푥1+푤2푥2+푏 ，则我们称y可以由 푥1、푥2 线性表示，y与 푥1、푥2 线性相关。

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)

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为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]])
#然后使用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.方程组求解与矩阵方程求解

在Lesson 1中，我们曾经利用损失函数的偏导函数方程组进行简单线性回归模型参数的求解：
∂ S S E L o s s ∂ ( w ) = 2 ( 2 − w − b ) ∗ ( − 1 ) + 2 ( 4 − 3 w − b ) ∗ ( − 3 ) = 20 w + 8 b − 28 = 0 ∂ S S E L o s s ∂ ( b ) = 2 ( 2 − w − b ) ∗ ( − 1 ) + 2 ( 4 − 3 w − b ) ∗ ( − 1 ) = 8 w + 4 b − 12 = 0 begin{aligned} frac{partial S S E L o s s}{partial(w)} &=2(2-w-b) *(-1)+2(4-3 w-b) *(-3) \ &=20 w+8 b-28 \ &=0 \ frac{partial S S E L o s s}{partial(b)} &=2(2-w-b) *(-1)+2(4-3 w-b) *(-1) \ &=8 w+4 b-12 \ &=0 end{aligned} ∂(w)∂SSELoss​∂(b)∂SSELoss​​=2(2−w−b)∗(−1)+2(4−3w−b)∗(−3)=20w+8b−28=0=2(2−w−b)∗(−1)+2(4−3w−b)∗(−1)=8w+4b−12=0​
尽管方程组求解有很多种方法，类似Lesson 1中所采用的，先通过方程变量相消法反解出一个变量（w=1），再将其带入其中任意一个式子求解出另一个变量（b=1），确实能够顺利求出方程组的解。但在如果想借助编程工具求解方程组，则需要将原先的方程组求解问题转化为矩阵方程的求解问题。例如：上述求解过程方程组为：
20 w + 8 b − 28 = 0 8 w + 4 b − 12 = 0 begin{array}{c} 20 w+8 b-28=0 \ 8 w+4 b-12=0 end{array} 20w+8b−28=08w+4b−12=0​
我们令：
A = [ 20 8 8 4 ] B = [ 28 12 ] X = [ w b ] begin{array}{c} A=left[begin{array}{cc} 20 & 8 \ 8 & 4 end{array}right] \ B=left[begin{array}{l} 28 \ 12 end{array}right] \ X=left[begin{array}{l} w \ b end{array}right] end{array} A=[208​84​]B=[2812​]X=[wb​]​
其中 푋 为参数向量。借助矩阵运算相关知识，上述方程组可等价表示为：
A ⋅ X − B = 0 A cdot X-B=0 A⋅X−B=0
即
A ⋅ X = B A cdot X=B A⋅X=B
至此我们就将方程组转化为了矩阵方程，并且，借助矩阵运算，我们可以直接在矩阵方程中对参数向量 푋 进行求解。首先我们利用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)
#array([[ 0.25, -0.5 ],
#       [-0.5 ,  1.25]])

然后在矩阵方程左右两端同时左乘其逆矩阵，即可解出X的取值
A − 1 A X = A − 1 B X = A − 1 B begin{array}{c} A^{-1} A X=A^{-1} B \ X=A^{-1} B end{array} A−1AX=A−1BX=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}{c} w \ b end{array}right]=left[begin{array}{l} 1 \ 1 end{array}right] X=[wb​]=[11​]
当然，此外，NumPy中还提供了一种求解矩阵方程的函数，类似于上述 퐴 ∗ 푋 = 퐵 퐴∗푋=퐵 A∗X=B的矩阵方程，我们还可以通过下式进行求解：

np.linalg.solve(A, B)
#array([[1.],
#       [1.]])

2.向量求导运算

由于在编程实践层面上更倾向于使用矩阵/向量而不是方程组的形式进行计算，因此包括最小二乘法在内的一系列优化方法和算法的理论讲解，我们也将采用矩阵/向量作为基本数据结构进行概念讲解和数学公式推导。在正式讲解最小二乘法的数学原理之前，我们需要先补充一些关于向量求导的相关知识。

2.1 向量求导基本方法

首先，我们先来看相对简单的向量求导方法，并借此探究对于具有一定结构化的变量进行求导运算的本质。

假设现有一个二元函数如下：
f ( x 1 , x 2 ) = 2 x 1 + x 2 fleft(x_{1}, x_{2}right)=2 x_{1}+x_{2} f(x1​,x2​)=2x1​+x2​
并且，我们可以分别对该函数中的两个变量 x 1 x_{1} x1​、 x 2 x_{2} x2​依次求偏导，可得：
∂ f ∂ x 1 = 2 ， ∂ f ∂ x 2 = 1 frac{partial f}{partial x_1} = 2 ， frac{partial f}{partial x_2} = 1 ∂x1​∂f​=2，∂x2​∂f​=1
现在我们考虑将上述求偏导的函数组改写为矩阵形式。则根据前述内容介绍，我们可以将函数中的两个变量依次排列，组成一个向量变元，即一个由多个变量所组成的向量，即
x = [ x 1 , x 2 ] T x = [x_1, x_2]^T x=[x1​,x2​]T
此时，如果我们按照向量变元内部的变量排列顺序，依次在每个变量位置填上该变量对应的偏导函数，则就构成了对于函数 푓 进行向量变元 푥 的向量求导的结果，即：
∂ f ( x ) ∂ x = [ 2 1 ] frac{partial f(x)}{partial x} = left [begin{array}{cccc} 2 \ 1 \ end{array}right] ∂x∂f(x)​=[21​]
其中，푥为向量变元。

至此，我们就完成了向量求导的基本过程。核心点在于，我们是依据向量变元中的变量排列顺序，依次填写了对应变量的偏导函数计算结果。不过，更进一步的来看，既然方程组需要改写成向量/矩阵形式，那么原始函数方程其实也同样需要改写成向量/矩阵形式。因此，原方程我们可以改写成：
f ( x ) = A T ⋅ x f(x) = A^T cdot x f(x)=AT⋅x
其中， A = [ 2 , 1 ] T ， x = [ x 1 , x 2 ] T A = [2, 1]^T ，x = [x_1, x_2]^T A=[2,1]T，x=[x1​,x2​]T
原方程为
y = 2 x 1 + x 2 y = 2x_1+x_2 y=2x1​+x2​
结合函数求导结果，我们不难发现， ∂ f ( x ) ∂ x frac{partial f(x)}{partial x} ∂x∂f(x)​最终计算结果就是 A A A，即
∂ f ( x ) ∂ x = ∂ ( A T ⋅ x ) ∂ x = A frac{partial f(x)}{partial x} = frac{partial(A^T cdot x)}{partial x} = A ∂x∂f(x)​=∂x∂(AT⋅x)​=A
其中 푥 是向量变元， 퐴 是列向量。当然，该结论也能推导至一般的情况，相关证明会在下一小节给出。此处给出向量变元的函数求导计算公式：

很多时候我们并不区分所谓向量方程和矩阵方程，一般所有自变量为向量或矩阵的方程，我们会统一称其为矩阵方程。包含向量或者矩阵的表达式，我们也会统一称其为矩阵表达式。

向量求导的定义法

设 푓(푥) 是一个关于 푥 的函数，其中 푥 是向量变元，并且
x = [ x 1 , x 2 , . . . , x n ] T x = [x_1, x_2,...,x_n]^T x=[x1​,x2​,...,xn​]T
则
∂ f ∂ x = [ ∂ f ∂ x 1 , ∂ f ∂ x 2 , . . . , ∂ f ∂ x n ] T frac{partial f}{partial x} = [frac{partial f}{partial x_1}, frac{partial f}{partial x_2}, ..., frac{partial f}{partial x_n}]^T ∂x∂f​=[∂x1​∂f​,∂x2​∂f​,...,∂xn​∂f​]T
而该表达式也被称为向量求导的梯度向量形式。
∇ x f ( x ) = ∂ f ∂ x = [ ∂ f ∂ x 1 , ∂ f ∂ x 2 , . . . , ∂ f ∂ x n ] T nabla _xf(x) = frac{partial f}{partial x} = [frac{partial f}{partial x_1}, frac{partial f}{partial x_2}, ..., frac{partial f}{partial x_n}]^T ∇x​f(x)=∂x∂f​=[∂x1​∂f​,∂x2​∂f​,...,∂xn​∂f​]T
通过求得函数的梯度向量求解向量导数的方法，也被称为定义法求解。

值得注意的是，多元函数是一定能够求得梯度向量的，但梯度向量或者说向量求导结果，能否由一些已经定义的向量解决表示，如 퐴 就是 푓(푥) 的向量求导结果，则不一定。

值得注意的是，多元函数是一定能够求得梯度向量的，但梯度向量或者说向量求导结果，能否由一些已经定义的向量解决表示，如 퐴 就是 푓(푥) 的向量求导结果，则不一定。

2.2 常见向量求导公式

在前期学习中，数学理论推导涉及到的求导以向量变元求导居多，因此，除了掌握基本的向量求导方法以外，我们还需要推导几个常用的向量求导公式，在这些公式中，向量求导结果都能通过一些已经定义的向量简洁表示。同样，此处我们假设x为包含n个变量的列向量， x = [ x 1 , x 2 , . . . , x n ] T x = [x_1, x_2,...,x_n]^T x=[x1​,x2​,...,xn​]T。

（1）
∂ a ∂ x = 0 frac{partial a}{partial x} = 0 ∂x∂a​=0
证明：
∂ a ∂ x = [ ∂ a ∂ x 1 , ∂ a ∂ x 2 , . . . , ∂ a ∂ x n ] T = [ 0 , 0 , . . . , 0 ] T frac{partial a}{partial x} = [frac{partial a}{partial x_1}, frac{partial a}{partial x_2}, ..., frac{partial a}{partial x_n}]^T = [0,0,...,0]^T ∂x∂a​=[∂x1​∂a​,∂x2​∂a​,...,∂xn​∂a​]T=[0,0,...,0]T（2）
∂ ( x T ⋅ A ) ∂ x = ∂ ( A T ⋅ x ) ∂ x = A frac{partial(x^T cdot A)}{partial x} = frac{partial(A^T cdot x)}{partial x} = A ∂x∂(xT⋅A)​=∂x∂(AT⋅x)​=A
证明：
  此时A为拥有n个分量的常数向量，设 A = [ a 1 , a 2 , . . . , a n ] T A = [a_1, a_2,...,a_n]^T A=[a1​,a2​,...,an​]T，则有
∂ ( x T ⋅ A ) ∂ x = ∂ ( A T ⋅ x ) ∂ x = ∂ ( a 1 ⋅ x 1 + a 2 ⋅ x 2 + . . . + a n ⋅ x n ) ∂ x = [ ∂ ( a 1 ⋅ x 1 + a 2 ⋅ x 2 + . . . + a n ⋅ x n ) ∂ x 1 ∂ ( a 1 ⋅ x 1 + a 2 ⋅ x 2 + . . . + a n ⋅ x n ) ∂ x 2 . . . ∂ ( a 1 ⋅ x 1 + a 2 ⋅ x 2 + . . . + a n ⋅ x n ) ∂ x n ] = [ a 1 a 2 . . . a n ] = A begin{aligned} frac{partial(x^T cdot A)}{partial x} & = frac{partial(A^T cdot x)}{partial x}\ & = frac{partial(a_1 cdot x_1 + a_2 cdot x_2 +...+ a_n cdot x_n)}{partial x}\ & = left [begin{array}{cccc} frac{partial(a_1 cdot x_1 + a_2 cdot x_2 +...+ a_n cdot x_n)}{partial x_1} \ frac{partial(a_1 cdot x_1 + a_2 cdot x_2 +...+ a_n cdot x_n)}{partial x_2} \ . \ . \ . \ frac{partial(a_1 cdot x_1 + a_2 cdot x_2 +...+ a_n cdot x_n)}{partial x_n} \ end{array}right] \ & =left [begin{array}{cccc} a_1 \ a_2 \ . \ . \ . \ a_n \ end{array}right] = A end{aligned} ∂x∂(xT⋅A)​​=∂x∂(AT⋅x)​=∂x∂(a1​⋅x1​+a2​⋅x2​+...+an​⋅xn​)​=⎣⎢⎢⎢⎢⎢⎢⎢⎡​∂x1​∂(a1​⋅x1​+a2​⋅x2​+...+an​⋅xn​)​∂x2​∂(a1​⋅x1​+a2​⋅x2​+...+an​⋅xn​)​...∂xn​∂(a1​⋅x1​+a2​⋅x2​+...+an​⋅xn​)​​⎦⎥⎥⎥⎥⎥⎥⎥⎤​=⎣⎢⎢⎢⎢⎢⎢⎡​a1​a2​...an​​⎦⎥⎥⎥⎥⎥⎥⎤​=A​（3）
∂ ( x T ⋅ x ) ∂ x = 2 x frac{partial (x^T cdot x)}{partial x} = 2x ∂x∂(xT⋅x)​=2x
证明：
∂ ( x T ⋅ x ) ∂ x = ∂ ( x 1 2 + x 2 2 + . . . + x n 2 ) ∂ x = [ ∂ ( x 1 2 + x 2 2 + . . . + x n 2 ) ∂ x 1 ∂ ( x 1 2 + x 2 2 + . . . + x n 2 ) ∂ x 2 . . . ∂ ( x 1 2 + x 2 2 + . . . + x n 2 ) ∂ x n ] = [ 2 x 1 2 x 2 . . . 2 x n ] = 2 x begin{aligned} frac{partial(x^T cdot x)}{partial x} & = frac{partial(x_1^2+x_2^2+...+x_n^2)}{partial x}\ & = left [begin{array}{cccc} frac{partial(x_1^2+x_2^2+...+x_n^2)}{partial x_1} \ frac{partial(x_1^2+x_2^2+...+x_n^2)}{partial x_2} \ . \ . \ . \ frac{partial(x_1^2+x_2^2+...+x_n^2)}{partial x_n} \ end{array}right] \ & =left [begin{array}{cccc} 2x_1 \ 2x_2 \ . \ . \ . \ 2x_n \ end{array}right] = 2x end{aligned} ∂x∂(xT⋅x)​​=∂x∂(x12​+x22​+...+xn2​)​=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎡​∂x1​∂(x12​+x22​+...+xn2​)​∂x2​∂(x12​+x22​+...+xn2​)​...∂xn​∂(x12​+x22​+...+xn2​)​​⎦⎥⎥⎥⎥⎥⎥⎥⎥⎤​=⎣⎢⎢⎢⎢⎢⎢⎡​2x1​2x2​...2xn​​⎦⎥⎥⎥⎥⎥⎥⎤​=2x​

此处 x T x x^Tx xTx也被称为向量的交叉乘积(crossprod)。

（4）
∂ ( x T A x ) x = A x + A T x frac{partial (x^T A x)}{x} = Ax + A^Tx x∂(xTAx)​=Ax+ATx
其中A是一个(n*n)的矩阵， A n ∗ n = ( a i j ) a = 1 , j = 1 n , n A_{n*n}=(a_{ij})_{a=1,j=1}^{n,n} An∗n​=(aij​)a=1,j=1n,n​
证明：
首先，
X T A X = [ x 1 , x 2 , . . . , x n ] ⋅ [ a 11 a 12 . . . a 1 n a 21 a 22 . . . a 2 n . . . . . . . . . . . . a n 1 a n 2 . . . a n n ] ⋅ [ x 1 , x 2 , . . . , x n ] T = [ x 1 a 11 + x 2 a 21 + . . . + x n a n 1 , x 1 a 12 + x 2 a 22 + . . . + x n a n 2 , . . . , x 1 a 1 n + x 2 a 2 n + . . . + x n a n n ] ⋅ [ x 1 x 2 . . . x n ] = x 1 ( x 1 a 11 + x 2 a 21 + . . . + x n a n 1 ) + x 2 ( x 1 a 12 + x 2 a 22 + . . . + x n a n 2 ) + . . . + x n ( x 1 a 1 n + x 2 a 2 n + . . . + x n a n n ) begin{aligned} X^TAX &= [x_1, x_2,...,x_n] cdot left [begin{array}{cccc} a_{11} &a_{12} &... &a_{1n}\ a_{21} &a_{22} &... &a_{2n}\ ... &... &... &... \ a_{n1} &a_{n2} &... &a_{nn}\ end{array}right] cdot [x_1, x_2,...,x_n]^T \ &=[x_1a_{11}+x_2a_{21}+...+x_na_{n1}, x_1a_{12}+x_2a_{22}+...+x_na_{n2},...,x_1a_{1n}+x_2a_{2n}+...+x_na_{nn}] cdot left [begin{array}{cccc} x_1 \ x_2 \ . \ . \ . \ x_n \ end{array}right] \ &=x_1(x_1a_{11}+x_2a_{21}+...+x_na_{n1})+x_2(x_1a_{12}+x_2a_{22}+...+x_na_{n2})+...+x_n(x_1a_{1n}+x_2a_{2n}+...+x_na_{nn}) end{aligned} XTAX​=[x1​,x2​,...,xn​]⋅⎣⎢⎢⎡​a11​a21​...an1​​a12​a22​...an2​​............​a1n​a2n​...ann​​⎦⎥⎥⎤​⋅[x1​,x2​,...,xn​]T=[x1​a11​+x2​a21​+...+xn​an1​,x1​a12​+x2​a22​+...+xn​an2​,...,x1​a1n​+x2​a2n​+...+xn​ann​]⋅⎣⎢⎢⎢⎢⎢⎢⎡​x1​x2​...xn​​⎦⎥⎥⎥⎥⎥⎥⎤​=x1​(x1​a11​+x2​a21​+...+xn​an1​)+x2​(x1​a12​+x2​a22​+...+xn​an2​)+...+xn​(x1​a1n​+x2​a2n​+...+xn​ann​)​
令 k ( x ) = x 1 ( x 1 a 11 + x 2 a 21 + . . . + x n a n 1 ) + x 2 ( x 1 a 12 + x 2 a 22 + . . . + x n a n 2 ) + . . . + x n ( x 1 a 1 n + x 2 a 2 n + . . . + x n a n n ) k(x) = x_1(x_1a_{11}+x_2a_{21}+...+x_na_{n1})+x_2(x_1a_{12}+x_2a_{22}+...+x_na_{n2})+...+x_n(x_1a_{1n}+x_2a_{2n}+...+x_na_{nn}) k(x)=x1​(x1​a11​+x2​a21​+...+xn​an1​)+x2​(x1​a12​+x2​a22​+...+xn​an2​)+...+xn​(x1​a1n​+x2​a2n​+...+xn​ann​)
则有
∂ k ( x ) ∂ x 1 = ( x 1 a 11 + x 2 a 21 + . . . + x n a n 1 ) + ( x 1 a 11 + x 2 a 12 + . . . + x n a 1 n ) frac{partial k(x)}{partial x_1} = (x_1a_{11}+x_2a_{21}+...+x_na_{n1})+ (x_1a_{11} + x_2a_{12}+...+x_na_{1n}) ∂x1​∂k(x)​=(x1​a11​+x2​a21​+...+xn​an1​)+(x1​a11​+x2​a12​+...+xn​a1n​)
类似可得：
令 k ( x ) = x 1 ( x 1 a 11 + x 2 a 21 + . . . + x n a n 1 ) + x 2 ( x 1 a 12 + x 2 a 22 + . . . + x n a n 2 ) + . . . + x n ( x 1 a 1 n + x 2 a 2 n + . . . + x n a n n ) k(x) = x_1(x_1a_{11}+x_2a_{21}+...+x_na_{n1})+x_2(x_1a_{12}+x_2a_{22}+...+x_na_{n2})+...+x_n(x_1a_{1n}+x_2a_{2n}+...+x_na_{nn}) k(x)=x1​(x1​a11​+x2​a21​+...+xn​an1​)+x2​(x1​a12​+x2​a22​+...+xn​an2​)+...+xn​(x1​a1n​+x2​a2n​+...+xn​ann​)

至此，我们完成相关向量求导常用公式的证明。不过从上面的证明不难看出，使用定义法进行公式证明往往会非常繁琐（尽管流程相对清晰）。因此，我们也会在后续补充除了定义法以外的向量乘法常用公式证明方法。

此外，对于矩阵来说，也有类似的求导方法。即如果变量以矩阵形式出现，则针对该矩阵的方程求导之际上也就是依照矩阵基本结构，在每个位置上对对应的变量分量求偏导函数。只不过相比向量，矩阵多了一个维度、结构更加复杂，因此求解过程也更加复杂。由于初期我们接触的大多是向量变元的方程，因此关于矩阵求导的常用公式推导，我们将在后续逐渐展开讨论。

最后，还需要简单进行一个概念辨析，那就是关于矩阵函数和矩阵方程二者概念的区别：

矩阵方程：指变量为矩阵的方程；矩阵函数：同函数矩阵，指自变量和因变量都是n阶矩阵的函数，也可以简单理解成由函数构成的矩阵，并且每个函数的变量都是矩阵。 三、最小二乘法的推导过程及使用方法

有了上述内容铺垫之后，接下来，我们从数学角度讨论最小二乘法的基本理论，并见尝试简单实现最小二乘法求解损失函数的一般过程。

1.模型及方程组的矩阵形式改写

首先，我们尝试对模型进行矩阵形式改写。

模型改写称矩阵表达式

首先，假设多元线性方程有如下形式
f ( x ) = w 1 x 1 + w 2 x 2 + . . . + w d x d + b f(x) = w_1x_1+w_2x_2+...+w_dx_d+b f(x)=w1​x1​+w2​x2​+...+wd​xd​+b
令 w = [ w 1 , w 2 , . . . w d ] T w = [w_1,w_2,...w_d]^T w=[w1​,w2​,...wd​]T， x = [ x 1 , x 2 , . . . x d ] T x = [x_1,x_2,...x_d]^T x=[x1​,x2​,...xd​]T，则上式可写为
f ( x ) = w T x + b f(x) = w^Tx+b f(x)=wTx+b

在机器学习领域，我们将线性回归自变量系数命名为w，其实是weight的简写，意为自变量的权重。

将带入数据后的方程组改写为矩阵方程
并且，假设现在总共有m条观测值， x ( i ) = [ x 1 ( i ) , x 2 ( i ) , . . . , x d ( i ) ] x^{(i)} = [x_1^{(i)}, x_2^{(i)},...,x_d^{(i)}] x(i)=[x1(i)​,x2(i)​,...,xd(i)​]，则带入模型可构成m个方程：
[ w 1 x 1 ( 1 ) + w 2 x 2 ( 1 ) + . . . + w d x d ( 1 ) + b w 1 x 1 ( 2 ) + w 2 x 2 ( 2 ) + . . . + w d x d ( 2 ) + b . . . w 1 x 1 ( m ) + w 2 x 2 ( m ) + . . . + w d x d ( m ) + b ] = [ y ^ 1 y ^ 2 . . . y ^ m ] left [begin{array}{cccc} w_1x_1^{(1)}+w_2x_2^{(1)}+...+w_dx_d^{(1)}+b \ w_1x_1^{(2)}+w_2x_2^{(2)}+...+w_dx_d^{(2)}+b \ . \ . \ . \ w_1x_1^{(m)}+w_2x_2^{(m)}+...+w_dx_d^{(m)}+b \ end{array}right] = left [begin{array}{cccc} hat y_1 \ hat y_2 \ . \ . \ . \ hat y_m \ end{array}right] ⎣⎢⎢⎢⎢⎢⎢⎢⎡​w1​x1(1)​+w2​x2(1)​+...+wd​xd(1)​+bw1​x1(2)​+w2​x2(2)​+...+wd​xd(2)​+b...w1​x1(m)​+w2​x2(m)​+...+wd​xd(m)​+b​⎦⎥⎥⎥⎥⎥⎥⎥⎤​=⎣⎢⎢⎢⎢⎢⎢⎡​y^​1​y^​2​...y^​m​​⎦⎥⎥⎥⎥⎥⎥⎤​
然后考虑如何将上述方程组进行改写，首先，我们可令
w ^ = [ w 1 , w 2 , . . . , w d , b ] T ， x ^ = [ x 1 , x 2 , . . . , x d , 1 ] T hat w = [w_1,w_2,...,w_d,b]^T，hat x = [x_1,x_2,...,x_d,1]^T w^=[w1​,w2​,...,wd​,b]T，x^=[x1​,x2​,...,xd​,1]T
X ^ = [ x 1 ( 1 ) x 2 ( 1 ) . . . x d ( 1 ) 1 x 1 ( 2 ) x 2 ( 2 ) . . . x d ( 2 ) 1 . . . . . . . . . . . . 1 x 1 ( m ) x 2 ( m ) . . . x d ( m ) 1 ] hat X = left [begin{array}{cccc} x_1^{(1)} &x_2^{(1)} &... &x_d^{(1)} &1 \ x_1^{(2)} &x_2^{(2)} &... &x_d^{(2)} &1 \ ... &... &... &... &1 \ x_1^{(m)} &x_2^{(m)} &... &x_d^{(m)} &1 \ end{array}right] X^=⎣⎢⎢⎢⎡​x1(1)​x1(2)​...x1(m)​​x2(1)​x2(2)​...x2(m)​​............​xd(1)​xd(2)​...xd(m)​​1111​⎦⎥⎥⎥⎤​
y = [ y 1 y 2 . . . y m ] y = left [begin{array}{cccc} y_1 \ y_2 \ . \ . \ . \ y_m \ end{array}right] y=⎣⎢⎢⎢⎢⎢⎢⎡​y1​y2​...ym​​⎦⎥⎥⎥⎥⎥⎥⎤​
y ^ = [ y ^ 1 y ^ 2 . . . y ^ m ] hat y = left [begin{array}{cccc} hat y_1 \ hat y_2 \ . \ . \ . \ hat y_m \ end{array}right] y^​=⎣⎢⎢⎢⎢⎢⎢⎡​y^​1​y^​2​...y^​m​​⎦⎥⎥⎥⎥⎥⎥⎤​
其中

w ^ hat w w^：方程系数所组成的向量，并且我们将自变量系数和截距放到了一个向量； x ^ hat x x^：方程自变量和1共同组成的向量； X ^ hat X X^：样本数据特征构成的矩阵，并在最后一列添加一个全为1的列； y y y：样本数据标签所构成的列向量; y ^ hat y y^​：预测值的列向量。

因此，上述方程组可表示为
X ^ ⋅ w ^ = y ^ hat X cdot hat w = hat y X^⋅w^=y^​

模型进一步改写

  在改写了 x ^ hat x x^和 w ^ hat w w^之后，线性模型也可以按照如下形式进行改写：
f ( x ^ ) = w ^ T ⋅ x ^ f(hat x) = hat w^T cdot hat x f(x^)=w^T⋅x^

2.构造损失函数

在方程组的矩阵表示基础上，我们可以以SSE作为损失函数基本计算流程构建关于 푤̂ 的损失函数： 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.最小二乘法求解损失函数的一般过程

在确定损失函数的矩阵表示形式之后，接下来即可利用最小二乘法进行求解。其基本求解思路仍然和Lesson 0中介绍的一样，先求导函数、再令导函数取值为零，进而解出参数取值。只不过此时求解的是矩阵方程。

在此之前，需要补充两点矩阵转置的运算规则：
( A − B ) T = A T − B T ， ( A B ) T = B T A T (A-B)^T=A^T-B^T，(AB)^T=B^TA^T (A−B)T=AT−BT，(AB)T=BTAT
接下来，对 푆푆퐸퐿표푠푠(푤) 求导并令其等于0：
S S E L o s s ( w ^ ) ∂ w ^ = ∂ ∣ ∣ y − X w ^ ∣ ∣ 2 2 ∂ w ^ = ∂ ( y − X w ^ ) T ( y − X w ^ ) ∂ w ^ = ∂ ( y T − w ^ T X T ) ( y − X w ^ ) ∂ w ^ = ∂ ( y T y − w ^ T X T y − y T X w ^ + w ^ T X T X w ^ ) ∂ w ^ = 0 − X T y − X T y + X T X w ^ + ( X T X ) T w ^ = 0 − X T y − X T y + 2 X T X w ^ = 2 ( X T X w ^ − X T y ) = 0 begin{aligned} frac{SSELoss(hat w)}{partial{boldsymbol{hat w}}} &= frac{partial{||boldsymbol{y} - boldsymbol{Xhat w}||_2}^2}{partial{boldsymbol{hat w}}} \ &= frac{partial(boldsymbol{y} - boldsymbol{Xhat w})^T(boldsymbol{y} - boldsymbol{Xhat w})}{partial{boldsymbol{hat w}}} \ & =frac{partial(boldsymbol{y}^T - boldsymbol{hat w^T X^T})(boldsymbol{y} - boldsymbol{Xhat w})}{partial{boldsymbol{hat w}}}\ &=frac{partial(boldsymbol{y}^Tboldsymbol{y} - boldsymbol{hat w^T X^Ty}-boldsymbol{y}^Tboldsymbol{X hat w} +boldsymbol{hat w^TX^T}boldsymbol{Xhat w})}{partial{boldsymbol{hat w}}}\ & = 0 - boldsymbol{X^Ty} - boldsymbol{X^Ty}+X^TXhat w+(X^TX)^That w \ &= 0 - boldsymbol{X^Ty} - boldsymbol{X^Ty} + 2boldsymbol{X^TXhat w}\ &= 2(boldsymbol{X^TXhat w} - boldsymbol{X^Ty}) = 0 end{aligned} ∂w^SSELoss(w^)​​=∂w^∂∣∣y−Xw^∣∣2​2​=∂w^∂(y−Xw^)T(y−Xw^)​=∂w^∂(yT−w^TXT)(y−Xw^)​=∂w^∂(yTy−w^TXTy−yTXw^+w^TXTXw^)​=0−XTy−XTy+XTXw^+(XTX)Tw^=0−XTy−XTy+2XTXw^=2(XTXw^−XTy)=0​
即 X T X w ^ = X T y X^TXhat w = X^Ty XTXw^=XTy
要使得此式有解，等价于 X T X X^TX XTX（也被称为矩阵的交叉乘积crossprod存在逆矩阵，若存在，则可解出
w ^ = ( X T X ) − 1 X T y hat w = (X^TX)^{-1}X^Ty w^=(XTX)−1XTy

4.最小二乘法的简单实现

接下来，我们回到最初的例子，简单尝试利用上述推导公式求解简单线性回归参数。原始数据如下：

因此利用矩阵表达式，可进行如下形式的改写：
特征矩阵： X ^ = [ 1 1 3 1 ] hat X = left [begin{array}{cccc} 1 &1 \ 3 &1 \ end{array}right] X^=[13​11​]
标签数组： y = [ 2 4 ] y = left [begin{array}{cccc} 2 \ 4 \ end{array}right] y=[24​]
参数向量： w ^ = [ w b ] hat w = left [begin{array}{cccc} w \ b \ end{array}right] w^=[wb​]
求解公式为：
w ^ = ( X T X ) − 1 X T y = ( [ 1 1 3 1 ] T [ 1 1 3 1 ] ) − 1 [ 1 1 3 1 ] T [ 2 4 ] begin{aligned} hat w &= (X^TX)^{-1}X^Ty \ &= (left [begin{array}{cccc} 1 &1 \ 3 &1 \ end{array}right]^{T} left [begin{array}{cccc} 1 &1 \ 3 &1 \ end{array}right])^{-1} left [begin{array}{cccc} 1 &1 \ 3 &1 \ end{array}right]^{T} left [begin{array}{cccc} 2 \ 4 \ end{array}right] \ end{aligned} w^​=(XTX)−1XTy=([13​11​]T[13​11​])−1[13​11​]T[24​]​
代码实现过程：

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

即解得푤=1,푏=1，即模型方程为푦=푥+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的奇异值。