高専プロコン2022競技部門・理論値解法

1回戦からエキシビションまでで参加した全問題で理論値（ $\stackrel{\text{def}}{=}$ 最も良い値）の得点を出して勝ち進んだので、「理論値解法」として紹介します。実装（特に並列化）が丁寧ではないので確かではないですが、乱択アルゴリズムではありません。

おことわり

これは実際に使用したアルゴリズムの詳細と考察です。

簡単のため、理論値に関係する部分のみ切り出して解説します。そうすると問題設定がシンプルで助かります。

音声解析は特に勉強しなかったのでその話はできません。

かなり競技プログラミング（ algorithm ）をする解法です。

問題概要

問題データ（ベクトル）を $\boldsymbol{R}_1,\boldsymbol{R}_2,\ldots ,\boldsymbol{R}_{88}$ とします。これらは英語と日本語の上毛かるた各 $44$ 枚分を読み上げた音声の波形を $48000[\text{Hz}]$ でサンプリングしたものを表します。

いわゆるベクトルを考えると都合が悪いので、ベクトルの添え字は任意の整数をとれるものとして考えます。

ベクトルに対して、以下のような操作を考えます。

• 足し算 ： ベクトルの足し算です。添え字が同じ要素をそれぞれ足します。 $\boldsymbol{x}+\boldsymbol{y}$ で表します。
• 内積 ： $\boldsymbol{x}\cdot\boldsymbol{y}=\sum_{i}\boldsymbol{x}_i\boldsymbol{y}_i$
• 絶対値の $2$ 乗 ： $\left\lvert \boldsymbol{z}\right\rvert ^2=\sum_{i}{\boldsymbol{z}_i}^2=\boldsymbol{z}\cdot\boldsymbol{z}$
• シフト ： ベクトル $\boldsymbol{x}$ と整数 $k$ から、 $\boldsymbol{x}$ の添え字に $k$ を加算して得られるベクトルを $\boldsymbol{x}\ll k$ と表します。つまり、これは音声の再生時刻を $k$ だけ遅らせる操作です。
• 反転 ： 添え字の符号を変えます。すなわち、 $\boldsymbol{x}_i$ と $\boldsymbol{x}_{-i}$ が入れ替わります。 $\boldsymbol{x}$ を反転して得られるベクトルを $\overline{\boldsymbol{x}}$ で表します。
• クリッピング ： ベクトルの添え字が $0$ 未満または $T$ 以上の要素を $0$ に変更し、さらに各要素が $-2^{15}$ 以上 $2^{15}-1$ 以下の最も近い整数に変更します。この操作は $\text{CLIP}(\boldsymbol{x},T)$ で表します。 WAVE 形式で保存するときなどに生じる制約のための操作です。
• 畳み込み ： $2$ つのベクトル $\boldsymbol{x}$ , $\boldsymbol{y}$ の畳み込みを $\boldsymbol{x}\ast\boldsymbol{y}$ と表します。 $\boldsymbol{z}=\boldsymbol{x}\ast\boldsymbol{y}$ は、次のように計算します。
  • $\boldsymbol{z}_k=\sum_{i+j=k}\boldsymbol{x}_i\boldsymbol{y}_j$
  • $\boldsymbol{z}_k=\boldsymbol{x}\cdot\left\lparen\overline{\boldsymbol{y}}\ll k\right\rparen$

$1$ つの問題は次のパラメータと制約から構築されます。

• 重ね合わせる読みデータの個数 $N$ $(3\leq N \leq 20)$
• 読みデータの番号 $I_i$ と再生時刻を表す整数 $W_i$ の組の集まり
• 問題データの長さ $T_1$ （実際の出題ではすべて $T_1=24000$ ）
• 制限時間 $\text{TL}$ （実際の出題ではすべて $\text{TL}=180[\text{sec}]$ ）
• 同じ札の読み上げは含まれない。

問題データの断片 $\boldsymbol{P}$ を次のように生成します。

• $\boldsymbol{P}=\text{CLIP}\left\lparen\sum_{i=1}^N \left\lparen\boldsymbol{R}_{I_i}\ll W_i\right\rparen ,T_1\right\rparen$

目標は、問題データ $\boldsymbol{P}$ が与えられてから時間 $\text{TL}$ 以内に $I_i$ をすべて求めることです。

閑話

この問題設定が、本当に運営が想定していたものと一致するのか、与えられた情報が表しているものであるのかには決定的な報告がありません。実際にこの記事のアルゴリズムを開発した人もあまり正確に把握できていなかったので説得力が小さいのですが、根拠づけとしては、配布された問題サンプルに添付されている「作問情報データ」の形式から導かれる問題設定として適していると言えるでしょう。

アプローチ

音声は重ね合わせると、たいてい大きい音になります。逆に、そのような重ね合わせを打ち消せば音は小さくなるはずです。すなわち、問題データの音が小さくなるように読みデータの波形を打ち消せば、それが問題データであった可能性が高いです。

まず、音声を表すベクトルに対して音の大きさを評価する指標 $f(\boldsymbol{x})$ を考え、それを多くの点で評価（点は、問題のパラメータ $(I_i,W_i)$ の候補）し、最善のものを選んで答えとします。

例えば $k$ 個の組 $(I_i,W_i)$ が判別したとすると、それらのみを用いて断片 $\boldsymbol{P}^{\prime}$ を作り、 $\boldsymbol{P}\leftarrow\boldsymbol{P}-\boldsymbol{P}^{\prime}$ （ $\leftarrow$ は代入操作）のように波形を打ち消して次の候補を探します。

計算の高速化

$f(\boldsymbol{x})$ は、

• 音の大きさをよく表していて、
• 多点評価が高速である

必要があります。これを満たすものとして考えられるのが、 $f(\boldsymbol{x})=\left\lvert \boldsymbol{x}\right\rvert ^2$ です。

音の大きさをよく表していることは、例えば、モデルとして $A_t=\sin(\omega t)$ と $B_t=\sin(2\omega t)$ を合成して積分で評価すると値がよく増加するのがわかります。

そして、実は多点評価が非常に高速です。

評価すべき式は $S_{i,w}=f\left\lparen\boldsymbol{x}-\text{CLIP}\left\lparen\boldsymbol{R}_i\ll w,T_1\right\rparen\right\rparen$ で、これを多くの点 $(i,w)$ について評価します。そのために、次の式変形を行います。

$$\begin{aligned} S_{i,w}&= {\left\lvert\boldsymbol{x}-\text{CLIP}\left\lparen\boldsymbol{R}_i\ll w,T_1\right\rparen\right\rvert } ^2 \\ &= {\left\lvert\boldsymbol{x}\right\rvert }^2-2\boldsymbol{x}\cdot\text{CLIP}\left\lparen\boldsymbol{R}_i\ll w,T_1\right\rparen + {\left\lvert\text{CLIP}\left\lparen\boldsymbol{R}_i\ll w,T_1\right\rparen\right\rvert }^2 \\ &= {\left\lvert\boldsymbol{x}\right\rvert }^2-2\boldsymbol{x}\cdot\left\lparen\boldsymbol{R}_i\ll w\right\rparen +{\left\lvert\text{CLIP}\left\lparen\boldsymbol{R}_i\ll w,T_1\right\rparen\right\rvert }^2 \\ &= {\left\lvert\boldsymbol{x}\right\rvert }^2-2\sum_{t}{\boldsymbol{x}_t(\boldsymbol{R}_i)_{t-w}}+\sum_{t=0}^{T_1-1}{(\boldsymbol{R}_i)_{t-w}}^2 \\ &= {\left\lvert\boldsymbol{x}\right\rvert }^2-2\sum_{t}{\boldsymbol{x}_t\overline{(\boldsymbol{R}_i)}_{w-t}}+\sum_{t=-w}^{T_1-w-1}{(\boldsymbol{R}_i)_{t}}^2 \\ &= {\left\lvert\boldsymbol{x}\right\rvert }^2-2{\left\lparen\boldsymbol{x}\ast\overline{(\boldsymbol{R}_i)}\right\rparen }_w+\sum_{t=-\infty}^{T_1-w-1}{(\boldsymbol{R}_i)_{t}}^2-\sum_{t=-\infty}^{-w-1}{(\boldsymbol{R}_i)_{t}}^2 \end{aligned}$$

$i$ を固定すると、重要な $w$ すべてに対する $S_{i,w}$ の値を高速に計算できることを示します。まず、 $\left\lvert\boldsymbol{x}\right\rvert$ はここでは定数です（第 $1$ 項）。次に、 $\boldsymbol{R}_i$ の要素数は有限なので、 $\sum_{t=-\infty}^{Z}{(\boldsymbol{R}_i)_{t}}^2$ をすべての $Z$ について計算しておくことができます（第 $3$ , $4$ 項）。

最後に、第 $2$ 項の畳み込みを計算します。これには数論変換を用いた準線形時間の畳み込みを使用しますが、値が $31$ ビットに収まらない可能性を考えて $2$ 個の除数で剰余を計算し、大きい除数に変換するアルゴリズムを併用します。この第 $2$ 項の部分は初めて原理から理解しようとするとまあまあ難しいので、外部ブログのリンクを貼っていきます。 modulo の世界がイメージしづらい場合は、代わりに FFT を学ぶとよいです。

• FFT と高速な畳み込みの解説（ by tatyam ） : https://trap.jp/post/1386/
• FFT から NTT へアルゴリズムを変形する解説（ by Sen ） ： https://sen-comp.hatenablog.com/entry/2021/02/06/180310
• この部分について、モンゴメリ乗算と SIMD （ AVX2 ）を用いてシングルスレッドでの計算を高速化します。この部分は Nyaan’s Library の実装を改造して使用しました。 https://nyaannyaan.github.io/library/ntt/ntt-avx2.hpp
• $M_0=998244353$ 、 $M_1=924844033$ で割った余りをそれぞれ計算したあと、答えに適切な係数を掛けてできた $2$ つを足し合わせると、 $\mod{M_0M_1}$ で正しい答えが得られます。「クリエイティヴなブログ」(https://www.creativ.xyz/ect-gcd-crt-garner-927/) が分かりやすいらしい？(https://twitter.com/red_black_gpgpu/status/1180451600499740672)
• この部分は自分でモンゴメリ乗算と SIMD で実装しました。参照した命令リスト https://www.intel.com/content/www/us/en/docs/intrinsics-guide/index.html
• CUDA でプログラムを書いて GPU を動かそうと思いましたが、開発途中からなぜかうまくいかなくなりました。

結局、ある $i$ に対して計算量 $O(L\log L)$ （ $L$ は $\boldsymbol{R}_i$ の要素数と $\boldsymbol{x}$ の要素数の和 ） で $S_{i,w}$ の値を列挙できることになります。 簡単のため、$L\leq 524288$ を仮定して実装しました。（断片の長さを $120000$ 以下に抑えることで $L\leq 524288$ が成り立ちます。たぶん。）

合成途中に 16bit 整数の範囲外の値がある場合について

実は上で計算した式の冒頭 $f\left\lparen\boldsymbol{x}-\text{CLIP}\left\lparen\boldsymbol{R}_i\ll w,T_1\right\rparen\right\rparen$ は、その意味合いと照らし合わせると正確ではありません。例えば、問題データ $\boldsymbol{x}$ の計算の過程でクリッピングにより値が変化した要素があると、合成する波形を $1$ つ減らした場合の問題データを単純な引き算で計算できません。このように、クリッピングの存在により、正確かつ高速に多点評価することは困難になります。

一方で、合成元の読みデータのサンプルの値の幅は $2^{15}$ と比べると小さく、また正負が混在するので、クリッピングが発生することはあれど全体に占める割合は小さいです。よって、これによる誤差はむしろ無視できます。

これを踏まえて、可能な対策として、「問題から引く（ $P^{\prime}=P-x_1-x_2-\ldots$ ）」のではなく、「構築しなおして差をとる（ $P^{\prime}=P-(x_1+x_2+\ldots )$ ）」計算順序をとりました。構築しなおした後に毎回クリッピングを入れることで、問題データのクリッピングによる誤差の累積が抑えられると考えられます。効果があったかはわかりません。

全体のアルゴリズム

読みデータの個数を $K$ とすると、すべての組 $(i,w)$ について $S_{i,w}$ の値を計算する際、実用上のボトルネックの計算量は $O(KL\log L)$ になります。これらの $S_{i,w}$ の最小値に従って、回答候補に $(i,w)$ を加えます。これを $N$ 回繰り返すことで回答候補が $N$ 組集まり、とりあえず問題の解答を作ることができます。全体の計算量（実用上のボトルネック）は $O(NKL \log L)$ です。 $i$ に関するループは並列化します。

もし、これでよい回答が得られなかった、すなわち回答の情報から問題データの再現がうまくいかなかった場合は、回答から組 $(i,w)$ を $1$ つ取り除き、さきほどと同様の方法で最善の組 $(i,w)$ を追加します。 $N$ 個それぞれについて取り除くのを試すため、 $1$ 周の計算量は $O(NKL \log L)$ です。

精度の改善

これだけでは日本語のデータがうまく処理できない場合があります。（エキシビションの問題はすべて日本語の読みデータの重ね合わせでした。）このため、少し関数の材料を増やします。

音の大きさの指標と同様の計算が可能な指標として、ベクトルの隣接要素の差が考えられます。ベクトル $\boldsymbol{x}$ の隣接要素の差 $(\boldsymbol{x}_{i+1}-\boldsymbol{x}_i)$ をとって得られるベクトルを $\Delta\boldsymbol{x}$ とします。（添え字があいまいですが、この解説ではそのままにします。）

$f(\boldsymbol{x})$ を次のように変更します。

$f(\boldsymbol{x})=\left\lvert \boldsymbol{x}\right\rvert ^2+40\left\lvert \Delta\boldsymbol{x}\right\rvert ^2$

$40$ 倍するのは、かなり綺麗な波形である読みデータの隣接サンプルの差はそのままでは小さすぎるからです。

実装しようとすると、添え字が $1$ 個ずれたり、ベクトルの長さがさっきと異なったりして処理が難しくなるので、妥協しながら実装すると、日本語の読みの重ね合わせでもかなりよい精度が得られます。

怪しい考察をすると、この差分をとる操作が波形の微分の働きをするので、強力なヒントとなるべき高周波の子音を強調していると考えられます。

実測

決勝で使用したのは 4 コアの CPU をもつノート PC で、重ね合わせ $1$ 個に対するイテレーション（計算量 $O(KL \log L)$ ）あたり $1$ 秒程度で処理できました。ちなみに $KL \log_2 L\approx 88\times 524288\times 19 \approx 8.8\times 10^8$ です。

高専プロコンで実際に解いた問題では $1$ 周で完全な復元ができたため、かかった時間はおよそ $N$ 秒（ $N$ は重ね合わせる読みデータの個数）でした。

実は

NAPROCK の問題は $2$ 周かかった。

ところであのパンフレット原稿は何だったのか

解法がかなりシンプルなので、書いたらバレるんですよ……