拉格朗日乘数法
昨天的模拟测考到了这样的一道题目：
给定一个 \(n\times n\) 的矩阵 \(w\)，构造一个有 \(n\) 个数的非负实数序列 \(x\)，需要保证其中所有数的和为 \(1\) ，求 \(\sum\limits^n_{i=1}\sum\limits^n_{j=i+1}w_{i, j}x_ix_j\) 的最大值。\(n\leq 10\)，保证 \(w_{i, i} = 0, w_{i, j} = w_{j, i}\)。
当时我非常兴奋，因为之前也考过一道非常相似的题目，只是 \(w\) 一定是 01 矩阵，且 \(n\) 的范围扩大到了 \(20\)。
之前的题目解法
不难发现这道题可以转化为，给定一个有 \(n\) 个点的无向图 \(G\)，连边给定，给每个点赋一个点权，且边权为所连接的两点点权之和，需要求出这个图的边权之和最大值 \(f(G)\)。
注意到在一个完全图的时候能够尽量应用到每个 \(w_{i, j}\)，所以 \(f(G)\) 一定有一个下界：设 \(G\) 的最大团大小为 \(k\)，那么可以令团内所有点权为 \(\frac{1}{k}\)，有\(f(G) \geq \frac{k-1}{2k}\)。
尝试找到比这更优的做法，发现这很难做到。有一个结论：\(f(G) = \frac{k-1}{2k}\)。
Proof
考虑归纳证明，当 \(n=1\) 时显然成立。
假设当图大小为当前 \(n-1\) 时成立。如果最优解中序列 \(x\) 存在一项为 \(0\)，那么去掉就回到了 \(n-1\)。若 \(G\) 是完全图，那么结论显然成立。接下来只需证明当 \(G\) 不是完全图且 \(x\) 为正实数序列时成立即可。
设得出的最优解为 \(x_1, \cdots, x_n\)，与 \(x\) 有连边的点集为 \(e_x\)。不妨设 \(1\notin e_2\)。由于 \(x\) 为最优解，因此 \(\sum\limits_{i\in{e_1}} x_i = \sum\limits_{j\in{e_2}} x_j\)。于是对 \(x_1\) 和 \(x_2\) 进行一定的偏移：新序列 \(x_1+\Delta, x_2 - \Delta, x_3, \cdots, x_n\) 的答案为最优解答案减去 \(\Delta (\sum\limits_{i\in{e_1}} x_i - \sum\limits_{j\in{e_2}} x_j))\)，并无变化。于是令 \(\Delta = -x_1\)，就可以把第一项消掉，仍是最优解。
因此只需要得出图的最大团大小了。
但是这道题的 \(w_{i, j}\) 值域改变，由于不好计算甚至难以定义 “带权最大团”，上道题的确定性解法无法再次应用。当然，两题依然可以用算法优化的模拟退火。但随机化不是好做法，我们寻求一个具有确定性的、可以稳定得出真正解的做法。
我们可以使用 “拉格朗日乘数法”。它适用于求多元函数 \(f(x_1, \cdots, x_n)\) 在另一个多元函数 \(g(x_1, \cdots, x_n) = 0\) 时的极值。只需要引入新变量 \(\lambda\)，求拉格朗日函数 \(F(x_1, \cdots, x_n, \lambda) = f(x_1, \cdots, x_n) + \lambda g(x_1, \cdots, x_n)\) 的极值即可。这点可以通过对每一元求偏导后结果定为 \(0\) 做到。
具体的，这里有一道数学题用以演示：已知 \(a^2+b^2=1\)，求 \((a+1)(b+2)\) 的最大值。
直接计算不好做。但发现这两个式子可以由两个多元函数表示：\(f(a, b) = (a+1)(b+2), \phi(a, b) = a^2+b^2-1 = 0\)。有 \(F(a, b, \lambda) = f(a, b) + \lambda\phi(a, b) = (a+1)(b+2) + \lambda(a^2+b^2-1)\)。
分别对 \(a\)、 \(b\) 与 \(\lambda\) 求偏导：\(F_a^{'}=2a\lambda+b+2=0, F_b^{'} = 2b\lambda +a+1=0, F_\lambda^{'}=a^2+b^2-1=0\)。解三元方程即可。
于是这道题可以用此方法做：\(f(x_1, \cdots, x_n)=\sum\limits^n_{i=1}\sum\limits^n_{j=i+1}w_{i, j}x_ix_j, \phi(x_1, \cdots, x_n) = (\sum x_i) - 1=0\)。\(F(x_1, \cdots, x_n, \lambda) = f(x_1, \cdots, x_n) + \lambda\phi(x_1, \cdots, x_n)\)。
分别对 \(x_i\) 与 \(\lambda\) 求偏导：\(F_{x_i}^{'}=\sum\limits_{j=1}^n w_{i, j}x_j + \lambda = 0, F_\lambda^{'}=(\sum x_i) - 1=0\)，\(n\) 元方程组暴力高斯消元即可。
但是有个 bug：直接做最优解可能会使 \(x_i<0\)，而这样是不合法的。我们只需要枚举哪些 \(x_i\) 是 \(0\)，剩下的都是正数，找出合法的去最大值即可。时间复杂度 \(O(n^32^n)\)。
Code
// STOOOOOOOOOOOOOOOOOOOOOOOOO hzt CCCCCCCCCCCCCCCCCCCCCCCORZ
#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <vector>
using namespace std;
using LL = long long;
using PII = pair<int, int>;
constexpr int kN = 11 + 1;
constexpr double kE = 1e-9;
int T, n;
int w[kN][kN];
double a[kN][kN], b[kN];
int p[kN], c;
bool G() {
for (int i = 1, o; i <= c; i++) {
o = i;
for (int j = i + 1; j <= c; j++) {
if (fabs(a[o][i]) < fabs(a[j][i])) {
o = j;
}
}
if (fabs(a[o][i]) < kE) {
return 0;
}
swap(a[i], a[o]), swap(b[i], b[o]);
for (int j = 1; j <= c; j++) {
if (i == j) {
continue;
}
double d = a[j][i] / a[i][i];
for (int k = 1; k <= c; k++) {
a[j][k] -= a[i][k] * d;
}
b[j] -= d * b[i];
}
}
for (int i = 1; i < c; i++) {
b[i] /= a[i][i];
if (b[i] < -kE) {
return 0;
}
}
return 1;
};
int main() {
cin.tie(0)->sync_with_stdio(0);
cin >> T;
cout << fixed << setprecision(7);
while (T--) {
cin >> n;
for (int i = 1; i <= n; i++) {
for (int j = 1; j <= n; j++) {
cin >> w[i][j];
}
}
double ans = 0;
for (int s = 1; s < 1 << n; s++) {
c = 0;
for (int i = 0; i < n; i++) {
if (s >> i & 1) {
p[++c] = i + 1;
}
}
for (int i = 1; i <= c; i++) {
for (int j = 1; j <= c; j++) {
a[i][j] = w[p[i]][p[j]];
}
a[i][c + 1] = 1, b[i] = 0;
}
c++, a[c][c] = 0, b[c] = 1;
for (int i = 1; i < c; i++) {
a[c][i] = 1;
}
if (!G()) {
continue;
}
double ret = 0;
for (int i = 1; i < c - 1; i++) {
for (int j = i + 1; j < c; j++) {
ret += b[i] * b[j] * w[p[i]][p[j]];
}
}
ans = max(ans, ret);
}
cout << ans << '\n';
}
return 0;
}
\(\\\)
练习题
T1
形式化题意：有 \(n\) 个路段，每个路段有三个实数描述：\(s_i, k_i, v_i^{'}\)，又给定一个实数 \(E_U\) 表示一个人的体力值。构造一个长度为 \(n\) 的序列 \(v\)，经过路段 \(i\) 所需体力值 \(e_i\) 为 \(k_i(v_i-v_i^{'})^2s_i\)，让序列满足 \((\sum\limits_{i=1}^n e_i) - E_U = 0\)，且 \(\sum\limits_{i=1}^n \frac{s_i}{v_i}\) 的值最小。\(n\leq 10^4\)。
Solution
由题有 \(f(v_1, \cdots, v_n)=\sum\limits_{i=1}^n \frac{s_i}{v_i}, \phi(v_1, \cdots, v_n) = \sum\limits_{i=1}^n k_i(v_i-v_i^{'})^2s_i - E_U = 0\)。\(F(v_1, \cdots, v_n, \lambda)=F(v_1, \cdots, v_n) + \lambda\phi(v_1, \cdots, v_n)\)。
分别对 \(v_i, \lambda\) 求偏导。\(f_{v_i}^{'} = -\frac{s_i}{v_i^2}+2\lambda k_i(v_i-v_i^{'})s_i=0, f_\lambda^{'}=\sum\limits_{i=1}^n k_i(v_i - v_i^{'})^2s_i - E_U = 0\)。但这题 \(n\) 是一万级别的了，无法再直接高斯消元。考虑分析 \(f_{v_i}^{'}\) 所得式子的性质。
至此，我们可以看出一些性质：首先显然 \(v_i^2(v_i-v_i^{'})\) 在 \(\lambda\) 一定时随 \(v_i\) 单调递增。其次，当 \(\lambda\) 上升时，\(\frac{1}{2\lambda k_i}\) 下降，\(v_i\) 下降，又有 \(e_i\) 上升，所以 \(e_i\) 在 \(v_i\) 一定时随 \(\lambda\) 单调递增。
有了这两个单调性，可以在内层对于当前给定 \(\lambda\) 求出 \(\sum\limits_{i=1}^n e_i\) 最小值，再在外层二分 \(\lambda\) 的值。时间复杂度 \(O(n\log^2 V)\)。
\(\\\)
Code
// STOOOOOOOOOOOOOOOOOOOOOOOOO hzt CCCCCCCCCCCCCCCCCCCCCCCORZ
#include <algorithm>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <vector>
using namespace std;
using LL = long long;
using PII = pair<int, int>;
constexpr int kN = 1e4 + 1;
constexpr double kE = 1e-12;
int n;
double eu;
double s[kN], k[kN], v0[kN];
double v[kN];
bool J(double h) {
double e = 0;
for (int i = 1; i <= n; i++) {
double l = max(v0[i], 0.0), r = 1e5;
while (l < r - kE) {
double m = (l + r) / 2;
if (2 * h * k[i] * m * m * (m - v0[i]) < 1.0 + kE) {
l = m;
} else {
r = m;
}
}
v[i] = l;
e += k[i] * (v[i] - v0[i]) * (v[i] - v0[i]) * s[i];
}
return e < eu + kE;
}
int main() {
cin.tie(0)->sync_with_stdio(0);
cin >> n >> eu;
for (int i = 1; i <= n; i++) {
cin >> s[i] >> k[i] >> v0[i];
}
double l = 0, r = 1e5;
while (l < r - kE) {
double m = (l + r) / 2;
if (J(m)) {
r = m;
} else {
l = m;
}
}
double ans = 0;
for (int i = 1; i <= n; i++) {
ans += s[i] / v[i];
}
cout << fixed << setprecision(7) << ans << '\n';
return 0;
}
T2
给定一个长度为 \(n\) 的序列 \(a\)，构造一个自然数序列 \(b\) 使第 \(i\) 个数都不超过 \(a_i\)，且每个数之和为 \(k\)，并使 \(\sum\limits_{i = 1}^nb_i(a_i-b_i^2)\) 最大。\(n\leq 10^5, 1\leq a_i\leq 10^9\)。
Solution
暂时忽略 \(b_i\leq a_i\) 的限制。有 \(f(b_1, \cdots, b_n)=\sum\limits_{i = 1}^nb_i(a_i-b_i^2), \phi(b_1, \cdots, b_n)=k-(\sum\limits_{i=1}^n b_i)\)。\(F(b_1, \cdots, b_n, \lambda)=f(b_1, \cdots, b_n) + \lambda\phi(b_1, \cdots, b_n)\)。
对 \(b_i\) 求偏导：\(F_{b_i}^{'}=a_i-3b_i^2-\lambda=0\)，\(\lambda=a_i-3b_i^2\)。而 \(b_i\) 因它是自然数而随 \(\lambda\) 单调递增。于是可以二分 \(\lambda\) 算出所有 \(b_i\)。
然后再考虑加入 \(b_i\leq a_i\)。显然只需要让 \(\lambda\) 的下界为 \(\max\limits_{i=1}^n a_i-3a_i^2\) 即可。但是 \(\lambda\) 二分后答案为其下界时，可能算出来的序列 \(b\) 不满足 \(\sum\limits_{i=1}^n b_i = k\)。
注意到最优解一定是 \(a\) 升序排序后的一段前缀所对应的 \(b\) 卡到了 \(b_i=a_i\)。可以在外层再二分前缀长度，就可以得到实数最优解 \(b\)。
感性理解一下自然数解 \(ans_i\) 一定是实数解 \(b_i\) 向上或向下取整。可以先默认向下取整，然后往序列 \(ans\) 中再找出 \(c = k-(\sum\limits^n_{i=1} ans_i)\) 个数 + 1。这点很容易做到，找出 + 1 后对答案贡献最大的前 \(c\) 个数加上即可。
时间复杂度 \(O(n\log n\log V)\)。
\(\\\)
Code
// STOOOOOOOOOOOOOOOOOOOOOOOOO hzt CCCCCCCCCCCCCCCCCCCCCCCORZ
#include <algorithm>
#include <cmath>
#include <iostream>
#include <numeric>
#include <vector>
using namespace std;
using LL = long long;
using LD = long double;
using PII = pair<int, int>;
constexpr int kN = 1e5 + 1;
constexpr double kE = 1;
#define w first
#define i second
int n, p[kN], a[kN], b[kN];
LL k;
PII q[kN];
vector<PII> v;
bool J(int x) {
double c = k - accumulate(a + 1, a + x + 1, 0ll);
if (c < 0) {
return 0;
}
bool t = 0;
LD l = -1e19, r = 1e19;
for (int i = x + 1; i <= n; i++) {
l = max(l, a[i] - LD(3) * a[i] * a[i]);
r = min(r, LD(1) * a[i]);
}
while (l < r - kE) {
LD m = (l + r) / 2, p = 0;
for (int i = x + 1; i <= n; i++) {
p += sqrt((a[i] - m) / 3);
}
if (p > c) {
l = m, t = 1;
} else {
r = m;
}
}
for (int i = x + 1; i <= n; i++) {
b[i] = sqrt((a[i] - l) / 3);
}
return !t;
}
int main() {
cin.tie(0)->sync_with_stdio(0);
cin >> n >> k;
for (int i = 1; i <= n; i++) {
cin >> a[i];
q[i] = {a[i], i};
}
sort(q + 1, q + n + 1);
for (int i = 1; i <= n; i++) {
p[q[i].i] = i;
}
sort(a + 1, a + n + 1);
int l = 0, r = n + 1;
while (l <= r) {
int m = (l + r) / 2;
if (J(m)) {
l = m + 1;
} else {
r = m - 1;
}
}
copy_n(a + 1, l, b + 1);
J(l);
for (int i = 1; i <= n; i++) {
k -= b[i];
if (b[i] < a[i]) {
v.emplace_back(a[i] - 3 * b[i] * b[i] - 3 * b[i], i);
}
}
sort(v.begin(), v.end(), greater<>());
for (int i = 1; i <= k; i++) {
b[v[i - 1].i]++;
}
for (int i = 1; i <= n; i++) {
cout << b[p[i]] << ' ';
}
return 0;
}
