Skilling-Hutchinson の跡推定量

Probability
Functional Analysis
Author

司馬博文

Published

8/20/2024

Modified

8/22/2024

概要
Skilling-Hutchinson の跡推定量は，跡の計算 \(O(d^2)\) を \(O(d)\) に落とすことができる Monte Carlo 法である．

1 命題

(Skilling, 1989)-(Hutchinson, 1990)

任意の正方行列 \(A\in M_n(\mathbb{R})\) と，\(\mathrm{V}[X]=I_n\) を満たす確率ベクトル \(X\in\mathcal{L}(\Omega;\mathbb{R}^n)\) について，

\[ \operatorname{Tr}(A)=\operatorname{E}[X^\top AX]. \]

\[ X^\top AX=\operatorname{Tr}(AXX^\top) \] に注意する．これは，一般に \(x,y\in\mathbb{R}^n\) に対して \[ yx^\top=\begin{pmatrix}y_1\\\vdots\\y_n\end{pmatrix}(x_1\;\cdots\;x_n)=\begin{pmatrix}y_1x_1&\cdots&y_1x_n\\\vdots&\ddots&\vdots\\y_nx_1&\cdots&y_nx_n\end{pmatrix} \] \[ \therefore\qquad\operatorname{Tr}(yx^\top)=x^\top y \] が成り立つためである．

よって，次のように計算できる：

\[\begin{align*} \operatorname{E}[X^\top AX]&=\operatorname{E}[\operatorname{Tr}(AXX^\top)]\\ &=\operatorname{Tr}(\operatorname{E}[AXX^\top])\\ &=\operatorname{Tr}(A\operatorname{E}[XX^\top])=\operatorname{Tr}(A). \end{align*}\]

(Hutchinson, 1990) では \(A\) を対称行列に，\(X\) を中心化された確率変数に限って示されている．

(Skilling, 1989) では (Hutchinson, 1990) のように命題の形では提示していないが，同様の推定量を提案しており，これと一般化跡 (generalized trace) と Chebyshev 多項式の議論を通じて，\(A\) のスペクトルのベイズ推定を議論している．

2 推定量の性質

実用上，\(X\) の分布は標準 Gauss や Rademacher 分布などが用いられる．

命題（推定量の分散）

\(A\in S_n(\mathbb{R})\) を対称行列とする．

  1. \(X\sim\mathrm{N}(0,I_n)\) のとき， \[ \mathrm{V}[X^\top AX]=2\operatorname{Tr}(A^2)=2\|A\|^2_\mathrm{HS}. \]

  2. \(X\sim\mathrm{Rad}^{\otimes n}\) のとき， \[ \mathrm{V}[X^\top AX]=2\sum_{i\ne j}a_{ij}^2. \]

  1. \(A\in S_n(\mathbb{R})\) が正定値対称であるとき，ある直交行列 \(U\in O_n(\mathbb{R})\) と対角行列 \(\Lambda=\mathrm{diag}(\lambda_1,\dots,\lambda_n)\) が存在して， \[ A=U\Lambda U^\top. \] \(Y:=U^\top X\) と定めるとやはり \(Y\sim\mathrm{N}(0,I_n)\) であり， \[ X^\top AX=X^\top U\Lambda U^\top X=Y^\top\Lambda Y, \] \[ \therefore\qquad\mathrm{V}[X^\top AX]=2\sum_{i=1}^n\lambda_i^2=2\operatorname{Tr}(\Lambda^2)=2\operatorname{Tr}(A^2). \]
  2. 一般の \(X\in\mathcal{L}(\Omega;\mathbb{R}^n)\) に関して， \[ \mathrm{V}[X^\top AX]=\sum_{i,j,k,l=1}^na_{ij}a_{kl}\biggr(\operatorname{E}[X_iX_jX_kX_l]-\operatorname{E}[X_iX_j]\operatorname{E}[X_kX_l]\biggl). \]
系

\(A\) が対称行列で \(\operatorname{E}[X]=0\) であるとき，\(X\) は Rademacher とした場合が最小分散不偏推定量を定める (Hutchinson, 1990, p. 命題1)．

3 応用

  • 残差フロー (residual flow) では Jacobian の推定が焦点になる．これに Skilling-Hutchinson の跡推定量を用いることができる．
  • Neural ODE において，Jacobian の跡 \(\operatorname{Tr}(J_{F_t}(x_t))\) の計算は Skilling-Hutchinson の跡推定量を用いれば \(O(d)\) で済む (Grathwohl et al., 2019)．
  • Sliced Score Matching の目的関数は，Skilling-Hutchinson の跡推定量により Jacobian \(Ds_\theta\) を推定したスコアマッチングと解釈できる．

4 文献紹介

(Adams et al., 2018) では (Skilling, 1989) の研究を踏襲し，大規模行列のスペクトル（密度）推定に向けて，Skilling-Hutchinson の跡推定量の拡張が議論されている．

(Meyer et al., n.d.) では Skilling-Hutchinson の跡推定量を改良したアルゴリズムが提案されている．

References

Adams, R. P., Pennington, J., Johnson, M. J., Smith, J., Ovadia, Y., Patton, B., and Saunderson, J. (2018). Estimating the spectral density of large implicit matrices.
Avron, H., and Toledo, S. (2011). Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix. J. ACM, 58(2).
Grathwohl, W., Chen, R. T. Q., Bettencourt, J., and Duvenaud, D. (2019). Scalable Reversible Generative Models with Free-form Continuous Dynamics. In International conference on learning representations.
Hutchinson, M. F. (1990). A Stochastic Estimator of the Trace of the Influence Matrix for Laplacian Smoothing Splines. Communications in Statistics - Simulation and Computation, 19(2), 433–450. doi: 10.1080/03610919008812866.
Meyer, R. A., Musco, C., Musco, C., and Woodruff, D. P. (n.d.). Hutch++: Optimal stochastic trace estimation. In 2021 symposium on simplicity in algorithms (SOSA), pages 142–155.
Skilling, J. (1989). The eigenvalues of mega-dimensional matrices. In J. Skilling, editor, Maximum entropy and bayesian methods: Cambridge, england, 1988, pages 455–466. Dordrecht: Springer Netherlands.

Footnotes

  1. (Hutchinson, 1990, p. 437), (Avron and Toledo, 2011, p. 補題9), (Adams et al., 2018, p. 命題4.2) も参照．↩︎