《转置原理的简单介绍》阅读随笔
目录
\(1.\) 概述
假设以下所有对单个元素的操作在一个数域 \(F\) 上进行。
\(\textbf{定义 1.1 }\text{(线性算法)}\)
输入一个 \(n\times 1\) 向量 \(\bm a\)，左乘一个 \(n\times n\) 的常数矩阵 \(A\) 后输出一个 \(n\times 1\) 的结果向量 \(\bm b = A\bm a\) 的算法，被称作线性算法。其中视矩阵 \(A\) 为常量，是不变的内容；其余内容为变量，由输入中得到。
容易发现 \(A\) 为方阵能包含所有可能的线性算法，因为可以适当补零。
\(\textbf{定义 1.1 }\text{(转置)}\)
对于形如 \(\bm b = A\bm a\) 的线性算法，称形如 \(\bm b' = A^{\mathsf T}\bm a'\) 的线性算法是原线性算法的转置。
转置原理断言，我们可以在保证时空复杂度不变的情况下，将一个线性算法改写为其的转置。这样我们可以首先观察解决原问题的线性算法的转置，从中得到优化方案，并将这个优化方案改写为解决原问题的算法。这常常可以简化问题，更容易导出结论。
接下来将讨论如何改写为转置。
\(2.\) 改写
我们首先需要明确的，是 \(A\) 矩阵表示着一系列对输入向量的基础操作的复合。这启发我们首先用矩阵表出基础操作，然后顺序地用这些矩阵操作输入向量，最后就能将矩阵等价地表为初等矩阵的复合。这将有利于观察 \(A^{\mathsf T}\) 的性质。
由线性代数的基础知识可以知道这复合定存在，而且根据高斯消元法，我们只需要三类初等矩阵。下面是对这三类矩阵的定义。
\(\textbf{定义 2.1 }\text{(初等矩阵)}\)
我们可以从单位矩阵 \(I\) 上做修改，得到以下三类矩阵：
- 将 \(I\) 的第 \(i\) 行和第 \(j\) 行交换。这矩阵左乘向量的效果是将向量的第 \(i\) 维和第 \(j\) 维值交换。
- 将 \(I\) 的第 \(i\) 行第 \(i\) 列元素由 \(1\) 改为数 \(a\)。这矩阵左乘向量的效果是将向量的第 \(i\) 维值乘 \(a\)。
- 将 \(I\) 的第 \(i\) 行第 \(i\) 列元素由 \(0\) 改为数 \(k\)。这矩阵左乘向量的效果是将向量的第 \(j\) 维值乘 \(k\) 得到的值加在第 \(i\) 维值上。
我们称这三类矩阵为初等矩阵。
随后可以将矩阵 \(A\) 表示为初等矩阵列 \(\langle E_i\rangle\) 的乘积：
不难得到
这样我们就可以得到矩阵对应操作在转置问题中的表示了：前两种操作原样不动，第三种是把贡献的位置调换。于是转置问题的算法就是倒序执行原问题的算法，并将形如 \(a_i \leftarrow a_i + ka_j\) 的操作改写为 \(a_j \leftarrow a_j + ka_i\)。具体实现时可能还需要将函数的输入和输出倒转。
容易发现复杂度是相同的。
这样我们就可以讨论转置问题的解法后无碍地应用在原问题上了。
这种做法即转置原理，或称特勒根原理(Tellegen's Principle)。
\(3.\) 例子
\(3.1\) \(\text{DFT}\)
考虑 \(\text{DFT}\) 的矩阵：
我们发现这个矩阵是对称矩阵，所以该算法的转置就是自身。\(\text{IDFT}\) 类似。
\(3.2\) 多项式乘法
设 \(\text M(n)\) 为求两个多项式的卷积的前 \(2n\) 位的时间。一般假定 \(\text M(n) = O(n\log n)\)，但需要注意常数问题。
假设 \(A, B\) 是两个多项式，\(C = A\times B\) 为这两个多项式作卷积的答案。考虑这个过程的转置。
我们视 \(A\) 为输入，\(B\) 为常量，则该线性算法 \(C = MA\) 的系数矩阵 \(M\) 由 \(M_{i, j} = [i\ge j]B[i - j]\) 确定。作转置，记转置问题为 \(C' = M^{\mathsf T} A'\)，可以得到
这就是减卷积的形式了。可以利用循环卷积自然溢出的性质做到 \(\frac{1}{2} \text M(n)\)。
我们记 \(\times^{\mathsf T}\) 为转置后的多项式乘法，转置符号右侧。
\(4.\) 一类特殊矩阵乘向量的快速算法
\(4.1\) 描述
考虑一个线性变换 \(A\bm\alpha = \bm\beta\)。考虑构造 \(n\times n\) 系数矩阵 \(A\) 的系数的二元生成函数
假设我们可以将 \(F(x, y)\) 表示为 \(u(x)v(y)f\left(g(x)h(y)\right)\) 的形式，其中 \(g, h\) 由常数个下方给出的简单函数复合而成，则我们就可以在 \(\widetilde O(\text M(n))\) 的时间复杂度内完成该矩阵（和他的逆）左乘向量。具体地，如果 \(g, h\) 中有 \(\exp, \log\)，则复杂度为 \(O(\text M(n)\log n)\)，反之为 \(O(\text M(n))\)。
为保证运算（主要是求逆）有意义，需要保证 \(gh\) 的常数项为 \(0\)，\(g, h\) 的一次项都非 \(0\)，\(u, v\) 常数项非 \(0\)，\(f\) 任一项非 \(0\)。
\(4.2\) 复合
\(\textbf{定义 4.1 }\text{(右复合)}\)
考虑 \(n - 1\) 次多项式 \(F\) 和常数多项式 \(G\)，以及它们的线性变换 \(F(G(x))\pmod{x^n}\)。我们称这是一次右复合，记作 \(F\circ G\)。
容易发现右复合满足结合律。
对于任意多项式的复合，我们并没有什么好方法，目前常用的复杂度仍然是 \(O(n^2)\)。然而对于一系列简单多项式，可以应用复合的结合律，每次将一个简单函数复合入其左侧所有函数的复合结果中，这样可以降低复杂度。除了少部分不封闭的函数需要右侧函数的信息外，复合的计算都可以每一步不丢失精度地截断到 \(x^n\) 次。
列出简单多项式：
-
加法 \(x + k\)：
可以简单地写为\[F(x + k) = \sum_{i = 0}^{n - 1} F[i](x + k)^i = \sum_{i = 0}^{n - 1}F[i] \sum_{j = 0}^i \binom{i}{j} x^j k^{i - j} = \sum_{j = 0}^{n - 1 }\left( \sum_{i = j}^{n - 1} i! \frac{k^{i - j} F[i]}{(i - j)!} \right) \frac{x^j}{j!} \]这可以在 \(O(\text M(n))\) 的复杂度内完成。
-
乘法 \(kx\)
代入即可，第 \(i\) 项系数乘 \(k^i\)。这可以在 \(O(n)\) 的复杂度内完成。 -
幂函数 \(x^k\)
直接变换下标即可，不进行任何运算。这可以在 \(O(n)\) 的复杂度内完成。 -
逆 \(x^{-1}\)
令 \(\text{rev}(f, n)\) 表示 \(n - 1\) 次多项式 \(f\) 翻转系数得到的多项式。有 \(F(x^{-1}) = \text{rev}(f, n)(x) \times x^{1 - n}\)。
这可以直接代入下一个多项式计算。计算多项式求逆和幂的复杂度为 \(O(\text M(n))\)，容易知道这是瓶颈。 -
根 \(x^{1/k}\)
这里要保证 \(x^{1 / k}\) 唯一。
如果接着复合，\(x\) 就代表了一个幂级数 \(g\)。假设 \(g = h^k\)，这里也需要讨论使得 \(h\) 唯一。我们将 \(F\) 的下标作带余除法 \(F_{ik + j}\)，按模 \(k\) 分组\[F_i(x) = \sum_{j} F_{jk + i} h^j \]则可以发现
\[F(h) = \sum_{i} F_{i}(g) h^i \]递归后拼合即可。求 \(h\) 需要 \(O(\text M(n))\) 的时间，随后需要 \(O(k\text M(n))\) 的时间求出各次幂，因此需要满足 \(k\) 是常数。
-
指数 \(\exp(x) - 1\)
减一是为了与对数互逆。
若不考虑减一，则有\[F(e^x) = \sum_{i = 0}^{n - 1} F[i] e^{ix} = \sum_{i = 0}^{n - 1} F[i] \sum_{j \ge 0} \frac{(ix)^j}{j !} = \sum_{j\ge 0} \left(\sum_{i=0}^{n-1} F[i]i^j \right) \frac{x^j}{j!} \]只需要考虑求
\[\sum_{j\ge 0} \left(\sum_{i=0}^{n-1} F[i]i^j \right) x^j \]这可以分治+多项式乘法解决。总时间复杂度 \(O(\text M(n)\log n)\)。观察内层求和号可以发现这类似于多项式 \(F\) 带入 \(i\in [0, n)\) 这些点值的情况，实际上不难发现这就是多点求值的转置问题。
-
对数 \(\log(x + 1)\)
不妨将其视作复合指数算法的逆。复合指数算法就是首先对 \(F\) 运行多点求值的转置算法，再在每一项乘 \(1 / j!\)。转置问题就是先在每一项乘 \(1 / j!\)，再对 \(F\) 运行多点求值。这个转置问题的逆显然为对 \(F\) 运行插值后在每一项乘 \(j!\)。再转置回去就能得到 \(O(\text M(n) \log n)\) 复杂度的算法。
多点求值和插值的做法将在下方详细说明。
\(4.3\) 计算
回到最开始的问题。我们假设 \(u = v = 1\)，可以得到
可以视作 \(g\) 右复合，点乘 \(f\)，转置 \(h\) 右复合。按 \(5.2\) 相关知识计算即可。\(u, v \neq 1\) 时相同，可以视作 \(u\) 乘原线性变换，再乘 \(v\) 的转置。
如果原问题满足限制，则这些操作都存在，且可以直接求逆，这就可以直接计算逆矩阵。
综上，我们以很低的复杂度解决了这些问题。
\(5.\) 范德蒙德矩阵的转置
\(5.1\) 定义与多点求值
\(\textbf{定义 5.1 }\text{(范德蒙德矩阵)}\)
我们定义一个范德蒙德矩阵 \(V_\alpha\) 是一个 \(n\times n\) 的方阵，其形如
\[V_\alpha = \begin{bmatrix} 1 & \alpha_0 & \alpha_0^2 & \cdots & \alpha_0^{n - 1} \\ 1 & \alpha_1 & \alpha_1^2 & \cdots & \alpha_1^{n - 1}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \alpha_{n - 1}^{1} & \alpha_{n - 1}^{2} & \cdots & \alpha_{n - 1}^{n - 1} \\ \end{bmatrix}\]其中 \(\langle \alpha \rangle\) 为一个长为 \(n\) 的序列。
我们可以通过范德蒙德矩阵来刻画一类操作。记 \(n - 1\) 次多项式 \(F\) 的系数组成的 \(1\times n\) 向量为 \(\bm f\)。
\(\textbf{定义 5.2 }\text{(多点求值)}\)
对于一个多项式 \(F\)，一类线性变换是具有重要意义的，其形如
\[V_\alpha \bm f^{\mathsf T} = \left(F(\alpha_0), F(\alpha_1), F(\alpha_2), \dots, F(\alpha_{n - 1}) \right)^{\mathsf T} \]此变换常被称为多点求值。
下面我们将应用转置原理大幅度减小多点求值的常数。首先介绍几个和该问题相关的问题。
\(5.2\) 多项式带余除法
给定 \(n - 1\) 次多项式 \(f\) 和 \(m - 1\) 次多项式 \(g(m \le n)\)，求 \(q, r\) 使得 \(f = qg + r\) 使得 \(deg(r) < m\)。\(m\) 一般为 \(n / 2\)。对于该问题，Sieveking-Kung's algorithm 给出了计算公式：
形式幂级数除法的复杂度是 \(\frac{7}{3}\text M(m) = \frac{7}{6} \text M(n)\)，可以查阅 negiizhao 的博客得到。随后算出 \(r = f - gh\) 还需要 \(\text M(m) = \frac{1}{2} \text M(n)\) 的复杂度，总时间复杂度为 \(\frac{5}{3} \text M(n)\)。
\(5.3\) 多个单项式求积
给定一系列 \(\langle a_i \rangle\)，求出 \(\prod_{i = 1}^n (1 + a_ix)\)。
这可以直接作分治。我们可以直接设 \(F(l, r) = \prod_{i = l}^r (1 + a_ix)\)，转移就是
其中 \(\text{mid} = \left\lfloor\dfrac{l + r}{2}\right\rfloor\)。
这样的时间复杂度为 \(T(n) = 2T(n / 2) + \text M(n / 2) = \frac{1}{2} \text M(n) \log_2 n + O(\text M(n))\)。这里需要注意 \(\text M(n / 2)\) 是得到前 \(n\) 项。
\(5.4\) 数列多幂次求和
给定常数列 \(\langle a_i\rangle\)，令 \(f(k) = \sum_{i = 1}^n a_i^k\)，求出 \(f(0), f(1), \dots, f(n - 1)\)。
直接构造 \(f\) 的生成函数 \(F\)，我们有
我们需要的就是 \(F(x)\pmod{x^n}\)。
不妨考虑将除法转化为乘法，也就是
这是 \(F\times A = B\) 的形式。我们只需要分别求出 \(A, B\) 在第 \(n\) 项截断的值即可应用多项式除法得到答案。\(A\) 是好求的，我们如上地设计 \(A(l, r)\) 即可得到。设计 \(B(l, r)\) 同 \(A(l, r)\)，可以发现若枚举缺口则得到
其中 \(\text{mid} = \left\lfloor\dfrac{l + r}{2}\right\rfloor\)。
分治的时间复杂度为 \(\frac{3}{2} \text M(n) \log_2 n + O(\text M(n))\)。
\(5.5\) 多点求值和插值
多点求值的做法传统上是应用分治求解。考虑我们需要让多项式 \(B\) 带入 \(\langle a_m \rangle\) 中的点值，当 \(deg(B) \ge m\) 时，我们可以将 \(B\) 对 \(\prod_{i = 0}^{m - 1} (x - a_i)\) 取模。这样 \(deg(B) < m\) 恒成立。这样我们每次可以将点值平均分为两个部分，分治求解。
这样的时间应用 \(4.\) 中的做法与 negiizhao 的实现，可以得到复杂度为 \(\frac{11}{3} \text M(n) \log_2 n + O(\text M(n))\)。对于流行的牛顿法求逆未优化的实现，这个速度还会慢一倍以上。
下面我们考虑范德蒙德矩阵的转置
考虑转置后的问题，即 \(\bm c = V_a^{\mathsf T}\bm b\)。可以发现
可以看作求多项式 \(G\) 满足
容易发现我们可以应用 \(5.3\) 中提到的分治做法求解。具体地，我们对分治区间 \([l, r)\) 维护分子 \(P_{[l, r)} (x) = \sum_{l\le i\le r} F_i\prod_{j \neq i, l \le i \le r} (1 - a_ix)\)，分母 \(Q_{[l, r)} (x) = \sum_{l\le i\le r} (1 - a_ix)\)。
可以写出和 \(5.3\) 相同的转移。
最后 \(\dfrac{P_{[0, n)}(x)}{Q_{[0, n)}(x)}\) 即为答案。
这个算法分为两步：
- 将分治树的每一层看作一个线性变换，按区间长度由短到长依次进行。
- 将 \(1/ Q_{[0, n)}(x)\) 看作由转移矩阵得到的常量，作乘法。
可以写出转置算法，具体如下。
\(Q_{[l,r)}(x)\) 与输入向量无关，是常量矩阵对应的值，可以预先采用分治法预处理。随后计算 \(P'_{[0, n)} = F\times^{\mathsf T} \dfrac{1}{ Q _{[0, n)}(x)}\)。
分析原分治树，将每一层看作一个线性变换。我们将 \(P_{[0, \text{mid})}\) 和 \(P_{[\text{mid}, r)}\) 看作输入，最终的 \(P_{[0, n)}(x)\) 看作输出。记 \(P_l = P_{[0, \text{mid})}, P_r = P_{[\text{mid}, r)}\)。
可以得到
可以得到该算法的转移矩阵 \(A_{k, i} = [0\le i < m] Q_r[k - i] + [m \le i < n] Q_l[k - i - m]\)。转置得到 \(A_{k, i}^{\mathsf T} = [0\le i < m] Q_r[i - k] + [m \le i < n] Q_l[i - k - m]\)
若转置算法的结果为 \(P'\)，则可以得到
不难发现这是转置乘法的形式，即
其实还可以直接导出这结论。
由线性性，我们有
直接作转置得到
最终叶子节点的 \(P\) 的常数项即为答案，正如转置中它们作为初始值一样。这样，我们就得到了一种常数较小且较好实现的做法。
插值的算法复杂度在于多点求值，总体可以自然过渡。
\(5.6\) 实践
作者（rushcheyo）的一份实现提交在此处。作者认为这种新方法码量小，常数好，一定程度上拓展了多点求值的实用性。
笔者也实现了一份，提交在此处。这里是经过封装的板子。核心代码如下：
多点求值与插值
namespace __multipoint_operation__ {
const int Len = 3e5 + 10;
#define ls (p << 1)
#define rs (p << 1 | 1)
poly __Q[Len << 2];
poly _E_Mul(poly A, poly B) {
int n = A.size(), m = B.size();
B.rev(); B = A * B;
for (int i = 0; i < n; ++ i)
A[i] = B[i + m - 1];
return A;
}
void _E_Init(int p, int l, int r, poly& a) {
if (l == r) {
__Q[p].resize(2);
__Q[p][0] = 1, __Q[p][1] = (a[l] ? mod - a[l] : a[l]);
return;
} int mid = l + r >> 1;
_E_Init(ls, l, mid, a), _E_Init(rs, mid + 1, r, a);
__Q[p] = __Q[ls] * __Q[rs];
}
void _E_Calc(int p, int l, int r, poly F, poly& g) {
F.resize(r - l + 1);
if (l == r) return void(g[l] = F[0]);
int mid = l + r >> 1;
_E_Calc(ls, l, mid, _E_Mul(F, __Q[rs]), g);
_E_Calc(rs, mid + 1, r, _E_Mul(F, __Q[ls]), g);
}
poly __P[Len << 2];
void _I_Init(int p, int l, int r, const poly& x) {
if (l == r) {
__P[p].resize(2), __P[p][0] = (x[l] ? mod - x[l] : 0), __P[p][1] = 1;
return;
} int mid = l + r >> 1;
_I_Init(ls, l, mid, x), _I_Init(rs, mid + 1, r, x);
__P[p] = __P[ls] * __P[rs];
}
poly _I_Calc(int p, int l, int r, const poly& t) {
if (l == r) return poly(1, t[l]);
int mid = l + r >> 1;
poly L(_I_Calc(ls, l, mid, t)), R(_I_Calc(rs, mid + 1, r, t));
L = L * __P[rs], R = R * __P[ls];
for (int i = 0; i < (int)R.size(); ++ i) {
L[i] = L[i] + R[i]; if (L[i] >= mod) L[i] -= mod;
} return L;
}
#undef ls
#undef rs
}
inline poly poly::eval(int n, poly a) const {
using namespace __multipoint_operation__;
n = max(n, size()); poly v(n);
poly F(f); F.resize(n + 1), a.resize(n);
_E_Init(1, 0, n-1, a); __Q[1].resize(n + 1);
_E_Calc(1, 0, n - 1, _E_Mul(F, __Q[1].inv()), v);
return v;
}
inline poly poly::intp(int n, const poly& x, const poly& y) {
using namespace __multipoint_operation__;
_I_Init(1, 0, n - 1, x);
__P[1] = __P[1].deri();
poly t = __P[1].eval(n, x);
for (int i = 0; i < n; ++ i)
t[i] = 1ll * y[i] * qp(t[i], mod - 2) % mod;
f = _I_Calc(1, 0, n - 1, t);
return *this;
}
可以发现这确实增加了多点求值的可实现性，希望能在未来看到这类算法更广泛的应用。
\(6.\) 应用
\(6.1\) 基变换
\(\textbf{定义 6.1 }\text{(多项式空间的基)}\)
考虑一个多项式集 \(\left\{ \varphi_0(x), \varphi_1(x), \dots, \varphi_n(x) \right\}\)。若任何 \(n - 1\) 次多项式都可以被他们线性表出，则称这是模 \(x^n\) 线性空间的一组基。OI 相关内容中常见单项式基和下降幂基。
\(6.2\) 例题
给定 \(n\) 个物品，第 \(i\) 个物品有属性 \(p_i(0 < p_i \le 1)\)。
主人公初始等级为 \(0\)。使用第 \(i\) 件物品有 \(p_i\) 的概率使等级加 \(1\)，\(1 - p_i\) 的概率不变。若最后等级为 \(k\)，则攻击力为 \(a_k\)。一件物品仅可使用一次。
设 \(f_i\) 表示使用除第 \(i\) 件物品外所有物品的情况下主人公攻击力的期望，求出 \(f_1, f_2, \dots, f_n\)。答案模 \(998244353\)，\(p_i\) 以分数 \(x_i / y_i\) 的方式给出，其中 \(1\le x_i\le y_i<998244353\)。
\(n\le 2\times 10^5, 0\le a_i \le a_{i + 1}\le 998244352\)。时间限制 \(2s\)。
容易发现，我们构造 \(C_i = (1 - p_i) + p_ix, C = \prod_{i = 1}^n C_i, A = \sum_{i = 0}^{n - 1} a_i x^i\)，则 \(f_i\) 就是 \(C / C_i\) 与 \(A\) 卷积的系数之和。
将 \(a\) 看作输入向量，\(f\) 看作输出向量，则转移矩阵 \(A_{i, j} = [x^j](C / C_i)\)。
作转置得到线性算法 \(g = A^{\mathsf T} a\)，可以写出 \(g = \sum_{i \ge 0} a_i C / C_i\)。
可以发现 \(g\) 形如上文中的 \(P_{[0, n)}\)，因此可以采用类似的方法得到答案，转置问题也如上文一样解决即可。
code
int n, pv[N], t1, t2;
poly a, seg[N], v;
#define ls (p << 1)
#define rs (p << 1 | 1)
void init(int p, int l, int r) {
if (l == r) {
seg[p].resize(2), seg[p][0] = (mod + 1 - pv[l]), seg[p][1] = pv[l];
return;
} int mid = l + r >> 1;
init(ls, l, mid); init(rs, mid + 1, r);
seg[p] = seg[ls] * seg[rs];
}
void calc(int p, int l, int r, poly F, poly& v) {
using namespace polynormial::__multipoint_operation__;
if (l == r) return void(v[l] = F[0]);
F.resize(r - l + 1); int mid = l + r >> 1;
calc(ls, l, mid, _E_Mul(F, seg[rs]), v);
calc(rs, mid + 1, r, _E_Mul(F, seg[ls]), v);
}
signed main() {
iot; cin >> n; a.resize(n), v.resize(n);
rep(i,0,n-1) cin >> a[i];
rep(i,0,n-1) cin >> t1 >> t2, pv[i] = 1ll * t1 * qp(t2, mod - 2) % mod;
init(1, 0, n - 1);
calc(1, 0, n - 1, a, v);
cout << v << '\n';
}
参考资料：
《转置原理的简单介绍》, rushcheyo,
转置原理小记, command_block,
关于优化形式幂级数计算的 Newton 法的常数, negiizhao.
以下是博客签名，与正文无关。
请按如下方式引用此页：
本文作者 joke3579，原文链接：https://www.cnblogs.com/joke3579/p/paperessay230119.html。
遵循 CC BY-NC-SA 4.0 协议。
请读者尽量不要在评论区发布与博客内文完全无关的评论，视情况可能删除。
