Information Theory | Kullback-Leibler Divergence

概述

KL 散度（Kullback-Leibler Divergence），也称为相对熵（Relative Entropy），是信息论和统计学中用于衡量两个概率分布差异的重要度量。它在机器学习、深度学习、变分推断等领域有广泛应用。
本笔记重点介绍 KL 散度的基础定义、数学推导、核心性质以及其在机器学习中的意义。

数学定义

离散概率分布

对于离散概率分布 $P$ 和 $Q$，KL 散度定义为：

$$D_{KL}(P||Q) = \sum_{x} P(x) \log\frac{P(x)}{Q(x)}$$

连续概率分布

对于连续概率分布，KL 散度定义为：

$$D_{KL}(P||Q) = \int p(x) \log\frac{p(x)}{q(x)} dx$$

其中：

• $P$ 是真实分布（或目标分布）
• $Q$ 是近似分布（或模型分布）
• $\log$ 通常使用自然对数（以 $e$ 为底），单位为 nats；使用以 2 为底的对数时，单位为 bits

核心性质

1. 非负性

$$D_{KL}(P||Q) \geq 0$$

当��仅当 $P = Q$ 时，等号成立（即 $D_{KL}(P||Q) = 0$）。

2. 非对称性

$$D_{KL}(P||Q) \neq D_{KL}(Q||P)$$

这意味着 KL 散度不是真正的距离度量，因为它不满足对称性。

3. 不满足三角不等式

KL 散度不满足三角不等式：

$$D_{KL}(P||R) \not\leq D_{KL}(P||Q) + D_{KL}(Q||R)$$

因此，KL 散度不是度量空间中的距离函数。

4. 凸性

KL 散度关于第一个参数 $P$ 是凸函数。

直观理解

信息论角度

KL 散度表示：用分布 $Q$ 来编码分布 $P$ 的样本时，相比用 $P$ 自身编码所需的额外信息量（以比特或 nats 为单位）。

• 如果 $Q$ 与 $P$ 完全相同，则不需要额外信息，$D_{KL}(P||Q) = 0$
• 如果 $Q$ 与 $P$ 差异很大，则需要更多额外信息来补偿编码效率的损失

统计学角度

KL 散度衡量分布 $Q$ 对分布 $P$ 的近似程度：

• 值越小，$Q$ 越接近 $P$
• 值越大，$Q$ 与 $P$ 的差异越大

机器学习角度

在机器学习中，通常：

• $P$ 是真实数据分布（未知但可以从数据中采样）
• $Q$ 是模型学习的分布（由参数 $\theta$ 决定）
• 最小化 $D_{KL}(P||Q)$ 等价于让模型分布尽可能接近真实分布

数学推导

1. 从信息熵推导

信息熵（Entropy）

分布 $P$ 的信息熵定义为：

$$H(P) = -\sum_{x} P(x) \log P(x)$$

它表示编码分布 $P$ 的样本所需的平均最小比特数。

交叉熵（Cross Entropy）

用分布 $Q$ 来编码分布 $P$ 的样本所需的平均比特数：

$$H(P, Q) = -\sum_{x} P(x) \log Q(x)$$

KL 散度 = 交叉熵 - 信息熵

$$\begin{align}
D_{KL}(P||Q) &= H(P, Q) - H(P) \
&= -\sum_{x} P(x) \log Q(x) + \sum_{x} P(x) \log P(x) \
&= \sum_{x} P(x) [\log P(x) - \log Q(x)] \
&= \sum_{x} P(x) \log\frac{P(x)}{Q(x)}
\end{align}$$

这表明 KL 散度是使用次优编码方案 $Q$ 相比最优编码方案 $P$ 所需的额外信息量。

2. 非负性证明（Gibbs 不等式）

我们使用 Jensen 不等式来证明 KL 散度的非负性。

Jensen 不等式

对于凸函数 $f(x)$ 和概率分布 $P$：

$$f\left(\sum_{x} P(x) \cdot x\right) \leq \sum_{x} P(x) \cdot f(x)$$

对于凹函数（如 $\log$），不等号反向。

证明过程

由于 $f(x) = -\log(x)$ 是凸函数，我们有：

$$\begin{align}
D_{KL}(P||Q) &= \sum_{x} P(x) \log\frac{P(x)}{Q(x)} \
&= -\sum_{x} P(x) \log\frac{Q(x)}{P(x)} \
&\geq -\log\left(\sum_{x} P(x) \cdot \frac{Q(x)}{P(x)}\right) \quad \text{[Jensen 不等式]} \
&= -\log\left(\sum_{x} Q(x)\right) \
&= -\log(1) \
&= 0
\end{align}$$

等号成立当且仅当 $\frac{Q(x)}{P(x)}$ 为常数，即 $P(x) = Q(x)$ 对所有 $x$ 成立。

结论：$D_{KL}(P||Q) \geq 0$，且仅当 $P = Q$ 时等号成立。

3. 与最大似然估计的关系

在机器学习中，我们通常希望找到参数 $\theta$ 使得模型分布 $P_\theta$ 尽可能接近真实数据分布 $P_{data}$。

最小化 KL 散度

$$\begin{align}
\min_\theta D_{KL}(P_{data}||P_\theta) &= \min_\theta \sum_{x} P_{data}(x) \log\frac{P_{data}(x)}{P_\theta(x)} \
&= \min_\theta \left[\sum_{x} P_{data}(x) \log P_{data}(x) - \sum_{x} P_{data}(x) \log P_\theta(x)\right] \
&= \min_\theta \left[-\sum_{x} P_{data}(x) \log P_\theta(x)\right] \quad \text{[第一项与 $\theta$ 无关]} \
&= \max_\theta \sum_{x} P_{data}(x) \log P_\theta(x)
\end{align}$$

经验分布近似

在实践中，我们无法直接访问 $P_{data}$，但可以从数据集 ${x_1, x_2, \ldots, x_N}$ 中采样。使用经验分布：
$$P_{data}(x) \approx \frac{1}{N} \sum_{i=1}^{N} \delta(x - x_i)$$

代入上式：

$$\max_\theta \sum_{x} P_{data}(x) \log P_\theta(x) \approx \max_\theta \frac{1}{N} \sum_{i=1}^{N} \log P_\theta(x_i)$$

这正是最大似然估计（Maximum Likelihood Estimation, MLE）的目标函数！
结论：最小化 KL 散度等价于最大化似然函数。

实际计算示例

示例 1：离散分布

假设真实分布 $P = [0.5, 0.3, 0.2]$，近似分布 $Q = [0.4, 0.4, 0.2]$。

计算 $D_{KL}(P||Q)$：

$$\begin{align}
D_{KL}(P||Q) &= 0.5 \cdot \log\frac{0.5}{0.4} + 0.3 \cdot \log\frac{0.3}{0.4} + 0.2 \cdot \log\frac{0.2}{0.2} \
&= 0.5 \cdot \log(1.25) + 0.3 \cdot \log(0.75) + 0.2 \cdot \log(1) \
&= 0.5 \times 0.2231 + 0.3 \times (-0.2877) + 0 \
&= 0.1116 - 0.0863 \
&\approx 0.0253 \text{ nats}
\end{align}$$

转换为 bits（除以 $\ln 2 \approx 0.693$）：

$$D_{KL}(P||Q) \approx 0.0365 \text{ bits}$$

示例 2：高斯分布

对于两个一维高斯分布：

• $P = \mathcal{N}(\mu_1, \sigma_1^2)$
• $Q = \mathcal{N}(\mu_2, \sigma_2^2)$

KL 散度的闭式解为：

$$D_{KL}(P||Q) = \log\frac{\sigma_2}{\sigma_1} + \frac{\sigma_1^2 + (\mu_1 - \mu_2)^2}{2\sigma_2^2} - \frac{1}{2}$$

特殊情况：当 $Q = \mathcal{N}(0, 1)$（标准正态分布）时：

$$D_{KL}(P||Q) = \frac{1}{2}\left(\mu_1^2 + \sigma_1^2 - \log\sigma_1^2 - 1\right)$$

这个公式在 VAE（变分自编码器）中被广泛使用。

机器学习中的应用

1. 变分推断（Variational Inference）

在贝叶斯推断中，我们希望计算后验分布 $P(\theta|X)$，但通常难以直接计算。变分推断使用一个简单的分布 $Q(\theta)$ 来近似：
$$\min_Q D_{KL}(Q(\theta)||P(\theta|X))$$

2. 变分自编码器（VAE）

VAE 的损失函数包含 KL 散度项，用于正则化潜在空间：

$$\mathcal{L} = \mathbb{E}{q(z|x)}[\log p(x|z)] - D{KL}(q(z|x)||p(z))$$

其中：

• 第一项是重建损失
• 第二项是 KL 散度，约束编码器输出接近先验分布

3. 策略优化（强化学习）

在强化学习中，KL 散度用于约束策略更新的幅度：
$$\max_\theta \mathbb{E}{\pi_\theta}[R] \quad \text{s.t.} \quad D{KL}(\pi_{old}||\pi_\theta) \leq \delta$$

这是 TRPO（Trust Region Policy Optimization）和 PPO（Proximal Policy Optimization）的核心思想。

4. 模型蒸馏（Knowledge Distillation）

使用 KL 散度让小模型（学生）学习大模型（教师）的输出分布：

$$\mathcal{L} = D_{KL}(P_{teacher}||P_{student})$$

5. 生成对抗网络（GAN）

虽然 GAN 不直接优化 KL 散度，但理论分析表明，GAN 的目标函数与 JS 散度（Jensen-Shannon Divergence）相关，而 JS 散度是基于 KL 散度定义的：

$$D_{JS}(P||Q) = \frac{1}{2}D_{KL}(P||M) + \frac{1}{2}D_{KL}(Q||M)$$

其中 $M = \frac{1}{2}(P + Q)$。

前向 KL 与反向 KL

前向 KL：$D_{KL}(P||Q)$

• 含义：用 $Q$ 近似 $P$
• 特点：
  • 当 $P(x) > 0$ 但 $Q(x) = 0$ 时，散度为无穷大
  • 要求 $Q$ 覆盖 $P$ 的所有支撑集（zero-avoiding）
  • $Q$ 倾向于覆盖 $P$ 的所有模式，可能过于分散

反向 KL：$D_{KL}(Q||P)$

• 含义：用 $Q$ 近似 $P$（但优化方向相反）
• 特点：
  • 当 $Q(x) > 0$ 但 $P(x) = 0$ 时，散度为无穷大
  • 要求 $Q$ 集中在 $P$ 的高概率区域（zero-forcing）
  • $Q$ 倾向于选择 $P$ 的某一个模式，可能过于集中

对比表格

与其他散度的关系

1. JS 散度（Jensen-Shannon Divergence）

JS 散度是 KL 散度的对称化版本：

$$D_{JS}(P||Q) = \frac{1}{2}D_{KL}(P||M) + \frac{1}{2}D_{KL}(Q||M)$$
其中 $M = \frac{1}{2}(P + Q)$。
性质：

• 对称：$D_{JS}(P||Q) = D_{JS}(Q||P)$
• 有界：$0 \leq D_{JS}(P||Q) \leq \log 2$
• 满足三角不等式的平方根

2. Wasserstein 距离

Wasserstein 距离（也称为 Earth Mover’s Distance）是另一种衡量分布差异的度量，在 GAN 中被广泛使用（WGAN）。
与 KL 散度相比：

• Wasserstein 距离是真正的度量（满足对称性和三角不等式）
• 即使两个分布没有重叠，Wasserstein 距离仍然有意义
• KL 散度在分布不重叠时可能为无穷大

3. $\chi^2$ 散度

$$D_{\chi^2}(P||Q) = \sum_{x} \frac{(P(x) - Q(x))^2}{Q(x)}$$
与 KL 散度相比，$\chi^2$ 散度对分布差异更敏感。

常见问题

问题 1：为什么 KL 散度不对称？

原因：KL 散度的定义中，$P$ 和 $Q$ 的角色不同：

$$D_{KL}(P||Q) = \sum_{x} P(x) \log\frac{P(x)}{Q(x)}$$

• $P$ 作为权重出现在求和中
• $Q$ 出现在对数的分母中

交换 $P$ 和 $Q$ 会得到不同的结果。

直观理解：

• $D_{KL}(P||Q)$ 衡量”用 $Q$ 编码 $P$ 的样本”的额外信息
• $D_{KL}(Q||P)$ 衡量”用 $P$ 编码 $Q$ 的样本”的额外信息
• 这两者通常不相等

问题 2：什么时候使用前向 KL，什么时候使用反向 KL？

**前向 KL $D_{KL}(P||Q)$**：

• 当你希望 $Q$ 覆盖 $P$ 的所有模式时
• 最大似然估计
• 监督学习
  **反向 KL $D_{KL}(Q||P)$**：
• 当你希望 $Q$ 集中在 $P$ 的高概率区域时
• 变分推断
• VAE、变分贝叶斯

问题 3：KL 散度可以为负吗？

不可以。根据 Gibbs 不等式，KL 散度始终非负：

$$D_{KL}(P||Q) \geq 0$$
当且仅当 $P = Q$ 时，$D_{KL}(P||Q) = 0$。

问题 4：如何处理 $Q(x) = 0$ 的情况？

当 $P(x) > 0$ 但 $Q(x) = 0$ 时，$\log\frac{P(x)}{Q(x)} = \infty$，导致 KL 散度为无穷大。

解决方法：

1. 拉普拉斯平滑：给 $Q(x)$ 加一个小的常数 $\epsilon$
2. 确保支撑集匹配：保证 $Q$ 的支撑集包含 $P$ 的支撑集
3. 使用其他散度：如 Wasserstein 距离，对不重叠分布更鲁棒

最佳实践

1. 数值稳定性

在实现 KL 散度时，避免直接计算 $\log\frac{P(x)}{Q(x)}$，而是使用：

1 2 3 4 5

# 不稳定的实现 kl = P * np.log(P / Q) # 稳定的实现 kl = P * (np.log(P) - np.log(Q))

### 2. 处理零概率

1 2 3

# 添加小的 epsilon 避免 log(0) epsilon = 1e-10 kl = P * (np.log(P + epsilon) - np.log(Q + epsilon))

### 3. 使用库函数
大多数深度学习框架提供了优化的 KL 散度实现：

1 2 3 4 5 6 7 8 9 10 11

# PyTorch import torch.nn.functional as F kl_loss = F.kl_div(Q.log(), P, reduction='batchmean') # TensorFlow import tensorflow as tf kl_loss = tf.keras.losses.KLDivergence()(P, Q) # SciPy from scipy.stats import entropy kl_div = entropy(P, Q)

参考资源

• Information Theory, Inference, and Learning Algorithms - David MacKay
• Pattern Recognition and Machine Learning - Christopher Bishop
• Deep Learning Book - Ian Goodfellow et al.
• Kullback-Leibler Divergence Explained

Project Proposal: 2D Point Cloud Registration using the Iterative Closest Point (ICP) Algorithm

Course: VV570/MATH6003J: Introduction to Engineering Numerical Analysis 1

1. Introduction & Problem Statement

In robotic perception, sensors like LiDAR capture the environment as a “point cloud”—a set of data points in space. To enable tasks like mapping or localization, a robot must align multiple point clouds taken from different positions. This project will address this fundamental robotics problem by implementing the Iterative Closest Point (ICP) algorithm, a numerical optimization method used to find the optimal alignment between two point clouds2.

2. Objective

The primary objective of this project is to develop a working 2D implementation of the ICP algorithm. We will apply numerical analysis techniques to find the rigid transformation (rotation and translation) that minimizes the distance between two point clouds, demonstrating a practical application of concepts from the course to the field of robotics.

3. Methodology

Our approach will be to apply a computational tool to solve this problem3. We will use Python with the NumPy library for numerical computations and Matplotlib for visualization. The core steps of our implementation will be:

• Data Association: For each point in a source cloud, find the closest corresponding point in the target cloud.

• Numerical Optimization: Calculate the optimal rotation and translation that minimizes the mean squared error between the corresponding point pairs. This will be solved using Singular Value Decomposition (SVD).

• Iterative Refinement: Apply the calculated transformation to the source cloud and repeat the process until the alignment error converges below a set threshold.

We will start with the simplification of working in 2D and assume a reasonable initial alignment between the point clouds4.

4. Expected Outcomes & Justification of Results

By the end of this three-week project5, we will deliver:

• A Python implementation of the 2D ICP algorithm.

• Visualizations demonstrating the algorithm’s effectiveness by showing the point clouds before and after alignment.

• A plot of the alignment error over iterations, which will serve to justify our results by showing the convergence of the numerical method6.

Ex 1

a

The result of hw4_1_a.m:

Vandermonde matrices are known to be ill-conditioned, with the condition number growing exponentially with $n$.

d

This tridiagonal matrix is well-conditioned. Its condition number grows polynomially with $n$ (like $O(n^2)$). The error $\epsilon_n$ should grow much more slowly than in the Vandermonde case, exhibiting polynomial rather than exponential growth. The residual $r_n$ will be a better, though not perfect, indicator of the true error compared to the previous case.

Comments on the obtained results:

• $\textbf{Error Growth:}$ The relative error $\epsilon_n$ for the tridiagonal system grows polynomially with $n$, which is significantly slower than the exponential growth observed for the Vandermonde matrix. This will be visible as a straight line on the loglog plot.
• $\textbf{Conditioning:}$ The dramatic difference in error behavior is due to the condition number. The Vandermonde matrix is severely ill-conditioned, whereas the tridiagonal matrix is well-conditioned.
• $\textbf{Residual as an Error Indicator:}$ For the well-conditioned tridiagonal matrix, the normalized residual $r_n$ and the relative error $\epsilon_n$ have much closer values. The residual is a far more reliable indicator of the true error in this case than it was for the Vandermonde system.

Ex 2

a

The task is to complete the files $jacobi.m$ and $gauss_seidel.m$. These files implement iterative methods based on a matrix splitting $A = P - (P-A)$, where $P$ is the preconditioner. The iterative update can be expressed as $x^{(k+1)} = x^{(k)} + z^{(k)}$, where $z^{(k)}$ is the solution to $Pz^{(k)} = r^{(k)}$, and $r^{(k)} = b - Ax^{(k)}$ is the residual.

For the $\textbf{Jacobi method}$, the preconditioner P is the diagonal of A. For the $\textbf{Gauss-Seidel method}$, P is the lower triangular part of A. The missing lines implement the calculation of the residual, the update step, and the storing of the residual norm at each iteration.

jacobi.m:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

function [x, iter, res]= jacobi(A,b,x0,nmax,tol)
% JACOBI iterative method
%   [X, Niter, Res] = JACOBI(A,B,X0,NMAX,TOL) attempts to solve the system of linear equations A*X=B
%   for X using the Jacobi method.

% Jacobi preconditioner
P=diag(diag(A)); %

% residual normalization
r0=norm(b);
if r0==0, r0=1; end

% first iteration
x=x0;
r=b-A*x; % The residual b - A*x
res(1)=norm(r); % The norm of the initial residual
iter=0;

while res(end)/r0 > tol && iter < nmax
    z=P\r; % Solve Pz=r
    x=x+z; % Update the solution
    r=b-A*x; % Recompute the residual
    iter=iter+1;
    res(iter+1)=norm(r); % Store the norm of the new residual
end

return

gauss_seidel.m:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

function [x, iter, res]= gauss_seidel(A,b,x0,nmax,tol)
% GAUSS-SEIDEL iterative method
%   [X, Niter, Res] = GAUSS_SEIDEL(A,B,X0,NMAX,TOL) attempts to solve the system of linear equations A*X=B
%   for X using the Gauss-Seidel method.

% Gauss-Seidel preconditioner
P=tril(A); %

% residual normalization
r0=norm(b);
if r0==0, r0=1; end

% first iteration
x=x0;
r=b-A*x; % The residual b - A*x
res(1)=norm(r); % The norm of the initial residual
iter=0;

while res(end)/r0 > tol && iter < nmax
    z=P\r; % Solve Pz=r using forward substitution
    x=x+z; % Update the solution
    r=b-A*x; % Recompute the residual
    iter=iter+1;
    res(iter+1)=norm(r); % Store the norm of the new residual
end

return

b

I solve the system for $n=4$ with the tridiagonal matrix $A$ having $2.5$ on the diagonal and $-1$ on the off-diagonals. The right-hand side is $b=(1.5,0.5,0.5,1.5)^{\top}$. The initial guess is $x^{(0)}=0$, tolerance is $10^{-10}$, and max iterations is $10^8$. I compare the number of iterations and computing time.

The matrix $A$ is strictly diagonally dominant since $|a_{ii}| = 2.5 > \sum_{j \neq i} |a_{ij}|$ for all rows. This property guarantees that both methods will converge.

For this well-conditioned matrix (as shown in part e), the number of iterations grows very slowly with $n$. The condition number remains small and constant. The relative error and the normalized residual are very close in value, indicating that the residual is a good estimator of the error.

d

The analysis from (c) is repeated for the matrix $A$ with 2 on the diagonal and -1 on the off-diagonals , and $b=(1,0,…,0,1)^{\top}$. The exact solution is again $x=(1,…,1)^{\top}$. This matrix is only weakly diagonally dominant, so I expect slower convergence.

The number of iterations will grow much more rapidly with $n$. This is because, as shown in part (e), the condition number of this matrix grows quadratically with $n$, making the problem progressively harder to solve.

e

I use the given formula for the eigenvalues of a symmetric tridiagonal matrix $M$ with $a$ on the diagonal and $b$ on the off-diagonals:
$$ \lambda_{i}(M)=a+2b\cos\left(\frac{\pi i}{n+1}\right), \quad i=1,…,n $$
The condition number is $\kappa_2(A) = |\lambda_{max}|/|\lambda_{min}|$.

Case 1: Matrix from part (c)
Here, $a=2.5$ and $b=-1$. The eigenvalues are $\lambda_i = 2.5 - 2\cos(\frac{\pi i}{n+1})$.

• As $n\to\infty$, $\lambda_{max} \to 2.5 - 2(-1) = 4.5$.
• As $n\to\infty$, $\lambda_{min} \to 2.5 - 2(1) = 0.5$.

The condition number $\kappa(A)$ approaches a constant, $\kappa(A) \to 4.5 / 0.5 = 9$. This is consistent with the numerical results from part (c), explaining the fast and stable convergence.

Case 2: Matrix from part (d)
Here, $a=2$ and $b=-1$. The eigenvalues are $\lambda_i = 2 - 2\cos(\frac{\pi i}{n+1}) = 4\sin^2(\frac{\pi i}{2(n+1)})$.

• As $n\to\infty$, $\lambda_{max} \to 4\sin^2(\pi/2) = 4$.
• For large $n$, $\lambda_{min} = 4\sin^2(\frac{\pi}{2(n+1)}) \approx 4(\frac{\pi}{2(n+1)})^2 = \frac{\pi^2}{(n+1)^2}$.
  The condition number $\kappa(A) \approx \frac{4}{\pi^2/(n+1)^2} = O(n^2)$. This quadratic growth in the condition number is consistent with the numerical results from part (d), explaining the significant increase in iterations with $n$.

Exercise 1

Problem Analysis

The objective is to find the coefficients $(\alpha, \beta, \gamma)$ for the finite difference formula:
$$
D_{h}f(\overline{x})=\frac{\alpha f(\overline{x})+\beta f(\overline{x}-h)+\gamma f(\overline{x}-2h)}{h} \quad
$$
The analysis begins by substituting the Taylor series expansions for $f(\overline{x}-h)$ and $f(\overline{x}-2h)$ around the point $\overline{x}$ into the formula.

The expansions are:
$$
f(\overline{x}-h) = f(\overline{x}) - hf’(\overline{x}) + \frac{h^2}{2}f’’(\overline{x}) - \frac{h^3}{6}f’’’(\overline{x}) + O(h^4)
$$
$$
f(\overline{x}-2h) = f(\overline{x}) - 2hf’(\overline{x}) + 2h^2f’’(\overline{x}) - \frac{4h^3}{3}f’’’(\overline{x}) + O(h^4)
$$
Substituting these into the formula and grouping terms by derivatives of $f(\overline{x})$ results in:
$$
D_{h}f(\overline{x}) = \frac{\alpha + \beta + \gamma}{h} f(\overline{x}) + (-\beta - 2\gamma) f’(\overline{x}) + \left( \frac{\beta}{2} + 2\gamma \right)h f’’(\overline{x}) + O(h^2)
$$
To approximate $f’(\overline{x})$, the coefficient of $f’(\overline{x})$ must be $1$, and the coefficient of $f(\overline{x})$ must be $0$. This yields two necessary equations:

1. $\alpha + \beta + \gamma = 0$
2. $-\beta - 2\gamma = 1$

a) Solution for First-Order Accuracy

For the formula to be of order 1, the two fundamental equations must be satisfied.

From equation (2), $\beta$ can be expressed in terms of $\gamma$:
$$
\beta = -1 - 2\gamma
$$
Substituting this into equation (1):
$$
\alpha + (-1 - 2\gamma) + \gamma = 0 \implies \alpha = 1 + \gamma
$$
Therefore, the family of coefficients $(\alpha, \beta, \gamma)$ that provides at least first-order accuracy is given by:
$$
(\alpha, \beta, \gamma) = (1+\gamma, -1-2\gamma, \gamma), \quad \forall \gamma \in \mathbb{R}
$$

b) Solution for Second-Order Accuracy

To achieve second-order accuracy, the truncation error must be $O(h^2)$. This requires the coefficient of the $h$ term in the error expansion to also be zero.

$$\frac{\beta}{2} + 2\gamma = 0$$

This creates a system of three linear equations to be solved for the unique values of $(\alpha, \beta, \gamma)$:

1. $\alpha + \beta + \gamma = 0$
2. $-\beta - 2\gamma = 1$
3. $\frac{\beta}{2} + 2\gamma = 0$

From equation (3), we find that $\beta = -4\gamma$. Substituting this into equation (2):
$$
-(-4\gamma) - 2\gamma = 1 \implies 2\gamma = 1 \implies \gamma = \frac{1}{2}
$$
With the value for $\gamma$, $\beta$ can be found:
$$
\beta = -4 \left( \frac{1}{2} \right) = -2
$$
Finally, substituting $\beta$ and $\gamma$ into equation (1) yields $\alpha$:
$$
\alpha + (-2) + \frac{1}{2} = 0 \implies \alpha = \frac{3}{2}
$$
The unique values which give a formula of order 2 are:
$$
\left(\alpha, \beta, \gamma\right) = \left(\frac{3}{2}, -2, \frac{1}{2}\right)
$$

Exercise 2

a) Midpoint Formula Calculation

The first task is to compute the integral of the function $f(x)=e^{-x}sin(x)$ on the interval $[a,b]=[0,2]$. The computation uses the composite midpoint formula with $M=10$ sub-intervals.

The resulting numerical approximation of the integral is 0.468592. The exact integral is 0.466630. So the error is 0.001962.

b) Convergence of the Midpoint Formula

This section analyzes the convergence of the composite midpoint formula for the same function and interval. The number of sub-intervals is varied as $M=10^{1},10^{2},10^{3},…10^{5}$.

The absolute error, $|I(f)-Q_{h}^{pm}(f)|$, is computed for each corresponding step size, $h=(b-a)/M$. This error is the difference between the numerical result and the exact value of the integral, $I=(1-e^{-2}(sin(2)+cos(2)))/2$.

A graph of the error versus the step size $h$ is then created using a logarithmic scale for both axes. In such a plot, the slope of the resulting line corresponds to the order of convergence.

The plot displays a straight line with a negative slope. A linear fit to the logarithmic data reveals a slope of approximately 2. This visually confirms that the method exhibits 2nd-order convergence, which aligns with the theoretical error bounds for the composite midpoint rule.

c) Convergence of Trapezoidal and Simpson’s Formulas

The analysis is repeated using the composite trapezoidal formula and the composite Simpson’s formula. The errors for both methods, $|I(f)-Q_{h}^{trap}(f)|$ and $|I(f)-Q_{h}^{simp}(f)|$, are plotted on the same graph in logarithmic scale.

The resulting graph contains two distinct straight lines, one for each method.

• Trapezoidal Rule: The error line for the trapezoidal rule has a slope of approximately 2. This confirms its theoretical 2nd-order convergence.
• Simpson’s Rule: The error line for Simpson’s rule is significantly steeper, with a slope of approximately 4. This demonstrates its theoretical 4th-order convergence.

A comparison of the results shows that for any given step size $h$, the error from Simpson’s rule is substantially smaller than the error from the trapezoidal rule, highlighting the superior accuracy of the higher-order method for smooth functions.

d) Analysis for a Non-Smooth Function

Finally, the convergence analysis is repeated for all three methods on a new function, $f(x)=\sqrt{|x|^{3}}$, over the interval $[a,b]=[-2,2]$. The exact value for this integral is $I=\frac{16}{5}\sqrt{2}$.

Comments on the Obtained Results

The function $f(x)=\sqrt{|x|^{3}}$ is not sufficiently smooth at the point $x=0$ within the integration interval. While the function and its first derivative are continuous, its second derivative, $f’’(x) \propto |x|^{-1/2}$, is unbounded at $x=0$.

The theoretical error estimates that predict 2nd and 4th-order convergence for these methods are based on the assumption that the function’s higher-order derivatives are continuous and bounded. Since this condition is violated, a degradation in the observed convergence rates is expected.

Book

Matrix Cook Book

Autonomous Driving

• Autonomous Driving
  • Perception
  • Decision-making
    • Architecture
      • Pre-trip decisions
      • In-trip decisions
      • Post-trip decision
  • Speed tracking
    • Dynamic system
      • Open-loop control
      • Closed-loop control
        • Linear controller
    • LTI system
      • Control of LTI system
  • Speed tracking
  • Trajectory tracking
    • Control policy

Perception

• Navigation tool
• Camera
• Radar
• LiDAR
• Infrared detector
• Speedometer
• veh-veh communication
• veh-infrastructure communication

Decision-making

• departure time
• good route
• good lane
• driving style
• way to accel/decel
• make a turn
• park location

Architecture

• Pre-trip

  trip generation, modal choice, routine choice

• In-trip

  Path planning, Maneuver regulation

• Post-trip

  Parking

Pre-trip decisions

• Centralized

  central command center send instructions
  Some win, some lose

• Decentralized

  plan trip independently, selfish
  little communication
  Sometime efficient, sometime inefficient

In-trip decisions

• offline info

  static, in form of map

• online info

  enable the route to be adapted to changing traffic conditions

Post-trip decision

Speed tracking

Dynamic system

• State Variable

  $v[t]$

• State space

  $R>= 0$

• Control input

  $u[t]$(accel)

• System dynamic

  $v[t+dt] = v[t]+u[t]*dt$

Open-loop control

Specify the control input $u[t]$ for all t in $[0,T]$
not a good idea for autonomous driving.

In practice, we need to consider the uncertainty of the system.
We should have:
$$
v[t+dt]=v[t]+u[t]dt+w(t)
$$
Where w(t) is the noise term capturing the uncertainty of the system.

Closed-loop control

Specify the control input $u[t]$ as a function of the state $v[t]$
able to adapt to the uncertainty of the system.
Compare the state $v[t]$ with the desired state $v_desired[t]$
Essentially, we need a mapping from state to control input
like map speed to acceleration
$$
v[t+dt]=v[t]+mu[v[t]]dt+w(t)
$$

Linear controller

the $u[t]$ has linear relationship with $v[t]$
$$
\mu[v] = k[v-v_{desired}]
$$

LTI system

linear time-invariant

Control of LTI system

exponential convergence is stronger than asymptotic convergent
if LTI system is asymptotic convergent, it is also exponential convergent
$$
x[t+1]=ax[t]+bu[t]
$$

Speed tracking

• given:

  $v_{desired}, \delta$

• determine:

  $u[t]$ $t = 0,1,2$

• state

  $v[t]$
  $v[t+1] = v[t] + u[t]\delta + w[t]$

select $u[t]$ so that $w[t]$ will not accumulate

$u[t]=\mu(v[t]) = kv[t]$

$$
v[t+1] = f(v[t],u[t])
$$

Remember the review question

Trajectory tracking

Overview:

• track means asymptotic convergent

• $x[t]$ and $x_{desired}[t]$

• $\lim_{t\rightarrow\infty}|x[t]-x_{desired}[t] | = 0$

• state of system:$[x[t], v[t]]^T$

• if successful:
  $$
  [x[t], v[t]]^T \rightarrow [x_{desired}[t], v_{desired}[t]]^T
  $$

Control policy

• A control policy is a function
• This function maps a state to a control input
  $$
  u[t] = f(x[t],v[t])
  $$

$$
u[t] = -k_1(x[t]-x_{desired}[t])-k_2(v[t]-v_{desired}[t])
$$

$$
k_1,k_2>0
$$