线性代数的代数应用

线性代数推导斐波那契递推公式

我们知道，斐波那契数列是一个经典的用递归定义的函数。

$$F(n) = \begin{cases} 0 & \text{if } n = 0 \\ 1 & \text{if } n = 1 \\ F(n-1) + F(n-2) & \text{if } n \geq 2 \end{cases}$$

不难发现这个公式的递推关系是线性的。

将初始状态用一个向量 $v_0$ 表示：

$$v_0 = \begin{bmatrix} F(1) \\ F(0) \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \end{bmatrix}$$

在向量的基础上进行递推，发现下一个向量 $v_1$ 是：

$$v_1 = \begin{bmatrix} F(2) \\ F(1) \end{bmatrix} = \begin{bmatrix} F(1) + F(0) \\ F(1) \end{bmatrix} = \begin{bmatrix} 1 + 0 \\ 1 \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \end{bmatrix}$$

继续递推，$v_2$ 是：

$$v_2 = \begin{bmatrix} F(3) \\ F(2) \end{bmatrix} = \begin{bmatrix} F(2) + F(1) \\ F(2) \end{bmatrix} = \begin{bmatrix} 1 + 1 \\ 1 \end{bmatrix} = \begin{bmatrix} 2 \\ 1 \end{bmatrix}$$

观察这个递推过程，我们可以发现一个线性变换矩阵 $A$，使得：

$$v_{n+1} = A \cdot v_n$$

其中矩阵 $A$ 满足：

$$A \cdot \begin{bmatrix} F(n) \\ F(n-1) \end{bmatrix} = \begin{bmatrix} F(n+1) \\ F(n) \end{bmatrix} = \begin{bmatrix} F(n) + F(n-1) \\ F(n) \end{bmatrix}$$

通过矩阵乘法，我们可以得到：

$$A = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}$$

因此，斐波那契数列的递推关系可以用矩阵形式表示为：

$$v_n = A^n \cdot v_0$$

其中 $v_0 = \begin{bmatrix} 1 \\ 0 \end{bmatrix}$ 是初始状态向量。

此时我们成功得到了斐波那契数列的通项公式，不过我们还不能直接知道 $A^n$ 在 n 取某个值的时候是多少，必须要算出来。

那么有没有什么方法可以让我们直接知道某个矩阵的 n 次方是多少呢？

答案就是特征值矩阵分解。

我们知道一个对称矩阵可以进行特征分解（论证这里不赘述）：

$$A = P \cdot D \cdot P^{-1}$$

其中 $P$ 是特征向量矩阵，$D$ 是对角矩阵，其对角线元素是特征值。

那么 $A^n$ 可以表示为：

$$A^n = (P \cdot D \cdot P^{-1})^n = P \cdot D^n \cdot P^{-1}$$

由于 $D$ 是对角矩阵，$D^n$ 的计算非常简单，只需要将每个对角线元素取 n 次方即可。

对于斐波那契数列的矩阵 $A = \begin{bmatrix} 1 & 1 \\ 1 & 0 \end{bmatrix}$，我们可以计算其特征值：

$$\det(A - \lambda I) = \det\begin{bmatrix} 1-\lambda & 1 \\ 1 & -\lambda \end{bmatrix} = (1-\lambda)(-\lambda) - 1 = \lambda^2 - \lambda - 1 = 0$$

解这个特征方程，得到特征值：

$$\lambda_1 = \frac{1 + \sqrt{5}}{2} = \phi \approx 1.618$$

$$\lambda_2 = \frac{1 - \sqrt{5}}{2} = 1 - \phi \approx -0.618$$

其中 $\phi$ 就是著名的黄金比例。

接下来计算对应的特征向量。对于特征值 $\lambda_1 = \phi$，解方程 $(A - \phi I)v = 0$：

$$\begin{bmatrix} 1-\phi & 1 \\ 1 & -\phi \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}$$

得到特征向量 $v_1 = \begin{bmatrix} \phi \\ 1 \end{bmatrix}$。

对于特征值 $\lambda_2 = 1-\phi$，解方程 $(A - (1-\phi) I)v = 0$：

$$\begin{bmatrix} \phi & 1 \\ 1 & \phi-1 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}$$

得到特征向量 $v_2 = \begin{bmatrix} 1-\phi \\ 1 \end{bmatrix}$。

因此特征向量矩阵 $P$ 为：

$$P = \begin{bmatrix} \phi & 1-\phi \\ 1 & 1 \end{bmatrix}$$

其逆矩阵 $P^{-1}$ 为：

$$P^{-1} = \frac{1}{\phi - (1-\phi)} \begin{bmatrix} 1 & -(1-\phi) \\ -1 & \phi \end{bmatrix} = \frac{1}{2\phi - 1} \begin{bmatrix} 1 & -(1-\phi) \\ -1 & \phi \end{bmatrix}$$

由于 $\phi = \frac{1+\sqrt{5}}{2}$，所以 $2\phi - 1 = \sqrt{5}$，因此：

$$P^{-1} = \frac{1}{\sqrt{5}} \begin{bmatrix} 1 & -(1-\phi) \\ -1 & \phi \end{bmatrix}$$

由此，我们最终得到 $A^n$ 为：

$$A^n = \begin{bmatrix} \phi & 1-\phi \\ 1 & 1 \end{bmatrix} \cdot \begin{bmatrix} \phi^n & 0 \\ 0 & (1-\phi)^n \end{bmatrix} \cdot \frac{1}{\sqrt{5}} \begin{bmatrix} 1 & -(1-\phi) \\ -1 & \phi \end{bmatrix}$$

进一步得到斐波那契公式的真正通项公式为：

$$F(n) = \frac{\phi^n - (1-\phi)^n}{\sqrt{5}} = \frac{1}{\sqrt{5}}\left[\left(\frac{1+\sqrt{5}}{2}\right)^n - \left(\frac{1-\sqrt{5}}{2}\right)^n\right]$$

这就是著名的斐波那契数列的闭式解（Binet公式）！

python复制
import numpy as np

# 代码实现斐波那契数列通项推导过程
A = np.array([[1, 1], [1, 0]])
v0 = np.array([1, 0])

# 计算特征值和特征向量
eigenvalues, eigenvectors = np.linalg.eig(A)

# 计算特征值矩阵
D = np.diag(eigenvalues)

# 计算特征向量矩阵
P = eigenvectors

def fibonacci(n):
    return (P @ np.diag(np.power(eigenvalues, n)) @ np.linalg.inv(P) @ v0)[0]

# 验证
for i in range(10):
    print(fibonacci(i))

plaintext复制
    0.9999999999999999
    1.0000000000000002
    1.9999999999999998
    3.0000000000000004
    5.0
    8.000000000000002
    13.000000000000004
    21.000000000000004
    34.000000000000014
    55.00000000000002

线性代数推导求导公式

现在我们来看一下多项式求导前后式子的关系

对于多项式 $f(x) = a_0 + a_1x + a_2x^2 + \cdots + a_{n-1}x^{n-1} + a_nx^n$，
其求导结果为：

$$f'(x) = a_1 + 2 a_2 x + 3 a_3 x^2 + \cdots + n a_n x^{n-1}$$

也就是说，多项式每一项 $a_k x^k$ 求导后变为 $k a_k x^{k-1}$。

假如我们用系数向量 $(a_0, a_1, a_2, \ldots, a_n)$ 来表示一个 $n$ 次多项式 $f(x)$，那么多项式求导的过程实际上就是对系数向量进行一次线性变换，这个变换可以用一个矩阵来表示。

具体来说，常数项 $a_0$ 求导后为 0，所以变换矩阵的第一列全为 0；一次项 $a_1 x$ 求导后变为 $a_1$，所以变换矩阵的第二列第一行为 1，其余为 0；二次项 $a_2 x^2$ 求导后变为 $2 a_2 x$，所以变换矩阵的第三列第二行为 2，其余为 0。依此类推。

因此，对于 $n$ 次多项式，求导的线性变换矩阵 $D$ 形式如下：

$$D = \begin{pmatrix} 0 & 1 & 0 & 0 & \cdots & 0 \\ 0 & 0 & 2 & 0 & \cdots & 0 \\ 0 & 0 & 0 & 3 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & \cdots & n \\ 0 & 0 & 0 & 0 & \cdots & 0 \\ \end{pmatrix}$$

python复制
import numpy as np
import sympy as sp
from IPython.display import Latex, display

# 🧮 多项式类，支持自定义打印，print时自动以LaTeX格式显示
class Poly:
    def __init__(self, coefficients):
        """
        初始化多项式，参数为系数列表
        :param coefficients: 多项式系数列表，低次幂在前
        """
        self.coefficients = coefficients

    def __call__(self, x):
        """
        计算多项式在x处的值
        :param x: 变量值
        :return: 多项式值
        """
        return sum(c * x**i for i, c in enumerate(self.coefficients))

    def __str__(self):
        """
        自定义print打印，自动以LaTeX格式显示多项式
        """
        x = sp.symbols('x')
        expr = sum(c * x**i for i, c in enumerate(self.coefficients))
        latex_expr = sp.latex(expr)
        display(Latex(f"${latex_expr}$"))
        return ''

    def diff(self):
        """
        返回多项式的导数（作为新的Poly对象）
        :return: 导数多项式
        """
        # 求导后最高次项的系数为0，长度与原多项式一致
        # return Poly([i * c for i, c in enumerate(self.coefficients[1:], start=1)] + [0])
        # 使用矩阵表示
        return Poly(np.diag(np.arange(1, len(self.coefficients)), k=1) @ self.coefficients)

# 原多项式
P = Poly([1, 2, 3, 4, 5])
print(P)  # print时自动以LaTeX格式显示
print(P(1))

# 求导后的多项式
D = P.diff()
print(D)  # print时自动以LaTeX格式显示
print(D(1))

$5 x^{4} + 4 x^{3} + 3 x^{2} + 2 x + 1$

15

$20 x^{3} + 12 x^{2} + 6 x + 2$

40

结果是正确的

我们还知道，任何连续函数都可以通过泰勒展开写成多项式的形式。

所以我们就可以通过这种方法求任意连续函数的导函数。（不过好像得到泰勒展开式本来就需要求导吧，所以这玩意运算上没有什么用）

不过他通过用矩阵的方式实现了求导这个算子，我们就可以通过对这个矩阵的研究来研究求导的性质

比如这个矩阵存在逆矩阵吗？如果存在的话是什么呢？🤔

在有限维的情况下，这个矩阵的行列式为0，因此它不可逆。

但是我们仍然可以得到一个接近于逆矩阵的东西——原矩阵转置后对非零元素求倒数。

$$D' = \begin{pmatrix} 0 & 0 & 0 & \cdots & 0 & 0 \\ 1 & 0 & 0 & \cdots & 0 & 0 \\ 0 & \frac{1}{2} & 0 & \cdots & 0 & 0 \\ 0 & 0 & \frac{1}{3} & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & \frac{1}{n} & 0 \\ \end{pmatrix}$$

如果把 D’ 作用在一个多项式系数向量上，会发现就相当于对这个多项式进行不完全积分（缺少常数项）。

把得到的这个矩阵与原矩阵相乘，得到：

$$D'D = \begin{pmatrix} 0 & 0 & 0 & \cdots & 0 \\ 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \\ \end{pmatrix}$$

$$DD' = \begin{pmatrix} 1 & 0 & 0 & \cdots & 0 \\ 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 0 \\ \end{pmatrix}$$

发现，当 n 是一个有限的整数时，两种矩阵相乘的结果都不是单位矩阵，所以这个矩阵在有限维空间中没有严格意义上的逆矩阵。

进一步，可以发现，先求导D再积分D’，相当于去掉原函数的常数项。
先积分D’再求导D，相当于把原函数的最高次项去掉。

如果，我们考虑n趋近于无穷大的情形，DD’舍去的最高次项就会变得越来越不重要，DD’的结果会无限逼近单位矩阵。

就可以从这个角度认识到，先求导再积分不可逆，但是先积分再求导在无限次函数中是近似可逆的。

线性代数求最小二乘法公式

先来看最基本的最小二乘法的问题：

已知有 $n$ 组已知的数据点 $(x_i, y_i)$，我们用两个 $n$ 维的向量 $X$ 和 $Y$ 来表示这些数据点。

假设这组数据满足一个带噪声的线性的关系 $y = kx + b + \epsilon$，其中 $k$ 和 $b$ 是未知的参数, $\epsilon \sim N(0, \sigma)$

我们想要用一条直线拟合这批数据，也就是计算出 $\beta = \begin{bmatrix} k \\ b \end{bmatrix}$，是的预测出的数据点 $\hat Y = X\beta$ 与 真实数据点 $Y$ 之间的 L2 范数最小。其中的 X 是原数据中的 X 向量在右边并上一列全为 1 的向量后得到的$(n, 2)$矩阵

由线性代数的知识，不难发现，$X$ 的列空间就是 n 维空间中的 2 维平面（前提是 n 组数据点的 x 坐标不全部相同）

那么 $\hat Y = X\beta$ 就是上述二维平面中的一个向量

我们的目标就是求出 $\beta$ 使得 $\|Y - \hat Y\|_2^2 = \|Y - X\beta\|_2^2$ 最小

由几何知识，当上述值最小时，$\hat Y$ 就是 $Y$ 在 $\beta$ 平面上的投影。我们知道，此时 $Y - \hat Y$ 必须垂直于 $X$ 的列空间

再由线性代数的知识，有 $X^T(Y-\hat Y) = 0$

即 $X^T(Y-X\beta) = 0$

最后得到 $X^TX\beta = X^TY$

$$\beta = (X^TX)^{-1}X^TY$$

这就是从线性代数的角度推导出的最小二乘法的求解公式

下面我们来实际运用一下

python复制
import numpy as np
import matplotlib.pyplot as plt

k = 3
b = 5

n = 100

sigma = 2

# 生成n组数据点
x = np.linspace(0, 10, n)
y = k * x + b + np.random.normal(0, sigma, n)

# 给 x 添加一列全1，表示截距项
x_b = np.vstack([x, np.ones(n)]).T
y = y.reshape(-1, 1)

# 使用最小二乘法拟合直线
beta = np.linalg.inv(x_b.T @ x_b) @ x_b.T @ y

print(f"Estimated parameters:\n {beta}")

# 画出拟合直线
plt.scatter(x, y, label='Data points')
plt.plot(x, x_b @ beta, color='red', label='Fitted line')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Linear Fit using Least Squares')
plt.legend()
plt.show()

Estimated parameters:
 [[3.01664099]
 [5.02909287]]

效果不错

不过，这个公式更厉害的点还在于他还能够适用于多项式的最小二乘法拟合的情况

python复制
# 多项式拟合

import numpy as np
import matplotlib.pyplot as plt

p = [0.2, 0, -2, 0, 3.6, 0]

# 原多项式
print(f"Original coefficients:\n {p}")

# 生成数据点
n = 200
sigma = 0.7
x = np.linspace(-3, 3, n)
y = np.polyval(p, x) + np.random.normal(0, sigma, n)  # 添加噪声

# 构造Vandermonde矩阵
X = np.vander(x, N=len(p), increasing=True)

# 使用最小二乘法拟合多项式
beta = np.linalg.inv(X.T @ X) @ X.T @ y

print(f"Estimated coefficients:\n {beta[::-1]}")

# 画出拟合曲线
plt.scatter(x, y, label='Data points')
plt.plot(x, X @ beta, color='red', label='Fitted line')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Polynomial Fit using Least Squares')
plt.legend()
plt.show()

Original coefficients:
 [0.2, 0, -2, 0, 3.6, 0]
Estimated coefficients:
 [ 0.20228931 -0.0107228  -2.00856267  0.0600106   3.53675279  0.0274673 ]