SciPy 中的鲁棒非线性回归¶

日期	2018-08-17（上次修改）

非线性最小二乘的主要应用之一是非线性回归或曲线拟合。也就是说，给定对 $\left\{ (t_i, y_i) \: i = 1, \ldots, n \right\}$，估计定义非线性函数 $\varphi(t; \mathbf{x})$ 的参数 $\mathbf{x}$，假设模型：\begin{equation} y_i = \varphi(t_i; \mathbf{x}) + \epsilon_i \end{equation}

其中 $\epsilon_i$ 是测量（观察）误差。在最小二乘估计中，我们将 $\mathbf{x}$ 作为以下优化问题的解：\begin{equation} \frac{1}{2} \sum_{i=1}^n (\varphi(t_i; \mathbf{x}) - y_i)^2 \rightarrow \min_\mathbf{x} \end{equation}

这种公式从数学角度来看是直观的和方便的。从概率的角度来看，最小二乘解被称为最大似然估计，前提是所有 $\epsilon_i$ 都是独立且正态分布的随机变量。

因此，从理论上讲，当 $\epsilon_i$ 的分布不是正态分布时，它不是最优的。尽管在工程实践中通常并不重要，即如果误差表现为一些具有零均值的合理随机变量，最小二乘估计的结果将是令人满意的。

当数据被异常值（完全错误的测量值）污染时，真正的问题就开始了。在这种情况下，最小二乘解可能会被显著地偏向，以避免异常值上的非常高的残差。为了定性地解释为什么会出现这种情况，让我们考虑最简单的最小二乘问题：\begin{align} \frac{1}{2} \sum_{i=1}^n (a - x_i)^2 \rightarrow \min_a \end{align} 并且解是平均值：\begin{align} a^* = \frac{1}{n} \sum_{i=1}^n x_i \end{align}

现在想象：$a = 1, n=4, x_1 = 0.9, x_2 = 1.05, x_3=0.95, x_4=10 \text{ (outlier) }$. 解 $a^* = 3.225$，即被单个异常值完全破坏。

众所周知的鲁棒估计器之一是 l1 估计器，其中最小化残差的绝对值的总和。为了演示，再次考虑最简单的问题：\begin{equation} \sum_{i=1}^n |a - x_i| \rightarrow \min_a \end{equation}

众所周知的解决方案是 $\{x_i\}$ 的中位数，这样少量异常值就不会影响解决方案。在我们的玩具问题中，我们有 $a^* = 1$。

l1 估计器的唯一缺点是产生的优化问题很困难，因为该函数在任何地方都不可微，这对于高效的非线性优化来说尤其麻烦。这意味着我们最好坚持可微问题，但以某种方式将鲁棒性纳入估计。为了实现这一点，我们引入了一个亚线性函数 $\rho(z)$（即其增长应该比线性慢）并制定了一个新的类似最小二乘的优化问题：\begin{equation} \frac{1}{2} \sum_{i=1}^n \rho \left((\varphi(t_i; x) - y_i)^2 \right) \rightarrow \min_x \end{equation}

事实证明，这个问题可以通过在每次迭代中修改残差向量和雅可比矩阵来简化为标准非线性最小二乘，这样计算出的梯度和海森矩阵近似值与目标函数的梯度和海森矩阵近似值相匹配。有关详细信息，请参阅论文。

实现的 $\rho$ 选择如下

线性函数，它给出了标准最小二乘：$\rho(z) = z$。
Huber 损失：$\rho(z) = \begin{cases} z & z \leq 1 \\ \sqrt{z} - 1 & z > 1 \end{cases}$
绝对值损失的平滑近似，“软 l1 损失”：$\rho(z) = 2 (\sqrt{1 + z} - 1)$
柯西损失：$\rho(z) = \ln(1 + z)$。
反正切损失：$\rho(z) = \arctan z$。

函数 2 和 3 相对温和，对于较大的残差，它们大约给出了绝对值损失。最后两个函数是强亚线性的，并且对异常值进行了显着的衰减。

上面的损失函数是在假设内点和异常值之间的软阈值等于 1.0 的情况下编写的。为了泛化，我们引入缩放参数 $C$ 并评估损失为 \begin{equation} \hat{\rho}(r^2) = C^2 \rho \left( \left(\frac{r}{C} \right)^2 \right) \end{equation}

请注意，如果残差很小，我们有 $\hat{\rho}(r^2) \approx \rho(r^2) \approx r^2$，对于上面定义的任何 $\rho(z)$。

为了说明引入函数的行为，我们构建了一个图

In [1]

import numpy as np

%matplotlib inline
import matplotlib.pyplot as plt

from matplotlib import rcParams
rcParams['figure.figsize'] = (10, 6)
rcParams['legend.fontsize'] = 16
rcParams['axes.labelsize'] = 16

In [2]

r = np.linspace(0, 5, 100)

linear = r**2

huber = r**2
huber[huber > 1] = 2 * r[huber > 1] - 1

soft_l1 = 2 * (np.sqrt(1 + r**2) - 1)

cauchy = np.log1p(r**2)

arctan = np.arctan(r**2)

In [3]

plt.plot(r, linear, label='linear')
plt.plot(r, huber, label='huber')
plt.plot(r, soft_l1, label='soft_l1')
plt.plot(r, cauchy, label='cauchy')
plt.plot(r, arctan, label='arctan')
plt.xlabel("$r$")
plt.ylabel(r"$\rho(r^2)$")
plt.legend(loc='upper left');

Out[3]

<matplotlib.legend.Legend at 0x1058c1ef0>

现在我们将展示鲁棒损失函数如何在模型示例中工作。我们将模型函数定义为 \begin{equation} f(t; A, \sigma, \omega) = A e^{-\sigma t} \sin (\omega t) \end{equation}

它可以模拟线性阻尼振荡器的观测位移。

定义数据生成器

In [4]

def generate_data(t, A, sigma, omega, noise=0, n_outliers=0, random_state=0):
    y = A * np.exp(-sigma * t) * np.sin(omega * t)
    rnd = np.random.RandomState(random_state)
    error = noise * rnd.randn(t.size)
    outliers = rnd.randint(0, t.size, n_outliers)
    error[outliers] *= 35
    return y + error

定义模型参数

In [5]

A = 2
sigma = 0.1
omega = 0.1 * 2 * np.pi
x_true = np.array([A, sigma, omega])

noise = 0.1

t_min = 0
t_max = 30

用于拟合参数的数据将包含 3 个异常值

In [6]

t_train = np.linspace(t_min, t_max, 30)
y_train = generate_data(t_train, A, sigma, omega, noise=noise, n_outliers=4)

定义用于最小二乘最小化的残差计算函数

在 [7]

def fun(x, t, y):
    return x[0] * np.exp(-x[1] * t) * np.sin(x[2] * t) - y

使用全为一的初始估计。

在 [8]

x0 = np.ones(3)

在 [9]

from scipy.optimize import least_squares

运行标准最小二乘

在 [10]

res_lsq = least_squares(fun, x0, args=(t_train, y_train))

运行鲁棒最小二乘，使用loss='soft_l1'，将f_scale设置为0.1，这意味着内点残差大约小于0.1。

在 [11]

res_robust = least_squares(fun, x0, loss='soft_l1', f_scale=0.1, args=(t_train, y_train))

定义数据以绘制完整曲线。

在 [12]

t_test = np.linspace(t_min, t_max, 300)
y_test = generate_data(t_test, A, sigma, omega)

使用找到的参数计算预测值

在 [13]

y_lsq = generate_data(t_test, *res_lsq.x)
y_robust = generate_data(t_test, *res_robust.x)

在 [14]

plt.plot(t_train, y_train, 'o', label='data')
plt.plot(t_test, y_test, label='true')
plt.plot(t_test, y_lsq, label='lsq')
plt.plot(t_test, y_robust, label='robust lsq')
plt.xlabel('$t$')
plt.ylabel('$y$')
plt.legend();

Out[14]

<matplotlib.legend.Legend at 0x107031898>

我们清楚地看到，标准最小二乘对异常值反应很大，并给出非常有偏差的解，而鲁棒最小二乘（使用软 l1 损失）给出的解非常接近真实模型。

章节作者：Nikolay Mayorov