休憩室

科学・技術・芸術

KaTeX数式+はてな+Markdown

せっかくBoostnoteで書いた記事をはてなブログに載せるときにレンダリングエラーが出まくって,挙句の果てには慎重に導入したつもりの数学的記法を変更するのが苦痛になってきたので,今後の運用方針を見直そうと思う.

具体的には

  • Boostnoteを用いない: 便利なだけに残念だが,致し方あるまい
  • (上に関連して)KaTeXをなるべく本番通りにレンダリングしてくれる環境を使う
  • よほどの理由がない限りインライン数式だけを用いる

あたりだろうか.ローカルでのKaTeXのレンダリングについては,次が分かりやすそうである.

7shi.hateblo.jp

これに関連して,七誌さんは次のような記事も書いている.参考にさせてもらいます.

7shi.hateblo.jp

これまでに公開した記事について

とりあえず,上記の方法でブログ全体の仕様を変更したので,ほとんどの記事でレンダリングが崩れています.そのうち修正します.

Web上で数式混じりの文章を書くことについて

ちょっとした引用程度なら手元のLaTeXitかなんかで画像化したものを貼ればいいと思うんですよ.

それを敢えてブラウザやサーバーサイドでレンダリングしようってのは,わざわざOverleafとかを使ってPDF化するほどでもない文章を,ずっと使い慣れている(僕みたいな)ほぼLaTeXネイティブの人間が公開するための方法なわけです.だから,サービスの説明に「LaTeX数式をサポートしています」みたいな文言が書かれていれば初手では信じますし,後から「この記法は使えない」みたいなことが発覚すると詐欺に遭った気分になります.

結城浩さんもBoostnoteについて何か書いてました.

"Latex Support" という部分に少し興味あり試してみる / “Boostnote” http://htn.to/QU4wh7 - 結城浩の連ツイ

本当に致命的なのは,MathjaxやKaTeXをサポートする(と謳う)様々なツールの間で記法・変換規則のコンセンサスが一切取れていないことだと思いますし,誰か統一したガイドラインを作ってくれないと時間を無駄にするばかりです.こういう時は普通,本家LaTeXの仕様とどう異なるかを明記するものなんじゃないですかね.

統計ノート: 確率モデルと識別関数(1)

確率・統計の基礎について,平井有三『はじめてのパターン認識』の第4章がまとまっていて分かりやすかったので,適度に情報を補いつつ自分で整理してみた.ちょくちょく

統計WEB - 統計学、調べる、学べる、BellCurve(ベルカーブ)

といったホームページを参照したので,まだ自分でもよく分かっていなかったんだなと実感.

本記事では,統計の基礎ということで,統計量・推定量・不変性・一致性などの概念を押さえていく.

以下,特に断りの無い場合は

  •  \boldsymbol{a},\boldsymbol{b},\ldots といった太字の小文字アルファベットはベクトル
  •  \boldsymbol{A},\boldsymbol{B},\ldots といった太字の大文字アルファベットは行列

を表し, \mathbb{R}^{d} と書いて  d 次元の実Euclid空間と解釈するものとする.

統計量と推定量

まず初めに,何らかの観測データ  \boldsymbol{x}\in\mathbb{R}^{d} が真の分布  p(\boldsymbol{x});\;p:\mathbb{R}^{d}\to[0,1] に従って生成されていたとする.実際に生成された観測データのことをここでは標本と呼び, \{\boldsymbol{x} _ {n}\} _ {n = 1}^{N} と表記する.まず,個別のサンプル  \boldsymbol{x}\in\{\boldsymbol{x} _ {n}\} _ {n = 1}^{N} と真の分布  p(\boldsymbol{x}) に対して,次の統計量を定義する.

  • 平均ベクトル   \boldsymbol{\mu}:=\int_{\mathbb{R}^{d}} d\boldsymbol{x}\;p(\boldsymbol{x})\boldsymbol{x},
  • 共分散行列   \boldsymbol{\Sigma}:=\int_{\mathbb{R}^{d}} d\boldsymbol{x}\;p(\boldsymbol{x})(\boldsymbol{x}-\boldsymbol{\mu})(\boldsymbol{x}-\boldsymbol{\mu})^{T}.

これらはまさしく平均値と平均値からのズレを表す基本的な統計量であるが,真の分布  p(\boldsymbol{x})積分核とする真の統計量であり,実際に計算できることは珍しい.多くの場合,我々が知りたいのは真の分布それ自身だからである.そこで,標本  \{\boldsymbol{x}_{n}\}_{n = 1}^{N} だけから

  • 標本平均ベクトル  \boldsymbol{\mu}_{N}:=\frac{1}{N}\sum_{n = 1}^{N}\boldsymbol{x}_{n},
  • 標本共分散行列  \boldsymbol{\Sigma}_{N}:=\frac{1}{N}\sum_{n = 1}^{N}(\boldsymbol{x}_{n}-\boldsymbol{\mu}_{N})(\boldsymbol{x}_{n}-\boldsymbol{\mu}_{N})^{T}.

も定義しておく.これらはそれぞれ,平均ベクトルと共分散行列の定量と呼ばれ,真の統計量とは区別される1

定量の不偏性と一致性

 (\boldsymbol{\mu}_{N}, \boldsymbol{\Sigma}_{N}) (\boldsymbol{\mu}, \boldsymbol{\Sigma}) の推定量であるとはどういうことかを考えたい.そこで,議論の見通しを良くするために,スカラー・ベクトル・行列を問わず観測データ  \boldsymbol{x} に依存する量  \hat{f}(\boldsymbol{x}) についての期待値標本平均操作的に次で定義する.

 E_{p}\{\hat{f}\} := \int d\boldsymbol{x} \;p(\boldsymbol{x})\hat{f}(\boldsymbol{x}),

 E_{N}\{\hat{f}\} := \frac{1}{N} \sum_{n = 1}^{N} \hat{f}(\boldsymbol{x}_{n}).

これを用いると, \boldsymbol{\mu}=E_{p}\{\boldsymbol{x}\} \boldsymbol{\mu}_{N}=E_{N}\{\boldsymbol{x}\} といった表現が可能になる.次に, N 個の標本の各々が真の分布から並列に生成されている状況を考え2,標本平均ベクトル・標本共分散行列のそれぞれの平均がどう振る舞うかを調べる.まず.標本平均ベクトルの期待値は

 E_{p}\{\boldsymbol{\mu}_{N}\} = \frac{1}{N}\sum_{n=1}^{N}E_{p}\{\boldsymbol{x}_{n}\} = \boldsymbol{\mu}

となり,自分自身の真の統計量に一致する.この性質は不偏性は呼ばれ.標本数  N に依存しないのが特徴である.一方,標本共分散行列の期待値については

 \boldsymbol{\Sigma}_{N} = E_{N}\{(\boldsymbol{x} - \boldsymbol{\mu}_{N})(\boldsymbol{x} - \boldsymbol{\mu}_{N})^{T}\} = E_{N}\{(\boldsymbol{x} - \boldsymbol{\mu})(\boldsymbol{x} - \boldsymbol{\mu})^{T}\} - E_{N}\{(\boldsymbol{\mu}_{N} - \boldsymbol{\mu})(\boldsymbol{\mu}_{N} - \boldsymbol{\mu})^{T}\}

という変形を事前に行うと

 E_{p}\{(\boldsymbol{x} - \boldsymbol{\mu})(\boldsymbol{x} - \boldsymbol{\mu})^{T}\} = \boldsymbol{\Sigma},

 E_{p}\{(\boldsymbol{\mu}_{N} - \boldsymbol{\mu})(\boldsymbol{\mu}_{N} - \boldsymbol{\mu})^{T}\} = \frac{1}{N} \boldsymbol{\Sigma}

が言えるので,結果的に

 E_{p}\{\boldsymbol{\Sigma}_{N}\} = E_{N}\left\{\boldsymbol{\Sigma} - \frac{1}{N} \boldsymbol{\Sigma}\right\} = \frac{N - 1}{N}\times E_{N}\left\{\boldsymbol{\Sigma}\right\} = \frac{N - 1}{N}\times\boldsymbol{\Sigma}

が成り立つ.すなわち,標本共分散行列は不偏性を満たさない推定量である.そこで,標本共分散行列の定義を

  \boldsymbol{\Sigma}_{N} \to \frac{N}{N - 1} \boldsymbol{\Sigma}_{N}

と変更し,不偏性を満たすようにすることもできる(不偏推定量).一方で,後で見るように,Gauss分布を仮定した最尤推定では元の定義のままの方が良い. ここで,標本数を無限大にする極限  N\to\infty を考えると,明らかに

 \lim_{N\to\infty}E_{p}\{\boldsymbol{\mu}_{N}\} = \boldsymbol{\mu},

 \lim_{N\to\infty}E_{p}\{\boldsymbol{\Sigma}_{N}\} = \boldsymbol{\Sigma}

が満たされていることが分かる.このような,標本数を増やした時に真の統計量に一致する性質のことを一致性と呼ぶ.

不偏推定量であればその統計量は明らかに一致性を持つが,逆が成り立つとは限らない.

補足: 式変形について

 E_{p}\{(\boldsymbol{\mu}_{N} - \boldsymbol{\mu})(\boldsymbol{\mu}_{N} - \boldsymbol{\mu})^{T}\} = \boldsymbol{\Sigma} / N が少し非自明かも知れないが,基本的には項別積分 E_{p}\{\bullet\} の持つ意味に注意すれば良くて,多変数関数の期待値を常に

 E_{p}\{F(\boldsymbol{x}_{1}, \boldsymbol{x}_{2}, \ldots, \boldsymbol{x}_{N})\} = \int\cdots\int_{\mathbb{R}^{d}}\prod_{n = 1}^{N}\left\{d\boldsymbol{x}_{n}\;p(\boldsymbol{x}_{n})\right\}F(\boldsymbol{x}_{1}, \boldsymbol{x}_{2}, \ldots, \boldsymbol{x}_{N})

という  (d\times N) 次元の多重積分で解釈すれば矛盾なく証明できる.


  1. 逆に,真の統計量である平均ベクトル・共分散行列のことを母平均母分散と呼ぶこともある

  2. このような仮定のことを独立同分布(independently and identically distributed; IID)という.

数値解析ノート: 連立1次方程式の定常反復解法(2). Gauss-Seidel法

はじめに

本記事では,$\mathbf{A}\mathbf{x} = \mathbf{b}$ 型の連立一次方程式を解くための反復法であるGauss-Seidel法の導出と,その性質について述べる.

Gauss-Seidel法は,Jacobi法と並んで定常反復法という大きな枠組みに分類される方法である.ここで,定常反復法というのは,方程式 $\mathbf{A}\mathbf{x} = \mathbf{b}$ を満たす解の候補としてある初期ベクトル $\mathbf{x} = \mathbf{x}^{(0)}$ を考え,ある反復行列 $\mathbf{B}$ および反復ベクトル $\mathbf{c}$ を用いて $$\mathbf{x}^{(k + 1)} = \mathbf{B} \mathbf{x}^{(k)} + \mathbf{c}; \;k = 0, 1, 2, \ldots$$ という漸化式にしたがって解ベクトル $\mathbf{x}$ の更新を行い,真の解に収束させていく方法であった.「定常反復法」というネーミングは,更新に用いられる $(\mathbf{B}, \mathbf{c})$ が変化しないことを指している1. 定常反復法に関する一般的については,前回の記事で取り扱ったのでそちらも適当に参照して欲しい.

数値解析ノート: 連立1次方程式の定常反復解法(1). Jacobiの反復法 - 休憩室

Jacobi法は,反復行列と反復ベクトルをそれぞれ $$\mathbf{B} = - \mathbf{D}^{-1} (\mathbf{L} + \mathbf{U}), \;\mathbf{c} = \mathbf{D}^{-1} \mathbf{b}$$ ととった場合に相当するが,Gauss-Seidel法ではこれが $$\mathbf{B} = - (\mathbf{D} + \mathbf{L})^{-1} \mathbf{U}, \;\mathbf{c} = (\mathbf{D} + \mathbf{L})^{-1}\mathbf{b}$$ に変更される.


定常反復法としてのGauss-Seidel法

Jacobi法は,更新したいベクトル $\mathbf{x}^{(k)}$ の各成分 ${x^{(k)}_{j}}_{j = 1, 2, \ldots, N}$ が,$\mathbf{x}^{(k)}$ から $\mathbf{x}^{(k + 1)}$ への更新の途中で互いに干渉しないという性質のおかげで,素直な並列化ができた.この辺に関しては,実装と数値実験の結果を示した次の記事を参照して欲しい.

数値解析ノート1. Jacobiの反復法 - Qiita

一方で,Gauss-Seidel法はこの並列化を見かけ上犠牲にする代わりに,ほんの少し収束性を良くする仕組みを備えている2.また,収束性のみならず実装が自然になるというメリットもある.この辺りの実用上の事情は次の記事に記しておいた. qiita.com


Gauss-Seidel法のアルゴリズム

Gauss-Seidel法は,成分表示で見ると

$$x^{(k + 1)}_{i} = (b_{i} - \sum_{j < i} a_{i, j} x^{(k + 1)}_{j} - \sum_{j > i} a_{i, j} x^{(k)}_{j}) / a_{i, i};\; k = 0, 1, 2, \ldots$$

という更新規則に基づくアルゴリズムである.一方,Jacobi法の更新規則は

$$x^{(k + 1)}_{i} = (b_{i} - \sum_{j < i} a_{i, j} x^{(k)}_{j} - \sum_{j > i} a_{i, j} x^{(k)}_{j}) / a_{i, i};\; k = 0, 1, 2, \ldots$$

で与えられており,これらはよく似ている.違いは更新したばかりのベクトルの成分を即座に用いるかどうかである.

このことを更に一般化すると,何らかの順列 ${p_{j}}_{j = 1, 2, \ldots, N}$ が与えられており,この順序に従って

$$x^{(k + 1)}_{p_{i}} = (b_{p_{i}} - \sum_{p_{j} < p_{i}} a_{p_{i}, p_{j}} x^{(k)}_{p_{j}} - \sum_{p_{j} > p_{i}} a_{p_{i}, p_{j}} x^{(k)}_{p_{j}}) / a_{p_{i}, p_{i}};\; k = 0, 1, 2, \ldots$$

という更新を施すのもまたGauss−Seidel法の資格を持つ.例えば本来の(素朴な)Gauss-Seidel法を逆順にした

$$x^{(k + 1)}_{i} = (b_{i} - \sum_{j < i} a_{i, j} x^{(k + 1)}_{j} - \sum_{j > i} a_{i, j} x^{(k)}_{j}) / a_{i, i};\; k = 0, 1, 2, \ldots$$

という更新規則はその例である.

ちなみに,Gauss-Seidel法は逐次的過緩和(SOR)法という比較的有名なアルゴリズムを導出する上で重要である.


  1. 一方で,共役勾配法(CG法)などに代表される非定常反復法では,これらの反復行列・反復ベクトルにもう少し複雑な情報をもたせることで収束性を良くする.一方,そのような場合には収束性に関する数学的な考察が難しくなる.

  2. この,見かけ上犠牲というのが微妙なステートメントになってしまったが,要するにJacobi法のようなstraightforwardな並列化はできないよ,という意味でしかない.実際,順序付けに関する工夫というのを行えば並列性が回復できるらしいが,その順序付けに関しては筆者も未だによく理解していないので,とりあえず大阪電通大の伊藤祥司氏の学位論文へのリンクを貼っておく: 大規模線形方程式に対するスーパーコンピュータ向きの前処理法 - 国立国会図書館デジタルコレクション(直でダウンロードはできないが,章毎にGoogle検索でPDFが引っかかり,その中の一つに赤黒順序付けという最も簡単な順序付けに関する解説がある)

数値解析ノート: 連立1次方程式の定常反復解法(1). Jacobiの反復法

はじめに

最近,朝倉書店の応用数値計算ライブラリというシリーズの一冊である『反復法Templates』に目を通した1

広い意味での数値解析に関する書物は,森正武・他『数値解析』(共立出版)や山本哲郎『数値解析入門』(サイエンス社)といった古典的名著に始まり,最近では齊藤宣一『数値解析入門』(東京大学出版会)のような画期的な教科書も出ているが,これらはいずれも,数値計算一般(偏微分方程式の解法・行列の対角化・etc.)に関する理論的洞察を与えてくれる.

一方で,反復法Templatesはその名の通り

$$ \mathbf{A} \mathbf{x} = \mathbf{b} $$

という形の連立方程式―適当な係数行列$ \mathbf{A} $および係数ベクトル$ \mathbf{b} $を与えると解ベクトル$ \mathbf{x} $を返すような問題―における反復法に話題を限定した数値解析の教科書だが,これがとてもいい本なのである.ここで,反復法というのは

あるベクトルの列 $ {\mathbf{x}_{n}}_{n=0}^{\infty} $($ \mathbf{x}_{0} $は初期値)を考え,その極限値 $ \mathbf{x}_{\infty}:=\lim_{n\to\infty}\mathbf{x}_{n} $ が真の解に一致するように,$( \mathbf{A}, \mathbf{b} )$ を用いて反復演算を行う

という手法の総称である.教科書なんかだと「LU分解に基づく直接法,そして反復法の2つがある」という感じで触れられるやつである.私はこれまで反復法をきちんと勉強したことがなく,個々の手法を知っていても体系的な理解には至らなかったが,この本を読んだことで随分と頭の中が整理された気がする2

というわけで,せっかく勉強したことを再度整理し直すことを目標に,本記事では

  • Jacobiの反復法

を反復法の観点から触れ,それに続く2つの記事で

  • Gauss-Seidelの反復法
  • 逐次的過緩和(SOR)法

を紹介しようと思う.これらは,反復法の中で更に特別な定常反復法というクラスに属し,アルゴリズムの導出も簡単,収束解析もきれいにできるという利点がある.この辺の一般論も以下である程度じっくり触れられればと思う.

なお,実機上での数値実験については,本稿から切り離して別の記事で扱うことにしたので,そちらも確認してみて欲しい.

qiita.com

以下,アルゴリズムの導出とちょっとした収束解析の話をまとめる.


Jacobiの反復法のアルゴリズム

Jacobi法は,$\mathbf{A}\mathbf{x} = \mathbf{b}$ 型の方程式に対するおそらくもっとも簡単な定常反復法である.まず,方程式を $$ \sum_{j = 1}^{n} a_{i, j}x_{j} = b_{i};\; i = 1, 2, \ldots, n \qquad\text{(1)} $$ と表記し,これを変形して $$ a_{i, i} x_{i} + \sum_{j \neq i}a_{i, j}x_{j} = b_{i} $$ $$ \Leftrightarrow x_{i} = (b_{i} - \sum_{j \neq i} a_{i, j} x_{j}) / a_{i, i} \qquad\text{(2)} $$ としたものがJacobi法の出発点である.

ここで,さきほど触れた「ベクトル列 ${\mathbf{x}_{\bullet}}$ が真の解ベクトルに収束する云々」という旨のことを思い出すと,方程式 $(2)$ を次のように適度に弱めても,(一定の条件下で)解ベクトルが収束すれば方程式 $(1)$ がきちんと満たされそうである3. $$x^{(k)}_{i} = (b_{i} - \sum_{j \neq i} a_{i, j} x^{(k - 1)}_{j}) / a_{i, i};\; k = 0, 1, 2, \ldots \qquad \text{(3)}$$ これが,Jacobi法の更新規則である.Jacobi法は異なる $i, j;\; i\neq j$ に対する $x^{(k)}_{i}$ と $x^{(k)}_{j}$ の更新が互いに依存しないので,一斉置き換え法(simultanious replacing method)とも呼ばれる4. Jacobi法に限らないことだが,式 $(3)$ のタイプの更新規則が定常反復法と呼ばれるのには理由がある.それは,このタイプの更新規則はいずれも $$\mathbf{x}^{(k)} = \mathbf{B} \mathbf{x}^{(k - 1)} + \mathbf{c}; \; k = 0, 1, 2, \ldots \qquad \text{(4)}$$ という格好で書けるからである.ここで登場した行列 $\mathbf{B}$ ・ベクトル $\mathbf{c}$ はそれぞれ反復行列反復ベクトルと呼ばれ,解ベクトルの情報には依存しないのがポイントである.特に反復行列 $\mathbf{B}$ を対角化したときの固有値分布からは,対応する定常反復法の収束に関する情報も得られて,数学的に性質の良いアルゴリズムの源泉になっている.

すなわち,定常反復法においては,反復行列 $\mathbf{B}$・反復ベクトル $\mathbf{c}$ の選び方がそのままアルゴリズムの選択になる.パフォーマンスの観点からはありえないが,こうやって定式化しておくと実際に $\mathbf{B}$・$\mathbf{c}$ をそのまま*axpy系の行列計算ルーチンで回すこともできる5.例えば,Jacobi法は $$\mathbf{B} = - \mathbf{D}^{-1} (\mathbf{L} + \mathbf{U}), \;\mathbf{c} = \mathbf{D}^{-1} \mathbf{b}$$ の場合に相当する.ただし,$\mathbf{D}, \;\mathbf{L}, \;\mathbf{U}$ はそれぞれ行列 $\mathbf{A}$ の対角成分,真下三角成分,そして真上三角成分である: $$\mathbf{A} = \mathbf{D} + \mathbf{L} + \mathbf{U}.$$


収束解析: 1次元Poisson方程式の場合

折角なので例を考える.例えばソース項 $f(x)$ を伴う,1次元領域 $(0, 1)\in \mathbb{R}$ 上の微分方程式 $$\mathcal{L} u(x) = - \partial_{x}^{2} u(x) = f(x);\; \mathcal{L} := - \partial_{x}^{2}$$ を考えると,これを $N$ 個の点で離散化した差分方程式として $$\hat{L}u(k / N) = u((k + 1) / N) + u((k - 1) / N) - 2 u(k / N) = N^{2}f(k / N) \qquad \text{(5)}$$ が得られる.ここで,$\hat{L}$ は $\mathcal{L}$ に対応する差分演算子であり,今の場合はそっくりそのまま係数行列 $\mathbf{A}$ を与える.まず,$u(x)$ の候補として $\hat{L}$ の固有関数 $$u_{n}(x) := \sin(n\pi x);\; n \in \mathbb{Z}$$ をとってくると,各 $n$ に対して固有値 $$\lambda_{n} := 4 \left(\sin\left[\frac{n\pi}{2N}\right]\right)^{2}$$ が付随する.ちなみに,関数 $u_{n}(x)$ は元の微分演算子 $\mathcal{L}$ の固有関数にもなっている.

固有値 $\lambda_{n}$ が固有関数 $u_{n}(x)$ の周波数成分を担うことに注意すると,次のことが言える.まず,差分方程式 $(5)$ に対応する係数行列 $\mathbf{A}$ は $$\mathbf{A} = \begin{pmatrix}2 & -1 & 0 & \cdots & \\ -1 & 2 & -1 & \cdots & \\ 0 & -1 & 2 & & \\ & & \cdots & 2 & -1 \\ & \cdots & 0 & -1 & 2\end{pmatrix}$$ という格好になっているはずである.これから直ちに反復行列 $\mathbf{B}$ を構成できる: $$\mathbf{B} = - \frac{1}{2} \begin{pmatrix}0 & -1 & 0 & \cdots & \\ -1 & 0 & -1 & \cdots & \\ 0 & -1 & 0 & & \\ & & \cdots & 0 & -1 \\ & \cdots & 0 & -1 & 0\end{pmatrix}.$$ ここで,反復行列 $\mathbf{B}$ の頭にかかる因子 $- 1 / 2$ は,係数行列 $\mathbf{A}$ の対角成分 $\mathbf{D}$ に,ゼロ埋めされた対角成分は上三角成分と下三角成分の和 $\mathbf{L} + \mathbf{U}$ にそれぞれ由来する.このことから,反復行列は $$\mathbf{B} = \mathbf{E} - \frac{1}{2} \mathbf{A}$$ とも書ける.すると,反復行列 $\mathbf{B}$ の(固有関数 $u_{n}(x)$ に付随する)固有値を係数行列 $\mathbf{A}$ の固有値を用いて表現することが許される: $$1 - \frac{1}{2} \lambda_{n} = 1 - 2 \left(\sin\left[\frac{n\pi}{2N}\right]\right)^{2}=:\lambda_{n}^{(B)}.$$ 反復行列 $\mathbf{B}$ の固有値 $\lambda_{n}^{(B)}$ は,区間 $0 \leq n < N$ で $n$ について単調減少なので,この区間では1次方程式の初期ベクトル $\mathbf{x}^{(0)}$ のうち大きな $n$ に対応するものから順に速く減衰するということが分かる.特に,十分小さな $n$ に対しては $$\lambda_{n} \sim 1 - \mathrm{const.} / N^{2}$$ が成り立ち,この $N$ に対する漸近挙動は今考えている1次元系の2階偏微分方程式のみならず,より高次元の系(2次元・3次元)においても普遍的に観測される.このことには,係数行列 $\mathbf{A}$ がバンド行列であるという事実が深く関わっている.ここで,バンド行列というのは

対角成分から非対角方向に次元 $N$ によらない定数の距離までしかノンゼロ成分を持たない行列

のことである.まさに今の場合は副対角成分のみがノンゼロなので,この資格を満たしている.

そもそも,何故バンド行列が出現するのかにもちょっと興味があるので,それに触れて本稿の締めくくりとする.少なくとも私は,ベースにある物理法則(支配方程式)が局所的にしか効かず,数学的にはそれが有限階の偏微分として表現されるのが全てのルーツであると考えている.よって,今の場合は「偏微分方程式の差分化」が行列の出処であったが,背景に同じような物理法則がある限り,差分法の代わりに境界要素法や有限要素法を考えたとしても,等しくバンド行列が出現するという風に理解している6. 以上のJacobi法を実装し,実際にコンパイル・実行可能なパッケージとして公開したので,良ければ参考にしていただきたい.

github.com


  1. 当該書物は既に絶版となっているが,研究所や大学等で過ごすことの多い方はぜひ手にとって眺めてみて欲しい.

  2. 今は地味に日本計算工学会の計算力学レクチャーコースというシリーズにも興味があり,頭の中が整理されているうちに手を出してみたいと思っている.

  3. この,「一定の条件下」というのがとても厄介で,実際の実行結果からも分かる通り,実は行列が優対角でないと発散してしまう.

  4. したがって,並列化も容易であると期待される.実際に並列化した結果についても,記事(まだ)に乗せておいた.

  5. もちろん,式 $(3)$ をそのままコードに翻訳するのが良い.システムのモジュールは切り離して可搬性を保つほど性能は落ちるという一般的なトレードオフの例である.どの程度切り離すかについては,経験則が支配しているように思える.ちなみに,とある計算科学系のハンズオン研究会で誰かが「ソフトウェア・エンジニアリングは,モジュール間のエンタングルメントエントロピーが最小となるように行われるべきである」と仰っていた記憶があり,言い得て妙な例えだと思った.

  6. 実は,同一の偏微分方程式のモデル化を異なる手法で行った時に,解くべき連立1次方程式が果たしてどこまで等価になるか,という問題についてはきちんと理解していない.ひょっとすると定数倍の不定性を覗いてexactに一致するのかもしれないし,システムサイズ $N$ についての漸近形が同じになるだけかも知れない.この辺りのことをご存じの方がいたらぜひ教えていただきたいです.とは言っても,実際に確かめるのはそんなに苦ではないかも知れない.

初等整数論ノート: 「剰余逆元と拡張Euclid互除法」についての補足

はじめに

以前投稿した初等整数論ノート: 剰余逆元と拡張Euclid互除法 - 休憩室の記事中のコードがちょっと微妙すぎたので,再投稿のつもりで補足をします.

詳細な問題背景については,全て当該記事に記載してありますので,そちらを参照して下さい.

trtn.hatenablog.com

具体的には,計算後に適当な整数(あるいは,ものによっては素数) $ m $ で割った余り(剰余)を算出するという設定で,

  • 四則演算
  • 累乗・二項係数
  • 階乗

を求める関数群をこちらの記事で提供しますので,ご興味があれば自由にお使い下さい.

前提条件とサブモジュール

言語は主にC++11以降です.C++11以前のC++C言語を用いる際には,上手くデータ構造を補って下さい.

また,プログラム冒頭(例えば #include <hogehoge> の直後)で

typedef long long ll;

という定義を行っている状態を仮定します.

そして,前回記事でも紹介した不完全Euclid互除法剰余逆元(の ll 版)についても,サブモジュールとして使える状態にもしておきます1,2:

std::pair<ll, ll> iegcd(const ll a, const ll b, const ll x = 1, const ll u = 0)
{
    if (b == 0)
    {
        return make_pair(a, x);
    }
    else
    {
        return iegcd(b, a % b, u, x - (a / b) * u);
    }
}

ll modinv(const ll a, const ll m = MOD)
{
    std::pair<ll, ll> res = iegcd(a, m);
    return res.second;
}

以上で,本記事中で紹介するライブラリは,g++などの標準的なコンパイラと実行環境であれば全て使えるようになるはずです.

剰余系ライブラリ1: 四則演算

四則演算のうち,途中でオーバーフローを起こす可能性のある計算は掛け算 $ \times $ であり,割り算 $ / $ は掛け算を使って定義されます.一方,剰余を取るための具体的な $ m $ の値は加減算 $ +,\; - $ でも必要なので,いっそ四則演算の全てを関数化してしまおうというのが今回の概略です.多少大げさですが,少なくともこうすることで

  • 間違いが少なくなる
  • 専用の整数クラスを作るのに便利

という恩恵を享受できるでしょう3. というわけで,一つ一つ紹介します.まずは,整数 $ a $ と整数 $ b $ の和の剰余(with $ m $)を返す関数です.流石に何の変哲もありません.

ll modpls(const ll a, const ll b, const ll m)
{
    return (a + b) % m;
}

次に,$ a $ と $ b $ の差の剰余(with $ m $)を返す関数.なぜ条件分岐をしているかというと,返り値を自然数全体 $ \mathbb{N} $ に制限したいからです.実際,入力値 $ a, \;b $ の範囲を自然数全体 $ \mathbb{N} $ に制限しただけでは,このことは保証されません.

ll modmin(const ll a, const ll b, const ll m)
{
    ll res = (a - b) % m;
    if (res < 0)
    {
        return res + m;
    }
    else
    {
        return res;
    }
}

そして,積の剰余.積の剰余について成り立つ一般的な関係を使っただけのコードです.

ll modmult(const ll a, const ll b, const ll m)
{
    return ((a % m) * (b % m)) % m;
}

最後に,商の剰余.前回記事とは打って変わって,積の剰余の実装に依存させてみました.もちろん,剰余逆元を用いています.

ll moddiv(const ll a, const ll b, const ll m)
{
    return modmult(a, modinv(b, m), m);
}

以上だけで,まず $ a, \;b $ を入力値とする四則演算の $ m $ による剰余が自由に計算できるようになりました.これを用いて,累乗や二項係数の計算が自由に行えるようになります.

剰余系ライブラリ2: 累乗・二項係数

まず,累乗.累乗は数学的には $ a^{n} $ の形で与えられるので,

$$ a^{n} \% m = ( (a^{n - 1} \% m) \times (a \% m) ) \% m = \cdots $$

という再帰的な分離を行っていけば良さそうですが,これを真面目に実装すると $ \mathcal{O}(n) $ のアルゴリズムになります4.ですが,これは工夫によって $ \mathcal{O}(\log(n)) $ に改善することができる計算です.そこで,繰り返し二乗法と呼ばれるアルゴリズムによって系統的な分離を行うことにします5.実装のミソとしては

  1. $ n_{0} := n = 2^{k_{1}} \times n_{1} $ の形にする
  2. $ n_{1} $ が奇数なら $ n_{1} = n_{2} + 1 $ と分離する
  3. $ n_{2} = 2^{k_{2}} \times n_{3} $ の形にする
  4. 以下略

みたいな感じです.この手続きによって必ず $ n_{\bullet}=1 $ に到達することが保証されるので,そこで改めて $ a \% m $ を評価し,その結果を上流に報告するようなイメージです.再帰関数の構造をフルに活かしたアルゴリズムになっています.

ll modpow(const ll a, const ll n, const ll m)
{
    if (n == 1)
    {
        return a % m;
    }
    else if (n % 2 == 1)
    {
        return modmult(a, modpow(a, n - 1, m), m);
    }
    else
    {
        ll t = modpow(a, n / 2, m);
        return modmult(t, t, m);
    }
}

また,二項係数の計算もある意味で四則演算を活かした高級なアルゴリズムですが,工夫のしどころは多少自明です.これも,前回記事のバージョンに比べてだいぶ計算量が減っています.特に,$ n $ が非常に大きいが $ \min(k, n - k) $ がそこまで大きくないような場合に大幅に速くなります6

ll modcombination(const ll n, const ll k, const ll m)
{
    ll k2 = min(n - k, k);

    ll num = 1;
    for (size_t i = n; i > n - k2; i--)
    {
        num = modmult(num, i, m);
    }

    ll den = 1;
    for (size_t i = k2; i > 1; i--)
    {
        den = modmult(den, i, m);
    }

    return moddiv(num, den, m);
}

剰余系ライブラリ3: 階乗

最後に,階乗を計算する関数の実装でこの記事を終えようと思います.これはほぼ再掲載で,前回記事から変わったポイントはほぼありません.また,$ n $ が非常に大きいときには計算量に目を瞑らねばならないので,標準ライブラリ扱いはちょっと厳しそうです.

ll modfactorial(const ll n, const ll m)
{
    ll res = 1;
    for (ll i = 2; i <= n; i++)
    {
        res = modmult(i, res, m);
    }
    return res;
}

まとめ

この手のアルゴリズムは,一度実装・実験するだけで結構勉強になることが多いです.とりあえず全体をRepl.itに貼り付けたものを公開しておきます.

repl.it

ちなみに,上記リンクの用途はAtCoder Beginner Contest 156D - Bouquetです.


  1. まず整数型 int の壁を取り払っています.こういうのは割と大事.

  2. 前回の記事では触れられませんでしたが,剰余逆元の計算はフェルマーの小定理と呼ばれるものがベースの,別のアルゴリズムを用いてもできます.

  3. 特に,$ (a - b) \% m $ を求める式を (a - b) % m と書いてしまう間違いはこれによって防げるはずです.

  4. 前回の記事の実装がまさにそんな感じで,$ n \sim 10^{10} $ で実行時間が数秒オーダーになり,競プロだと間に合わなくなります.

  5. 繰り返し二乗法については,【競プロ】繰り返し二乗法と再帰関数 | なかけんの数学ノートを大いに参考しました.

  6. 実は,前回の記事のバージョンはこの後で触れる階乗関数に依存した形になっていますが,これが大幅な遅延のもとになります.

初等整数論ノート: 剰余逆元と拡張Euclid互除法

2020/02/25(火)更新: 注意事項

本記事の内容は,数学的には何とかなっているものの,掲載したコードのパフォーマンスが微妙なため,直接利用するのはお勧め致しません.詳しくは,次の記事をご参照下さい.

trtn.hatenablog.com

はじめに

AtCoder Beginner Contest 154 - AtCoder のF問題 - Many Many Pathsの解説動画を見て,自分も剰余系の問題をテキパキ捌けるようになりたいと思ったのがきっかけで,初等整数論っぽい内容の記事を書きました.競技プログラミングに限らず,整数論はそれ自体が数学としても面白く,人を惹き付ける魅力を持っています.この記事が,誰かの何らかのお役に立てば光栄です.

剰余逆元を拡張ユークリッドの互除法を用いて計算する方法

整数が四則演算で閉じるようにするには,適当な自然数 $ m\in\mathbb{N} $ を考え,次のように演算規則を定めれば良い.

  • 加減乗法に関しては,結果を $ m $ で割った剰余(modular)を採用する
  • 除法に関しては,乗法における対応関係から定める

これは,整数環を適当なパラメータ $ m $ のもとに縮小してとして扱っていることに他ならない1.そして,このようなことが常にできるのは実は $ m $ が素数の場合に限られる. 具体的には剰余逆元(modular multipricative inverse)というものを考えて,それとの掛け算を用いて割り算を定義する.ここで,整数 $a$ の(法 $ m $ についての)剰余逆元を,任意の整数 $ a \in \mathbb{Z} $ と自然数 $ m\in \mathbb{N} $ が与えられた時 $$ a x \equiv 1\;(\mathrm{mod}\; m) $$ が成り立つような $ x\in \mathbb{Z} $ と定義する.この式における両辺の掛け算は通常の意味での整数同士の掛け算であるが,$ a $ の剰余逆元が存在する場合には $$ x \equiv a^{-1}\;(\mathrm{mod}\; m) $$ という書き方をよくする,この右辺は額面通りの(実数としての)逆数ではないというのがポイントである2.このような剰余逆元は常に存在するわけではなく,存在するための必要十分条件は「$ a $ と $ m $ が互いに素である」ことである.これを認めると,任意の $ a $ と互いに素であるような $ m $ は素数に限られるので,以下では $ m $ が素数の場合だけを考えていく.

手計算

例えば $ (a, m) = (3, 11) $ の場合を考えてみると3 $$ 3 x \equiv 1\; (\mathrm{mod}\; 11) $$ を満たす $x$ は次のように求まる.$ 3x \equiv 1 \;(\mathrm{mod}\;11) $ ということは適当な整数 $ k \in \mathbb{Z} $ が存在して $$ 3x = 11 \times k + 1 $$ が成立するはずである.そこで,$ k $ の値を当てずっぽうに $ k = 1 $ としてみるとあら不思議, $$ 3x = 11 \times 1 + 1 = 12 \Leftrightarrow x = 4 $$ が成り立つではありませんか.というわけで,こうして発見法的に求めた $x$ の特解を使って,今度は $$ x \mapsto x'_{l} := x + 11 \times l; (l \in\mathbb{Z}) $$ という一般解を構成できる.今は $ \mathrm{mod}\;11 $ の世界で考えているので,これが解空間を包む極大集合であることは明らかである.この一般解を上の定義式に代入すると, $$ 3(x + 11 \times l) = 11 \times k + 1 \Leftrightarrow 3x = 11 \times(k - 3\times l) + 1 $$ となり,これを満たす $ x $ の範囲はやはり変わらないことが分かる($ k $ が $ k - 3\times l $に変わっただけで動ける範囲は同じである).ゆえに,こうして求めた $ x $ の解空間 $ \{ 4 + 11 \times k \mid k \in \mathbb{Z} \} $ はexactに正しいことが分かった.手計算では,こういう感じで剰余逆元を求めるのが普通である.

計算機を使った計算

計算機を使えば,上記の手続きをもう少し系統的に追うことができる.そこで役に立つのが,拡張Euclid互除法と呼ばれるアルゴリズムである.まず,Bézoutの等式と呼ばれるものを認めると,任意の整数のペア $ (m, n)\in \mathbb{Z} \times \mathbb{Z} $ について $$ mx + ny = \mathrm{gcd}(m, n) $$ を満たす解 $ (x, y) \in \mathbb{Z}\times \mathbb{Z} $ が必ず存在することが分かる.剰余逆元の定義が,与えられた $ (a, m) $ について $$ a x \equiv 1\;(\mathrm{mod}\; m) $$ を満たす $ x $ であったことを思い出すと,この条件は変形できて,適当な整数 $ q\in\mathbb{Z} $ に対して $$ ax - 1 = mq \Leftrightarrow ax - mq = 1 $$ を満たす $ x $ という風にも言い換えられる.この新しい定義は,見れば分かる通りBézoutの等式の例になっている.というのも,$ m $ が素数の場合は任意の整数 $ q $ について $ \mathrm{gcd}(a, m) = 1 $ が成り立つからである.逆に,「$ m $ が素数でなければ任意の $ a $ に対する(法 $ m $ についての)剰余逆元が定義できない」のは,それがBézoutの等式の保証対象外だからである.拡張Euclid互除法は,Bézoutの等式が成り立つ範囲で次のように解 $ (x, y) $ を求める4

アルゴリズム

まず,本来の素朴なEuclid互除法を思い出すと,与えられた整数のペア $ (m, n) $ に対して

  • 初期値: $ r_{0} = m, \;r_{1} = n $
  • 漸化式1: $ r_{j} = r_{j - 2} \% k_{j - 2};\; j \geq 2 $
  • 漸化式2: $ k_{j} = r_{j} / r_{j + 1};\;j \geq 0 $
  • 停止条件: $ r_{h+1} = 0 $

をベースに整数列 $ {r_{j}} $ を出力するものであった.これを行列表示すると $$ \begin{pmatrix} r_{j} \\ r_{j + 1} \end{pmatrix} = \begin{pmatrix} k_{j} & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} r_{j + 1} \\ r_{j + 2} \end{pmatrix}; \; 0\leq j \leq h - 1 $$ となる.更に,これを繰り返し用いると $$ \begin{pmatrix} r_{0} \\ r_{1} \end{pmatrix} = \begin{pmatrix} k_{0} & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} k_{1} & 1 \\ 1 & 0 \end{pmatrix} \cdots \begin{pmatrix} k_{h - 1} & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} r_{h} \\ 0 \end{pmatrix} $$ を得る.このとき,行列 $ \begin{pmatrix} k_{j} & 1 \ 1 & 0 \end{pmatrix}; \;{0\leq j\leq h - 1} $ は必ず逆行列を持つ.ゆえに,上の式を反転させて $$ \begin{pmatrix} r_{h} \\ 0 \end{pmatrix} = \begin{pmatrix} k_{h - 1} & 1 \\ 1 & 0 \end{pmatrix}^{-1} \begin{pmatrix} k_{h - 2} & 1 \\ 1 & 0 \end{pmatrix}^{-1} \cdots \begin{pmatrix} k_{0} & 1 \\ 1 & 0 \end{pmatrix}^{-1} \begin{pmatrix} r_{0} \\ r_{1} \end{pmatrix} $$ と書くと,アルゴリズムの順序通りに $ r_{0} = m, \;r_{1} = n $ を初期値として $ r_{h} = \mathrm{gcd}(m, n) $ を求める公式になっていることが分かる.ここで, $$ \begin{pmatrix} x & y \\ u & v \end{pmatrix} := \begin{pmatrix} k_{h - 1} & 1 \\ 1 & 0 \end{pmatrix}^{-1} \begin{pmatrix} k_{h - 2} & 1 \\ 1 & 0 \end{pmatrix}^{-1} \cdots \begin{pmatrix} k_{0} & 1 \\ 1 & 0 \end{pmatrix}^{-1} $$ とおくと $ \mathrm{gcd}(m, n) = mx + ny, \; 0 = ux + vy $ ,すなわち行列の要素がBézoutの等式の解になっていることが分かる.これが拡張Euclidの互除法の概略である.

実装方法

拡張Euclid互除法は,プログラム上では次のように実装すればよい.本来の素朴なEuclid互除法だと,入力値 $ (a, b) $ に対して $$ (a, b) \mapsto (b, a \% b) \mapsto (a \% b, b / (a \% b)) \mapsto \cdots \mapsto (\mathrm{gcd}(a, b), 0) $$ という列が生成されるが,これに付随する行列の要素も保持するのが拡張Euclid互除法であった.行列は掛け算によって更新されるので初期値を $$ \begin{pmatrix} x & y \\ u & v\end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 0 & 1\end{pmatrix} $$ という風に単位行列で与れば,次の更新規則を適用したときの停止時の $ (x, y) $ がBézoutの等式の解になる. $$ \begin{pmatrix} x & \\ u & v\end{pmatrix} \mapsto \begin{pmatrix} u & v \\ x - (a / b) u & y - (a / b) v\end{pmatrix}. $$

実際に,C++で拡張Euclidの互除法を実装すると,次のような感じになる.

#include <tuple>

std::tuple<int, int, int> egcd(int a, int b, int x = 1, int y = 0, int u = 0, int v = 1) {
  if (b == 0) {
    return std::make_tuple(a, x, y);
  } else {
    return egcd(b, a % b, u, v, x - (a / b) * u, y - (a / b) * v);
  }
}

/* a simple gcd */
int gcd(int a, int b) {
  if (b == 0) {
    return a;
  } else {
    return gcd(b, a % b);
  }
}

素朴なEuclid互除法の実装と見比べれば分かると思うが,引数にx, y, u, v,返り値にx, yを補うだけで拡張Euclid互除法にできるところがポイントである.また,この実装ではC++11の機能であるstd::tupleを用いているが,他にもx, y, u, vを参照渡しにするなどの方法が考えられる.今は整数型としてintを用いているが,一斉にlong longなどに変更すれば適当にサイズも大きくできる.

もちろん,このプログラムは一般のBézoutの等式を解くのには役に立つが,与えられた整数 $ a $ に対する(法 $ m $ についての)剰余逆元 $ a^{-1} $ を求めるという目的にはいくらか不要な箇所がある.上で与えた行列の更新規則 $$ \begin{pmatrix} x & y \\ u & v\end{pmatrix} \mapsto \begin{pmatrix} u & v \\ x - (a / b) u & y - (a / b) v\end{pmatrix} $$ をよく眺めてみると,実は $ (x, u) $ と $ (y, v) $ が独立に更新されていることが分かる.そこで,思い切って $ (y, v) $ の更新部分をカットした不完全な拡張ユークリッドの互除法を用意し,これを用いて剰余逆元を求める関数も実装したのが次のコード例である.

pair<int, int> iegcd(int a, int b, int x = 1, int u = 0) {
  if (b == 0) {
    return make_pair(a, x);
  } else {
    return iegcd(b, a % b, u, x - (a / b) * u);
  }
}

int modinv(int a, int m) {
  /* m should be a prime */
  std::pair<int, int> res = iegcd(a, m);
  return res.second;
}

剰余逆元が計算できると何が嬉しいか

これだけでも初等整数論の実験素材としては面白いのだが,世の中には奇妙な需要というものがいくらでもあり,

途中でオーバーフローするかも知れない計算を経て得られた(それ自体もオーバーフローするかも知れない)値を適当な素数,例えば $ 1000000007 $($ 10 $億$ 7 $)で割った値を求める

というタスクが存在する.具体的には,オーバーフローすれすれの $ 2 $ つの大きな整数の組 $ (N, K) $ に対して $$ {}_{N}C_{K} \% 1000000007 $$ を計算しなくてはいけない場面が存在する5.一見すると無茶苦茶で,何も考えずに $ {}_{\bullet}C_{\circ} $ に含まれる階乗を計算しようとすると呆気なく死ぬ.ところが,合同式に関する次のルールを知っていれば,この値を涼しい顔で計算できるのである6. $$ (a \times b) \% m = ( (a \% m) \times (b \% m) ) \% m. $$ 計算タスク自体は

  1. 組み合わせ $ {}_{N}C_{K} = N! / (K!\times (N - K)!) $ の各因子の剰余 $ (N!)\%m, \;(K!)\%m, \;((N - K)!)\%m $ を見積もる
  2. 上の公式を用いて,見積もった $ 3 $ つの因子を順次掛け算・割り算する

という手続きに尽きる.この時,手続き 1. における逆数の剰余というのが曲者で,これの計算に上で与えた modinv が必要になる.それは,剰余の世界ではある数で割ることは,その数の剰余逆元を掛けることと等価だからである.というわけで,上のルールを加味した剰余付きの掛け算・割り算・階乗・組み合わせの計算プログラムを実装すると次のようになる7

int modmult(int a, int b, int m) {
  return ((a % m) * (b % m)) % m;
}

int moddiv(int a, int b, int m) {
  return ((a % m) * modinv(b, m)) % m;
}

int modfactorial(int n, int m) {
  /* it should be that n >= 1 */
  int res = 1;
  for (int i = 2; i <= n; i++) {
    res = modmult(i, res, m);
  }
  return res;
}

int modcombination(int n, int k, int m) {
  /* it should be that n >= k >= 1 */
  int res = modfactorial(n, m);
  res = moddiv(res, modfactorial(k, m), m);
  res = moddiv(res, modfactorial(n - k, m), m);
  return res;
}

以上のコードは,repl.itコンパイル・実行可能な状態で上げておいたので,興味があればご確認いただきたい.

repl.it


  1. 整数同士の割り算は,割る数が $ 0 $ でないときに割られる数・割る数の両方を実数と同一視することで実数同士の割り算に翻訳することはできるが,今考えているのは飽くまで整数に対して定まる $ \mathbb{Z} \times \mathbb{Z} \to \mathbb{Z} $ という格好の演算である.これを良い感じに保証することで,整数を体として扱えるようになる.

  2. 混乱の元になるので,個人的には $ [a^{-1}]_{m} $ という書き方をなるべくならしたい.実際,それに近い表記をしている書物も存在する: 筑波大学・坂井公さんの計算機数学の講義ノート.

  3. このとても素晴らしい例は,剰余の逆元 (modinv) の手計算という記事からお借りした.

  4. このアルゴリズム拡張Euclid互除法と呼ばれるゆえんは,何らかの整数のペア $ (m, n) $ に対して $ \mathrm{gcd}(m, n) $ を求めるのに加えて,Bézoutの等式を解いて具体的に $ (x, y) $ を出力するからである.すなわち,$ \mathrm{gcd}(m, n) $ が予め分かっている必要はないのである.一方,今の場合は既に $ \mathrm{gcd}(m, n) = 1 $ が保証されているシステムなので,拡張Euclid互除法は多少オーバーキルである.これを分かっていると,後で見るように具体的な実装を簡便化できる.

  5. まあ競技プログラミングが代表例です.drken氏も頑張って解説してくださっています: 「1000000007 で割ったあまり」の求め方を総特集! 〜 逆元から離散対数まで 〜 - Qiita

  6. $ (a \% m) = 0 $ や $ (b \% m) = 0 $ の場合は自明なので,それ以外の場合を考えると多分簡単に証明できる.

  7. 実は,階乗計算にかかる $ n $ 回のループは $ \log n $ 回に落とせるので,これが最速のプログラムには決してならない.飽くまで簡単のためである.一方,$ 1000000007 $ というのは典型的な競技プログラミングにおけるTLEの基準値 $ \mathcal{O}(10^{10}) $ に届くか届かないかのギリギリのラインなので,実用上はやはりもう一段階の高速化が必要である.

AtCoder Beginners Contest 155 - D. Pairs

マッチングアプリと同じ名前の問題です(はい

二分探索法を二重にして解くような問題を知らなかったのでメモがてら.

問題のリンクは以下から. atcoder.jp

端的に言えば,適当な順序集合 $ \mathcal{A} := \{ A_{i} \mid 0 \leq i < N \} $ を与えた時に,集合 $ \mathcal{A} $ から次のようにして新しい集合 $ \mathcal{B} $ への変換を定めることができる.

$$ \mathcal{A} \mapsto \mathcal{B} := \{ A_{i} \times A_{j} \mid 0\leq i < j < N; \; A_{i} \in \mathcal{A} \} $$.

この時,$ A_{i} $ が今考えている問題(ABC155)のように整数ならば,$ \mathcal{B} $ もまた順序集合となるが,

$ \mathcal{B} $ の $ K $ 番目に小さい元を求めよ

というのが問題の趣旨である.ABC155の制約を見てみると

$ 2 \leq N \leq 2 \times 10^{5} $

とあり,集合 $ \mathcal{B} $ のサイズは $ N(N - 1) / 2 $ なので,$ \mathcal{B} $ の全ての元を列挙すると間に合わないことが分かる1. したがって,陽に $ \mathcal{B} $ を走査しないようなアルゴリズムが必要になるが,$ \mathcal{O}(\log |\mathcal{B}|) $ の時間計算量であれば間に合いそうであり,$ \log $ 系のアルゴリズムといえば二分探索法である2. というわけで,二分探索法による解法が有効であることをこれから見ていく3. 結局のところ,必要なのは問題文の言い換えであり,下から $ K $ 番目の $ \mathcal{B} $ の元というのは

全ての元が $ B $ 未満であるような $ \mathcal{B} $ の部分集合のサイズを $ k_{B} $ としたとき,$ k_{B} = K $ となるような $ B $

であることに気付けば後は実装するだけである.とは言っても,筆者もこれに気付くことはできなかったので,とりあえず簡単な例を扱うことにする.

仮に,$ \mathcal{A} = \{ 1, 2, 3, 4 \} $ とすると,$ \mathcal{B} = \{ 2, 3, 4, 6, 8, 12 \} $ である.この時,下から $ K = 5 $ 番目の元は $ 8 $ である.一方,(空集合も含めた)$ \mathcal{B} $ の部分集合をその制約ごとに分類すると

  • $ | \{ B \in \mathcal{B} \mid B < B^{*}; \; -\infty < B^{*} \leq 2 \} | = 0 $
  • $ | \{ B \in \mathcal{B} \mid B < 3 \} | = 1 $
  • $ | \{ B \in \mathcal{B} \mid B < 4 \} | = 2 $
  • $ | \{ B \in \mathcal{B} \mid B < 5 \} | = 3 $
  • $ | \{ B \in \mathcal{B} \mid B < 6 \} | = 3 $
  • $ | \{ B \in \mathcal{B} \mid B < 7 \} | = 4 $
  • $ | \{ B \in \mathcal{B} \mid B < 8 \} | = 4 $
  • $ | \{ B \in \mathcal{B} \mid B < 9 \} | = 5 $
  • $ | \{ B \in \mathcal{B} \mid B < 10 \} | = 5 $
  • $ | \{ B \in \mathcal{B} \mid B < 11 \} | = 5 $
  • $ | \{ B \in \mathcal{B} \mid B < 12 \} | = 5 $
  • $ | \{ B \in \mathcal{B} \mid B < B^{*}; \; 13 \leq B^{*} < \infty \} | = 6 $

という上下に有界で階段状の構造を持つことが分かる(確率論における累積分布関数のようなもの).このうち,全ての元が $ 8 $ 未満となるような $ \mathcal{B} $ の部分集合は,上から7つ目の $ | \{ B \in \mathcal{B} \mid B < 8 \} | = 4 $ であり,部分集合のサイズは $ K - 1 = 4 $ となっている.

こう書くと,ただ数学的な主張を逆にしただけに見えるが,全ての元が $ 8 $ 未満となるような $ \mathcal{B} $ の部分集合は,サイズが$ 4 $(すなわち$ 5 $未満)であるもののうち最大であることに気付けば勝ちである4.この最大のポイントを探り当てるのが二分探索法であり,実装すると次のようになる.

#include <bits/stdc++.h>
using namespace std;

const long long INF = static_cast<long long>(1e18) + 1LL;

// Generated by 1.1.6 https://github.com/kyuridenamida/atcoder-tools  (tips: You use the default template now. You can remove this line by using your custom template)
int main()
{
    long long N = 0, K = 0;
    cin >> N >> K;

    vector<long long> A(N, 0);
    for (size_t i = 0; i < N; i++)
    {
        cin >> A[i];
    }

    sort(A.begin(), A.end());

    long long res = 0;

    long long left = -INF, right = INF;
    while (left + 1 < right)
    {
        long long center = (left + right) / 2, x = 0;
        for (size_t i = 0; i < N; i++)
        {
            if (A[i] < 0)
            {
                long long l = -1, r = N;
                while (l + 1 < r)
                {
                    long long c = (l + r) / 2;
                    if (A[i] * A[c] < center)
                    {
                        r = c;
                    }
                    else
                    {
                        l = c;
                    }
                }
                x += (N - r);
            }
            else
            {
                long long l = -1, r = N;
                while (l + 1 < r)
                {
                    long long c = (l + r) / 2;
                    if (A[i] * A[c] < center)
                    {
                        l = c;
                    }
                    else
                    {
                        r = c;
                    }
                }
                x += r;
            }

            if (A[i] * A[i] < center)
            {
                x--;
            }
        }
        x /= 2;

        if (x < K)
        {
            left = center;
        }
        else
        {
            right = center;
        }
    }

    cout << left << endl;

    return 0;
}

まず,大きな括りの部分で,上からどんな値で部分集合の元を抑えるべきかという変数の範囲leftrightを狭めている(外殻部分の二分探索).そして,その中で各値 $ A_{i} \in \mathcal{A};\; i = 0, 1, 2, \ldots, N - 1 $ について,自身との積がcenter = left + right未満に抑えられるような $ A_{j};\; j\neq i $ の個数をカウントしている5. このカウントの過程が内殻部分の二分探索である.そして,カウントされた個数がK未満となるようなギリギリのところに収束すれば,それを満たすrightを出力せよ,という算段である.

以上から,ナイーブにソート・ランダムアクセスを行った場合には時間計算量 $ \mathcal{O}(N^{2}) $ であったものが,$ \mathcal{O}(\log (N\log N)) $ に改善されていることが分かる.

こうやって執筆している段階でも,未だに腑に落ちた間隔がしないのが二分探索法を使う系統の難しさである(個人の感想).慣れていくしかなさそうです.


  1. 大まかに見積もって,比較的実行速度の早いC++を用いる場合でも,$ 10^{10} $ 回程度の操作を行うとなるとTLEになる可能性が高い.

  2. 別にそういう思いつき方をして解けたわけでもなく,実際のところ本番ではC問題に時間を使いすぎて時間を溶かしてしまった.無念である.虚無じゃのう.

  3. snukeさんの解説動画(https://www.youtube.com/watch?v=SG60Cp9pSog)を多いに参考にした.本当にありがたい.

  4. この手の数学的事実については,それ自体を式として導出・理解するのは大変なので,整数変数・整数値関数のグラフを適当に書いて図式的に理解するのが早い.早そうである.

  5. 今は $ A_{i} $ の区間が負の整数まで拡張されたことで場合分けを伴い,実装が多少複雑になっている.