QuadraticForm - 题解

标签与难度

标签: 线性代数, 矩阵, 高斯消元, 拉格朗日乘数法, 最优化问题, 数学 难度: 2200

题目大意喵~

主人你好呀，这道题是说，我们有一个 $n \times n$ 的正定矩阵 $A$ 和一个 $n$ 维向量 $b$。我们需要找到一个 $n$ 维实数向量 $x = (x_1, x_2, \dots, x_n)$，它需要满足一个不等式约束：

$$\sum_{i=1}^n \sum_{j=1}^n A_{i, j} x_i x_j \le 1$$

同时，要让另一个线性表达式的值最大化：

$$\sum_{i=1}^n b_i x_i$$

题目告诉我们，这个最大值的平方，也就是 $\left(\sum_{i=1}^n b_i x_i\right)^2$，一定可以表示成一个分数 $\frac{P}{Q}$。我们的任务就是计算 $P \cdot Q^{-1} \pmod{998244353}$ 的值，喵~

用矩阵的语言来描述一下就是：

• 约束条件：$x^T A x \le 1$
• 优化目标：最大化 $b^T x$

我们要找到这个最大值的平方，并在模意义下输出结果，呐。

解题思路分析

这道题看起来充满了复杂的数学公式，但别怕，让我一步步带你解开它的神秘面纱，喵~

第一步：理解问题本质

这是一个经典的约束最优化问题。我们想在由 $x^T A x \le 1$ 定义的一个 $n$ 维空间中的“椭球体”内部或边界上，找到一个点 $x$，使得它在向量 $b$ 方向上的投影最远（也就是 $b^T x$ 最大）。

直觉上，这个最大值点肯定是在椭球的边界上取得的，对吧？如果在内部，我们总可以再向 $b$ 的方向挪动一小步，让 $b^T x$ 的值变得更大，同时还不会跑出椭球。所以，我们的约束条件可以收紧为等式：

$$x^T A x = 1$$

第二步：处理二次型与矩阵对称性

表达式 $x^T A x$ 叫做二次型。这里有一个小技巧哦！对于任何方阵 $A$，二次型的值都只和它的对称部分有关。我们可以把 $A$ 分解成一个对称矩阵 $S$ 和一个反对称矩阵 $K$ 的和，$A = S+K$。

$$S = \frac{A + A^T}{2} \quad (\text{对称部分})$$

$$K = \frac{A - A^T}{2} \quad (\text{反对称部分})$$

其中 $A^T$ 是 $A$ 的转置矩阵。对于任何向量 $x$，都有 $x^T K x = 0$。所以：

$$x^T A x = x^T (S+K) x = x^T S x + x^T K x = x^T S x$$

这意味着，约束 $x^T A x = 1$ 其实等价于 $x^T S x = 1$。$S$ 是一个对称矩阵，处理起来会方便很多！

所以，我们的问题转化为了：

• 约束条件：$x^T S x = 1$，其中 $S = \frac{A+A^T}{2}$
• 优化目标：最大化 $b^T x$

第三步：拉格朗日乘数法登场！

对于这种带等式约束的优化问题，拉格朗日乘数法是我们的超级武器，喵！我们构造一个拉格朗日函数 $\mathcal{L}$：

$$\mathcal{L}(x, \lambda) = b^T x - \lambda (x^T S x - 1)$$

其中 $\lambda$ 是拉格朗日乘子。为了找到极值点，我们让 $\mathcal{L}$ 对 $x$ 的梯度（所有偏导数组成的向量）为零向量：

$$\nabla_x \mathcal{L} = \nabla_x(b^T x) - \lambda \nabla_x(x^T S x - 1) = 0$$

计算梯度：

• $\nabla_x(b^T x) = b$
• $\nabla_x(x^T S x) = 2Sx$ (因为 $S$ 是对称的)

所以梯度方程是：

$$b - 2\lambda S x = 0 \implies 2\lambda S x = b$$

因为 $A$ 是正定的，所以它的对称部分 $S$ 也是正定的，这意味着 $S$ 是可逆的。我们可以解出 $x$:

$$x = \frac{1}{2\lambda} S^{-1} b$$

第四步：求解与化简

我们得到了 $x$ 的表达式，但它还依赖于未知的 $\lambda$。别急，我们还有约束条件 $x^T S x = 1$ 没用呢！把 $x$ 的表达式代入约束中：

$$\left( \frac{1}{2\lambda} S^{-1} b \right)^T S \left( \frac{1}{2\lambda} S^{-1} b \right) = 1$$

$$\frac{1}{4\lambda^2} (S^{-1} b)^T S (S^{-1} b) = 1$$

因为 $(S^{-1})^T = (S^T)^{-1}$ 且 $S$ 对称 ($S=S^T$)，所以 $(S^{-1})^T = S^{-1}$。

$$\frac{1}{4\lambda^2} b^T S^{-1} S S^{-1} b = 1$$

$$\frac{1}{4\lambda^2} b^T S^{-1} b = 1 \implies 4\lambda^2 = b^T S^{-1} b$$

现在我们知道了 $4\lambda^2$ 的值。我们再来看看我们的目标函数 $b^T x$：

$$b^T x = b^T \left( \frac{1}{2\lambda} S^{-1} b \right) = \frac{1}{2\lambda} (b^T S^{-1} b)$$

我们要求的是 $(b^T x)^2$，那就平方一下：

$$(b^T x)^2 = \left( \frac{1}{2\lambda} (b^T S^{-1} b) \right)^2 = \frac{(b^T S^{-1} b)^2}{4\lambda^2}$$

把刚才求出的 $4\lambda^2 = b^T S^{-1} b$ 代入分母：

$$(b^T x)^2 = \frac{(b^T S^{-1} b)^2}{b^T S^{-1} b} = b^T S^{-1} b$$

哇！经过一番推导，我们得到了一个超级简洁的结论！我们要求的答案就是 $b^T S^{-1} b$。

第五步：制定算法

现在，任务就很明确了，就是计算 $b^T S^{-1} b \pmod{998244353}$。

1. 构造矩阵 $S$: $S = \frac{A+A^T}{2}$。在模运算下，除以2等于乘以2的逆元。所以 $S_{ij} = (A_{ij} + A_{ji}) \cdot (2^{-1}) \pmod{M}$。
2. 计算 $S^{-1} b$: 我们不想直接求出 $S^{-1}$ 这个大矩阵，太慢啦。我们可以把 $y = S^{-1} b$ 看作是线性方程组 $Sy = b$ 的解。
3. 求解线性方程组: 我们可以用高斯消元法来解方程组 $Sy = b$，得到向量 $y$。
4. 计算最终结果: 有了 $y = S^{-1} b$ 之后，我们要求的 $b^T S^{-1} b$ 就是 $b^T y$。

   $$b^T y = \sum_{i=1}^n b_i y_i$$

   计算这个点积，然后取模，就是最终答案啦！

总结一下我们的最终算法：

1. 根据输入的 $A$ 构造矩阵 $C = A+A^T$。
2. 构造增广矩阵 $[C | 2b]$。为什么是 $2b$？因为我们想解的方程是 $Sy=b$，即 $\frac{A+A^T}{2}y=b$, 这等价于 $(A+A^T)y = 2b$。
3. 使用高斯消元法求解线性方程组 $Cy = 2b$，得到解向量 $y$。
4. 计算最终答案 $\sum_{i=1}^n b_i y_i \pmod M$。

这样是不是清晰多啦？喵~ 让我们把这个思路变成代码吧！

代码实现

这是我根据上面的思路，精心重构的代码哦！注释超级详细，希望能帮到你，喵~

#include <iostream>
#include <vector>
#include <numeric>

using namespace std;

// 使用 long long 防止计算过程中溢出
using ll = long long;

const int MOD = 998244353;

// 快速幂函数，用于计算 a^b % MOD
// 在这里主要用来求逆元，即 a^(MOD-2) % MOD
ll power(ll base, ll exp) {
    ll res = 1;
    base %= MOD;
    while (exp > 0) {
        if (exp % 2 == 1) res = (res * base) % MOD;
        base = (base * base) % MOD;
        exp /= 2;
    }
    return res;
}

// 模块逆元函数
ll modInverse(ll n) {
    return power(n, MOD - 2);
}

// 高斯消元求解线性方程组 Ax = b
// A 是 n x n 矩阵, b 是 n 维向量
// 返回解向量 x
vector<ll> gauss_elimination(vector<vector<ll>>& a, vector<ll>& b) {
    int n = a.size();

    // 构造增广矩阵 [A | b]
    vector<vector<ll>> augmented_matrix(n, vector<ll>(n + 1));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            augmented_matrix[i][j] = a[i][j];
        }
        augmented_matrix[i][n] = b[i];
    }

    // --- 正向消元 ---
    for (int i = 0; i < n; ++i) {
        // 找到第 i 列绝对值最大的元素作为主元，以增加数值稳定性
        // 在模运算下，我们只需要找到一个非零元素即可
        int pivot_row = i;
        for (int k = i + 1; k < n; ++k) {
            if (augmented_matrix[k][i] != 0) {
                pivot_row = k;
                break;
            }
        }
        swap(augmented_matrix[i], augmented_matrix[pivot_row]);

        // 如果主元是0，说明矩阵奇异，可能无解或有无穷解
        // 在本题中，因为 S 是正定的，所以必有唯一解
        ll pivot_inv = modInverse(augmented_matrix[i][i]);

        // 将主元所在行归一化，使得 a[i][i] = 1
        for (int j = i; j <= n; ++j) {
            augmented_matrix[i][j] = (augmented_matrix[i][j] * pivot_inv) % MOD;
        }

        // 将其他行的第 i 列元素消为 0
        for (int k = 0; k < n; ++k) {
            if (k != i) {
                ll factor = augmented_matrix[k][i];
                for (int j = i; j <= n; ++j) {
                    ll term = (factor * augmented_matrix[i][j]) % MOD;
                    augmented_matrix[k][j] = (augmented_matrix[k][j] - term + MOD) % MOD;
                }
            }
        }
    }

    // --- 回代得到解 ---
    // 经过上面的操作，矩阵 A 部分已经变成了单位矩阵
    // 所以增广矩阵的最后一列就是解向量 x
    vector<ll> x(n);
    for (int i = 0; i < n; ++i) {
        x[i] = augmented_matrix[i][n];
    }

    return x;
}


int main() {
    // 加速输入输出，喵~
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);

    int n;
    while (cin >> n) {
        vector<vector<ll>> A(n, vector<ll>(n));
        vector<ll> b(n);

        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                cin >> A[i][j];
            }
        }
        for (int i = 0; i < n; ++i) {
            cin >> b[i];
        }

        // 1. 构造矩阵 C = A + A^T
        vector<vector<ll>> C(n, vector<ll>(n));
        for (int i = 0; i < n; ++i) {
            for (int j = 0; j < n; ++j) {
                C[i][j] = (A[i][j] + A[j][i] + MOD) % MOD;
            }
        }

        // 2. 构造向量 2b
        vector<ll> two_b(n);
        for (int i = 0; i < n; ++i) {
            two_b[i] = (2 * b[i] % MOD + MOD) % MOD;
        }
        
        // 3. 求解线性方程组 Cy = 2b, 得到 y
        // y 就是我们推导中的 S^{-1}b
        vector<ll> y = gauss_elimination(C, two_b);

        // 4. 计算最终答案 b^T * y
        ll ans = 0;
        for (int i = 0; i < n; ++i) {
            ans = (ans + (b[i] * y[i]) % MOD + MOD) % MOD;
        }

        cout << ans << endl;
    }

    return 0;
}

复杂度分析

• 时间复杂度: $O(N^3)$

  • 构造矩阵 $C=A+A^T$ 需要 $O(N^2)$ 的时间。
  • 求解 $N$ 元线性方程组的高斯消元法，主要计算量在三重循环，复杂度是 $O(N^3)$。
  • 最后计算点积 $b^T y$ 需要 $O(N)$ 的时间。
  • 所以总的时间复杂度由高斯消元主导，为 $O(N^3)$，喵~

• 空间复杂度: $O(N^2)$

  • 我们需要存储矩阵 $A$, $C$ 和增广矩阵，以及向量 $b$ 和 $y$。其中矩阵占用的空间是主要的，为 $O(N^2)$。

知识点总结

这道题是数学和算法的美妙结合，涉及到的知识点还不少呢！

1. 二次型与对称矩阵: 理解 $x^T A x = x^T \left(\frac{A+A^T}{2}\right) x$ 是关键，它让我们能把问题标准化，用对称矩阵来处理。
2. 约束最优化: 问题的本质是在一个约束条件下求函数的极值。
3. 拉格朗日乘数法: 解决带等式约束优化问题的强大数学工具，是推导出最终简洁公式的核心。
4. 线性方程组求解: 理论推导的最后，问题转化为了求解一个线性方程组 $Sy=b$。
5. 高斯消元法: 求解线性方程组的标准算法，也是本题在编程实现上的核心。
6. 模运算: 别忘了所有计算都要在模 $998244353$ 下进行，特别是除法要用乘以逆元的方式来代替。

希望这篇题解能帮助你完全理解这道题目！如果还有不明白的地方，随时可以再来问我哦，喵~ >w<