线性回归是有监督学习中的经典问题,其核心在于找到样本的多个特征与标签值之间的线性关系。样本集中的第j个样本可被表示为:
- 特征向量:
- 标签值:
而线性回归系统给出权重向量:
使得该样本的预测值为:
当所有样本的预测值与标签值的误差最小时,即代表该线性回归系统找到了最优的拟合曲线
本文采用了梯度下降法以解决线性回归问题。梯度下降法与现代控制系统相似,现代控制系统实现的梯度下降法如下:
该系统将线性回归问题的权重向量
作为状态变量,以损失函数
反向传播的梯度
作为
,并通过
对权重向量进行更新 (其中的
可视为学习率),使得所有样本的误差越来越小
状态空间描述
在前文的讨论中,样本的预测值
表示为“特征向量×权重向量+偏置”的形式。为了后续计算的整洁,将样本的特征向量表示为:
则样本的预测值可被重写为:
以样本的预测值与标签值的差值作为误差
:
误差
的值域为
,最优值为 0,显然不可直接作为梯度下降法的损失函数。故以其原函数 (MSE) 作为损失函数:
由梯度下降法可知,以
的方式对权重向量进行更新,将使得损失函数
逐渐减小。故令
,以
对权重向量进行更新。综上所述,该系统的状态空间描述为:
- 状态方程:
- 输出方程:
其中,状态方程描述了权重向量
随时间变化的速度,输出方程描述了该线性回归系统对每一个样本的误差
实验案例
给定定义域
,在下式所示的曲线上均匀地选取样本点 (上图中蓝色散点),并在纵坐标 (标签值) 上添加噪声:
横坐标为
的样本的特征向量被定义为:
将样本的特征向量与标签值输入线性回归系统,该系统将用曲线
对样本进行拟合。运行 Python 程序(见附录)后,还得到状态空间表达式中的四个矩阵
最终拟合曲线为上图中的橙色曲线,该拟合曲线的权重向量为:
该拟合曲线的方程为:
代码实现
上述系统主要通过 Auto_Ctrl_System 类实现,该类的实例属性有:
- A, B, C, D:系统矩阵 ,输入矩阵 ,输出矩阵 ,直接传递矩阵
- w:权重向量
- eig_coef:矩阵 的特征多项式的系数
- dstate, output:状态方程、输出方程的结果
实例方法有:
- adjust:根据状态方程,调整权重向量;同时根据 fine_grit 参数对误差值之和进行采样
- predict:利用当前的权重向量,对样本 (可另给定) 给出预测值
- ctrl:转换系统为能控标准型
- jordan:转换系统为 Jordan 标准型 (该系统中的矩阵 A 为对称矩阵,可对角化), 由矩阵 的特征向量构成
- state_tran:利用矩阵 可被对角化这一点,对该系统的状态转移矩阵进行求解
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()