定数整列線形繰返し最適化マトリクス高速べき乗
23261 ワード
一般行列高速べき乗の形式:f(n)=Σi=1 k 1 a i f(n−i)+Σi=1 k 2 b i g(n−i)+cf(n)=sum_{i=1}^{k_1} a_if(n-i)+\sum_{i=1}^{k_2}b_ig(n−i)+cf(n)=i=1Σk 1 ai f(n−i)+i=1Σk 2 bi g(n−i)+cは、k 3 log
遷移f(n)=Σi=1 k a i f(n−i)f(n)=sum_{i=1}^k a_if(n−i)f(n)=i=1Σkaif(n−i)これが定数線形繰返しの基本形態である.一般的にはk logk logn klog klog n klogklognができますが、比較的簡単な書き方でもk 2 logn k^2log n k 2 lognができますが、最適化方法について説明します.
##特徴多項式前置技能:行列のランク
フィーチャー値、フィーチャーベクトル:定数がある場合λ\lambda λ,ベクトルv⃗vec v,満足λ v ⃗ = A v ⃗\lambda\vec v=A\vec v λv=Avは、ベクトルv⃗vec v vと呼ばれるマトリクスA A A Aの特徴ベクトルのセットであり、$lambdaはマトリクスA$の特徴値のセットである.ランクkのマトリクスには、k群の線形非相関特徴ベクトルがある.
フィーチャー多項式は関係式を変換します:(λ I − A ) v ⃗ = 0 (\lambda I-A)\vec v=0 (λI−A)v =0. この式の非ゼロ解の要件はd e t(λ I − A ) = 0 det(\lambda I-A)=0 det(λI−A)=0、すなわち行列がランク未満である.k k k次多項式を得ることができ,この多項式の値が0に等しい場合,この方程式を行列A A Aの特徴方程式と呼ぶ.この多項式を行列A A Aの特徴多項式と呼ぶ.k個の解はもちろんk k個の特徴値が辛い(解く必要はなく、どう処理すれば後で言う)!
計算数の基本定理に基づいて、最終的な多項式はk個の解があって、f(x)=を書くことができますΠ i ( λ i − x ) f(x)=\Pi_i(\lambda_i-x) f(x)=Πi(λi−x)であり,A−A行列自体に持ち込むと0 0 0行列(Cayley−Hamilton定理)が得られる.
証明:(以下の証明はブロガーが以前OIを学んだときに書いた比較的複雑な証明であり,実際には伴随行列で直接証明できる)f(A)=Π i ( λ i I − A ) f(A)=\Pi_i(\lambda_i I -A) f(A)=Πi(λi I−A)は、f(A)f(A)f(A)f(A)から得られる行列をそれぞれ特徴ベクトルに乗じることを考慮し、最後に0行列が得られた場合、これらの特徴ベクトルが線形に相関しないため、f(A)f(A)f(A)f(A)f(A)に任意のベクトルを乗じて0 0 0行列が得られ、A A A Aが0 0 0 0 0行列であることが証明される.
問題は、f(A)f(A)f(A)f(A)が任意の特徴ベクトル値に乗じて0 0行列を得ることを証明することに変換される.まず証明(λ i I − A ) ( λ j I − A ) = ( λ j I − A ) ( λ i I − A ) (\lambda_i I-A)(\lambda_j I - A)=(\lambda_j I - A)(\lambda_i I-A) (λiI−A)(λjI−A)=(λjI−A)(λi I−A)行列が加算交換則を満たすので,これは直接展開すればよい.証明が終わる.
現在乗算されている特徴ベクトルがva v a,あり:f(A)v⃗a=Π i ( λ i I − A ) ∗ v ⃗ a = Π i ! = a ( λ i I − A ) ∗ ( λ a I − A ) ∗ v ⃗ a\begin{aligned}f(A)\vec v_a&=\Pi_{i}(\lambda_iI-A)*\vec v_a\\&=\Pi_{i!=a}(\lambda_iI-A)*(\lambda_aI-A)*\vec v_a &&\end{aligned} f(A)v a=Πi(λiI−A)∗v a=Πi!=a(λiI−A)∗(λaI−A)∗v a
行列乗算は結合則を満たし,後半では,(λ a I − A ) ∗ v ⃗ a = λ a I v ⃗ a − A v ⃗ a (\lambda_aI -A)*\vec v_a=\lambda_aI\vec v_a-A\vec v_a (λaI−A)∗v a=λaIv a−Av a.
定義に基づいて、λ a I v ⃗ a − A v ⃗ a = 0\lambda_aI\vec v_a-A\vec v_a=0 λa Iv a−Av a=0であれば、代入還元式は0 0 0行列を得ることができる.
同時にA A Aに関するk k k k次多項式である.では、A w A^w Awを求めると、多項式の型取りを行い、余分な0行列を型抜きしてAに関するk−1 k−1次多項式を得ることができる:g(A)=A w%f(A)=Σi=0 k−1 a i A i g(A)=A^w%f(A)=sum_{i=0}^{k-1}a_iA^i g(A)=Aw%f(A)=i=0∑k−1aiAi
初期繰返し行列をh=h k,h k−1,h k−2,.,h 1 h={h_k,h_{k-1},h_{k-2},...,h_1} h=hk,hk−1,hk−2,...,h 1では、Σa i∗[(h∗A i)]11sum a_i*[(h*A^i)]_Σai∗[(h∗Ai)]11,h∗A i h*A^i h∗Aiの第1項はh i+1 h_{i+1}hi+1であれば,直接O(n)O(n)O(n)計算が可能になる.多項式型取りは高速フーリエ変換でk log便宜上、多項式はk 2 k^2 k 2アルゴリズムを模写することができ、総複雑度はO(k 2 logn)O(k^2log n)O(k 2 log n)である.
##特徴多項式を解く多項式の型取りの前提は多項式を知ることである.明らかにk k k個の特徴値を解くのは現実的ではなく、常係数の整列繰返しの特徴を考慮して、行列の特徴多項式を迅速に求めることができる:まず、転移方程式をリストする:
( a 1 a 2 a 3 ⋯ a k − 2 a k − 1 a k 1 0 0 ⋯ 0 0 0 0 1 0 ⋯ 0 0 0 0 0 1 ⋯ 0 0 0 ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ 0 0 0 ⋯ 1 0 0 0 0 0 ⋯ 0 1 0 ) ( h k − 1 h k − 2 h k − 3 ⋮ h 2 h 1 h 0 ) = ( h k h k − 1 h k − 2 ⋮ h 3 h 2 h 1 )\begin{pmatrix} a_1 & a_2 & a_3 &\cdots & a_{k - 2} & a_{k - 1} & a_k\\1 & 0 & 0 &\cdots & 0 & 0 & 0\\0 & 1 & 0 &\cdots & 0 & 0 & 0\\0 & 0 & 1 &\cdots & 0 & 0 & 0\\\vdots &\vdots &\vdots &\ddots &\vdots &\vdots &\vdots\\0 & 0 & 0 &\cdots & 1 & 0 & 0\\0 & 0& 0 &\cdots & 0 & 1 & 0\end{pmatrix}\begin{pmatrix}h_{k - 1}\\h_{k - 2}\\h_{k - 3}\\\vdots\\h_2\\h_1\\h_0\end{pmatrix} =\begin{pmatrix}h_{k}\\h_{k - 1}\\h_{k - 2}\\\vdots\\h_3\\h_2\\h_1\end{pmatrix} ⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛a1100⋮00a2010⋮00a3001⋮00⋯⋯⋯⋯⋱⋯⋯ak−2000⋮10ak−1000⋮01ak000⋮00⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛hk−1hk−2hk−3⋮h2h1h0⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛hkhk−1hk−2⋮h3h2h1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
フィーチャー多項式:f(λ ) = ∣ λ I − A ∣ = ( λ − a 1 − a 2 − a 3 ⋯ − a k − 2 − a k − 1 − a k − 1 λ 0 ⋯ 0 0 0 0 − 1 λ ⋯ 0 0 0 0 0 − 1 ⋯ 0 0 0 ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ 0 0 0 ⋯ − 1 λ 0 0 0 0 ⋯ 0 − 1 λ ) f(\lambda) = |\lambda I - A| =\begin{pmatrix}\lambda - a_1 & -a_2 & -a_3 &\cdots & -a_{k - 2} & -a_{k - 1} & -a_k\\-1 &\lambda & 0 &\cdots & 0 & 0 & 0\\0 & -1 &\lambda &\cdots & 0 & 0 & 0\\0 & 0 & -1 &\cdots & 0 & 0 & 0\\\vdots &\vdots &\vdots &\ddots &\vdots &\vdots &\vdots\\0 & 0 & 0 &\cdots & -1 &\lambda & 0\\0 & 0& 0 &\cdots & 0 & -1 &\lambda\end{pmatrix} f(λ)=∣λI−A∣=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛λ−a1−100⋮00−a2λ−10⋮00−a30λ−1⋮00⋯⋯⋯⋯⋱⋯⋯−ak−2000⋮−10−ak−1000⋮λ−1−ak000⋮0λ⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
最初の行の展開に対してそれぞれ和を求めることを考慮して、あります:
f ( λ ) = ( λ − a 1 ) A 11 + ( − a 2 ) A 12 + ⋯ + ( − a k ) A 1 n = λ k − a 1 λ k − 1 − a 2 λ k − 2 − ⋯ − a k f(\lambda) = (\lambda - a_1)A_{11} + (-a_2)A_{12} +\cdots + (-a_k)A_{1n} =\lambda ^ k - a_1\lambda ^ {k - 1} - a_2\lambda ^ {k - 2} -\cdots - a_k f(λ)=(λ−a1)A11+(−a2)A12+⋯+(−ak)A1n=λk−a1λk−1−a2λk−2−⋯−ak
そのうちA 1 i A_{1 i}A 1 iは、A A Aの代数余子式である.
これにより,定数整列線形繰返しも行列の高速べき乗を最適化できるようになった.
##小最適化1.結果を計算するのは面倒だと思いますが、もう一つ乗算しなければなりませんか?いいえ!A n−k A^{n−k}An−kではなく、A n−A^{n−k}An−kを直接計算すると、最後の結果はベクトルの最後のビットになる.このときA i(i≦k)A^{i}(ile k)Ai(i≦k)の最下位は原繰返し式の1つであり,直接求めることができる.
2.注意k 2 k^2 k 2は、i,j i,j i,j iの順で、i+j i+j i+j i+jに、i,j i,j i,j i,j,a i=a j∗a−j a_ではなく、i+j i+j i+j i+j i+jに加算される.i=a_j*a_{i-j} ai=aj∗ai−j.
BZOJ 4161のコードを添付します.これは定数係数が順次線形に繰り返される裸の問題です.k 2 logn k^2log n k 2 lognの方法で过ごすことができます.
遷移f(n)=Σi=1 k a i f(n−i)f(n)=sum_{i=1}^k a_if(n−i)f(n)=i=1Σkaif(n−i)これが定数線形繰返しの基本形態である.一般的にはk logk logn klog klog n klogklognができますが、比較的簡単な書き方でもk 2 logn k^2log n k 2 lognができますが、最適化方法について説明します.
##特徴多項式前置技能:行列のランク
フィーチャー値、フィーチャーベクトル:定数がある場合λ\lambda λ,ベクトルv⃗vec v,満足λ v ⃗ = A v ⃗\lambda\vec v=A\vec v λv=Avは、ベクトルv⃗vec v vと呼ばれるマトリクスA A A Aの特徴ベクトルのセットであり、$lambdaはマトリクスA$の特徴値のセットである.ランクkのマトリクスには、k群の線形非相関特徴ベクトルがある.
フィーチャー多項式は関係式を変換します:(λ I − A ) v ⃗ = 0 (\lambda I-A)\vec v=0 (λI−A)v =0. この式の非ゼロ解の要件はd e t(λ I − A ) = 0 det(\lambda I-A)=0 det(λI−A)=0、すなわち行列がランク未満である.k k k次多項式を得ることができ,この多項式の値が0に等しい場合,この方程式を行列A A Aの特徴方程式と呼ぶ.この多項式を行列A A Aの特徴多項式と呼ぶ.k個の解はもちろんk k個の特徴値が辛い(解く必要はなく、どう処理すれば後で言う)!
計算数の基本定理に基づいて、最終的な多項式はk個の解があって、f(x)=を書くことができますΠ i ( λ i − x ) f(x)=\Pi_i(\lambda_i-x) f(x)=Πi(λi−x)であり,A−A行列自体に持ち込むと0 0 0行列(Cayley−Hamilton定理)が得られる.
証明:(以下の証明はブロガーが以前OIを学んだときに書いた比較的複雑な証明であり,実際には伴随行列で直接証明できる)f(A)=Π i ( λ i I − A ) f(A)=\Pi_i(\lambda_i I -A) f(A)=Πi(λi I−A)は、f(A)f(A)f(A)f(A)から得られる行列をそれぞれ特徴ベクトルに乗じることを考慮し、最後に0行列が得られた場合、これらの特徴ベクトルが線形に相関しないため、f(A)f(A)f(A)f(A)f(A)に任意のベクトルを乗じて0 0 0行列が得られ、A A A Aが0 0 0 0 0行列であることが証明される.
問題は、f(A)f(A)f(A)f(A)が任意の特徴ベクトル値に乗じて0 0行列を得ることを証明することに変換される.まず証明(λ i I − A ) ( λ j I − A ) = ( λ j I − A ) ( λ i I − A ) (\lambda_i I-A)(\lambda_j I - A)=(\lambda_j I - A)(\lambda_i I-A) (λiI−A)(λjI−A)=(λjI−A)(λi I−A)行列が加算交換則を満たすので,これは直接展開すればよい.証明が終わる.
現在乗算されている特徴ベクトルがva v a,あり:f(A)v⃗a=Π i ( λ i I − A ) ∗ v ⃗ a = Π i ! = a ( λ i I − A ) ∗ ( λ a I − A ) ∗ v ⃗ a\begin{aligned}f(A)\vec v_a&=\Pi_{i}(\lambda_iI-A)*\vec v_a\\&=\Pi_{i!=a}(\lambda_iI-A)*(\lambda_aI-A)*\vec v_a &&\end{aligned} f(A)v a=Πi(λiI−A)∗v a=Πi!=a(λiI−A)∗(λaI−A)∗v a
行列乗算は結合則を満たし,後半では,(λ a I − A ) ∗ v ⃗ a = λ a I v ⃗ a − A v ⃗ a (\lambda_aI -A)*\vec v_a=\lambda_aI\vec v_a-A\vec v_a (λaI−A)∗v a=λaIv a−Av a.
定義に基づいて、λ a I v ⃗ a − A v ⃗ a = 0\lambda_aI\vec v_a-A\vec v_a=0 λa Iv a−Av a=0であれば、代入還元式は0 0 0行列を得ることができる.
同時にA A Aに関するk k k k次多項式である.では、A w A^w Awを求めると、多項式の型取りを行い、余分な0行列を型抜きしてAに関するk−1 k−1次多項式を得ることができる:g(A)=A w%f(A)=Σi=0 k−1 a i A i g(A)=A^w%f(A)=sum_{i=0}^{k-1}a_iA^i g(A)=Aw%f(A)=i=0∑k−1aiAi
初期繰返し行列をh=h k,h k−1,h k−2,.,h 1 h={h_k,h_{k-1},h_{k-2},...,h_1} h=hk,hk−1,hk−2,...,h 1では、Σa i∗[(h∗A i)]11sum a_i*[(h*A^i)]_Σai∗[(h∗Ai)]11,h∗A i h*A^i h∗Aiの第1項はh i+1 h_{i+1}hi+1であれば,直接O(n)O(n)O(n)計算が可能になる.多項式型取りは高速フーリエ変換でk log便宜上、多項式はk 2 k^2 k 2アルゴリズムを模写することができ、総複雑度はO(k 2 logn)O(k^2log n)O(k 2 log n)である.
##特徴多項式を解く多項式の型取りの前提は多項式を知ることである.明らかにk k k個の特徴値を解くのは現実的ではなく、常係数の整列繰返しの特徴を考慮して、行列の特徴多項式を迅速に求めることができる:まず、転移方程式をリストする:
( a 1 a 2 a 3 ⋯ a k − 2 a k − 1 a k 1 0 0 ⋯ 0 0 0 0 1 0 ⋯ 0 0 0 0 0 1 ⋯ 0 0 0 ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ 0 0 0 ⋯ 1 0 0 0 0 0 ⋯ 0 1 0 ) ( h k − 1 h k − 2 h k − 3 ⋮ h 2 h 1 h 0 ) = ( h k h k − 1 h k − 2 ⋮ h 3 h 2 h 1 )\begin{pmatrix} a_1 & a_2 & a_3 &\cdots & a_{k - 2} & a_{k - 1} & a_k\\1 & 0 & 0 &\cdots & 0 & 0 & 0\\0 & 1 & 0 &\cdots & 0 & 0 & 0\\0 & 0 & 1 &\cdots & 0 & 0 & 0\\\vdots &\vdots &\vdots &\ddots &\vdots &\vdots &\vdots\\0 & 0 & 0 &\cdots & 1 & 0 & 0\\0 & 0& 0 &\cdots & 0 & 1 & 0\end{pmatrix}\begin{pmatrix}h_{k - 1}\\h_{k - 2}\\h_{k - 3}\\\vdots\\h_2\\h_1\\h_0\end{pmatrix} =\begin{pmatrix}h_{k}\\h_{k - 1}\\h_{k - 2}\\\vdots\\h_3\\h_2\\h_1\end{pmatrix} ⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛a1100⋮00a2010⋮00a3001⋮00⋯⋯⋯⋯⋱⋯⋯ak−2000⋮10ak−1000⋮01ak000⋮00⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛hk−1hk−2hk−3⋮h2h1h0⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛hkhk−1hk−2⋮h3h2h1⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
フィーチャー多項式:f(λ ) = ∣ λ I − A ∣ = ( λ − a 1 − a 2 − a 3 ⋯ − a k − 2 − a k − 1 − a k − 1 λ 0 ⋯ 0 0 0 0 − 1 λ ⋯ 0 0 0 0 0 − 1 ⋯ 0 0 0 ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ ⋮ 0 0 0 ⋯ − 1 λ 0 0 0 0 ⋯ 0 − 1 λ ) f(\lambda) = |\lambda I - A| =\begin{pmatrix}\lambda - a_1 & -a_2 & -a_3 &\cdots & -a_{k - 2} & -a_{k - 1} & -a_k\\-1 &\lambda & 0 &\cdots & 0 & 0 & 0\\0 & -1 &\lambda &\cdots & 0 & 0 & 0\\0 & 0 & -1 &\cdots & 0 & 0 & 0\\\vdots &\vdots &\vdots &\ddots &\vdots &\vdots &\vdots\\0 & 0 & 0 &\cdots & -1 &\lambda & 0\\0 & 0& 0 &\cdots & 0 & -1 &\lambda\end{pmatrix} f(λ)=∣λI−A∣=⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎛λ−a1−100⋮00−a2λ−10⋮00−a30λ−1⋮00⋯⋯⋯⋯⋱⋯⋯−ak−2000⋮−10−ak−1000⋮λ−1−ak000⋮0λ⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎞
最初の行の展開に対してそれぞれ和を求めることを考慮して、あります:
f ( λ ) = ( λ − a 1 ) A 11 + ( − a 2 ) A 12 + ⋯ + ( − a k ) A 1 n = λ k − a 1 λ k − 1 − a 2 λ k − 2 − ⋯ − a k f(\lambda) = (\lambda - a_1)A_{11} + (-a_2)A_{12} +\cdots + (-a_k)A_{1n} =\lambda ^ k - a_1\lambda ^ {k - 1} - a_2\lambda ^ {k - 2} -\cdots - a_k f(λ)=(λ−a1)A11+(−a2)A12+⋯+(−ak)A1n=λk−a1λk−1−a2λk−2−⋯−ak
そのうちA 1 i A_{1 i}A 1 iは、A A Aの代数余子式である.
これにより,定数整列線形繰返しも行列の高速べき乗を最適化できるようになった.
##小最適化1.結果を計算するのは面倒だと思いますが、もう一つ乗算しなければなりませんか?いいえ!A n−k A^{n−k}An−kではなく、A n−A^{n−k}An−kを直接計算すると、最後の結果はベクトルの最後のビットになる.このときA i(i≦k)A^{i}(ile k)Ai(i≦k)の最下位は原繰返し式の1つであり,直接求めることができる.
2.注意k 2 k^2 k 2は、i,j i,j i,j iの順で、i+j i+j i+j i+jに、i,j i,j i,j i,j,a i=a j∗a−j a_ではなく、i+j i+j i+j i+j i+jに加算される.i=a_j*a_{i-j} ai=aj∗ai−j.
BZOJ 4161のコードを添付します.これは定数係数が順次線形に繰り返される裸の問題です.k 2 logn k^2log n k 2 lognの方法で过ごすことができます.
#include
using namespace std;
inline int read(){
char ch=getchar();int i=0,f=1;
while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
while(isdigit(ch)){i=(i<<1)+(i<<3)+ch-'0';ch=getchar();}
return i*f;
}
const int Maxn=2e3+20,Mod=1000000007;
#define rg register
int n,k,a[Maxn<<1],b[Maxn<<1],c[Maxn<<1],h[Maxn],tp[Maxn<<1],ans,lim;
inline void add(int &x,int y){
x+=y;(x>Mod?x-=Mod:(x<0?x+=Mod:0));
}
inline void mul(int *A,int *B){
for(rg int i=0;i<=lim;i++)tp[i]=0;
for(rg int i=0;i<=k;i++)
for(rg int j=0;j<=k;j++)
tp[i+j]=(tp[i+j]+1ll*A[i]*B[j])%Mod;
for(rg int i=lim;i>=k;i--){
for(rg int j=0;j<k;j++)
tp[i-k+j]=(tp[i-k+j]+1ll*tp[i]*a[k-j])%Mod;
tp[i]=0;
}
for(int i=0;i<=lim;i++)A[i]=tp[i];
}
inline void power(int *A,int B,int *C){
for(;B;B>>=1,mul(A,A))if(B&1)mul(C,A);
}
int main(){
n=read(),k=read();lim=k<<1;
for(int i=1;i<=k;i++)a[i]=read();
for(int i=1;i<=k;i++)h[i]=(read()%Mod+Mod)%Mod;
if(n<k){printf("%d
",h[n+1]);return 0;}
b[0]=1;c[1]=1;power(c,n,b);
for(int i=0;i<k;i++)add(ans,1ll*b[i]*h[i+1]%Mod);
printf("%d
",ans);
}