文章目录
- 一、积和式的定义
- 二、Ryser算法
- 三、代码实现
一、积和式的定义
积和式(permanent)是一种和行列式长得很像的矩阵函数。在介绍积和式之前,我们先看看行列式(determinant)的定义。
首先需要引入“排列”(permutation)的概念。对于集合S={1,2,⋯,n}S=\{1,2,\cdots,n\}S={1,2,⋯,n},它的一个排列σ\sigmaσ就是对SSS中元素的一个重排。σ\sigmaσ的第iii个元素记作σi\sigma_iσi。例如,对于n=5n=5n=5,我们令σ={2,5,1,4,3}\sigma=\{2,5,1,4,3\}σ={2,5,1,4,3},则σ3=1\sigma_3=1σ3=1,σ5=3\sigma_5=3σ5=3。
排列的逆序对就是aaa在bbb前面但σa>σb\sigma_a>\sigma_bσa>σb的情况。例如σ={2,1,3,5,4}\sigma=\{2,1,3,5,4\}σ={2,1,3,5,4},有两个逆序对:(σ1,σ2)=(2,1)(\sigma_1,\sigma_2)=(2,1)(σ1,σ2)=(2,1)和(σ4,σ5)=(5,4)(\sigma_4,\sigma_5)=(5,4)(σ4,σ5)=(5,4)。一个排列σ\sigmaσ中逆序对的个数记作τ(σ)\tau(\sigma)τ(σ)。令sgn(σ)=(−1)τ(σ)\mathrm{sgn}(\sigma)=(-1)^{\tau(\sigma)}sgn(σ)=(−1)τ(σ)。对于一个排列σ\sigmaσ,如果你把其中的两个数互换,则sgn(σ)\mathrm{sgn}(\sigma)sgn(σ)会变号。所有nnn个元素的排列的集合记作SnS_nSn。例如,S3={(123),(132),(213),(231),(312),(321)}S_3=\{(1\ 2\ 3),(1\ 3\ 2),(2\ 1\ 3),(2\ 3\ 1),(3\ 1\ 2),(3\ 2\ 1)\}S3={(1 2 3),(1 3 2),(2 1 3),(2 3 1),(3 1 2),(3 2 1)}。
给定一个n×nn\times nn×n的矩阵A=(aij)n×nA=(a_{ij})_{n\times n}A=(aij)n×n,它的行列式为det(A)=∑σ∈Sn(sgn(σ)∏i=1nai,σi)\det(A)=\sum\limits_{\sigma\in S_n}\left(\mathrm{sgn}(\sigma)\prod\limits_{i=1}^{n}a_{i,\sigma_{i}}\right) det(A)=σ∈Sn∑(sgn(σ)i=1∏nai,σi)例如,当n=3n=3n=3时,设A=[abcdefghi]A=\begin{bmatrix}a&b&c\\d&e&f\\g&h&i\end{bmatrix}A=adgbehcfi,则det(A)=aei−afh+bfg−bdi+cdh−ceg\det(A)=aei-afh+bfg-bdi+cdh-ceg det(A)=aei−afh+bfg−bdi+cdh−ceg而积和式的定义就是在行列式中把sgn(σ)\mathrm{sgn}(\sigma)sgn(σ)去掉:perm(A)=∑σ∈Sn(∏i=1nai,σi)\mathrm{perm}(A)=\sum\limits_{\sigma\in S_n}\left(\prod\limits_{i=1}^{n}a_{i,\sigma_{i}}\right) perm(A)=σ∈Sn∑(i=1∏nai,σi)可以理解为:在矩阵中每行选取一个元素,且要求这些元素的列各不相同;将这些元素乘起来,得到一个乘积,积和式就是所有可能的选法对应的乘积之和。例如,当n=3n=3n=3时,设A=[abcdefghi]A=\begin{bmatrix}a&b&c\\d&e&f\\g&h&i\end{bmatrix}A=adgbehcfi,则perm(A)=aei+afh+bfg+bdi+cdh+ceg\mathrm{perm}(A)=aei+afh+bfg+bdi+cdh+ceg perm(A)=aei+afh+bfg+bdi+cdh+ceg积和式在量子场论、图论等领域中有应用。
积和式与行列式看起来只是某些项的符号不同,而且积和式看起来更简单了(没有sgn(σ)\mathrm{sgn}(\sigma)sgn(σ)),那是不是比行列式好算呢?答案是:大错特错!行列式可以用高斯消元法在O(n3)O(n^3)O(n3)的时间内算出来,而积和式目前最快的算法需要指数级的时间。事实上,1979年,Leslie G. Valiant证明了积和式的计算是#P\mathsf{\# P}#P完全问题,如果发现积和式有多项式时间的算法,那么将意味着FP=#P\mathsf{FP}=\mathsf{\#P}FP=#P,这是比P=NP\mathsf{P}=\mathsf{NP}P=NP还要强的命题。而大多数计算机科学家认为P≠NP\mathsf{P}\ne\mathsf{NP}P=NP,所以积和式大概率没有多项式时间的算法。我们要介绍的Ryser算法就是O(n2n)O(n 2^n)O(n2n)时间的。
二、Ryser算法
Ryser算法的核心思想就是容斥原理。我们还是先考察一下n=3n=3n=3的情况:令A=[abcdefghi]A=\begin{bmatrix}a&b&c\\d&e&f\\g&h&i\end{bmatrix}A=adgbehcfi,则perm(A)=aei+afh+bfg+bdi+cdh+ceg\mathrm{perm}(A)=aei+afh+bfg+bdi+cdh+ceg perm(A)=aei+afh+bfg+bdi+cdh+ceg观察式子T=(a+b+c)(d+e+f)(g+h+i)T=(a+b+c)(d+e+f)(g+h+i)T=(a+b+c)(d+e+f)(g+h+i),你会发现它的展开式中包含积和式的666个项(用蓝色标出):T=adg+adh+adi+aeg+aeh+aei+afg+afh+afi+bdg+bdh+bdi+beg+beh+bei+bfg+bfh+bfi+cdg+cdh+cdi+ceg+ceh+cei+cfg+cfh+cfi\begin{aligned} T&=a d g + a d h + a d i + a e g + a e h + \textcolor{blue}{a e i} + a f g + \textcolor{blue}{a f h} + a f i\\ &+b d g + b d h + \textcolor{blue}{b d i} + b e g + b e h + b e i + \textcolor{blue}{b f g} + b f h + b f i\\ &+c d g + \textcolor{blue}{c d h} + c d i + \textcolor{blue}{c e g} + c e h + c e i + c f g + c f h + c f i \end{aligned}T=adg+adh+adi+aeg+aeh+aei+afg+afh+afi+bdg+bdh+bdi+beg+beh+bei+bfg+bfh+bfi+cdg+cdh+cdi+ceg+ceh+cei+cfg+cfh+cfi于是,我们只需要在TTT的展开式中剔除不属于积和式的项就可以了。不属于积和式的项,也就是选取的某两个元素在同一列的项。这些项的特点是:元素的列组成的集合大小不超过222。比如adhadhadh一项,它只涉及第一和第二列,而没有涉及第三列,所以它不是积和式中的项。同样,cficficfi只涉及第三列,它也不是积和式中的项。我们可以枚举元素的列组成的集合(集合的大小为222),将对应的项剔除出去。
- 只涉及第一、二列的项:H12=(a+b)(d+e)(g+h)=adg+adh+aeg+aeh+bdg+bdh+beg+behH_{12}=(a+b)(d+e)(g+h)=a d g + a d h + a e g + a e h + b d g + b d h + b e g + b e hH12=(a+b)(d+e)(g+h)=adg+adh+aeg+aeh+bdg+bdh+beg+beh
- 只涉及第二、三列的项:H23=(b+c)(e+f)(h+i)=beh+bei+bfh+bfi+ceh+cei+cfh+cfiH_{23}=(b+c)(e+f)(h+i)=b e h + b e i + b f h + b f i + c e h + c e i + c f h + c f iH23=(b+c)(e+f)(h+i)=beh+bei+bfh+bfi+ceh+cei+cfh+cfi
- 只涉及第一、三列的项:H13=(a+c)(d+f)(g+i)=adg+adi+afg+afi+cdg+cdi+cfg+cfiH_{13}=(a+c)(d+f)(g+i)=a d g + a d i + a f g + a f i + c d g + c d i + c f g + c f iH13=(a+c)(d+f)(g+i)=adg+adi+afg+afi+cdg+cdi+cfg+cfi
只需要从TTT中把这些项剔除出去就可以了。但答案是perm(A)=T−H12−H23−H13\mathrm{perm}(A)=T-H_{12}-H_{23}-H_{13}perm(A)=T−H12−H23−H13吗?非也,因为H12H_{12}H12、H23H_{23}H23、H13H_{13}H13之间还有重叠部分,我们减的时候把重叠部分减了两次,还得加回来。H12H_{12}H12和H23H_{23}H23的重叠部分,就是只涉及第二列的项:behbehbeh。H12H_{12}H12和H13H_{13}H13的重叠部分则是只涉及第一列的项:adgadgadg。同理,H23H_{23}H23和H13H_{13}H13的重叠部分就是只涉及第三列的项——cficficfi了。
这样,我们得到计算三阶矩阵积和式的公式为:perm(A)=T−H12−H23−H13+adg+beh+cfi=(a+b+c)(d+e+f)(g+h+i)−(a+b)(d+e)(g+h)−(b+c)(e+f)(h+i)−(a+c)(d+f)(g+i)+adg+beh+cfi\begin{aligned} \mathrm{perm}(A)&=T-H_{12}-H_{23}-H_{13}+adg+beh+cfi\\ &=(a+b+c)(d+e+f)(g+h+i)-(a+b)(d+e)(g+h)-(b+c)(e+f)(h+i)-(a+c)(d+f)(g+i)+adg+beh+cfi \end{aligned}perm(A)=T−H12−H23−H13+adg+beh+cfi=(a+b+c)(d+e+f)(g+h+i)−(a+b)(d+e)(g+h)−(b+c)(e+f)(h+i)−(a+c)(d+f)(g+i)+adg+beh+cfi我们可以把这种容斥原理的思想推广到nnn阶矩阵的积和式。计算nnn阶矩阵的积和式的Ryser公式如下:perm(An×n)=(−1)n∑S⊆{1,2,⋯,n}[(−1)∣S∣∏i=1n(∑j∈Saij)]\mathrm{perm}(A_{n\times n})={(-1)}^{n} \sum\limits_{S\subseteq \{1,2,\cdots,n\}}\left[{(-1)}^{|S|}\prod\limits_{i=1}^{n}\left(\sum\limits_{j\in S}a_{ij}\right)\right] perm(An×n)=(−1)nS⊆{1,2,⋯,n}∑(−1)∣S∣i=1∏nj∈S∑aij这个公式可以这么理解:我们把AAA的行和之积展开,里面一定包含我们要求的积和式;然后减去涉及n−1n-1n−1列的项,加上涉及n−2n-2n−2列的项,减去涉及n−3n-3n−3列的项,……式中SSS就是涉及的列的集合,(−1)∣S∣(-1)^{|S|}(−1)∣S∣用于计算是加还是减;前面的(−1)n{(-1)}^{n}(−1)n是修正项,用于解决当nnn是奇数时,S={1,2,⋯,n}S=\{1,2,\cdots,n\}S={1,2,⋯,n}的情况下(−1)∣S∣{(-1)}^{|S|}(−1)∣S∣是负数的问题。
三、代码实现
理论上讲,如果我们按照格雷码顺序枚举SSS,那么时间复杂度可以降到O(n2n)O(n2^n)O(n2n)。但在这里我们为了方便起见就递归枚举SSS,对于每个SSS,计算各行的、列号为SSS的元素之和的乘积即可。下面给出一个时间复杂度为O(n22n)O(n^2 2^n)O(n22n)的C++实现:
#include <cstdint>typedef std::int64_t num;num recursion(int i, bool* b, int n, num** A)// 枚举S
{if(i == n) // 递归终点,已经得到一个S{num prod = 1;for(int row = 0; row < n; row++){num sum = 0;for(int col = 0; col < n; col++){if(b[col]){sum += A[row][col];}}prod *= sum;}int S_size = 0; // |S|for(int col = 0; col < n; col++){if(b[col]){S_size++;}}if(S_size % 2 == 1) // (-1)^|S|{prod = -prod;}return prod;}num result = 0;b[i] = true; // 选第i列result += recursion(i + 1, b, n, A);b[i] = false; // 不选第i列result += recursion(i + 1, b, n, A);return result;
}num ryser(int n, num** A)// 计算n x n矩阵A的积和式
{bool* b = new bool[n]; // S中是否含有第i列num result = recursion(0, b, n, A);delete []b;if(n % 2 == 1){result = -result; // (-1)^n}return result;
}