Python 现代控制理论 —— 梯度下降法实现的线性回归系统

荷碧TongZJ

已于 2023-06-05 17:44:25 修改

阅读量1k

点赞数 3

分类专栏：课程设计文章标签： python 线性回归

于 2022-12-27 22:21:03 首次发布

本文链接：https://blog.csdn.net/qq_55745968/article/details/128461481

版权

课程设计专栏收录该内容

6 篇文章 5 订阅

订阅专栏

线性回归是有监督学习中的经典问题，其核心在于找到样本的多个特征与标签值之间的线性关系。样本集中的第j个样本可被表示为：

特征向量： $X_j=\begin{bmatrix} x_{j1} &x_{j2} & \hdots & x_{jn} \end{bmatrix}$
标签值： $u_j$

而线性回归系统给出权重向量：

$w=\begin{bmatrix} w_0 & w_1 & w_2 & \hdots & w_n \end{bmatrix}^T$

使得该样本的预测值为：

$y_j=w_0+X_jw=w_0+\sum^n_{i=1}x_{ji}w_i$

当所有样本的预测值与标签值的误差最小时，即代表该线性回归系统找到了最优的拟合曲线

本文采用了梯度下降法以解决线性回归问题。梯度下降法与现代控制系统相似，现代控制系统实现的梯度下降法如下：

该系统将线性回归问题的权重向量 $w$ 作为状态变量，以损失函数 $L(w)$ 反向传播的梯度 $-\frac{dL(w)}{dw}$ 作为 $\frac{dw}{dt}$ ，并通过 $w= w +\frac{dw}{dt} dt$ 对权重向量进行更新 (其中的 $dt$ 可视为学习率)，使得所有样本的误差越来越小

状态空间描述

在前文的讨论中，样本的预测值 $y_j$ 表示为“特征向量×权重向量+偏置”的形式。为了后续计算的整洁，将样本的特征向量表示为：

$X_j=\begin{bmatrix} 1 & x_{j1} &x_{j2} & \hdots & x_{jn} \end{bmatrix}$

则样本的预测值可被重写为：

$y_j=X_jw$

以样本的预测值与标签值的差值作为误差 $e$ ：

$e=y-u=\begin{bmatrix} 1& x_{11} & x_{12} & \hdots & x_{1n}\\ 1& x_{21}&x_{22} & \hdots & x_{2n}\\ \vdots &\vdots &\vdots & \ddots &\vdots \\ 1& x_{m1}& x_{m2}& \hdots & x_{mn}\end{bmatrix} w - \begin{bmatrix} u_1\\ u_2\\ \hdots\\ u_m \end{bmatrix}$

误差 $e$ 的值域为 $(- \infty, \infty)$ ，最优值为 0，显然不可直接作为梯度下降法的损失函数。故以其原函数 (MSE) 作为损失函数：

$L(w)=\frac{\overline{e^2}}{2} \in [0, \infty)$

$\frac{\partial L(w)}{\partial w_i}=\overline{(e\frac{\partial e}{\partial w_i})} = \overline{(x_{.i}X)}w - \overline{(x_{.i}u)},\ x_{.0}=1$

$\frac{dL(w)}{dw}=\begin{bmatrix} 1 & \overline{x_{.1}} & \overline{x_{.2}} & \hdots&\overline{x_{.n}} \\ \overline{x_{.1}} & \overline{x_{.1}^2} & \overline{x_{.1}x_{.2}} & \hdots & \overline{x_{.1}x_{.n}}\\ \vdots &\vdots & \vdots &\ddots &\vdots \\ \overline{x_{.n}} & \overline{x_{.n}x_{.1}} & \overline{x_{.n}x_{.2}} & \hdots & \overline{x_{.n}^2} \end{bmatrix} w - \begin{bmatrix} \overline{u}\\ \overline{(x_{.1}u)}\\ \vdots\\ \overline{(x_{.n}u)} \end{bmatrix}$

由梯度下降法可知，以 $w = w - \frac{dL(w)}{dw} p$ 的方式对权重向量进行更新，将使得损失函数 $L(w)$ 逐渐减小。故令 $\frac{dw}{dt}=-\frac{dL(w)}{dw},\ dt=p$ ，以 $w=w+\frac{dw}{dt}dt$ 对权重向量进行更新。综上所述，该系统的状态空间描述为：

状态方程： $\dot{w}=-\begin{bmatrix} 1 & \overline{x_{.1}} & \overline{x_{.2}} & \hdots&\overline{x_{.n}} \\ \overline{x_{.1}} & \overline{x_{.1}^2} & \overline{x_{.1}x_{.2}} & \hdots & \overline{x_{.1}x_{.n}}\\ \vdots &\vdots & \vdots &\ddots &\vdots \\ \overline{x_{.n}} & \overline{x_{.n}x_{.1}} & \overline{x_{.n}x_{.2}} & \hdots & \overline{x_{.n}^2} \end{bmatrix} w + \begin{bmatrix} \overline{u}\\ \overline{(x_{.1}u)}\\ \vdots\\ \overline{(x_{.n}u)} \end{bmatrix} 1(t)$
输出方程： $e=\begin{bmatrix} 1& x_{11} & x_{12} & \hdots & x_{1n}\\ 1& x_{21}&x_{22} & \hdots & x_{2n}\\ \vdots &\vdots &\vdots & \ddots &\vdots \\ 1& x_{m1}& x_{m2}& \hdots & x_{mn}\end{bmatrix} w - \begin{bmatrix} u_1\\ u_2\\ \hdots\\ u_m \end{bmatrix} 1(t)$

其中，状态方程描述了权重向量 $w$ 随时间变化的速度，输出方程描述了该线性回归系统对每一个样本的误差

实验案例

给定定义域 $x \in [-3, 3]$ ，在下式所示的曲线上均匀地选取样本点 (上图中蓝色散点)，并在纵坐标 (标签值) 上添加噪声：

$y=2+x+0.3x^2-0.5x^3+4\sin{x}$

横坐标为 $x_j$ 的样本的特征向量被定义为：

$X_j = \begin{bmatrix} 1 &x_j & e^{-x_j} & e^{x_j} \end{bmatrix}$

将样本的特征向量与标签值输入线性回归系统，该系统将用曲线 $y=w_0 + w_1x +w_2 e^{-x} + w_3 e^x$ 对样本进行拟合。运行 Python 程序（见附录)后，还得到状态空间表达式中的四个矩阵

最终拟合曲线为上图中的橙色曲线，该拟合曲线的权重向量为：

$w=\begin{bmatrix} 2.02 & 6.86 & 1.73 & -1.46 \end{bmatrix}^T$

该拟合曲线的方程为：

$y=2.02 + 6.86 x +1.73 e^{-x} - 1.46 e^x$

代码实现

上述系统主要通过 Auto_Ctrl_System 类实现，该类的实例属性有：

A, B, C, D：系统矩阵 $A$ ，输入矩阵 $B$ ，输出矩阵 $C$ ，直接传递矩阵 $D$
w：权重向量
eig_coef：矩阵 $A$ 的特征多项式的系数
dstate, output：状态方程、输出方程的结果

实例方法有：

adjust：根据状态方程，调整权重向量；同时根据 fine_grit 参数对误差值之和进行采样
predict：利用当前的权重向量，对样本 (可另给定) 给出预测值
ctrl：转换系统为能控标准型
jordan：转换系统为 Jordan 标准型 (该系统中的矩阵 A 为对称矩阵，可对角化)， $T$ 由矩阵 $A$ 的特征向量构成
state_tran：利用矩阵 $A$ 可被对角化这一点，对该系统的状态转移矩阵进行求解

import control as ct
import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
from tqdm import trange

np.set_printoptions(precision=3, suppress=True)
INFO = lambda *args: print(*args, sep='\n')


def simplify(x, eps=1e-5):
    ''' 针对 exp 多项式的简化函数, 不一定普适'''
    if isinstance(x, sp.Add):
        x = list(x.args)
        for i in range(len(x)):
            coef, exp = x[i].args
            coef, exp = map(float, [coef, exp.args[0].args[0]])
            # exp 系数小于阈值, exp 指数项为负数
            if abs(coef) < eps and exp < 0: x[i] = 0
        x = sum(x)
    return x


class Auto_Ctrl_System:
    n = property(fget=lambda self: self.w.size)
    # 状态方程, 输出方程: 输入为阶跃函数
    dstate = property(fget=lambda self: self.A @ self.w + self.B)
    output = property(fget=lambda self: self.C @ self.w + self.D)
    # 对输入的矩阵预处理 (添加偏置)
    _process = lambda self, x: np.concatenate([np.ones([x.shape[0], 1]), x], axis=1) if self.bias else x

    def __init__(self, x, y, bias=True):
        x, y = map(lambda i: np.float16(i), [x, y])
        n = x.shape[1] + bias
        self.bias = bias
        # 输出方程的矩阵 C, D
        self.C = self._process(x)
        self.D = - y[:, None]
        # 状态方程的矩阵 A, B
        self.A = - (self.C[..., None] * self.C[:, None]).mean(axis=0)
        self.B = - (self.C * self.D).mean(axis=0, keepdims=True).T
        # 初始化权重向量
        self.w = np.random.random([n, 1])
        # Lyapunov 矩阵 P, Q = I
        p = ct.lyap(self.A, np.eye(self.n))
        stab = ' (正定)' * (np.linalg.eig(p)[0].min() > 0).item()
        INFO(f'Lyapunov 矩阵 P{stab}:', p)
        # 特征多项式的系数: a0, a1 ...
        tran = ct.ctrb(self.A, self.B)
        tran_inv = self._inv(tran)
        A = tran_inv @ self.A @ tran
        self.eig_coef = - A[:, -1]
        # 不能控的状态数
        print(f'不能控状态数: {self.n - sp.Matrix(tran).T.rank()}')

    def adjust(self, epochs=2e4, dt=1e-3, fine_grit=100):
        qbar = trange(round(epochs))
        # 采样间隔, 采样结果
        sample_t = int(epochs / fine_grit)
        sum_error = []
        for i in qbar:
            self.w += self.dstate * dt
            output = self.output
            loss = np.square(output).mean()
            qbar.set_description(f'MSE {loss:.2f}')
            # 对输出采样
            if not i % sample_t: sum_error.append(output.sum())
        t = np.linspace(0, epochs * dt, fine_grit)
        return t, np.array(sum_error)

    def predict(self, x=None):
        x = self._process(x) if isinstance(x, np.ndarray) else self.C
        return (x @ self.w).flatten()

    def ctrl(self, i=1):
        string = 'ⅠⅡ'
        assert i in (1, 2), f'只支持能控标准{"/".join(string)}型'
        tran = ct.ctrb(self.A, self.B)
        # 能控标准Ⅰ型
        if i == 1:
            coef = self.eig_coef * (np.arange(self.n) > 0)
            tran = tran[:, ::-1]
            tran_a = np.stack(np.meshgrid(*[np.arange(self.n) for _ in range(2)]),
                              -1).sum(axis=-1)[::-1] + 1
            tran_a = coef[tran_a * (tran_a < self.n)] + np.eye(self.n)
            tran = tran @ tran_a
        self._ns_tran(tran)
        INFO(f'能控标准{string[i - 1]}型:', self)
        return tran

    def jordan(self):
        tran = np.linalg.eig(self.A)[1]
        self._ns_tran(tran)
        INFO('Jordan 标准型:', self)
        return tran

    def state_tran(self, response=False):
        t = sp.symbols('t')
        # 对角化的变换矩阵
        tran = np.linalg.eig(self.A)[1]
        tran_inv = self._inv(tran)
        # 状态转移矩阵
        diag = sp.Matrix(np.eye(self.n) * (tran_inv @ self.A @ tran) * t).exp()
        state_tran = tran @ diag @ tran_inv
        # 对 exp 多项式进行简化
        for i in range(self.n):
            for j in range(self.n):
                state_tran[i, j] = simplify(state_tran[i, j])
        INFO('状态转移矩阵:', np.array(state_tran))
        # fixme: 阶跃响应 (针对性求解, 勿用, 需根据状态转移矩阵设计)
        if response:
            fc, f, gc, g, a, b = sp.symbols('f_c, f, g_c, g, a, b')
            # 利用新变量替代状态转移矩阵中的 exp 多项式
            state_tran = sp.Matrix([[fc, 0, f, f],
                                    [0, gc, -g, g],
                                    [f, -g, a + b, a - b],
                                    [f, g, a - b, a + b]])
            # 零输入响应, 零状态响应
            zero_input = state_tran @ self.w
            zero_state = self._inv(self.A) @ (state_tran - sp.eye(self.n)) @ self.B
            return zero_input, zero_state
        return state_tran

    def _ns_tran(self, tran):
        tran_inv = self._inv(tran)
        INFO('非奇异变换矩阵:', tran)
        self.A = tran_inv @ self.A @ tran
        self.B = tran_inv @ self.B
        self.C = self.C @ tran
        self.w = tran_inv @ self.w

    def _inv(self, arr):
        inv = np.linalg.inv(arr)
        assert np.all(np.isfinite(inv))
        return inv

    def __repr__(self):
        return str(np.concatenate([self.A, self.B], axis=1))
    
    
def standard_coord(*args):
    ''' args: subplot 参数'''
    fig = plt.subplot(*args)
    for key in 'right', 'top':
        fig.spines[key].set_color('None')
    return fig


if __name__ == '__main__':
    x = np.linspace(-3, 3, 50)
    z = np.stack([x, np.exp(-x), np.exp(x)], axis=1)

    np.random.seed(10)
    # 原函数: x + 0.3 x^2 - 0.5 x^3 + 4 sin(x) + 噪声
    y = 2 + x + 0.3 * np.power(x, 2) - 0.5 * np.power(x, 3) + 4 * np.sin(x) + 5 * (np.random.random(len(x)) - 0.5)
    plt.scatter(x, y, c='deepskyblue', label='true')

    acs = Auto_Ctrl_System(z, y, bias=True)
    # 转化为 Jordan 标准型
    acs.jordan()
    t, y_ = acs.adjust(2e4)
    # 可视化拟合结果
    plt.plot(x, acs.predict(), c='orange', label='pred')
    plt.legend(frameon=False)
    plt.show()

    plt.rcParams['figure.figsize'] = [3.0, 2.0]
    # 绘制阶跃响应: 误差值之和
    standard_coord()
    plt.plot(t, y_, c='orange')
    plt.plot([0, t[-1]], [0, 0], c='gray', linestyle='--')
    plt.tight_layout()
    plt.show()

荷碧TongZJ

关注

3
点赞
踩
7

收藏

觉得还不错? 一键收藏
打赏
0
评论
Python 现代控制理论 —— 梯度下降法实现的线性回归系统

线性回归是有监督学习中的经典问题，其核心在于找到样本的多个特征与标签值之间的线性关系。而线性回归系统给出权重向量：使得该样本的预测值为：当所有样本的预测值与标签值的误差最小时，即代表该线性回归系统找到了最优的拟合曲线本文采用了梯度下降法以解决线性回归问题。梯度下降法与现代控制系统相似，现代控制系统实现的梯度下降法如下：该系统将线性回归问题的权重向量作为状态变量，以损失函数反向传播的梯度作为，并通过对权重向量进行更新(其中的可视为学习率)，使得所有样本的误差越来越小。
复制链接

扫一扫