第十三章：线性因子模型

1. 线性因子模型概述

线性因子模型是统计机器学习和深度学习中一类基本的潜变量模型。这些模型假设观测数据 \(\mathbf{x} \in \mathbb{R}^D\) 是由一组潜因子 \(\mathbf{h} \in \mathbb{R}^H\) 通过线性变换生成的，并叠加了加性高斯噪声。其核心思想是：高维数据的复杂性可以通过相对较少的潜因子有效捕捉，每个潜因子代表数据变异的一个有意义的潜在原因。

线性因子模型的一般形式将观测数据在给定潜因子条件下的条件分布指定为：

\[ p(\mathbf{x} | \mathbf{h}) = \mathcal{N}\left(\mathbf{W}\mathbf{h} + \mathbf{b}, \sigma^2 \mathbf{I}\right) \]

其中： - \(\mathbf{W} \in \mathbb{R}^{D \times H}\) 是因子加载矩阵，负责将潜因子映射到观测空间 - \(\mathbf{b} \in \mathbb{R}^D\) 是偏置向量 - \(\sigma^2\) 是各向同性高斯噪声的方差 - \(\mathbf{I}\) 是 \(D \times D\) 的单位矩阵

潜因子上的先验分布通常假设为高斯分布：

\[ p(\mathbf{h}) = \mathcal{N}(\mathbf{0}, \mathbf{I}) \]

这一 formulation 为降维、特征提取和数据表示学习提供了一种有原则的概率框架。线性因子模型是更复杂的深度生成模型的重要构建基础，包括变分自编码器和深度信念网络。

---

2. 概率主成分分析（PPCA）

2.1 模型 formulation

概率主成分分析（Probabilistic PCA，简称 PPCA）由 Tipping 和 Bishop 于 1999 年提出，它为经典 PCA 提供了概率解释，同时扩展了处理缺失数据和概率推理的能力。模型假设观测数据 \(\mathbf{x}_i \in \mathbb{R}^D\) 由 \(H\) 维潜编码 \(\mathbf{h}_i\) 通过以下方式生成：

\[ \mathbf{x}_i = \mathbf{W} \mathbf{h}_i + \boldsymbol{\mu} + \boldsymbol{\epsilon}_i \]

其中： - \(\mathbf{h}_i \sim \mathcal{N}(\mathbf{0}, \mathbf{I}_H)\) 是潜因子 - \(\boldsymbol{\epsilon}_i \sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I}_D)\) 是各向同性噪声，且与 \(\mathbf{h}_i\) 独立 - \(\mathbf{W} \in \mathbb{R}^{D \times H}\) 是加载矩阵，其列向量代表主成分 - \(\boldsymbol{\mu} \in \mathbb{R}^D\) 是均值向量

\(\mathbf{x}_i\) 的边缘分布为：

\[ p(\mathbf{x}_i) = \mathcal{N}\left(\boldsymbol{\mu}, \mathbf{C}\right), \quad \mathbf{C} = \mathbf{W}\mathbf{W}^\top + \sigma^2 \mathbf{I}_D \]

2.2 最大似然估计

观测数据 \(\mathcal{X} = \{\mathbf{x}_i\}_{i=1}^N\) 的对数似然为：

\[ \log p(\mathcal{X} | \mathbf{W}, \boldsymbol{\mu}, \sigma^2) = -\frac{ND}{2}\log(2\pi) - \frac{N}{2}\log|\mathbf{C}| - \frac{1}{2}\sum_{i=1}^N (\mathbf{x}_i - \boldsymbol{\mu})^\top \mathbf{C}^{-1} (\mathbf{x}_i - \boldsymbol{\mu}) \]

最大似然估计的解具有以下关系。最优均值 \(\boldsymbol{\mu}\) 就是样本均值 \(\hat{\boldsymbol{\mu}} = \frac{1}{N}\sum_{i=1}^N \mathbf{x}_i\)。最优加载矩阵 \(\hat{\mathbf{W}}\) 的列向量与样本协方差矩阵 \(\mathbf{S} = \frac{1}{N}\sum_{i=1}^N (\mathbf{x}_i - \hat{\boldsymbol{\mu}})(\mathbf{x}_i - \hat{\boldsymbol{\mu}})^\top\) 的主特征向量成比例：

\[ \hat{\mathbf{W}} = \mathbf{U}_H (\boldsymbol{\Lambda}_H - \sigma^2 \mathbf{I})^{1/2} \mathbf{R} \]

其中 \(\mathbf{U}_H\) 包含 \(\mathbf{S}\) 的前 \(H\) 个特征向量，\(\boldsymbol{\Lambda}_H\) 是对应特征值组成的对角矩阵，\(\mathbf{R} \in \mathbb{R}^{H \times H}\) 是任意正交旋转矩阵。噪声方差的估计为：

\[ \hat{\sigma}^2 = \frac{1}{D - H} \sum_{j=H+1}^D \lambda_j \]

这代表了被丢弃维度的平均未解释方差。

2.3 潜因子的后验分布

给定观测数据后，潜编码的后验分布是高斯分布：

\[ p(\mathbf{h} | \mathbf{x}) = \mathcal{N}\left(\mathbf{M}^{-1} \mathbf{W}^\top (\mathbf{x} - \boldsymbol{\mu}), \sigma^2 \mathbf{M}^{-1}\right) \]

其中 \(\mathbf{M} = \mathbf{W}^\top \mathbf{W} + \sigma^2 \mathbf{I}_H\)。后验均值提供了降维的最优表示：

\[ \mathbb{E}[\mathbf{h} | \mathbf{x}] = \mathbf{M}^{-1} \mathbf{W}^\top (\mathbf{x} - \boldsymbol{\mu}) \]

---

3. 因子分析

3.1 模型 specification

因子分析（Factor Analysis，简称 FA）通过允许异构噪声扩展了 PPCA。与各向同性高斯噪声不同，因子分析假设对角但维度特定的噪声：

\[ p(\mathbf{x} | \mathbf{h}) = \mathcal{N}\left(\mathbf{W}\mathbf{h} + \boldsymbol{\mu}, \boldsymbol{\Psi}\right) \]

其中 \(\boldsymbol{\Psi} = \text{diag}(\psi_1, \psi_2, \ldots, \psi_D)\) 是一个对角协方差矩阵，允许每个观测维度具有其自身的噪声方差。这一 formulation 在不同测量模态或传感器具有不同可靠性时特别适用。

3.2 可辨识性与旋转歧义

因子分析中的一个基本挑战是模型可辨识性。似然函数对潜因子的正交旋转是不变的：

\[ \mathbf{W} \rightarrow \mathbf{W}\mathbf{R}, \quad \mathbf{h} \rightarrow \mathbf{R}^\top \mathbf{h} \]

其中 \(\mathbf{R} \in \mathbb{R}^{H \times H}\) 是任意正交矩阵。这种旋转不确定性意味着潜因子本身不是唯一确定的——只有它们张成的子空间是可辨识的。为了解决这个问题，研究者通常施加旋转准则（例如 varimax 旋转）以获得更具可解释性的因子载荷。

3.3 因子分析与 PPCA 的关系

因子分析与 PPCA 之间的关键区别在于噪声协方差结构：

方面	PPCA	因子分析
噪声模型	各向同性 \(\sigma^2 \mathbf{I}\)	对角 \(\boldsymbol{\Psi} = \text{diag}(\psi_1, \ldots, \psi_D)\)
参数数量	\(DH + D + 1\)	\(DH + 2D\)
计算复杂度	\(\mathcal{O}(ND^2 + DH^2)\)	\(\mathcal{O}(ND^2 + DH^2)\)（迭代）
最优旋转	仅正交	任意可逆线性变换

当所有 \(\psi_d = \sigma^2\) 被约束相等时，因子分析退化为 PPCA。当 \(H = D\) 时，模型与满秩因子分析等价，本质上是直接建模数据协方差。

---

4. 独立成分分析（ICA）

4.1 信息最大化原则

独立成分分析（Independent Component Analysis，简称 ICA）旨在将观测数据 \(\mathbf{x} \in \mathbb{R}^D\) 分解为统计独立的潜源 \(\mathbf{s} \in \mathbb{R}^D\) 的线性组合：

\[ \mathbf{x} = \mathbf{A}\mathbf{s} \]

其中 \(\mathbf{A} \in \mathbb{R}^{D \times D}\) 是混合矩阵。信息最大化原则（Infomax principle）由 Linsker 于 1988 年、Bell 和 Sejnowski 于 1995 年提出，通过最大化编码器网络输入与输出之间的互信息来驱动 ICA。对于单调非线性函数 \(g(\cdot)\)，信息最大化目标为：

\[ \mathcal{L}_{\text{Infomax}} = \sum_{i=1}^N \mathbb{E}\left[\log g'(\mathbf{w}_i^\top \mathbf{x} + b_i)\right] \]

解混矩阵 \(\mathbf{W} = \mathbf{A}^{-1}\) 的梯度更新为：

\[ \Delta \mathbf{W} \propto \left(\mathbf{I} - \tanh(\mathbf{W}\mathbf{x})\mathbf{x}^\top\right)\mathbf{W} \]

4.2 FastICA 算法

FastICA（Hyvärinen 和 Oja，1997）使用定点迭代为信息最大化提供了一种高效替代方案。对于非高斯源信号，主成分分析识别最大方差的方向，而 FastICA 最大化峰度或负熵以找到统计独立的成分。

基于峰度的 FastICA 对于单个成分 \(\mathbf{w}^\top \mathbf{x}\)：

1. 随机初始化 \(\mathbf{w}\)，满足 \(\|\mathbf{w}\| = 1\)
2. 更新：\(\mathbf{w}^+ = \mathbb{E}\left[\mathbf{x}(\mathbf{w}^\top \mathbf{x})^3\right] - 3\mathbf{w}\)
3. 与现有成分正交化（Gram-Schmidt 过程）
4. 归一化：\(\mathbf{w} = \mathbf{w}^+ / \|\mathbf{w}^+ \|\)
5. 重复直至收敛

基于负熵的 FastICA 最大化负熵的近似：

\[ J(\mathbf{w}) \propto \left(\mathbb{E}[G(\mathbf{w}^\top \mathbf{x})] - \mathbb{E}[G(\nu)]\right)^2 \]

其中 \(G\) 是非二次对比函数（例如 \(G(u) = \log\cosh(au)\)，\(1 \leq a \leq 2\)），\(\nu \sim \mathcal{N}(0,1)\)。

4.3 ICA 的应用

ICA 已被广泛应用于： - 盲源分离：分离混合音频信号（鸡尾酒会问题） - 脑电图/脑磁图神经成像：分离神经信号和伪迹信号 - 金融时间序列：发现驱动资产收益的独立因子 - 图像处理：去除伪迹和提取特征

---

5. 稀疏编码与字典学习

5.1 稀疏编码模型

稀疏编码假设每个数据点 \(\mathbf{x}_i\) 可以表示为来自过完备字典 \(\mathbf{D} \in \mathbb{R}^{D \times K}\)（其中 \(K > D\)）中少量基向量（原子）的线性组合：

\[ \mathbf{x} \approx \mathbf{D}\boldsymbol{\alpha} \]

其中系数向量 \(\boldsymbol{\alpha} \in \mathbb{R}^K\) 是稀疏的，意味着其大多数元素为零。稀疏编码目标结合了重构误差与系数的 \(\ell_1\) 惩罚：

\[ \min_{\mathbf{D}, \boldsymbol{\alpha}_i} \sum_{i=1}^N \|\mathbf{x}_i - \mathbf{D}\boldsymbol{\alpha}_i\|_2^2 + \lambda \|\boldsymbol{\alpha}_i\|_1 \]

约束条件为 \(\mathbf{D}\) 的所有列 \(\mathbf{d}_k\) 满足 \(\|\mathbf{d}_k\|_2 \leq 1\)。

5.2 字典学习算法

最优方向法（MOD）：在稀疏编码（使用 LARS 或 OMP）与字典更新之间交替：

\[ \mathbf{D} = \mathbf{X}\mathbf{A}^\top (\mathbf{A}\mathbf{A}^\top)^{-1} \]

K-SVD 算法：一次更新一个字典原子，以及使用该原子的系数，实现更精确的字典更新：

1. 稀疏编码步骤：对每个 \(\mathbf{x}_i\)，计算稀疏系数 \(\boldsymbol{\alpha}_i\)
2. 字典更新步骤：对每个原子 \(\mathbf{d}_k\)，使用残差的 SVD 进行更新

5.3 与神经机制的关联

稀疏编码为生物神经系统中的信号处理提供了一个令人信服的模型。视觉皮层中的 V1 简单细胞对边缘状特征产生响应，其响应表现出稀疏且具有选择性的活动模式——与稀疏表示假设一致。这一关联促使将稀疏编码原理纳入深度学习架构。

---

6. 慢特征分析（SFA）

6.1 时间缓变性原则

慢特征分析（Wiskott 和 Sejnowski，2002）利用了环境中的潜在因果因子随时间缓慢变化的假设，而高维感官观测可能快速变化。对于时间序列 \(\mathbf{x}(t)\)，SFA 提取变化缓慢的特征：

\[ \min_{f} \frac{\langle\dot{z}^2\rangle}{\langle z^2 \rangle} = \frac{\mathbb{E}[\dot{z}^2]}{\mathbb{E}[z^2] - \mathbb{E}[z]^2} \]

其中 \(z = f(\mathbf{x})\) 是提取的特征，\(\dot{z} = \frac{dz}{dt}\) 是其时间导数。

6.2 线性 SFA

对于线性特征 \(z = \mathbf{w}^\top \mathbf{x}\)，目标变为：

\[ \min_{\mathbf{w}} \frac{\mathbf{w}^\top \langle \ddot{\mathbf{x}}\ddot{\mathbf{x}}^\top \rangle \mathbf{w}}{\mathbf{w}^\top \langle \tilde{\mathbf{x}}\tilde{\mathbf{x}}^\top \rangle \mathbf{w}} \]

其中 \(\tilde{\mathbf{x}} = \mathbf{x} - \langle \mathbf{x} \rangle\) 是中心化的，\(\ddot{\mathbf{x}}\) 表示二阶时间导数（或者等价地，时间差分的协方差）。解涉及求解广义特征值问题：

\[ \mathbf{C}_{\dot{x}\dot{x}} \mathbf{w} = \lambda \mathbf{C}_{\tilde{x}\tilde{x}} \mathbf{w} \]

6.3 分层 SFA

深度 SFA 堆叠多个 SFA 层以提取越来越抽象、变化缓慢的特征。每一层执行： 1. 线性展开（例如，二次基展开） 2. SFA 寻找慢特征 3. 非线性压缩（例如，阈值处理）

这种分层方法已成功应用于从视频数据学习不变视觉表示。

---

7. 状态空间模型：卡尔曼滤波与隐马尔可夫模型

7.1 卡尔曼滤波

卡尔曼滤波为线性高斯状态空间模型提供最优递归推理：

\[ \mathbf{h}_t = \mathbf{A}\mathbf{h}_{t-1} + \mathbf{w}_t, \quad \mathbf{w}_t \sim \mathcal{N}(\mathbf{0}, \mathbf{Q}) \]

\[ \mathbf{x}_t = \mathbf{C}\mathbf{h}_t + \mathbf{v}_t, \quad \mathbf{v}_t \sim \mathcal{N}(\mathbf{0}, \mathbf{R}) \]

其中 \(\mathbf{h}_t\) 是隐藏状态，\(\mathbf{x}_t\) 是观测值，\(\mathbf{A}\) 是状态转移矩阵，\(\mathbf{C}\) 是观测矩阵。卡尔曼滤波在以下步骤之间交替：

预测步骤：

\[ \hat{\mathbf{h}}_{t|t-1} = \mathbf{A}\hat{\mathbf{h}}_{t-1|t-1} \]

\[ \mathbf{P}_{t|t-1} = \mathbf{A}\mathbf{P}_{t-1|t-1}\mathbf{A}^\top + \mathbf{Q} \]

更新步骤：

\[ \mathbf{K}_t = \mathbf{P}_{t|t-1}\mathbf{C}^\top (\mathbf{C}\mathbf{P}_{t|t-1}\mathbf{C}^\top + \mathbf{R})^{-1} \]

\[ \hat{\mathbf{h}}_{t|t} = \hat{\mathbf{h}}_{t|t-1} + \mathbf{K}_t(\mathbf{x}_t - \mathbf{C}\hat{\mathbf{h}}_{t|t-1}) \]

\[ \mathbf{P}_{t|t} = (\mathbf{I} - \mathbf{K}_t\mathbf{C})\mathbf{P}_{t|t-1} \]

7.2 隐马尔可夫模型（HMM）

隐马尔可夫模型将线性因子模型扩展到离散潜状态。HMM 由以下特征刻画： - 初始状态分布 \(\boldsymbol{\pi} = P(h_1)\) - 转移矩阵 \(\mathbf{A}\)，其中 \(A_{ij} = P(h_t = j | h_{t-1} = i)\) - 发射分布 \(p(\mathbf{x}_t | h_t = k)\)

HMM 的三个基本问题是： 1. 评估：使用前向算法计算 \(P(\mathbf{x} | \boldsymbol{\lambda})\) 2. 解码：使用 Viterbi 算法找到最可能的状态序列 3. 学习：通过 Baum-Welch 算法（一种特殊的 EM 算法）估计参数

7.3 线性因子模型的 EM 算法

期望最大化（EM）算法为训练具有潜变量的线性因子模型提供了统一框架。对于模型 \(p(\mathbf{x}, \mathbf{h} | \boldsymbol{\theta})\)，EM 在以下步骤之间交替：

E 步骤： 计算潜变量上的后验分布：

\[ Q(\mathbf{h}^{(i)}) = p\left(\mathbf{h}^{(i)} | \mathbf{x}^{(i)}, \boldsymbol{\theta}^{\text{old}}\right) \]

M 步骤： 最大化期望完整数据对数似然：

\[ \boldsymbol{\theta}^{\text{new}} = \arg\max_{\boldsymbol{\theta}} \sum_i \mathbb{E}_{\mathbf{h}^{(i)} \sim Q}\left[\log p\left(\mathbf{x}^{(i)}, \mathbf{h}^{(i)} | \boldsymbol{\theta}\right)\right] \]

对于 PPCA，EM 算法为直接最大似然估计提供了一种计算高效替代方案，特别是在 \(H \ll D\) 时：

M 步骤更新：

\[ \mathbf{W}^{\text{new}} = \left(\sum_i (\mathbf{x}_i - \boldsymbol{\mu})(\mathbb{E}[\mathbf{h}_i]^\top)\right) \left(\sum_i \mathbb{E}[\mathbf{h}_i\mathbf{h}_i^\top]\right)^{-1} \]

\[ \boldsymbol{\mu}^{\text{new}} = \bar{\mathbf{x}} - \mathbf{W}^{\text{new}}\mathbb{E}[\mathbf{h}] \]

\[ (\sigma^2)^{\text{new}} = \frac{1}{ND}\sum_i \left(\|\mathbf{x}_i - \boldsymbol{\mu}^{\text{new}}\|^2 - 2\mathbb{E}[\mathbf{h}_i]^\top (\mathbf{W}^{\text{new}})^\top (\mathbf{x}_i - \boldsymbol{\mu}^{\text{new}}) + \text{tr}(\mathbb{E}[\mathbf{h}_i\mathbf{h}_i^\top])\right) \]

其中 \(\mathbb{E}[\mathbf{h}_i]\) 和 \(\mathbb{E}[\mathbf{h}_i\mathbf{h}_i^\top]\) 使用第 2.3 节中的后验公式计算。

---

8. 作为深度生成模型的构建基础

线性因子模型在许多深度生成模型中作为基础组件：

8.1 深度信念网络（DBN）

深度信念网络由堆叠的受限玻尔兹曼机（RBM）组成，其中第 \(k\) 层的 RBM 隐藏层作为第 \(k-1\) 层可见层的潜因子。顶部两层形成无向 RBM，而较低层形成有向自顶向下连接。DBN 可以使用唤醒-睡眠算法或对比散度逐层贪婪训练。

8.2 变分自编码器（VAE）

VAE 将线性因子模型（具体为高斯潜变量模型）与深度神经网络结合用于编码器和解码器。编码器参数化 \(q_\phi(\mathbf{h}|\mathbf{x}) = \mathcal{N}(\boldsymbol{\mu}_\phi(\mathbf{x}), \boldsymbol{\sigma}_\phi^2(\mathbf{x})\mathbf{I})\)，而解码器定义 \(p_\theta(\mathbf{x}|\mathbf{h})\)。VAE 目标为：

\[ \mathcal{L} = \mathbb{E}_{q_\phi(\mathbf{h}|\mathbf{x})}[\log p_\theta(\mathbf{x}|\mathbf{h})] - D_{\text{KL}}(q_\phi(\mathbf{h}|\mathbf{x}) \| p(\mathbf{h})) \]

8.3 Helmholtz 机器与唤醒-睡眠算法

Helmholtz 机器通过使用任意神经网络参数化识别模型 \(q(\mathbf{h}|\mathbf{x})\) 和生成模型 \(p(\mathbf{x}|\mathbf{h})\) 来推广线性因子模型。唤醒-睡眠算法通过以下交替过程训练这些网络： - 唤醒阶段：固定识别网络，使用来自后验的样本更新生成网络 - 睡眠阶段：固定生成网络，使用来自生成先验的样本更新识别网络

---

9. 公式汇总表

| 模型 | 似然 \(p(\mathbf{x}|\mathbf{h})\) | 先验 \(p(\mathbf{h})\) | 关键参数 | |------|--------------------------------|---------------------|---------| | PPCA | \(\mathcal{N}(\mathbf{W}\mathbf{h} + \boldsymbol{\mu}, \sigma^2\mathbf{I})\) | \(\mathcal{N}(\mathbf{0}, \mathbf{I})\) | \(\mathbf{W} \in \mathbb{R}^{D \times H}, \sigma^2\) | | 因子分析 | \(\mathcal{N}(\mathbf{W}\mathbf{h} + \boldsymbol{\mu}, \boldsymbol{\Psi})\) | \(\mathcal{N}(\mathbf{0}, \mathbf{I})\) | \(\mathbf{W} \in \mathbb{R}^{D \times H}, \boldsymbol{\Psi} = \text{diag}(\psi_1, \ldots, \psi_D)\) | | ICA | \(\mathbf{x} = \mathbf{A}\mathbf{s}\)（确定性） | 边缘分布的乘积 \(p(s_1)\cdots p(s_D)\) | 混合矩阵 \(\mathbf{A}\) | | 稀疏编码 | \(\mathcal{N}(\mathbf{D}\boldsymbol{\alpha}, \sigma^2\mathbf{I})\) | Laplace 先验 \(p(\alpha_k) \propto \exp(-\lambda|\alpha_k|)\) | 字典 \(\mathbf{D} \in \mathbb{R}^{D \times K}\) | | SFA | \(p(\mathbf{x}|z) = \mathcal{N}(\mathbf{W}_z z, \boldsymbol{\Psi}_z)\) | 因实现而异 | 线性投影 \(\mathbf{w}\) | | 卡尔曼滤波 | \(\mathcal{N}(\mathbf{C}\mathbf{h}_t, \mathbf{R})\) | \(\mathcal{N}(\mathbf{A}\mathbf{h}_{t-1}, \mathbf{Q})\) | \(\{\mathbf{A}, \mathbf{C}, \mathbf{Q}, \mathbf{R}\}\) | | HMM | 类别或高斯发射 | 马尔可夫链 \(\{\boldsymbol{\pi}, \mathbf{A}\}\) | \(\{\boldsymbol{\pi}, \mathbf{A}, \mathbf{B}\}\) |

---

10. 结论

线性因子模型为理解高维数据中的潜变量依赖关系提供了丰富的数学框架。从经典的 PCA 到现代的深度生成模型，线性因子建模的原则——通过潜因子进行降维、通过 EM 进行概率推理、以及从噪声中分离信号——继续支撑着表示学习和生成建模的当代研究。

线性因子模型的数学优雅性，结合其可处理的推理过程和可解释的潜表示，确保了它们作为深度学习和统计机器学习中基础概念的持续相关性。

---

参考文献

• Bishop, C. M. (1999). Training with Noise is Equivalent to Tikhonov Regularization. Neural Computation.
• Hyvärinen, A., & Oja, E. (1997). Independent Component Analysis: Algorithms and Applications. Neural Networks.
• Tipping, M. E., & Bishop, C. M. (1999). Probabilistic Principal Component Analysis. Journal of the Royal Statistical Society.
• Wiskott, L., & Sejnowski, T. J. (2002). Slow Feature Analysis: A Theoretical Analysis of Optimal Free Responses. Neural Computation.
• Murphy, K. P. (2012). Machine Learning: A Probabilistic Perspective. MIT Press.
• Goodfellow, I., Bengio, Y., & Courville, A. (2016). Deep Learning. MIT Press.