演算法學習筆記(3)：與斜率最佳化共舞

敲鍵盤的大鯊魚發表於2023-02-02
演算法筆記

斜率最佳化 DP

對於像如下這樣的 dp 方程，我們可以使用斜率最佳化解決。

\[dp_i = \min / \max \{ a_i \times b_j + c_i + d_j \} \]

例題：P2900 Land Acquisition G

顯然，如果 \(w_i \ge w_j\) 且 \(l_i \ge l_j\)，那麼土地 \(j\) 可以直接被土地 \(i\) 併購。

考慮將所有土地按 \(w_i\) 降序排序，\(w_i\) 相同的按 \(l_i\) 降序排序，再透過雙指標保留不能被其它土地併購的土地。

注意到此時被保留下的土地，\(w_i\) 滿足 單調遞減，\(l_i\) 滿足 單調遞增。

令 \(dp_i\) 表示將第 \(1 \sim i\) 塊土地併購的最小代價，不難列出 dp 方程：

\[dp_i = \min_{0 \le j < i} \{ dp_j + w_{j + 1} \times l_{i} \} \]

上面的方程複雜度為 \(O(n^2)\)，然而 \(n \leq 5 \times 10^4\)，我們需要想辦法最佳化。

先將 \(\min\) 去掉，再把 \(w_{j + 1} \times l_{i}\) 移到左邊，得：

\[-w_{j + 1} \times l_i + dp_i = dp_j \]

將所有能夠轉移到 \(i\) 的 \(j\) 視為一個點 \((-w_{j + 1}, dp_j)\)，那麼問題就轉化成了：

  • 有一條斜率為 \(l_i\) 的直線，這條直線需要經過上述的這些點中的一個或多個，並且希望它的截距最小。

上面這句話有點抽象，不妨將 原方程 與 一次函式 做一個對比，發現它們確實十分相似：

\[\textcolor{red}{-w_{j + 1}} \times \textcolor{blue}{l_i} + \textcolor{green}{dp_i} = \textcolor{orange}{dp_j} \]

\[\newline \textcolor{red}{k}\textcolor{blue}{x} + \textcolor{green}{b} = \textcolor{orange}{y} \]

本質就是我們拿著一條斜率為 \(l_i\) 的直線，從下往上靠（增加它的截距），直到某一時刻這條線經過了我們維護的一個或多個點，停下來，此時截距（\(dp_i\)）一定是最小的。

仔細觀察可以發現，無論斜率 (\(l_i\)) 是多少，有些點一定不會被第一次經過（如下圖灰色點）。將這些點去除，剩餘的點恰好形成了一個右下凸殼：

不妨先來維護這個凸殼，假設現在加入點 \(C(-w_{i + 1}, dp_i)\)，考慮凸殼上最新的兩個點 \(A\) 和 \(B\)，只可能有以下兩種情況：

  • \(C\) 點可以直接加入凸殼。

記點 \(X\) 與點 \(Y\) 所連成的直線斜率為 \(\text{slope}(X, Y)\)，則 \(C\) 能直接加入凸殼當且僅當 \(\text{slope}(A, B) \leq \text{slope}(B, C)\)。

  • \(C\) 點加入後無法形成凸殼。

這時候一定有 \(\text{slope}(A, B) \geq \text{slope}(B, C)\)，故需要將 \(B\) 彈出凸殼，在拿剩餘凸殼最新的兩個點與 \(C\) 作比較。

由於這些點的 \(x\) 單調遞增，可以使用類似單調棧的方法維護，時間複雜度 \(O(n)\)。

接下來就是要在維護出的凸殼上找到最優決策點，隨便畫一張圖，由於這些直線的斜率一定不斷增加，所以最優決策點一定不斷向凸殼上方移動：

進一步觀察，可以發現，如果凸殼上相鄰的兩點 \(A\) 和 \(B\) 滿足 \(\text{slope}(A, B) \leq l_i\)，那麼點 \(B\) 永遠在 \(i\) 以後一定 不會成為最優決策點。不妨把原來維護凸殼的單調棧換成單調佇列，透過凸殼最前面兩點的斜率判斷是否彈出隊頭。最終，隊頭一定是當前的最優決策點。

時間複雜度 \(O(n)\)（本題由於要排序所以總複雜度應該為 \(O(n \log n)\)，這裡只計算了 dp 的複雜度）。程式碼如下：

#include <bits/stdc++.h>
using namespace std;
int n, pos, q[50005], h, t;
long long dp[50005];
struct land
{
  int w, l;
  bool operator<(const land t) const
  {
    return w == t.w ? l > t.l : w > t.w;
  }
} p[50005];

double X(int j)
{
  return -p[j + 1].w;
}

double Y(int j)
{
  return dp[j];
}

double slope(int x, int y)
{
  return (Y(x) - Y(y)) / (X(x) - X(y));
}

int main()
{
  scanf("%d", &n);
  for (int i = 1; i <= n; i++)
    scanf("%d%d", &p[i].w, &p[i].l);
  sort(p + 1, p + n + 1);
  for (int i = 1; i <= n; i++)
    if (p[i].l > p[pos].l)
      p[++pos] = p[i];
  for (int i = 1; i <= pos; i++)
  {
    while (h < t && slope(q[h], q[h + 1]) <= p[i].l)
      h++;
    dp[i] = dp[q[h]] + p[q[h] + 1].w * 1LL * p[i].l;
    while (h < t && slope(q[t - 1], q[t]) >= slope(q[t], i))
      t--;
    q[++t] = i;
  }
  printf("%lld", dp[pos]);
  return 0;
}

例題：P3195 玩具裝箱

記 \(s_i = \sum_{j=1}^i C_i\)，則可以列出 dp 方程：

\[dp_i = \min_{0 \le j < i} \{ dp_j + (s_i - s_j + i - j - 1 - L)^2 \} \]

令 \(a_i = s_i + i - 1 - L, b_i = s_j + j\)，原方程變為

\[dp_i = \min_{0 \le j < i} \{ dp_j + (a_i - b_j)^2 \} \]

即

\[dp_i = \min_{0 \le j < i} \{ dp_j - 2a_ib_j + b_j^2 \} + a_i^2 \]

將 \(\min\) 及其以外得東西去掉，得

\[dp_i = dp_j - 2a_ib_j + b_j^2 \]

化為直線的形式

\[2a_ib_j + dp_i = dp_j + b_j^2 \]

類似與之前的形式，把每個 \(j\) 視作點 \((b_j, dp_j + b_j^2)\)，每次直線的斜率為 \(2a_i\)，做斜率最佳化即可。

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
int n, L, C, h, t, q[50005];
ll s[50005], dp[50005];

double X(int j)
{
  return s[j] + j;
}

double Y(int j)
{
  return (s[j] + j) * (s[j] + j) + dp[j];
}

double slope(int i, int j)
{
  return (Y(i) - Y(j)) / (X(i) - X(j));
}

int main()
{
  scanf("%d%d", &n, &L);
  for (int i = 1; i <= n; i++)
  {
    scanf("%d", &C);
    s[i] = s[i - 1] + C;
  }
  h = t = 1;
  for (int i = 1; i <= n; i++)
  {
    while (h < t && slope(q[h], q[h + 1]) <= 2 * (s[i] + i - L - 1))
      h++;
    dp[i] = dp[q[h]] + sq(s[i] - s[q[h]] + i - q[h] - 1 - L);
    while (h < t && slope(q[t - 1], q[t]) >= slope(q[t], i))
      t--;
    q[++t] = i;
  }
  printf("%lld", dp[n]);
  return 0;
}

例題：P5785 任務安排（弱化版 P2365 任務安排）

至於方程是如何提前計算貢獻的我就不細說了，記 \(\text{sumT}_i = \sum_{j=1}^i T_i,\text{sumC}_i = \sum_{j=1}^i C_i\)，可列出方程：

\[dp_i = \min \{ dp_j + \text{sumT}_i \times (\text{sumC}_i - \text{sumC}_j) + s \times (\text{sumC}_n - \text{sumC}_j) \} \]

展開整理成直線形式，即

\[(\text{sumT}_i + s) \times \text{sumC}_j + dp_i = dp_j \]

每個 \(j\) 所代表的點：\((\text{sumC}_i, dp_j)\)，直線斜率 \(\text{sumT}_i + s\)。

正準備敲板子的你突然發現了這一條限制：

\[|T_i| \leq 2^8 \]

也就是說

\[-2^8 \le T_i \le 2^8 \]

那麼 \(\text{sumT}_i + s\) 就不單調了，就不能一味地彈出隊頭了！

可凸殼上相鄰兩點地斜率還是單調的啊！那麼只需要在凸殼上二分，找決策點即可。

時間複雜度 \(O(n \log n)\)。

#include <bits/stdc++.h>
using namespace std;
long long n, s, sumT[300005], sumC[300005], dp[300005], top, stk[300005];

long long X(int j)
{
  return sumC[j];
}

long long long Y(int j)
{
  return dp[j];
}

int main()
{
  scanf("%lld%lld", &n, &s);
  for (int i = 1; i <= n; i++)
  {
    long long t, f;
    scanf("%lld%lld", &t, &f);
    sumT[i] = sumT[i - 1] + t;
    sumC[i] = sumC[i - 1] + f;
  }
  stk[++top] = 0;
  for (int i = 1; i <= n; i++)
  {
    int l = 1, r = top - 1, pos = stk[top];
    while (l <= r)
    {
      int mid = (l + r) >> 1;
      if (Y(stk[mid + 1]) - Y(stk[mid]) > (sumT[i] + s) * (X(stk[mid + 1]) - X(stk[mid])))
        pos = stk[mid], r = mid - 1;
      else
        l = mid + 1;
    }
    dp[i] = dp[pos] + sumT[i] * (sumC[i] - sumC[pos]) + s * (sumC[n] - sumC[pos]);
    while (top > 1 && (Y(stk[top]) - Y(stk[top - 1])) * (X(i) - X(stk[top])) >= (Y(i) - Y(stk[top])) * (X(stk[top]) - X(stk[top - 1])))
      top--;
    stk[++top] = i;
  }
  printf("%lld", dp[n]);
  return 0;
}

注：斜率最佳化比較時建議移項，把除法化為乘法；如果使用 double 型別的 slope 很可能會產生精度誤差（本題就卡了。

斜率最佳化小結：

  • 列出方程，先透過一些簡單的代換化簡成直線形式。

  • 透過 \(\min / \max\) 以及 \(X\) 和 \(Y\) 座標的單調性判斷凸殼方向。

  • 如果直線斜率不單調，使用二分維護。

  • 如果 \(Y\) 座標不單調，可以證明對最終凸殼沒有影響。

  • 如果 \(X\) 座標不單調，可以考慮 \(\text{CDQ}\) 分治解決，時間複雜度 \(O(n \log^2 n)\)。

相關文章