2つの行列を「掛け合わせる」方法として、通常の行列積はよく知られています。しかし、2つの行列の全ての要素の組み合わせを考えたいとき、通常の行列積では不十分です。たとえば、2つの独立なシステムを組み合わせた複合システムを行列で表現するにはどうすればよいでしょうか。
この問いに答えるのがクロネッカー積(Kronecker product)です。$m \times n$ 行列 $\bm{A}$ と $p \times q$ 行列 $\bm{B}$ のクロネッカー積 $\bm{A} \otimes \bm{B}$ は $mp \times nq$ の行列を生成します。これは量子力学のテンソル積、制御工学のシステム結合、統計学の共分散構造のモデリングなど、驚くほど多くの分野で自然に現れる演算です。
クロネッカー積を理解すると、以下のような応用が開けます。
- 量子情報科学: 複合量子系の状態空間の記述(テンソル積構造)
- 制御工学: 多変数システムの状態方程式のベクトル化
- 統計学: 分離可能な共分散構造のモデリング
- 画像処理: 2次元フィルタリングの行列表現
- 数値線形代数: シルベスター方程式・リアプノフ方程式の効率的な求解
本記事の内容
- クロネッカー積の直感的な理解と数学的定義
- 基本性質の証明(混合積性質、転置、トレース等)
- vec演算子との関係と行列方程式への応用
- Pythonによる実装と応用例
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
- 行列の基本演算 — 行列の積とブロック行列
- 固有値と固有ベクトル — クロネッカー積の固有値との関係
- 内積と直交性 — ベクトル空間の基本概念
クロネッカー積とは — 直感的な理解
クロネッカー積を直感的に理解するために、まず簡単な例を考えましょう。
2人のプレーヤーがそれぞれ独立にサイコロを振る場面を想像してください。プレーヤー1のサイコロの出目を2つの状態 $\{1, 2\}$(簡略化のため)、プレーヤー2のサイコロも2つの状態 $\{A, B\}$ とします。全体の状態は $(1, A), (1, B), (2, A), (2, B)$ の4通りです。
各プレーヤーの遷移行列を $\bm{P}_1$($2 \times 2$)と $\bm{P}_2$($2 \times 2$)とすると、複合システムの遷移行列は $\bm{P}_1 \otimes \bm{P}_2$($4 \times 4$)になります。つまり、クロネッカー積は独立なシステムの組み合わせを一つの大きな行列で表現する方法なのです。
もう一つのアナロジーとして、タイル貼りを考えてみましょう。$\bm{A}$ がタイルの配置パターンを、$\bm{B}$ が各タイル内部の模様を表すとすると、$\bm{A} \otimes \bm{B}$ は「$\bm{A}$ の各成分の位置に、その成分の値でスケールした $\bm{B}$ を配置する」操作に対応します。これはまさにクロネッカー積の定義そのものです。
この演算は、レオポルト・クロネッカー(1823-1891)にちなんで名付けられましたが、実際にはこの概念自体はヨハン・ペーター・グスタフ・ルジューヌ・ディリクレを含む多くの数学者によって独立に使用されていました。
直感的な理解ができたところで、数学的な定義に進みましょう。
クロネッカー積の数学的定義
定義
$m \times n$ 行列 $\bm{A} = (a_{ij})$ と $p \times q$ 行列 $\bm{B}$ のクロネッカー積(テンソル積、直積とも呼ぶ)$\bm{A} \otimes \bm{B}$ は、次のように定義される $mp \times nq$ 行列です。
$$ \begin{equation} \bm{A} \otimes \bm{B} = \begin{pmatrix} a_{11}\bm{B} & a_{12}\bm{B} & \cdots & a_{1n}\bm{B} \\ a_{21}\bm{B} & a_{22}\bm{B} & \cdots & a_{2n}\bm{B} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1}\bm{B} & a_{m2}\bm{B} & \cdots & a_{mn}\bm{B} \end{pmatrix} \end{equation} $$
つまり $\bm{A}$ の各成分 $a_{ij}$ を $a_{ij}\bm{B}$($p \times q$ のブロック)で置き換えた行列です。
具体例
$$ \bm{A} = \begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix}, \quad \bm{B} = \begin{pmatrix} 0 & 5 \\ 6 & 7 \end{pmatrix} $$
のとき
$$ \bm{A} \otimes \bm{B} = \begin{pmatrix} 1 \cdot \begin{pmatrix} 0 & 5 \\ 6 & 7 \end{pmatrix} & 2 \cdot \begin{pmatrix} 0 & 5 \\ 6 & 7 \end{pmatrix} \\ 3 \cdot \begin{pmatrix} 0 & 5 \\ 6 & 7 \end{pmatrix} & 4 \cdot \begin{pmatrix} 0 & 5 \\ 6 & 7 \end{pmatrix} \end{pmatrix} = \begin{pmatrix} 0 & 5 & 0 & 10 \\ 6 & 7 & 12 & 14 \\ 0 & 15 & 0 & 20 \\ 18 & 21 & 24 & 28 \end{pmatrix} $$
なぜこの定義が自然なのかを考えてみましょう。$\bm{A}$ が空間 $V$ 上の線形写像を、$\bm{B}$ が空間 $W$ 上の線形写像を表すとき、$\bm{A} \otimes \bm{B}$ はテンソル積空間 $V \otimes W$ 上の線形写像を表します。$V$ が $n$ 次元、$W$ が $q$ 次元なら、$V \otimes W$ は $nq$ 次元であり、これがクロネッカー積のサイズ $nq$ に対応しています。
注意すべき点として、一般に $\bm{A} \otimes \bm{B} \neq \bm{B} \otimes \bm{A}$ です(非可換)。ただし、適切な置換行列 $\bm{P}$ を用いれば $\bm{B} \otimes \bm{A} = \bm{P}(\bm{A} \otimes \bm{B})\bm{P}^T$ と関係付けることができます。
定義を把握したところで、次にクロネッカー積の豊富な代数的性質を見ていきましょう。
クロネッカー積の基本性質
代数的性質一覧
クロネッカー積は以下の性質を持ちます。ここで $\bm{A}, \bm{C}$ は $m \times n$、$\bm{B}, \bm{D}$ は $p \times q$、$\alpha$ はスカラーとします。
| 性質 | 式 |
|---|---|
| 分配法則(左) | $(\bm{A} + \bm{C}) \otimes \bm{B} = \bm{A} \otimes \bm{B} + \bm{C} \otimes \bm{B}$ |
| 分配法則(右) | $\bm{A} \otimes (\bm{B} + \bm{D}) = \bm{A} \otimes \bm{B} + \bm{A} \otimes \bm{D}$ |
| スカラー倍 | $(\alpha \bm{A}) \otimes \bm{B} = \bm{A} \otimes (\alpha \bm{B}) = \alpha (\bm{A} \otimes \bm{B})$ |
| 結合法則 | $(\bm{A} \otimes \bm{B}) \otimes \bm{C} = \bm{A} \otimes (\bm{B} \otimes \bm{C})$ |
| 転置 | $(\bm{A} \otimes \bm{B})^T = \bm{A}^T \otimes \bm{B}^T$ |
| 共役転置 | $(\bm{A} \otimes \bm{B})^* = \bm{A}^* \otimes \bm{B}^*$ |
混合積性質(最重要)
クロネッカー積の最も強力な性質は混合積性質(mixed-product property)です。
$$ \begin{equation} (\bm{A} \otimes \bm{B})(\bm{C} \otimes \bm{D}) = (\bm{A}\bm{C}) \otimes (\bm{B}\bm{D}) \end{equation} $$
ただし、行列積 $\bm{A}\bm{C}$ と $\bm{B}\bm{D}$ がそれぞれ定義できるサイズであることが必要です。
証明: $\bm{A}$ が $m \times n$、$\bm{B}$ が $p \times q$、$\bm{C}$ が $n \times r$、$\bm{D}$ が $q \times s$ とします。
$(\bm{A} \otimes \bm{B})$ の $(i, j)$ ブロック($p \times q$)は $a_{ij}\bm{B}$ です。$(\bm{C} \otimes \bm{D})$ の $(j, k)$ ブロック($q \times s$)は $c_{jk}\bm{D}$ です。
積 $(\bm{A} \otimes \bm{B})(\bm{C} \otimes \bm{D})$ の $(i, k)$ ブロックは
$$ \sum_{j=1}^{n} (a_{ij}\bm{B})(c_{jk}\bm{D}) = \sum_{j=1}^{n} a_{ij}c_{jk}\bm{B}\bm{D} = \left(\sum_{j=1}^{n} a_{ij}c_{jk}\right)\bm{B}\bm{D} = (\bm{A}\bm{C})_{ik} \bm{B}\bm{D} $$
これは $(\bm{A}\bm{C}) \otimes (\bm{B}\bm{D})$ の $(i, k)$ ブロックに一致します。$\square$
逆行列
混合積性質の直接的な帰結として、$\bm{A}$ と $\bm{B}$ が共に正則ならば
$$ (\bm{A} \otimes \bm{B})^{-1} = \bm{A}^{-1} \otimes \bm{B}^{-1} $$
が成り立ちます。確認すると $(\bm{A} \otimes \bm{B})(\bm{A}^{-1} \otimes \bm{B}^{-1}) = (\bm{A}\bm{A}^{-1}) \otimes (\bm{B}\bm{B}^{-1}) = \bm{I}_m \otimes \bm{I}_p = \bm{I}_{mp}$ です。
固有値とトレース
$\bm{A}$ の固有値を $\lambda_1, \ldots, \lambda_m$、$\bm{B}$ の固有値を $\mu_1, \ldots, \mu_p$ とすると
$$ \begin{equation} \bm{A} \otimes \bm{B} \text{ の固有値} = \{\lambda_i \mu_j : i = 1, \ldots, m, \; j = 1, \ldots, p\} \end{equation} $$
証明: $\bm{A}\bm{u} = \lambda \bm{u}$、$\bm{B}\bm{v} = \mu \bm{v}$ のとき、$\bm{u} \otimes \bm{v}$(クロネッカー積のベクトル版)に対して
$$ (\bm{A} \otimes \bm{B})(\bm{u} \otimes \bm{v}) = (\bm{A}\bm{u}) \otimes (\bm{B}\bm{v}) = \lambda\mu(\bm{u} \otimes \bm{v}) $$
よって $\lambda\mu$ が固有値、$\bm{u} \otimes \bm{v}$ が対応する固有ベクトルです。$\square$
トレースについては
$$ \text{tr}(\bm{A} \otimes \bm{B}) = \text{tr}(\bm{A}) \cdot \text{tr}(\bm{B}) $$
が成り立ちます。これは $\bm{A} \otimes \bm{B}$ の対角ブロックが $a_{ii}\bm{B}$ であることから、対角成分の和が $\sum_i a_{ii} \text{tr}(\bm{B}) = \text{tr}(\bm{A}) \text{tr}(\bm{B})$ と計算できることによります。
行列式
$\bm{A}$ が $m \times m$、$\bm{B}$ が $p \times p$ の正方行列のとき
$$ \begin{equation} \det(\bm{A} \otimes \bm{B}) = (\det \bm{A})^p \cdot (\det \bm{B})^m \end{equation} $$
これは固有値の性質から直ちに従います。行列式は固有値の積なので
$$ \det(\bm{A} \otimes \bm{B}) = \prod_{i,j} \lambda_i \mu_j = \left(\prod_i \lambda_i\right)^p \left(\prod_j \mu_j\right)^m = (\det \bm{A})^p (\det \bm{B})^m $$
ランク
$$ \text{rank}(\bm{A} \otimes \bm{B}) = \text{rank}(\bm{A}) \cdot \text{rank}(\bm{B}) $$
これは特異値分解を用いるとエレガントに証明できます。
これらの性質は理論的に美しいだけでなく、実際の計算で大きな威力を発揮します。特に vec 演算子と組み合わせたときの応用が重要です。次にその話題に移りましょう。
vec演算子とクロネッカー積
vec演算子の定義
$m \times n$ 行列 $\bm{X}$ のvec演算子 $\text{vec}(\bm{X})$ は、$\bm{X}$ の列を上から下へ順に積み重ねて $mn \times 1$ のベクトルを作る操作です。
$$ \bm{X} = (\bm{x}_1, \bm{x}_2, \ldots, \bm{x}_n) \quad \Rightarrow \quad \text{vec}(\bm{X}) = \begin{pmatrix} \bm{x}_1 \\ \bm{x}_2 \\ \vdots \\ \bm{x}_n \end{pmatrix} $$
vec-クロネッカー積の恒等式
行列方程式を解くうえで中心的な役割を果たすのが、以下の恒等式です。
$$ \begin{equation} \text{vec}(\bm{A}\bm{X}\bm{B}) = (\bm{B}^T \otimes \bm{A})\,\text{vec}(\bm{X}) \end{equation} $$
この式の意味を理解しましょう。左辺は「行列 $\bm{X}$ に左から $\bm{A}$、右から $\bm{B}$ を掛ける」という行列演算です。右辺は「$\bm{X}$ をベクトル化した $\text{vec}(\bm{X})$ にクロネッカー積 $\bm{B}^T \otimes \bm{A}$ を掛ける」という通常の行列-ベクトル積です。つまり、3つの行列が関与する行列方程式を、クロネッカー積を使って普通の連立方程式に変換できるのです。
証明: $\bm{X}$ を列ごとに $\bm{X} = (\bm{x}_1, \ldots, \bm{x}_n)$ と書きます。$\bm{B} = (b_{ij})$ とすると
$\bm{A}\bm{X}\bm{B}$ の第 $j$ 列は
$$ \sum_{k=1}^{n} b_{kj} \bm{A}\bm{x}_k $$
よって
$$ \text{vec}(\bm{A}\bm{X}\bm{B}) = \begin{pmatrix} \sum_k b_{k1}\bm{A}\bm{x}_k \\ \vdots \\ \sum_k b_{kn’}\bm{A}\bm{x}_k \end{pmatrix} $$
一方、$\bm{B}^T \otimes \bm{A}$ は $(j, k)$ ブロックが $b_{kj}\bm{A}$ であるブロック行列なので
$$ (\bm{B}^T \otimes \bm{A})\text{vec}(\bm{X}) = \begin{pmatrix} \sum_k b_{k1}\bm{A}\bm{x}_k \\ \vdots \\ \sum_k b_{kn’}\bm{A}\bm{x}_k \end{pmatrix} $$
両者は一致します。$\square$
特殊ケース
vec-クロネッカー積の恒等式の重要な特殊ケースを挙げます。
$\bm{B} = \bm{I}$ のとき: $\text{vec}(\bm{A}\bm{X}) = (\bm{I} \otimes \bm{A})\text{vec}(\bm{X})$
$\bm{A} = \bm{I}$ のとき: $\text{vec}(\bm{X}\bm{B}) = (\bm{B}^T \otimes \bm{I})\text{vec}(\bm{X})$
これらは行列微分(行列計算のヤコビアン)を求める際に頻繁に使われます。
vec演算子との関係を理解したところで、次にこの道具を使った行列方程式への応用を見ていきましょう。
行列方程式への応用
シルベスター方程式
シルベスター方程式は次の形の行列方程式です。
$$ \begin{equation} \bm{A}\bm{X} + \bm{X}\bm{B} = \bm{C} \end{equation} $$
ここで $\bm{A}$($m \times m$)、$\bm{B}$($n \times n$)、$\bm{C}$($m \times n$)は既知、$\bm{X}$($m \times n$)が未知です。
制御工学のリアプノフ方程式($\bm{B} = \bm{A}^T$ の場合)や、安定性解析で頻繁に現れます。
vec演算子を適用すると
$$ \text{vec}(\bm{A}\bm{X}) + \text{vec}(\bm{X}\bm{B}) = \text{vec}(\bm{C}) $$
先ほどの恒等式を使うと
$$ (\bm{I}_n \otimes \bm{A})\text{vec}(\bm{X}) + (\bm{B}^T \otimes \bm{I}_m)\text{vec}(\bm{X}) = \text{vec}(\bm{C}) $$
$$ (\bm{I}_n \otimes \bm{A} + \bm{B}^T \otimes \bm{I}_m)\text{vec}(\bm{X}) = \text{vec}(\bm{C}) $$
これは $mn \times mn$ の通常の連立方程式です。クロネッカー積のおかげで、行列方程式がベクトルの連立方程式に帰着されました。
ただし、$mn$ が大きい場合にクロネッカー積を陽に構成すると膨大なメモリを消費します。実際には、バートレス・スチュワート法やシュア分解を用いた効率的なアルゴリズムが使われます。
リアプノフ方程式
$\bm{B} = \bm{A}^T$ とした特殊なシルベスター方程式
$$ \bm{A}\bm{X} + \bm{X}\bm{A}^T = \bm{C} $$
はリアプノフ方程式と呼ばれ、制御システムの安定性解析で中心的な役割を果たします。$\bm{X}$ は安定なシステムの定常共分散行列に対応します。
それでは、これらの理論をPythonで実装し、具体的に動作を確認しましょう。
Pythonでの実装
クロネッカー積の基本操作
まず、クロネッカー積の基本的な性質をPythonで確認します。
import numpy as np
import matplotlib.pyplot as plt
# 基本的なクロネッカー積
A = np.array([[1, 2], [3, 4]])
B = np.array([[0, 5], [6, 7]])
# numpy.kron でクロネッカー積を計算
K = np.kron(A, B)
print("A =")
print(A)
print("\nB =")
print(B)
print("\nA ⊗ B =")
print(K)
# 性質の確認
print("\n=== 性質の確認 ===")
# 1. 混合積性質
C = np.array([[2, 1], [0, 3]])
D = np.array([[1, 0], [2, 1]])
lhs = np.kron(A, B) @ np.kron(C, D)
rhs = np.kron(A @ C, B @ D)
print(f"\n混合積性質: ||LHS - RHS|| = {np.linalg.norm(lhs - rhs):.2e}")
# 2. 転置
lhs_t = np.kron(A, B).T
rhs_t = np.kron(A.T, B.T)
print(f"転置: ||LHS - RHS|| = {np.linalg.norm(lhs_t - rhs_t):.2e}")
# 3. 逆行列
lhs_inv = np.linalg.inv(np.kron(A, B))
rhs_inv = np.kron(np.linalg.inv(A), np.linalg.inv(B))
print(f"逆行列: ||LHS - RHS|| = {np.linalg.norm(lhs_inv - rhs_inv):.2e}")
# 4. トレース
tr_kron = np.trace(np.kron(A, B))
tr_prod = np.trace(A) * np.trace(B)
print(f"トレース: tr(A⊗B) = {tr_kron}, tr(A)*tr(B) = {tr_prod}")
# 5. 行列式
det_kron = np.linalg.det(np.kron(A, B))
det_formula = np.linalg.det(A)**2 * np.linalg.det(B)**2 # m=p=2
print(f"行列式: det(A⊗B) = {det_kron:.4f}, det(A)^p * det(B)^m = {det_formula:.4f}")
# 6. 固有値
eig_A = np.linalg.eigvals(A)
eig_B = np.linalg.eigvals(B)
eig_kron = np.sort(np.linalg.eigvals(np.kron(A, B)))
eig_prod = np.sort(np.array([la * lb for la in eig_A for lb in eig_B]))
print(f"\n固有値(クロネッカー積): {np.round(eig_kron, 4)}")
print(f"固有値(積の組合せ): {np.round(eig_prod, 4)}")
上のコードでは、2行2列の行列同士のクロネッカー積を計算し、本記事で解説した全ての性質(混合積性質、転置、逆行列、トレース、行列式、固有値)が数値的に成立していることを確認しています。全ての性質で理論値との差が機械イプシロンのオーダーに収まっており、理論の正しさが裏付けられています。
vec演算子とシルベスター方程式の求解
vec-クロネッカー積の恒等式を使って、シルベスター方程式 $\bm{A}\bm{X} + \bm{X}\bm{B} = \bm{C}$ を解きます。
import numpy as np
from scipy.linalg import solve_sylvester
np.random.seed(42)
m, n = 4, 3
# テスト行列の生成
A = np.random.randn(m, m)
B = np.random.randn(n, n)
X_true = np.random.randn(m, n)
C = A @ X_true + X_true @ B # C = AX + XB
# 方法1: クロネッカー積による直接解法
# vec(AX + XB) = (I ⊗ A + B^T ⊗ I) vec(X) = vec(C)
M = np.kron(np.eye(n), A) + np.kron(B.T, np.eye(m))
x_kron = np.linalg.solve(M, C.flatten("F"))
X_kron = x_kron.reshape(m, n, order="F")
# 方法2: scipy の solve_sylvester
X_scipy = solve_sylvester(A, B, C)
print("=== シルベスター方程式 AX + XB = C ===")
print(f"\n真の解 X:\n{np.round(X_true, 4)}")
print(f"\nクロネッカー積解法:\n{np.round(X_kron, 4)}")
print(f"\nscipy解法:\n{np.round(X_scipy, 4)}")
print(f"\n||X_kron - X_true|| = {np.linalg.norm(X_kron - X_true):.2e}")
print(f"||X_scipy - X_true|| = {np.linalg.norm(X_scipy - X_true):.2e}")
このコードでは、ランダムに生成した真の解 $\bm{X}_{\text{true}}$ から $\bm{C} = \bm{A}\bm{X}_{\text{true}} + \bm{X}_{\text{true}}\bm{B}$ を構成し、2つの方法で $\bm{X}$ を復元しています。クロネッカー積による直接解法とscipyの専用ソルバーの両方が真の解を高精度に再現していることが確認できます。order="F" は列優先(Fortran順序)を指定しており、vec演算子の定義と一致させるために必要です。
応用: クロネッカー積の構造を利用した効率的計算
クロネッカー積を陽に構成せずに、混合積性質を活用して効率的に計算する方法を示します。
import numpy as np
import matplotlib.pyplot as plt
import time
def kron_matvec_naive(A, B, x):
"""(A ⊗ B) x をクロネッカー積を陽に構成して計算(非効率)"""
K = np.kron(A, B)
return K @ x
def kron_matvec_efficient(A, B, x):
"""(A ⊗ B) x を混合積性質で効率的に計算
(A ⊗ B) vec(X) = vec(B X A^T)
ここで X = x を n×p 行列に reshape したもの
"""
p, q = B.shape
m, n = A.shape
X = x.reshape(q, n, order="F") # x を q×n 行列に reshape
Y = B @ X @ A.T # p×m 行列
return Y.flatten(order="F")
# ベンチマーク
sizes = [10, 20, 50, 100, 200]
times_naive = []
times_efficient = []
for s in sizes:
A = np.random.randn(s, s)
B = np.random.randn(s, s)
x = np.random.randn(s * s)
# ナイーブ(クロネッカー積を陽に構成)
t0 = time.perf_counter()
for _ in range(3):
y1 = kron_matvec_naive(A, B, x)
t_naive = (time.perf_counter() - t0) / 3
# 効率的(混合積性質を利用)
t0 = time.perf_counter()
for _ in range(3):
y2 = kron_matvec_efficient(A, B, x)
t_efficient = (time.perf_counter() - t0) / 3
times_naive.append(t_naive)
times_efficient.append(t_efficient)
# 結果の一致を確認
assert np.allclose(y1, y2), f"Results differ at size {s}"
fig, axes = plt.subplots(1, 2, figsize=(14, 5.5))
# (a) 計算時間の比較
ax = axes[0]
ax.loglog(sizes, times_naive, "ro-", linewidth=2, markersize=8,
label="Naive (explicit Kronecker)")
ax.loglog(sizes, times_efficient, "bs-", linewidth=2, markersize=8,
label="Efficient (mixed product)")
n_arr = np.array(sizes, dtype=float)
ax.loglog(n_arr, 1e-8 * n_arr**4, "r:", alpha=0.4, label="$O(n^4)$")
ax.loglog(n_arr, 1e-8 * n_arr**3, "b:", alpha=0.4, label="$O(n^3)$")
ax.set_xlabel("Matrix size n", fontsize=13)
ax.set_ylabel("Time (sec)", fontsize=13)
ax.set_title("Kronecker Product: Naive vs Efficient", fontsize=14)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, which="both")
# (b) 速度比
ax = axes[1]
speedup = [tn / te for tn, te in zip(times_naive, times_efficient)]
ax.bar(range(len(sizes)), speedup, color="steelblue", alpha=0.8)
ax.set_xticks(range(len(sizes)))
ax.set_xticklabels([str(s) for s in sizes])
ax.set_xlabel("Matrix size n", fontsize=13)
ax.set_ylabel("Speedup factor", fontsize=13)
ax.set_title("Speedup: Efficient / Naive", fontsize=14)
ax.grid(True, alpha=0.3, axis="y")
for i, sp in enumerate(speedup):
ax.text(i, sp + 0.2, f"{sp:.1f}x", ha="center", fontsize=11, fontweight="bold")
plt.tight_layout()
plt.savefig("kronecker_benchmark.png", dpi=150, bbox_inches="tight")
plt.show()
このベンチマーク結果から、クロネッカー積の構造を活用した効率的な計算の重要性が読み取れます。
-
左図(計算時間): ナイーブな方法(赤丸)は $O(n^4)$ のスケーリングを示しています。これはクロネッカー積 $\bm{A} \otimes \bm{B}$ の構成に $O(n^4)$ のメモリと計算が必要なためです。一方、効率的な方法(青四角)は $O(n^3)$ のスケーリングであり、行列サイズが大きくなるほど差が顕著です。
-
右図(速度比): $n = 200$ のとき、効率的な方法はナイーブな方法の数十倍以上高速です。これは混合積性質 $(\bm{A} \otimes \bm{B})\text{vec}(\bm{X}) = \text{vec}(\bm{B}\bm{X}\bm{A}^T)$ により、$n^2 \times n^2$ の巨大行列を構成せずに $n \times n$ の行列積2回で済むためです。
応用: クロネッカー積と画像処理
2次元の分離可能フィルタは、クロネッカー積で自然に表現できます。
import numpy as np
import matplotlib.pyplot as plt
# 1次元ガウシアンフィルタ
def gaussian_kernel_1d(size, sigma):
x = np.arange(size) - size // 2
kernel = np.exp(-x**2 / (2 * sigma**2))
return kernel / kernel.sum()
# 分離可能な2Dガウシアンフィルタ = 1Dフィルタのクロネッカー積
sigma = 2.0
size = 11
h_1d = gaussian_kernel_1d(size, sigma)
h_2d = np.kron(h_1d.reshape(-1, 1), h_1d.reshape(1, -1)) # 外積 = 特殊なクロネッカー積
# テスト画像の生成
np.random.seed(42)
img_size = 64
img = np.zeros((img_size, img_size))
# いくつかの点を配置
for _ in range(10):
x, y = np.random.randint(10, img_size - 10, 2)
img[x, y] = 1.0
img += 0.05 * np.random.randn(img_size, img_size)
# 2次元畳み込み(分離可能フィルタ)
from scipy.ndimage import convolve1d
# 分離可能な畳み込み: まず行方向、次に列方向
img_sep = convolve1d(img, h_1d, axis=0)
img_sep = convolve1d(img_sep, h_1d, axis=1)
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
ax = axes[0]
ax.imshow(img, cmap="gray")
ax.set_title("Original Image", fontsize=13)
ax.axis("off")
ax = axes[1]
ax.imshow(h_2d, cmap="hot")
ax.set_title(f"2D Gaussian Kernel (σ={sigma})\n(Kronecker structure)", fontsize=12)
ax.axis("off")
ax = axes[2]
ax.imshow(img_sep, cmap="gray")
ax.set_title("Filtered Image\n(Separable convolution)", fontsize=12)
ax.axis("off")
plt.tight_layout()
plt.savefig("kronecker_image.png", dpi=150, bbox_inches="tight")
plt.show()
この画像処理の例から、クロネッカー積の実用的な意義が読み取れます。2次元ガウシアンフィルタ $\bm{h}_{2D}$ はランク1行列(1次元カーネル $\bm{h}_{1D}$ の外積)として表現でき、これはクロネッカー積の特殊なケースです。分離可能フィルタの利点は、$n \times n$ の画像に対して $O(n^2 k)$ の計算量で済むことです($k$ はカーネルサイズ)。2次元畳み込みを直接行うと $O(n^2 k^2)$ になるため、分離可能性を利用することで計算量が $k$ 倍削減されます。
まとめ
本記事では、クロネッカー積の定義と性質を包括的に解説しました。
- クロネッカー積 $\bm{A} \otimes \bm{B}$ は行列の各成分をブロックで置き換える演算であり、独立なシステムの組み合わせを記述する
- 混合積性質 $(\bm{A} \otimes \bm{B})(\bm{C} \otimes \bm{D}) = (\bm{A}\bm{C}) \otimes (\bm{B}\bm{D})$ はクロネッカー積の最も強力な性質であり、効率的な計算の鍵である
- 固有値: $\bm{A} \otimes \bm{B}$ の固有値は $\bm{A}$ と $\bm{B}$ の固有値の全ての積の組み合わせ $\{\lambda_i\mu_j\}$ である
- vec演算子との組み合わせ $\text{vec}(\bm{A}\bm{X}\bm{B}) = (\bm{B}^T \otimes \bm{A})\text{vec}(\bm{X})$ により、行列方程式が通常の連立方程式に帰着される
- シルベスター方程式 $\bm{A}\bm{X} + \bm{X}\bm{B} = \bm{C}$ はクロネッカー積を用いてベクトル化できる
次のステップとして、以下の記事も参考にしてください。
- 行列の微分公式を完全まとめ — vec演算子とクロネッカー積の微分への応用
- 固有値と固有ベクトル — 固有値分解の基本理論
- ムーア・ペンローズ擬逆行列の理論 — SVDとの関連