转置原理小记
command_block
2021-03-23 08:22:08
**参考资料** : 候选队论文 2020，陈宇 ：《转置原理的简单介绍》
# 转置原理 ： 定义 & 概述
- **线性算法**
给定过程矩阵 $A$ 以及输入向量 $a$ ，求解输出向量 $b=Aa$ 的算法，被称为线性算法。（又称线性变换）
写成更显式的求和，即 ：
$$b_k=\sum\limits_{i=0}A_{k,i}a_i$$
这里我们认为输入向量是变量，矩阵 $A$ 是常量。可以将 $A_{k,i}$ 理解为“该算法中 $a_i$ 向 $b_k$ 的贡献系数”。
为了方便处理，约定矩阵 $A$ 为方阵，否则也只需适当补零，不影响理论推导。
- **线性代数**
以下内容是线性代数常识，已熟悉的同学可以略过。
- **转置定理** ：对于矩阵 $C=AB$ ，有 $C^T=B^TA^T$。
- **初等矩阵**
- ① 将 $I$ 的第 $i,j$ 行交换的矩阵。该矩阵左乘一个向量的效果为，将该向量的第 $i,j$ 位交换。
- ② 将 $I_{i,i}$ 从 $1$ 改为 $c$ 的矩阵。该矩阵左乘一个向量的效果为，将该向量的第 $i$ 位乘以 $c$。
- ③ 将 $I_{i,j}$ 从 $0$ 改为 $c$ 的矩阵。该矩阵左乘一个向量的效果为，将该向量的第 $i$ 位加上第 $j$ 位的 $c$ 倍。
初等矩阵的转置仍然是初等矩阵。
不难发现，上述三种矩阵对应着高斯消元中的三种基本操作。这说明，所有非奇异矩阵都能表示为若干初等矩阵的乘积。
- **线性算法的转置**
形如 $b=Aa$ 的算法，与形如 $b'=A^Ta'$ 的算法互为转置。
- **转置原理**
转置原理给出了互为转置的两个线性算法之间的转化方法。
对于矩阵 $A$ ，将其表示为若干初等矩阵的乘积，得 ：
$$A=E_1E_2...E_m$$
在对某个向量执行线性算法时，只需将表示中的初等矩阵逐个乘上去。
根据上文中的转置定理，则有 ：
$$A^T=E_m^T...E_2^TE_1^T$$
这说明，根据 $A$ 的初等矩阵表示法，容易得到 $A^T$ 的初等矩阵表示法。
倒序逐乘初等矩阵的转置即可以同样的时间复杂度完成转置算法。
当然，在实际操作中，为了方便，我们不会真的以初等矩阵为单位来拆分并转置，而是以整合后的矩阵为单位描述算法。
在下文中，我们将介绍一系列经典线性算法及其转置。
# 线性算法的例子
- ### 离散傅里叶变换
观察 $\rm DFT$ 的矩阵 ：
$$
\begin{bmatrix}
1&1&1&...&1\\
1&\omega_n^1&\omega_n^2&\dots&\omega_n^{n-1}\\
1&\omega_n^2&\omega_n^4&\dots&\omega_n^{2(n-1)}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
1&\omega_n^{n-1}&\omega_n^{2(n-1)}&\dots&\omega_n^{(n-1)^2}\\
\end{bmatrix}
$$
该矩阵关于主对角线对称，故 $\rm DFT$ 算法的转置恰是其本身。
对于 $\rm IDFT$ 算法也有类似结论。
- ### 加法卷积
若将 $H=F*G$ 中的 $F$ 看成输入，$G$ 看成常量，则有 ：
$$H_k=\sum\limits_{i=0}G_{k-i}F_i$$
可知该线性算法 $H=AF$ 的矩阵为 $A_{i,j}=[i\geq j]G_{i-j}$
将其转置，记 $H'=A^TF$ ，则有 :
$$H_k'=\sum\limits_{i=k}G_{i-k}F_i=\sum\limits_{i=0}G_iF_{i+k}$$
即为减法卷积。
下面用 $\times^T$ 或 $*^T$ 表示转置后的多项式乘法。（将符号右边进行转置）
- ### 多点求值
考虑多项式 $F$ 和序列 $X$ ，求 $F$ 在每个 $X_i$ 处的取值。记答案为序列 $Y$。
可写出 ：
$$Y_k=F(X_k)=\sum\limits_{i=0}F_iX_k^i$$
可知该线性算法 $Y=AF$ 的矩阵为 $A_{i,j}=X_k^i$。
将其转置，记 $Y'=A^TF$ ，则有 :
$$Y_k'=\sum\limits_{i=0}F_iX_i^k$$
可以看做求一个多项式 $G$
$$
\begin{aligned}
G(x)&=\sum\limits_{k=0}x^k\sum\limits_{i=0}F_iX_i^k\\
&=\sum\limits_{i=0}F_i\sum\limits_{k=0}x^kX_i^k\\
&=\sum\limits_{i=0}\dfrac{F_i}{1-xX_i}
\end{aligned}
$$
接下来，我们要找出一个求解上式的高效线性算法。一个经典的算法是 ： 分治 `FFT`，复杂度为 $O(n\log^2n)$。
对于分治区间 $[l,r)$ ，维护分母 $Q_{[l,r)}(x)=\prod_{i=l}^r(1-xX_i)$ 与分子 $P_{[l,r)}(x)=\sum_{k=l}^rF_k\prod_{i=l,\ i\neq k}^r(1-xX_i)$。
在利用 $[l,m),[m,r)$ 的信息计算 $[l,r)$ 时，有 ：
$$
\begin{aligned}
Q_{[l,r)}&=Q_{[l,m)}Q_{[m,r)}\\
P_{[l,r)}&=P_{[l,m)}Q_{[m,r)}+P_{[m,r)}Q_{[l,m)}
\end{aligned}
$$
最终，计算 $P_{[0,n)}/Q_{[0,n)}$ 即可得到答案。
这个线性算法分为两步。
1. 将分治树的每一层看做一个线性算法，从低到高依次进行（利用多项式乘法）。
2. 最终，将 $\frac{1}{Q_{[0,n)}}$ 看做由转移矩阵得到的常量，进行多项式乘法。
所以，转置后的算法可以这样实现。
1. $Q_{[l,r)}$ 与输入向量无关，只和算法的矩阵有关，故为常量。先使用分治 FFT 预处理。
2. 计算 $P'_{[0,n)}=F\times^{T}\frac{1}{Q_{[0,n)}}$。
3. 分析原分治的转置，将每一层的转移看做一个线性算法。
将 $P_{[0,m)},P_{[m,n)}$ 看做输入，最终得到的 $P_{[0,n)}$ 看做输出。
简记 $Q^L=Q_{[0,m)},Q^R=Q_{[m,n)}$。
$$
\begin{aligned}
P_{[0,n)}&=P_{[0,m)}Q_R+P_{[m,n)}Q_L\\
\Rightarrow P_k&=\sum\limits_{i=0}^{\min(m-1,k)}P^*_iQ^R_{k-i}+\sum\limits_{i=0}^{\min(n-m-1,k)}P^*_{i+m}Q^L_{k-i}
\end{aligned}
$$
可以观察到该算法的转移矩阵为 $A_{k,i}=[0\leq i<m]Q^R_{k-i}+[m\leq i<n]Q^L_{k-i-m}$
故转置后的矩阵为 $A^T_{k,i}=[0\leq k<m]Q^R_{i-k}+[m\leq k<n]Q^L_{i-k-m}$
即转置算法的结果为 $P'$ ，写作和式，可得 ：
$$
\begin{aligned}
P_k'&=\sum\limits_{i=k}^{m-1}P^*_iQ^R_{i-k}&(0\leq k<m)\\
P_k'&=\sum\limits_{i=k+m}^{n-1}P^*_iQ^L_{i-k-m}&(m\leq k<n)
\end{aligned}
$$
不难发现，两者实际上是转置乘法。即
$$
\begin{aligned}
P'_{[0,m)}=P_{[0,n)}\times ^TQ^R\\
P'_{[m,n)}=P_{[0,n)}\times ^TQ^L
\end{aligned}
$$
有没有更简洁的方法来得到这一结论呢？
$P_{[0,m)},P_{[m,n)}$ 对 $P_{[0,n)}$ 的贡献是互相独立的（线性组合），可以写出贡献图：
$$P_{[0,m)}\xrightarrow{\times Q^R}P_{[0,n)}\xleftarrow{\times Q^L}P_{[m,n)}$$
故将贡献流向（以及边内部）取反后，有贡献图 ：
$$P_{[0,m)}\xleftarrow{\times^TQ^R}P_{[0,n)}\xrightarrow{\times^TQ^L}P_{[m,n)}$$
此时转移就一目了然了。
最终，分治树叶节点的 $P$ 的常数项即为答案。正如在转置算法中他们作为初始值。
至此，我们得到了一个常数较小，且易于实现（不需要多项式取模）的多点求值算法。
```cpp
#include<algorithm>
#include<cstring>
#include<cstdio>
#include<vector>
#define ll long long
#define ull unsigned ll
#define clr(f,n) memset(f,0,sizeof(int)*(n))
#define cpy(f,g,n) memcpy(f,g,sizeof(int)*(n))
#define pb push_back
const int _G=3,mod=998244353,Maxn=1<<16|500;
using namespace std;
ll powM(ll a,ll t=mod-2){
ll ans=1;
while(t){
if(t&1)ans=ans*a%mod;
a=a*a%mod;t>>=1;
}return ans;
}
const int invG=powM(_G);
int tr[Maxn<<1],tf;
void tpre(int n){
if (tf==n)return ;tf=n;
for(int i=0;i<n;i++)
tr[i]=(tr[i>>1]>>1)|((i&1)?n>>1:0);
}
void NTT(int *g,bool op,int n)
{
tpre(n);
static ull f[Maxn<<1],w[Maxn];w[0]=1;
for (int i=0;i<n;i++)f[i]=(((ll)mod<<4)+g[tr[i]])%mod;
for(int l=1;l<n;l<<=1){
ull tG=powM(op?_G:invG,(mod-1)/(l+l));
for (int i=1;i<l;i++)w[i]=w[i-1]*tG%mod;
for(int k=0;k<n;k+=l+l)
for(int p=0;p<l;p++){
int tt=w[p]*f[k|l|p]%mod;
f[k|l|p]=f[k|p]+mod-tt;
f[k|p]+=tt;
}
}if (!op){
ull invn=powM(n);
for(int i=0;i<n;++i)
g[i]=f[i]%mod*invn%mod;
}else for(int i=0;i<n;++i)g[i]=f[i]%mod;
}
void px(int *f,int *g,int n)
{for(int i=0;i<n;++i)f[i]=1ll*f[i]*g[i]%mod;}
#define Poly vector<int>
Poly operator + (const Poly &A,const Poly &B)
{
Poly C=A;C.resize(max(A.size(),B.size()));
for (int i=0;i<B.size();i++)C[i]=(C[i]+B[i])%mod;
return C;
}
Poly operator - (const Poly &A,const Poly &B)
{
Poly C=A;C.resize(max(A.size(),B.size()));
for (int i=0;i<B.size();i++)C[i]=(C[i]+mod-B[i])%mod;
return C;
}
Poly operator * (int c,const Poly &A)
{
Poly C;C.resize(A.size());
for (int i=0;i<A.size();i++)C[i]=1ll*c*A[i]%mod;
return C;
}
Poly operator * (const Poly &A,const Poly &B)
{
static int a[Maxn<<1],b[Maxn<<1];
if (min(A.size(),B.size())<=40){
Poly C;C.resize(A.size()+B.size()-1);
for (int i=0;i<A.size();i++)
for (int j=0;j<B.size();j++)
C[i+j]=(C[i+j]+1ll*A[i]*B[j])%mod;
return C;
}
cpy(a,&A[0],A.size());
cpy(b,&B[0],B.size());
Poly C;C.resize(A.size()+B.size()-1);
int n=1;for(n;n<C.size();n<<=1);
NTT(a,1,n);NTT(b,1,n);
px(a,b,n);NTT(a,0,n);
cpy(&C[0],a,C.size());
clr(a,n);clr(b,n);
return C;
}
Poly mulT(Poly &A,const Poly &B)
{
reverse(A.begin(),A.end());
Poly C=A*B;C.resize(A.size());
reverse(C.begin(),C.end());reverse(A.begin(),A.end());
return C;
}
void inv(const Poly &A,Poly &B,int n)
{
if (n==1)B.push_back(powM(A[0]));
else if (n&1){
inv(A,B,--n);
int sav=0;
for (int i=0;i<n;i++)sav=(sav+1ll*B[i]*A[n-i])%mod;
B.push_back(1ll*sav*powM(mod-A[0])%mod);
}else {
inv(A,B,n/2);
Poly sA;sA.resize(n);
cpy(&sA[0],&A[0],n);
B=2*B-B*B*sA;
B.resize(n);
}
}
Poly inv(const Poly &A)
{Poly C;inv(A,C,A.size());return C;}
Poly Q[Maxn<<1],P[22];
int sx[Maxn];
void solve1(int l,int r,int u)
{
if (l==r){
Q[u].pb(1);
Q[u].pb(mod-sx[l]);
return ;
}int mid=(l+r)>>1;
solve1(l,mid,u<<1);
solve1(mid+1,r,u<<1|1);
Q[u]=Q[u<<1]*Q[u<<1|1];
}
int sy[Maxn];
void solve2(int l,int r,int u,int dep)
{
if (l==r){sy[l]=P[dep][0];return ;}
int mid=(l+r)>>1;
P[dep].resize(r-l+1);
P[dep+1]=mulT(P[dep],Q[u<<1|1]);
solve2(l,mid,u<<1,dep+1);
P[dep+1]=mulT(P[dep],Q[u<<1]);
solve2(mid+1,r,u<<1|1,dep+1);
}
int n,m;Poly F;
int main()
{
scanf("%d%d",&n,&m);
F.resize(++n);
for (int i=0;i<n;i++)scanf("%d",&F[i]);
for (int i=0;i<m;i++)scanf("%d",&sx[i]);
solve1(0,m-1,1);Q[1].resize(max(n,m));
P[0]=mulT(F,inv(Q[1]));
solve2(0,m-1,1,0);
for (int i=0;i<m;i++)
printf("%d\n",sy[i]);
return 0;
}
```
# 例题
- [Comet#1179 [Contest #2]真实无妄她们的人生之路](https://www.cometoj.com/problem/1179)
**题意** ： 有 $n$ 件物品，第 $i$ 件物品属性为 $p_i$。
主人公等级初始为 $0$，使用第 $i$ 件物品会有 $p_i$ 的概率让等级加一，$1-p_i$ 的概率不变。
若最后等级为 $j$，则会具有 $a_j$ 的攻击力。
令 $f_i$ 表示使用除 $i$ 外的 $n-1$ 个物品后的期望攻击力。求出 $f_1,f_2, . . . ,f_n$。答案对 $998244353$ 取模。
记 $F_i(x)=(1-p_i)+p_ix,\ F(x)=\prod_{i=1}^nF_i$ ，则 $F/F_i$ 与 $a$ 点积后的系数和即为 $f_1$。
将 $a$ 看做输入向量， $f$ 看做输出向量，求解的过程是一个线性变换 $f=Aa$ ，其中 $A_{k,i}=[x^i](F/F_k)$。
考虑转置问题 $g=A^Ta$ ，则有 $g_k=\sum\limits_{i}[x^k](F/F_i)$ ，于是问题变为了求 $\sum\limits_{i=1}^nF/F_i$。
这显然也可以利用分治 FFT 求解 ：
对于分治区间 $[l,r)$ ，记 $Q_{[l,r)}(x)=\prod\limits_{i\in[l,r)}F_i,\ P_{[l,r)}(x)=\sum\limits_{i\in[l,r)}Q_{[l,r)}/F_i$
转移和上文中的多点求值类似。故转置后的形式也类似。
[评测记录](https://www.cometoj.com/status/bpgiy46veuwfksiealt0t5zhw6aqp7b4)
- [[数学记录]P7438~7440 「KrOI2021」Feux Follets 三连](https://www.luogu.com.cn/blog/command-block/post-shuo-xue-ji-lu-p7438-geng-jian-dan-di-pai-lie-ji-shuo)