2020牛客暑期多校训练营（第一场）D

D. Quadratic Form

题目描述

Bobo has a positive-definite $n \times n$ matrix $A$ and an $n$ -dimension vector $b$ . He would like to find $x_1, x_2, \dots, x_n$ where $x_1, x_2, \dots, x_n \in \mathbb{R}$ , $\sum\limits_{i = 1}^n \sum\limits_{j = 1}^n a_{ij} x_i x_j \leqslant 1$ and $\sum\limits_{i = 1}^n b_i x_i$ is maximum.
It can be shown that $\left(\sum\limits_{i = 1}^n b_i x_i\right)^2 = \frac{P}{Q}$ , which is rational.
Find the value of $P \cdot Q^{-1} \bmod 998244353$ .

输入描述

The input consists of several test cases terminated by end-of-file.
The first line of each test case contains an integer $n$ .
The $i$ -th of the following $n$ lines contains $n$ integers $a_{i1}, a_{i2}, \cdots, a_{i n}$ .
The last line contains $n$ integers $b_1, b_2, \cdots, b_n$ .
$1 \leqslant n \leqslant 200$ , $0 \leqslant |a_{i j}|, |b_{i}| \leqslant 10^9$ , $a_{ij} =a_{ji}$ .
$\sum\limits_{i = 1}^n \sum\limits_{j = 1}^n a_{i j} x_i x_j > 0$ , for any $x \in \mathbb{R}^n$ .
$\mathrm{det}(A) \not\equiv 0 \pmod {998244353}$ .
The sum of $n$ does not exceed $10^4$ .

输出描述

For each test case, print an integer which denotes the result.

示例1

输入

2
1 0
0 1
0 0
2
1 0
0 1
1 1
2
8 4
4 3
1 1

输出

0
2
623902721

分析

将题目转化为约束条件下的最值问题。此处只讨论极值，不讨论边界可能存在的最值。
考虑到约束条件 $\sum\limits_{i = 1}^n \sum\limits_{j = 1}^n a_{ij} x_i x_j \leqslant 1$ 为不等式约束，可以添加 $\mathrm{KKT}$ 条件，利用拉格朗日乘数法求解。
首先需要讨论矩阵 $A$ 的一些性质。由于 $\mathrm{det}(A) \not\equiv 0 \pmod {998244353}$ ，故 $A$ 可逆。根据 $a_{ij}\in\mathbb Z$ 且 $a_{ij}=a_{ji}$ ，得 $A$ 为实对称矩阵，继而 $A^{-1}$ 也为实对称矩阵。又有 $\sum\limits_{i = 1}^n \sum\limits_{j = 1}^n a_{i j} x_i x_j > 0$ ，故 $A$ 为正定矩阵，继而 $A^{-1}$ 也为正定矩阵。上述性质都会在接下来的推导中涉及。
定义： $x=\begin{pmatrix}x_1\\x_2\\\vdots\\x_n\end{pmatrix}$ ， $b=\begin{pmatrix}b_1\\b_2\\\vdots\\b_n\end{pmatrix}$ 。则有约束条件 $x^TAx\leqslant1$ ，要求最大值的表达式为 $x^Tb$ 或 $b^Tx$ 。
令 $L(x_1,x_2,\cdots,x_n,\lambda)=\sum\limits_{i=1}^n b_ix_i+\lambda\left(\sum\limits_{i=1}^n\sum\limits_{j=1}^{n}a_{ij}x_ix_j-1\right)$ ；有 $\mathrm{KKT}$ 方程组如下。

$\begin{cases}\frac{\partial L}{\partial x_i}=0\quad i=1,2\cdots,n\\x^TAx\leqslant 1\\\lambda\geqslant 0\\\lambda (x^TAx-1)=0\end{cases}$

将 $\mathrm{KKT}$ 方程组的第一式展开，得到如下方程组。

$\begin{cases} b_1+\lambda(2a_{11}x_1+2a_{12}x_2+\cdots+2a_{1n}x_n)=0\\ b_2+\lambda(2a_{21}x_1+2a_{22}x_2+\cdots+2a_{2n}x_n)=0\\\quad\quad\quad\quad\quad\quad\quad\quad\quad\vdots\\b_n+\lambda(2a_{n1}x_1+2a_{n2}x_2+\cdots+2a_{nn}x_n)=0\end{cases}$

将上述方程组写成矩阵形式，有 $b+2\lambda Ax=0$ 。即 $\lambda x=-\frac{1}{2}A^{-1}b\ (1)$ ，等式两边取转置，得 $\lambda x^T=-\frac{1}{2}b^TA^{-1}\ (2)$ ； $(2)\times A\times(1)$ ，得 $\lambda^2x^TAx=\frac{1}{4}b^TA^{-1}b$ 。代入约束条件，移项后有 $x^TAx\le 1$ ， $\lambda^2\ge\frac{1}{4}b^TA^{-1}b$ ；由于 $A^{-1}$ 为正定矩阵，故 $\lambda\ge\frac{1}{2}\sqrt{b^TA^{-1}b}>0$ ；所以约束条件是有效的，即 $x^TAx=1\ (3)$ 。根据式 $(1),(2)$ ，可得到最优解为 $x=-\frac{1}{2\lambda}A^{-1}b,x^T=-\frac{1}{2\lambda}b^TA^{-1}$ ；将其带入式 $(3)$ ，有 $\frac{1}{4\lambda^2}b^TA^{-1}b=1\ (\ast)$ 。
由于 $\sum\limits_{i=1}^nb_ix_i=b^Tx$ ，代入最优解 $x=-\frac{1}{2\lambda}A^{-1}b$ ，得 $\sum\limits_{i=1}^nb_ix_i$ 的最大值为 $-\frac{1}{2\lambda}b^TA^{-1}b$ 。于是， $\left(\sum\limits_{i=1}^nb_ix_i\right)^2$ 最大值为 $\frac{1}{4\lambda^2}\left(b^TA^{-1}b\right)^2$ ；将式 $(\ast)$ 代入，有 $\frac{1}{4\lambda^2}\left(b^TA^{-1}b\right)^2=b^TA^{-1}b$ 。
综上所述，最终答案为 $b^TA^{-1}b\bmod 998244353$ 。套用矩阵求逆模板和矩阵乘法公式即可。

代码

/******************************************************************
Copyright: 11D_Beyonder All Rights Reserved
Author: 11D_Beyonder
Problem ID: 2020牛客暑期多校训练营（第一场） Problem D
Date: 7/22/2020 
Description: 
    Use lagrange multiplier method
    Compute the inverse matrix 
*******************************************************************/
#include<algorithm>
#include<iostream>
#include<cstdio>
using namespace std;
typedef long long ll;
const int N=202;
const ll mod=998244353;
ll a[N][N<<1],b[N],c[N];
int n;
void inverse();
ll qpow_mod(ll a,int b);
int main()
{
    while(~scanf("%d",&n))
    {
        int i,j,k;
        for(i=1;i<=n;i++)
        {
            for(j=1;j<=n;j++)
            {
                scanf("%lld",&a[i][j]);
                a[i][n+j]=0;//初始化
                a[i][j]%=mod;
            }
        }
        for(i=1;i<=n;i++)
        {
            scanf("%lld",b+i);
            b[i]%=mod;
        }
        inverse();
        //=========================================
        //a[1][1+n]*b[1]+a[2][1+n]*b[2]...
        //a[2][1+n]*b[1]+a[2][1+n]*b[2]...
        //...
        for(j=1;j<=n;j++)
        {
            ll res=0;
            for(i=1;i<=n;i++)
            {
                res+=a[i][j+n]*b[i];
                res%=mod;
            }
            c[j]=res%mod;
        }
        //c存b的转置由乘A的结果
        //==========================================
        //计算c右乘B
        ll ans=0;
        for(i=1;i<=n;i++) 
        {
            ans+=b[i]*c[i];
            ans%=mod;
        }
        //==========================================
        printf("%lld\n",ans);
    }
    return 0;
}
void inverse()//矩阵求逆模板
{
    int m=n+n;
    int i,j,k;
    for(i=1;i<=n;i++) a[i][i+n]=1;
    for(i=1;i<=n;i++)
    {
        for(j=i;j<=n;j++)
        {
            if(a[j][i])
            {
                for(k=1;k<=m;k++)
                {
                    swap(a[i][k],a[j][k]);
                }
            }
        }
        ll r=qpow_mod(a[i][i],mod-2);
        for(j=i;j<=m;j++) a[i][j]=r*a[i][j]%mod;
        for(j=1;j<=n;j++)
        {
            if(i==j) continue;
            ll rate=a[j][i]*a[i][i]%mod;
            for(k=i;k<=m;k++) 
            {
                a[j][k]=a[j][k]-rate*a[i][k]%mod;
                a[j][k]=(a[j][k]%mod+mod)%mod;
            }
        }
    }
}
ll qpow_mod(ll a,int b)
{
    ll res=1;
    while(b)
    {
        if(b&1) res=res*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return res;
}