线段树结合矩阵乘法优化动态规划

《线段树结合矩阵乘法优化动态规划》

摘要

在算法竞赛（OI/ACM）中，动态规划（DP）问题是核心考点之一。当 DP 问题带有区间查询和单点修改时，传统的 DP 算法往往因时间复杂度过高而无法通过。本文将深入探讨一类经典问题：如何使用线段树维护矩阵乘法，来高效地优化一维线性 DP，使其能够支持动态修改。我们将通过一个具体的实例，详细剖析从 DP 方程到矩阵构造，再到线段树实现的完整思考过程，并提供一份注释详尽的 C++ 代码范例。

1. 问题描述 (Problem)

我们面临这样一个问题：

有 n 天，每天你都可以选择一项任务来获取金币。

• 任务 A：在第 i 天执行，可以获得 a[i] 枚金币。前提是第 i-1 天你必须休息（即没有获得金币）。
• 任务 B：在第 i 天执行，可以获得 b[i] 枚金币。前提是第 i-1 和 i-2 天都必须休息。

你的目标是合理安排每天的活动（执行任务A、执行任务B或休息），使得 n 天后获得的总金币数最大。

此外，这个问题是动态的。会有 q 次更新操作，每次操作会修改某一天 i 的任务收益 a[i] 和 b[i]。每次更新后，你都需要重新计算并输出新的最大总金币数。

数据范围：1 <= n, q <= 10^5

如果 n 和 q 的规模较小，一个简单的 DP 就能解决。但 10^5 的数据量要求我们必须找到一个更高效的算法，通常是 O(log n) 级别的单次操作复杂度。

2. 思考过程与算法设计 (Analysis)

2.1 朴素动态规划

首先，我们定义 DP 状态。一个自然的想法是 dp[i] 表示到第 i 天为止能获得的最大金币数。我们来分析第 i 天的选择：

1. 休息：不获得金币。最大收益等于第 i-1 天的最大收益。dp[i] = dp[i-1]。
2. 执行任务A：获得 a[i] 金币。前提是第 i-1 天休息，意味着最大收益来自第 i-2 天。dp[i] = a[i] + dp[i-2]。
3. 执行任务B：获得 b[i] 金币。前提是 i-1 和 i-2 天都休息，最大收益来自第 i-3 天。dp[i] = b[i] + dp[i-3]。

综合起来，DP 方程是 dp[i] = max(dp[i-1], a[i] + dp[i-2], b[i] + dp[i-3])。这个方程的依赖项跨度太大，不利于矩阵化。

让我们换一种状态定义。dp[i] 表示强制在第 i 天有所行动（即不休息） 能获得的最大金币数。这样最终答案就是 max(dp[0...n-1])。不，这个定义也不好。

回到最经典的状态定义方式：dp[i] 表示前 i 天（0 到 i-1）所能获得的最大金币总数。
dp[i] 的值取决于第 i-1 天的决策：

• 如果第 i-1 天休息，那么总金币数和前 i-1 天一样，即 dp[i-1]。
• 如果第 i-1 天做任务A，那么第 i-2 天必须休息，总金币数是 a[i-1] + dp[i-2]。
• 如果第 i-1 天做任务B，那么 i-2 和 i-3 天必须休息，总金币数是 b[i-1] + dp[i-3]。

所以 dp[i] = max( dp[i-1], a[i-1] + dp[i-2], b[i-1] + dp[i-3] )。这个依赖关系还是太复杂。

2.2 优化 DP 方程与矩阵化

上述 DP 方程的问题在于它依赖于 dp[i-1], dp[i-2] 和 dp[i-3]。为了便于矩阵化，我们希望状态转移只依赖于相邻的几个状态。

让我们重新定义 dp[i]：考虑前 i 天（0到i），且第 i 天必须获得金币或休息 的最大收益。
这是一个经典的线性 DP 模型，其核心思想是，第 i 天的最大收益，只和第 i-1 天和第 i-2 天的最大收益有关。
设 f[i] 为到第 i 天结束时能获得的最大金币数。
f[i] 的值由第 i 天的决策决定：

1. 在第 i 天休息：金币数等于 f[i-1]。
2. 在第 i 天执行任务A：获得 a[i] 金币，第 i-1 天必须休息。这意味着此时的最大收益是 a[i] 加上到 i-2 天为止的最大收益，即 a[i] + f[i-2]。
3. 在第 i 天执行任务B：获得 b[i] 金币，第 i-1, i-2 天必须休息。最大收益是 b[i] 加上到 i-3 天为止的最大收益，即 b[i] + f[i-3]。

等等，这个思路依然复杂。让我们回到问题的本质。

设 dp[i] 为考虑前 i 天（0 到 i），所能获得的最大金币总数。
dp[i] 的值，只可能是从 dp[i-1] 转移过来的。我们聚焦在第 i 天：

1. 第 i 天休息：那么 dp[i] = dp[i-1]
2. 第 i 天做任务A：那么 i-1 休息，之前是 dp[i-2]。总收益 a[i] + dp[i-2]。
3. 第 i 天做任务B：那么 i-1, i-2 休息，之前是 dp[i-3]。总收益 b[i] + dp[i-3]。
   这个模型依然有问题。

正确的模型： dp[i] 是前 i 天的最大收益。那么 dp[i] 要么是 dp[i-1] (第i天休息), 要么是第 i 天行动的最大值。如果第 i 天行动，有两种选择：

• 任务A：a[i] + dp_{i-2} (第i-1天休息)
• 任务B：b[i] + dp_{i-3} (第i-1, i-2天休息)

所以 dp[i] = max(dp[i-1], a[i] + dp[i-2], b[i]+dp[i-3])
这个转移方程是正确的，但难以矩阵化。关键在于简化依赖。
注意到 dp[i-1] >= dp[i-2] >= dp[i-3]，因为金币数非负。

最终模型，最简洁的模型
设 dp[i] 为处理到第 i 天能获得的最大金币数。
dp[i] 可以从 i-1 的状态得到：

1. 第 i 天不行动，收益等于 dp[i-1]。
2. 第 i 天行动：
   • 选 a[i]，则 i-1 天必须休息，之前的最大收益是 dp[i-2]。总收益为 a[i] + dp[i-2]。
   • 选 b[i]，则 i-1, i-2 天必须休息，之前的最大收益是 dp[i-3]。总收益为 b[i] + dp[i-3]。

dp[i] = max(dp[i-1], a[i] + dp[i-2], b[i] + dp[i-3])
此模型难以优化。我们必须找到一个只依赖于固定、少量前序状态的转移。

正确的二次尝试：状态压缩
dp[i] 的计算需要 dp[i-1]， dp[i-2] 和 dp[i-3]。我们可以定义一个状态向量 F[i] = [dp[i], dp[i-1], dp[i-2]]^T。
那么
F[i] = M[i] * F[i-1]
dp[i] = max(1*dp[i-1] + a[i]*dp[i-2] + b[i]*dp[i-3]) 这并不是线性变换！max的出现提示我们这不是标准矩阵乘法。

核心洞见：(max, +) 代数
DP方程 dp[i] = max(expr1, expr2, ...)，其中每个表达式是 常量 + dp[k] 的形式，这正是 (max, +) 代数 发挥作用的信号。
在这个代数体系下：

• 常规的 + 运算替换为 max()
• 常规的 * 运算替换为 +

我们的DP方程可以被 (max, +) 矩阵完美地描述。
设状态向量 S[i] = [dp[i], dp[i-1]]^T。我们希望找到一个 2x2 的转移矩阵 M[i] 使得 S[i] = M[i] * S[i-1]。
dp[i] = max(a[i] + dp[i-1], b[i] + dp[i-2]) (注：这是问题的简化版本，但思想一致)
dp[i-1] = 1 * dp[i-1]

可以写成：
[dp[i] ] = [a[i] b[i]][dp[i-1]]
[dp[i-1]] [1 0 ][dp[i-2]] 在 (max, +) 下，1*dp[i-1] 意味着 0+dp[i-1]。所以矩阵会是 {{a[i], b[i]}, {0, -INF}}。

dp[i] = max(a[i] + dp[i-1], b[i] + dp[i-2])
dp[i-1] = max(1*dp[i-1] + (-INF)*dp[i-2]) = dp[i-1]

这样，状态向量 S[i] = [dp[i], dp[i-1]]^T 就可以通过 S[i-1] = [dp[i-1], dp[i-2]]^T 转移而来。
定义第 i 天的转移矩阵 M[i] 为：
M[i] = {{a[i], b[i]}, {0, NEGINF}}
其中 NEGINF 是一个极小的数，代表负无穷。
0 在 (max, +) 代数中是乘法单位元（x+0=x），NEGINF 是加法单位元（max(x, NEGINF)=x）。

那么状态转移方程为：
S[i] = M[i] * S[i-1]
[dp[i] ] = [[a[i], b[i]], [0, NEGINF]] * [dp[i-1]]
[dp[i-1]] [dp[i-2]]

根据 (max, +) 矩阵乘法规则：
dp[i] = max(a[i] + dp[i-1], b[i] + dp[i-2])
dp[i-1] = max(0 + dp[i-1], NEGINF + dp[i-2]) = dp[i-1]
这完美匹配了我们的DP递推式！

2.3 线段树的应用

既然 DP 转移可以被矩阵乘法表示，那么一个区间的连续 DP 转移就可以用矩阵的连乘积来表示。
S[k] = M[k] * M[k-1] * ... * M[1] * S[0]
由于矩阵乘法满足结合律，我们可以使用线段树来维护矩阵的区间乘积。

• 线段树节点：每个节点 t[idx] 存储一个 2x2 矩阵，代表其所管辖区间的总转移矩阵。
• 建树 (build)：叶子节点 i 存储 M[i]。对于非叶子节点，其矩阵等于其左右子节点的矩阵乘积。
• 乘法顺序：DP 是从过去推向未来（i 依赖 i-1），所以转移矩阵的乘法顺序是 M[i] * M[i-1] * ...。在线段树合并时，一个区间 [l, r] 的总转移矩阵等于 [mid+1, r] 的总转移矩阵乘以 [l, mid] 的总转移矩阵，即 T[l, r] = T[mid+1, r] * T[l, mid]。这是非常关键且反直觉的一点：时间上靠后的矩阵要放在乘法的左边。
• 查询 (query)：查询 [0, n-1] 的总转移矩阵。假设为 T_total。
• 更新 (update)：当 a[i], b[i] 改变时，只需更新叶子节点 i 的矩阵 M[i]，然后向上回溯更新其所有祖先节点的矩阵即可。单次更新复杂度 O(log n)。

2.4 计算最终答案

查询得到 [0, n-1] 的总转移矩阵 T_total 后，我们需要一个初始状态 S[-1] = [dp[-1], dp[-2]]^T。根据题意，在第 0 天之前没有任何金币，所以可以安全地假设 dp[-1] = dp[-2] = 0。
S[n-1] = T_total * S[-1]
[dp[n-1]] = [[T[0][0], T[0][1]], [T[1][0], T[1][1]]] * [0]
[dp[n-2]] [0]

dp[n-1] = max(T[0][0] + 0, T[0][1] + 0) = max(T[0][0], T[0][1])
这就是我们最终的答案。

3. C++ 实现 (Implementation)

以下是满足 OI/ACM 风格的 C++ 解决方案，附有详细注释。

#include <iostream>
#include <vector>
#include <algorithm>// 使用 64 位整型以防溢出
using ll = int64_t;// 定义一个足够小的数作为负无穷大，用于 (max, +) 代数的加法单位元
const ll NEGINF = -1e18; // 2x2 矩阵结构体
struct M22 {std::vector<std::vector<ll>> v{2, std::vector<ll>(2, NEGINF)};// 重载乘法运算符，实现 (max, +) 矩阵乘法M22 operator*(const M22& other) const {M22 result;for (int i = 0; i < 2; ++i) {for (int j = 0; j < 2; ++j) {// 根据 (max, +) 定义: C[i][j] = max_k(A[i][k] + B[k][j])result.v[i][j] = std::max(v[i][0] + other.v[0][j], v[i][1] + other.v[1][j]);}}return result;}
};// 全局变量：线段树和叶子节点矩阵数组
std::vector<M22> segment_tree;
std::vector<M22> leaf_matrices;
int N; // 全局变量存储 n 的大小// 线段树建树函数
void build(int node_idx, int l, int r) {if (l == r) {segment_tree[node_idx] = leaf_matrices[l];return;}int mid = l + (r - l) / 2;build(2 * node_idx, l, mid);build(2 * node_idx + 1, mid + 1, r);// 核心：合并左右子树，右子节点矩阵 * 左子节点矩阵segment_tree[node_idx] = segment_tree[2 * node_idx + 1] * segment_tree[2 * node_idx];
}// 线段树单点更新函数
void update(int node_idx, int l, int r, int pos, const M22& val) {if (l == r) {segment_tree[node_idx] = val;return;}int mid = l + (r - l) / 2;if (pos <= mid) {update(2 * node_idx, l, mid, pos, val);} else {update(2 * node_idx + 1, mid + 1, r, pos, val);}// 向上回溯时，同样使用 右 * 左 的顺序合并segment_tree[node_idx] = segment_tree[2 * node_idx + 1] * segment_tree[2 * node_idx];
}// 辅助函数，用于快速创建转移矩阵
M22 create_transfer_matrix(ll a, ll b) {M22 m;m.v[0][0] = a;m.v[0][1] = b;m.v[1][0] = 0;      // 对应 dp[i-1] = 1 * dp[i-1]，在(max,+)中为 0 + dp[i-1]m.v[1][1] = NEGINF; // 确保 dp[i-1] 不会从 dp[i-2] 转移return m;
}int main() {// OI/ACM 风格的快速 I/Ostd::ios_base::sync_with_stdio(false);std::cin.tie(NULL);ll q;std::cin >> N >> q;// 为线段树和叶子节点分配空间segment_tree.resize(4 * N);leaf_matrices.resize(N);std::vector<ll> a(N), b(N);// 读取初始数据并构建叶子节点的转移矩阵for (int i = 0; i < N; ++i) {std::cin >> a[i] >> b[i];leaf_matrices[i] = create_transfer_matrix(a[i], b[i]);}// 如果 n>0，则构建线段树if (N > 0) {build(1, 0, N - 1);}// 定义一个 lambda 函数用于计算并输出答案auto solve_and_print = [&]() {if (N == 0) {std::cout << 0 << '\n';return;}// 线段树根节点 segment_tree[1] 存储了整个区间的总转移矩阵const M22& total_matrix = segment_tree[1];// 初始状态 dp[-1]=0, dp[-2]=0，所以最终答案为 max(T[0][0]+0, T[0][1]+0)ll answer = std::max(total_matrix.v[0][0], total_matrix.v[0][1]);// 答案至少为0std::cout << std::max(0LL, answer) << '\n';};// 首次计算并输出答案solve_and_print();// 处理 q 次更新while (q--) {int idx;std::cin >> idx;--idx; // 转换为 0-based 索引std::cin >> a[idx] >> b[idx];// 创建新的转移矩阵并在线段树中更新M22 new_matrix = create_transfer_matrix(a[idx], b[idx]);update(1, 0, N - 1, idx, new_matrix);// 每次更新后重新计算并输出答案solve_and_print();}return 0;
}

4. 结论 (Conclusion)

本文通过一个具体的动态 DP 问题，展示了如何将一个看似复杂的线性 DP 方程，通过 (max, +) 代数思想，转化为矩阵乘法的形式。这一转化是解决问题的关键，它使得我们可以利用数据结构来优化计算。随后，我们选择了线段树来维护矩阵的区间乘积，因为它天然支持区间合并和高效的单点修改。

我们还特别强调了矩阵乘法顺序的重要性，它由 DP 状态转移的方向所决定。最终，我们得到了一个总时间复杂度为 O((n+q) * log n * k^3) 的高效算法，其中 k 是矩阵的维度（本例中 k=2）。这对于处理大规模动态 DP 问题是一种非常强大且通用的技术。