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

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

import numpy as np
import pandas as pd

一、NumPy 矩阵运算基础

• 在机器学习基础阶段，需要掌握的矩阵及线性代数基本理论包括：
• （1） 矩阵的形变及特殊矩阵的构造方法：包括矩阵的转置、对角矩阵的创建、单位矩阵的创建、上/下三角矩阵的创建等；
• （2） 矩阵的基本运算：包括矩阵乘法、向量内积、矩阵和向量的乘法等；
• （3） 矩阵的线性代数运算：包括矩阵的迹、矩阵的秩、逆矩阵的求解、伴随矩阵和广义逆矩阵等；
• （4） 矩阵分解运算：特征分解、奇异值分解和 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

• 关于两种对象类型的选取，此处进行简单说明：
• （1） NumPy 中的 matrix 类型对象和 MATLAB 中的 matrix 类型等价，和 NumPy 中数组类型对象底层基本结构不同；
• （2） 在 NumPy 中，针对大规模数据，数组类型对象的计算速度要快于矩阵类型对象；
• （3） 矩阵类型对象可以通过运算符直接进行矩阵乘法，而二维数组要进行矩阵乘法（及其他矩阵运算），则必须要使用包括 linalg（线性代数运算）模块在内的相关函数。

AM * AM
#matrix([[3, 4],
#        [2, 3]])

A.dot(A)
#array([[3, 4],
#       [2, 3]])

# 新版NumPy也支持使用符号进行矩阵乘法
A @ A
#array([[3, 4],
#       [2, 3]])

• 在上述代码当中，A.dot(A) 表示矩阵 A 与矩阵 A 乘积。
• 为了执行更高效的计算、以及确保代码整体基本对象类型统一，如无说明，将统一使用二维数组表示矩阵。

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.arange 函数共有三个参数，第一个参数是 start 表示起点，第二个参数是 stop 表示终点，第三个参数是 step 表示步长。
• （1） 一个参数时，参数值为终点，起点取默认值 0，步长取默认值 1。
• （2） 两个参数时，第一个参数为起点，第二个参数为终点，步长取默认值 1。
• （3） 当有三个参数时，第一个参数为起点，第二个参数为终点，第三个参数为步长，步长可以是小数。
• np.reshape 函数用于矩阵规格变换，将矩阵转换为特定的行和列的矩阵。
• （1） 若是 np.reshape(x, -1) 则是将矩阵变成行数为 x，列数不规定的矩阵，具体列数按照总元素个数除行数，均分得到。
• （2） 若是 np.reshape(-1, x) 则是将矩阵变成列数为 x，行数不规定的矩阵，具体行数按照总元素个数除列数，均分得到。

# 创建单位矩阵
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、@

• * ：逐元素相乘

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 并不是内积。
• 矩阵乘法，在 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]])
       
np.trace(B)
#6

• 矩阵的秩（rank），是指矩阵中行或列的极大线性无关数，且矩阵中行、列极大无关数总是相同的，任何矩阵的秩都是唯一值，满秩指的是方阵（行数和列数相同的矩阵）中行数、列数和秩相同，满秩矩阵有线性唯一解等重要特性，而其他矩阵也能通过求解秩来降维，同时，秩也是奇异值分解等运算中涉及到的重要概念。
• 线性相关，其实也就是线性表示，如果 y=wx+b，我们则称 y 可以由 x 线性表示，二者线性相关，反之则线性无关。类似，如果 y=w₁x₁w₂x₂+b，则我们称 y 可以由 x₁、x₂ 线性表示，y与 x₁、x₂ 线性相关。
• 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₋₁。
• 逆矩阵的性质是相互的，我们也可称 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]])

• 当然，对于逆矩阵，还有很多其他理解角度。
• 例如，从方程组求解角度来看，逆矩阵的存在就代表着方程组存在唯一解，并且逆矩阵本身也是方程组求解的关键；从矩阵分解角度来看，逆矩阵是一种最为基础的矩阵分解的形式。

二、矩阵方程与向量求导方法

• 在铺垫了基本矩阵和线性代数相关知识后，接下来，我们尝试将 Lesson 1 中的方程组表示形式转化为矩阵表示形式，并借助矩阵方法进行相关方程的求解。
• 在 Lesson 1 中，我们已经简单讨论了最小二乘法这一优化算法的基本思路，最小二乘法一个最基础的优化算法。
• 因此，本节开始，我们先从矩阵方程入手，先进行矩阵运算的相关方法的回顾、以及进行矩阵求导方法的讲解，讨论最小二乘法的基本原理。

1. 方程组求解与矩阵方程求解

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

• 然后在矩阵方程左右两端同时左乘其逆矩阵，即可解出 X 的取值。
  $$A^{-1}AX=A^{-1}B$$$$X=A^{-1}B$$

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

# 也可以使用dot方法，对于二维数组，dot就是执行矩阵乘法
np.linalg.inv(A).dot(B)
#array([[1.],
#       [1.]]

• 即$$X=\left[\begin{array}{c}w\\ b\end{array}\right]=\left[\begin{array}{c}1\\ 1\end{array}\right]$$
• 当然，此外，NumPy 中还提供了一种求解矩阵方程的函数，类似于上述 A∗XT=B 的矩阵方程，我们还可以通过下式进行求解。

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

2. 向量求导运算

2.1 向量求导基本方法

• 首先，我们先来看相对简单的向量求导方法。假设现有一个二元函数如下。
  $$f(x_1,x_2)=2x_1+x_2$$
• 并且，我们可以分别对该函数中的两个变量 x₁、x₂ 依次求偏导，可得。
  $$\frac{\partial f}{\partial x_1}=2$$ $$\frac{\partial f}{\partial x_2}=1$$
• 现在我们考虑将上述求偏导的函数组改写为矩阵形式。则根据前述内容介绍，我们可以将函数中的两个变量依次排列，组成一个向量变元，即一个由多个变量所组成的向量。
  $$x=[x_1,x_2{]}^T$$
• 此时，如果我们按照向量变元内部的变量排列顺序，依次在每个变量位置填上该变量对应的偏导函数，则就构成了对于函数 f 进行向量变元 x 的向量求导的结果。
  $$\frac{\partial f(x)}{\partial x}=\left[\begin{array}{c}2\\ 1\end{array}\right]$$
• 其中， x 为向量变元。
• 至此，我们就完成了向量求导的基本过程。核心点在于，我们是依据向量变元中的变量排列顺序，依次填写了对应变量的偏导函数计算结果。
• 不过，更进一步的来看，既然方程组需要改写成向量/矩阵形式，那么原始函数方程其实也同样需要改写成向量/矩阵形式。因此，原方程我们可以改写成。
  $$f(x)=A^T\cdot x$$
• 其中。
  $$A=[2,1{]}^T$$
  $$x=[x_1,x_2{]}^T$$
• 原方程为。
  $$y=2x_1+x_2$$
• 结合函数求导结果，我们不难发现，$\frac{\partial f(x)}{\partial x}$最终计算结果就是$A$。
  $$\frac{\partial f(x)}{\partial x}=\frac{\partial (A^T\cdot x)}{\partial x}=A$$
• 其中 x 是向量变元， A 是列向量。
• 很多时候我们并不区分所谓向量方程和矩阵方程，一般所有自变量为向量或矩阵的方程，我们会统一称其为矩阵方程。包含向量或者矩阵的表达式，我们也会统一称其为矩阵表达式。
• 向量求导的定义法。
• 设 $f(x)$ 是一个关于 $x$ 的函数，其中 $x$ 是向量变元，并且 $$x=[x_1,x_2,...,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$$
• 而该表达式也被称为向量求导的梯度向量形式。
  $${\nabla }_{x}f(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$$
• 通过求得函数的梯度向量求解向量导数的方法，也被称为定义法求解。
• 值得注意的是，多元函数是一定能够求得梯度向量的，但梯度向量或者说向量求导结果，能否由一些已经定义的向量解决表示，如 $A$ 就是 $f(x)$ 的向量求导结果，则不一定。

2.2 常见向量求导公式

• 在前期学习中，数学理论推导涉及到的求导以向量变元求导居多，因此，除了掌握基本的向量求导方法以外，我们还需要推导几个常用的向量求导公式，在这些公式中，向量求导结果都能通过一些已经定义的向量简洁表示。
• 同样，此处我们假设 x 为包含 n 各变量的列向量，$x=[x_1,x_2,...,x_n{]}^T$。
• （1）$$\frac{\partial a}{\partial x}=0$$证明：
  $$\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$$
• （2）
  $$\frac{\partial (x^T\cdot A)}{\partial x}=\frac{\partial (A^T\cdot x)}{\partial x}=A$$
• 证明，此时 A 为拥有 n 个分量的常数向量，设 $A=[a_1,a_2,...,a_n{]}^T$，则有
  $$\begin{array}{rl}\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}{c}\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}{c}{a}_1\\ a_2\\ .\\ .\\ .\\ a_n\end{array}\right]=A\end{array}$$
• （3）
  $$\frac{\partial (x^T\cdot x)}{\partial x}=2x$$
• 证明：
  $$\begin{array}{rl}\frac{\partial (x^T\cdot x)}{\partial x}& =\frac{\partial (x_1^2+x_2^2+...+x_n^2)}{\partial x}\\ & =\left[\begin{array}{c}\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}{c}2x_1\\ 2x_2\\ .\\ .\\ .\\ 2x_n\end{array}\right]=2x\end{array}$$
• 此处 $x^{T}x$ 也被称为向量的交叉乘积（crossprod）。
• （4）
  $$\frac{\partial (x^{T}Ax)}{x}=Ax+A^{T}x$$
• 其中 A 是一个（n*n）的矩阵，$A_{n\ast n}=(a_{ij}{)}_{a=1,j=1}^{n,n}$
• 证明：
• 首先，
  $$\begin{array}{rl}{X}^{T}AX& =[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_1{a}_{11}+x_2{a}_{21}+...+x_n{a}_{n1},x_1{a}_{12}+x_2{a}_{22}+...+x_n{a}_{n2},...,x_1{a}_{1n}+x_2{a}_{2n}+...+x_n{a}_{nn}]\cdot \left[\begin{array}{c}{x}_1\\ x_2\\ .\\ .\\ .\\ x_n\end{array}\right]\\ & =x_1(x_1{a}_{11}+x_2{a}_{21}+...+x_n{a}_{n1})+x_2(x_1{a}_{12}+x_2{a}_{22}+...+x_n{a}_{n2})+...+x_n(x_1{a}_{1n}+x_2{a}_{2n}+...+x_n{a}_{nn})\end{array}$$
• 令$$k(x)=x_1(x_1{a}_{11}+x_2{a}_{21}+...+x_n{a}_{n1})+x_2(x_1{a}_{12}+x_2{a}_{22}+...+x_n{a}_{n2})+...+x_n(x_1{a}_{1n}+x_2{a}_{2n}+...+x_n{a}_{nn})$$
• 则有
  $$\frac{\partial k(x)}{\partial x_1}=(x_1{a}_{11}+x_2{a}_{21}+...+x_n{a}_{n1})+(x_1{a}_{11}+x_2{a}_{12}+...+x_n{a}_{1n})$$
• 类似可得：
  $$\begin{array}{rl}\frac{\partial k(x)}{\partial x}& =\left[\begin{array}{c}\frac{\partial k(x)}{\partial x_1}\\ \frac{\partial k(x)}{\partial x_2}\\ .\\ .\\ .\\ \frac{\partial k(x)}{\partial x_n}\end{array}\right]\\ & =\left[\begin{array}{c}(x_1{a}_{11}+x_2{a}_{21}+...+x_n{a}_{n1})+(x_1{a}_{11}+x_2{a}_{12}+...+x_n{a}_{1n})\\ (x_1{a}_{12}+x_2{a}_{22}+...+x_n{a}_{n2})+(x_1{a}_{21}+x_2{a}_{22}+...+x_n{a}_{2n})\\ .\\ .\\ .\\ (x_1{a}_{1n}+x_2{a}_{2n}+...+x_n{a}_{nn})+(x_1{a}_{n1}+x_2{a}_{n2}+...+x_n{a}_{nn})\end{array}\right]\\ & =\left[\begin{array}{c}(x_1{a}_{11}+x_2{a}_{21}+...+x_n{a}_{n1})\\ (x_1{a}_{12}+x_2{a}_{22}+...+x_n{a}_{n2})\\ .\\ .\\ .\\ (x_1{a}_{1n}+x_2{a}_{2n}+...+x_n{a}_{nn})\end{array}\right]+\left[\begin{array}{c}(x_1{a}_{11}+x_2{a}_{12}+...+x_n{a}_{1n})\\ (x_1{a}_{21}+x_2{a}_{22}+...+x_n{a}_{2n})\\ .\\ .\\ .\\ (x_1{a}_{n1}+x_2{a}_{n2}+...+x_n{a}_{nn})\end{array}\right]\\ & =\left[\begin{array}{cccc}{a}_{11}& a_{21}& ...& a_{n1}\\ a_{12}& a_{22}& ...& a_{n2}\\ ...& ...& ...& ...\\ a_{1n}& a_{2n}& ...& a_{nn}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ .\\ .\\ .\\ x_n\end{array}\right]+\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]\left[\begin{array}{c}{x}_1\\ x_2\\ .\\ .\\ .\\ x_n\end{array}\right]\\ & =A^{T}x+Ax\end{array}$$
• 至此，我们完成相关向量求导常用公式的证明。不过从上面的证明不难看出，使用定义法进行公式证明往往会非常繁琐（尽管流程相对清晰）。
• 此外，对于矩阵来说，也有类似的求导方法。即如果变量以矩阵形式出现，则针对该矩阵的方程求导之际上也就是依照矩阵基本结构，在每个位置上对对应的变量分量求偏导函数。只不过相比向量，矩阵多了一个维度、结构更加复杂，因此求解过程也更加复杂。
• 最后，还需要简单进行一个概念辨析，那就是关于矩阵函数和矩阵方程二者概念的区别。
• 矩阵方程：指变量为矩阵的方程；
• 矩阵函数：同函数矩阵，指自变量和因变量都是n阶矩阵的函数，也可以简单理解成由函数构成的矩阵，并且每个函数的变量都是矩阵。

三、最小二乘法的推导过程及使用方法

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

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

• 首先，我们尝试对模型进行矩阵形式改写。
• 模型改写成矩阵表达式。
• 假设多元线性方程有如下形式 $$f(x)=w_1{x}_1+w_2{x}_2+...+w_d{x}_d+b$$
• 令 $w=[w_1,w_2,...w_d{]}^T$，$x=[x_1,x_2,...x_d{]}^T$，则上式可写为 $$f(x)=w^{T}x+b$$
• 在机器学习领域，我们将线性回归自变量系数命名为 w，其实是 weight 的简写，意为自变量的权重。
• 将带入数据后的方程组改写为矩阵方程。
• 并且，假设现在总共有m条观测值，$x^{(i)}=[x_1^{(i)},x_2^{(i)},...,x_d^{(i)}]$，则带入模型可构成 m 个方程：
  $$\left[\begin{array}{c}{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\end{array}\right]=\left[\begin{array}{c}{\hat{y}}_1\\ {\hat{y}}_2\\ .\\ .\\ .\\ {\hat{y}}_m\end{array}\right]$$
• 然后考虑如何将上述方程组进行改写，我们可令
  $$\hat{w}=[w_1,w_2,...,w_d,b{]}^T$$$$\hat{x}=[x_1,x_2,...,x_d,1{]}^T$$$$\hat{X}=\left[\begin{array}{ccccc}{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]$$$$y=\left[\begin{array}{c}{y}_1\\ y_2\\ .\\ .\\ .\\ y_m\end{array}\right]$$$$\hat{y}=\left[\begin{array}{c}{\hat{y}}_1\\ {\hat{y}}_2\\ .\\ .\\ .\\ {\hat{y}}_m\end{array}\right]$$
• 其中
• $\hat{w}$：方程系数所组成的向量，并且我们将自变量系数和截距放到了一个向量；
• $\hat{x}$：方程自变量和1共同组成的向量；
• $\hat{X}$：样本数据特征构成的矩阵，并在最后一列添加一个全为1的列；
• $y$：样本数据标签所构成的列向量;
• $\hat{y}$：预测值的列向量。
• 因此，上述方程组可表示为
  $$\hat{X}\cdot \hat{w}=\hat{y}$$
• 在改写了 $\hat{x}$ 和 $\hat{w}$ 之后，线性模型也可以按照如下形式进行改写：
  $$f(\hat{x})={\hat{w}}^T\cdot \hat{x}$$

2. 构造损失函数

• 在方程组的矩阵表示基础上，我们可以以 SSE 作为损失函数基本计算流程构建关于 $\hat{w}$ 的损失函数：$$SSELoss(\hat{w})=||y-X\hat{w}|{|}_2^2=(y-X\hat{w}{)}^T(y-X\hat{w})$$
• 需要补充两点基础知识：
• 向量的 2- 范数计算公式
  上式中，$||y-X{\hat{w}}^T|{|}_2$为向量的 2- 范数的计算表达式。向量的 2- 范数计算过程为各分量求平方和再进行开平方。例如$a=[1,-1,]$，则$||a|{|}_2=\sqrt{1^2+(-1{)}^2}=\sqrt{2}$。
• 向量的 1- 范数为各分量绝对值之和。值得注意的是，矩阵也有范数的概念，不过矩阵的范数计算要比向量复杂得多。
• 2- 范数计算转化为内积运算
  向量的 2- 范数计算结果其实就是向量（将其是做矩阵）的交叉乘积计算结果后开平方。例如，$a=[1,-1]$，则$a$的交叉乘积为$a\cdot a^T=[1,-1]\cdot \left[\begin{array}{c}1\\ -1\end{array}\right]=2$，开平方后等于其 2- 范数计算结果。

3. 最小二乘法求解损失函数的一般过程

• 在确定损失函数的矩阵表示形式之后，接下来即可利用最小二乘法进行求解。其基本求解思路，先求导函数、再令导函数取值为零，进而解出参数取值。只不过此时求解的是矩阵方程。
• 在此之前，需要补充两点矩阵转置的运算规则：
  $$(A-B{)}^T=A^T-B^T$$$$(AB{)}^T=B^T{A}^T$$
• 接下来，对 $SSELoss(w)$ 求导并令其等于 0。
  $$\begin{array}{rl}\frac{SSELoss(\hat{w})}{\partial \hat{\mathbf{w}}}& =\frac{\partial {||\mathbf{y}-\mathbf{X}\hat{\mathbf{w}}|{|}_2}^2}{\partial \hat{\mathbf{w}}}\\ & =\frac{\partial (\mathbf{y}-\mathbf{X}\hat{\mathbf{w}}{)}^T(\mathbf{y}-\mathbf{X}\hat{\mathbf{w}})}{\partial \hat{\mathbf{w}}}\\ & =\frac{\partial ({\mathbf{y}}^T-{\hat{\mathbf{w}}}^{\mathbf{T}}{\mathbf{X}}^{\mathbf{T}})(\mathbf{y}-\mathbf{X}\hat{\mathbf{w}})}{\partial \hat{\mathbf{w}}}\\ & =\frac{\partial ({\mathbf{y}}^T\mathbf{y}-{\hat{\mathbf{w}}}^{\mathbf{T}}{\mathbf{X}}^{\mathbf{T}}\mathbf{y}-{\mathbf{y}}^T\mathbf{X}\hat{\mathbf{w}}+{\hat{\mathbf{w}}}^{\mathbf{T}}{\mathbf{X}}^{\mathbf{T}}\mathbf{X}\hat{\mathbf{w}})}{\partial \hat{\mathbf{w}}}\\ & =0-{\mathbf{X}}^{\mathbf{T}}\mathbf{y}-{\mathbf{X}}^{\mathbf{T}}\mathbf{y}+X^{T}X\hat{w}+(X^{T}X{)}^T\hat{w}\\ & =0-{\mathbf{X}}^{\mathbf{T}}\mathbf{y}-{\mathbf{X}}^{\mathbf{T}}\mathbf{y}+2{\mathbf{X}}^{\mathbf{T}}\mathbf{X}\hat{\mathbf{w}}\\ & =2({\mathbf{X}}^{\mathbf{T}}\mathbf{X}\hat{\mathbf{w}}-{\mathbf{X}}^{\mathbf{T}}\mathbf{y})=0\end{array}$$
• 即$$X^{T}X\hat{w}=X^{T}y$$
• 要使得此式有解，等价于 $X^{T}X$（也被称为矩阵的交叉乘积 crossprod 存在逆矩阵，若存在，则可解出$$\hat{w}=(X^{T}X{)}^{-1}{X}^{T}y$$

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

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

Whole weight	Rings
1	2
3	4

• 因此利用矩阵表达式，可进行如下形式的改写。
• 特征矩阵：$$\hat{X}=\left[\begin{array}{cc}1& 1\\ 3& 1\end{array}\right]$$
• 标签数组：$$y=\left[\begin{array}{c}2\\ 4\end{array}\right]$$
• 参数向量：$$\hat{w}=\left[\begin{array}{c}w\\ b\end{array}\right]$$
• 求解公式为：
  $$\begin{array}{rl}\hat{w}& =(X^{T}X{)}^{-1}{X}^{T}y\\ & =({\left[\begin{array}{cc}1& 1\\ 3& 1\end{array}\right]}^T\left[\begin{array}{cc}1& 1\\ 3& 1\end{array}\right]{)}^{-1}{\left[\begin{array}{cc}1& 1\\ 3& 1\end{array}\right]}^T\left[\begin{array}{c}2\\ 4\end{array}\right]\end{array}$$
• 代码实现过程：

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

• 即解得 $w=1，b=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 的奇异值。