Rでクロネッカー積

クロネッカー積とは

クロネッカー積(Kronecker product)は、2つの行列から、より大きな1つの行列を作る演算であり、記号は \(\otimes\) で表されます。

\(m\) 行 \(n\) 列の行列 \(A = (a_{ij})\) と、\(p\) 行 \(q\) 列の行列 \(B = (b_{kl})\) に対して、クロネッカー積 \(A \otimes B\) は、\(A\) の各要素 \(a_{ij}\) を「\(a_{ij}\) 倍した \(B\)」に置き換えたブロック行列として定義されます。

\[
A \otimes B =
\begin{pmatrix}
a_{11} B & a_{12} B & \cdots & a_{1n} B \\
a_{21} B & a_{22} B & \cdots & a_{2n} B \\
\vdots & \vdots & \ddots & \vdots \\
a_{m1} B & a_{m2} B & \cdots & a_{mn} B
\end{pmatrix}
\]

ここで、\(a_{ij}\) は \(A\) の第 \(i\) 行第 \(j\) 列の要素、\(B\) は \(p \times q\) 行列です。

結果の行列は \(mp\) 行 \(nq\) 列になります。

要素ごとに書きますと、次のとおりになります。

\[
\begin{align}
&(A \otimes B)_{(i-1)p+k,\;(j-1)q+l} = a_{ij} \, b_{kl} \\
&i = 1, \ldots, m \quad j = 1, \ldots, n \quad k = 1, \ldots, p \quad l = 1, \ldots, q
\end{align}
\]

ここで、\((i-1)p+k\) は結果の行番号、\((j-1)q+l\) は結果の列番号です。

つまり、\(A\) の位置 \((i, j)\) が「ブロックの位置」を決め、\(B\) の位置 \((k, l)\) が「ブロックの中の位置」を決めています。

主な性質

混合積の性質

\(A\) は \(m \times n\)、\(B\) は \(p \times q\)、\(C\) は \(n \times r\)、\(D\) は \(q \times s\) の行列。

\[(A \otimes B)(C \otimes D) = (AC) \otimes (BD)\]

転置の性質

\[(A \otimes B)^{\top} = A^{\top} \otimes B^{\top}\]

行列式の性質

\(A\) は \(m \times m\)、\(B\) は \(p \times p\) の正方行列。

\[\det(A \otimes B) = \det(A)^{p} \, \det(B)^{m}\]

トレースの性質

\(A\) と \(B\) はともに正方行列。

\[\mathrm{tr}(A \otimes B) = \mathrm{tr}(A) \, \mathrm{tr}(B)\]

vec作用素の性質

\(X\) は \(n \times p\) 行列、\(B\) は \(p \times q\) 行列であり、\(\mathrm{vec}(\cdot)\) は行列を列ごとに縦に積み重ねてベクトルにする操作を表します。

\[\mathrm{vec}(AXB) = (B^{\top} \otimes A) \, \mathrm{vec}(X)\]

その他の性質

  • \(A\) の固有値を \(\lambda_i\) (\(i = 1, \ldots, m\))、\(B\) の固有値を \(\mu_j\) (\(j = 1, \ldots, p\)) とすると、\(A \otimes B\) の固有値は、すべての組合せの積 \(\lambda_i \mu_j\) になります。
  • クロネッカー積は交換できません。一般に \(A \otimes B \neq B \otimes A\) となります。ただし、両者は要素の並びが異なるだけで、同じ要素から成ります。

クロネッカー積は、多変量回帰や見かけ上無関係な回帰(SUR)、時系列モデルの行列表現など、行列を「ベクトルにして」扱う場面でよく使われます。

Rコード

Rでは kronecker() 関数、または演算子 %x% でクロネッカー積を求められます。

以下では、定義に沿った自作関数と結果を比べながら、構造と性質を順に確認します。

小さな例で結果を確認

\(A\) を \(2 \times 2\) 行列、\(B\) を \(2 \times 3\) 行列とします。結果は \(4 \times 6\) 行列になるはずです。

# 行列A(2行2列)を定義する
mat_a <- matrix(c(
  1, 2,
  3, 4
), nrow = 2, byrow = TRUE)

# 行列B(2行3列)を定義する
mat_b <- matrix(c(
  1, 0, 2,
  0, 1, 3
), nrow = 2, byrow = TRUE)

# クロネッカー積を求める
kron_ab <- kronecker(mat_a, mat_b)

mat_a
     [,1] [,2]
[1,]    1    2
[2,]    3    4
mat_b
     [,1] [,2] [,3]
[1,]    1    0    2
[2,]    0    1    3
kron_ab
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    0    2    2    0    4
[2,]    0    1    3    0    2    6
[3,]    3    0    6    4    0    8
[4,]    0    3    9    0    4   12
# 大きさを確認する(mp行nq列になっているはず)
dim(kron_ab)
[1] 4 6

%x% でも同じ結果になります。

# 演算子 %x% は kronecker() と同じ
all.equal(mat_a %x% mat_b, kron_ab)
[1] TRUE

ブロック構造を確認

結果の左上のブロックは \(a_{11} B\)、右上のブロックは \(a_{12} B\) と等しくなります。

ブロックの位置と、元の行列の要素との対応を確かめます。

# Bの行数と列数を取り出す
nrow_b <- nrow(mat_b)
ncol_b <- ncol(mat_b)

# A[1, 2] のブロック(第1行・第2列のブロック)を取り出す
rows_12 <- (1 - 1) * nrow_b + seq_len(nrow_b)
cols_12 <- (2 - 1) * ncol_b + seq_len(ncol_b)
block_12 <- kron_ab[rows_12, cols_12]

block_12
     [,1] [,2] [,3]
[1,]    2    0    4
[2,]    0    2    6
mat_a[1, 2] * mat_b
     [,1] [,2] [,3]
[1,]    2    0    4
[2,]    0    2    6
# ブロックが a_12 * B に一致するか
all.equal(block_12, mat_a[1, 2] * mat_b)
[1] TRUE

要素の対応式 \((A \otimes B)_{(i-1)p+k,\,(j-1)q+l} = a_{ij} b_{kl}\) も、1つの要素で確かめておきます。

# i = 2, j = 1, k = 1, l = 3 の場合
i <- 2
j <- 1
k <- 1
l <- 3

# 結果の行番号と列番号を計算する
row_idx <- (i - 1) * nrow(mat_b) + k
col_idx <- (j - 1) * ncol(mat_b) + l

kron_ab[row_idx, col_idx]
[1] 6
mat_a[i, j] * mat_b[k, l]
[1] 6

定義どおりに自作する

定義(「\(A\) の各要素を \(a_{ij} B\) に置き換える」)をそのままコードにした関数を作り、kronecker() と比べます。

# 定義に沿ったクロネッカー積の自作関数
kron_manual <- function(mat_a, mat_b) {
  nrow_a <- nrow(mat_a)
  ncol_a <- ncol(mat_a)
  nrow_b <- nrow(mat_b)
  ncol_b <- ncol(mat_b)

  # 結果を入れる行列(mp行nq列)を用意する
  result <- matrix(0, nrow = nrow_a * nrow_b, ncol = ncol_a * ncol_b)

  for (i in seq_len(nrow_a)) {
    for (j in seq_len(ncol_a)) {
      # ブロック(i, j)が占める行と列の範囲を求める
      rows <- (i - 1) * nrow_b + seq_len(nrow_b)
      cols <- (j - 1) * ncol_b + seq_len(ncol_b)

      # ブロックに a_ij * B を入れる
      result[rows, cols] <- mat_a[i, j] * mat_b
    }
  }

  result
}

kron_manual(mat_a, mat_b)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    0    2    2    0    4
[2,]    0    1    3    0    2    6
[3,]    3    0    6    4    0    8
[4,]    0    3    9    0    4   12
# 組み込み関数と一致するか
all.equal(kron_manual(mat_a, mat_b), kronecker(mat_a, mat_b))
[1] TRUE

性質を数値で確かめる

ここでは、\(A\) と \(B\) を対称行列(\(A\) は \(2 \times 2\)、\(B\) は \(3 \times 3\))、\(C\) と \(D\) を対称でない行列(\(C\) は \(2 \times 2\)、\(D\) は \(3 \times 3\))とします。

# 対称行列 A (2行2列)
sq_a <- matrix(c(
  2, 1,
  1, 3
), nrow = 2, byrow = TRUE)

# 対称行列 B (3行3列)
sq_b <- matrix(c(
  2, 0, 1,
  0, 3, 0,
  1, 0, 4
), nrow = 3, byrow = TRUE)

# 対称でない行列 C (2行2列)
sq_c <- matrix(c(
  1, 2,
  0, 1
), nrow = 2, byrow = TRUE)

# 対称でない行列 D (3行3列)
sq_d <- matrix(c(
  1, 0, 1,
  2, 1, 0,
  0, 1, 1
), nrow = 3, byrow = TRUE)

混合積の性質 \((A \otimes B)(C \otimes D) = (AC) \otimes (BD)\)

lhs <- kronecker(sq_a, sq_b) %*% kronecker(sq_c, sq_d)
rhs <- kronecker(sq_a %*% sq_c, sq_b %*% sq_d)

all.equal(lhs, rhs)
[1] TRUE

転置の性質 \((A \otimes B)^{\top} = A^{\top} \otimes B^{\top}\)

# 対称行列だと転置しても変わらないため、対称でないCとDで確かめます
all.equal(t(kronecker(sq_c, sq_d)), kronecker(t(sq_c), t(sq_d)))
[1] TRUE

行列式の性質 \(\det(A \otimes B) = \det(A)^{p} \det(B)^{m}\)(ここでは \(m = 2\)、\(p = 3\))

m <- nrow(sq_a)
p <- nrow(sq_b)

det(kronecker(sq_a, sq_b))
[1] 55125
det(sq_a)^p * det(sq_b)^m
[1] 55125

トレースの性質 \(\mathrm{tr}(A \otimes B) = \mathrm{tr}(A)\,\mathrm{tr}(B)\)

sum(diag(kronecker(sq_a, sq_b)))
[1] 45
sum(diag(sq_a)) * sum(diag(sq_b))
[1] 45

固有値の性質 \(A \otimes B\) の固有値は \(\lambda_i \mu_j\) になる。

# 対称行列ですので、固有値はすべて実数になる
eigen_a <- eigen(sq_a, symmetric = TRUE)$values
eigen_b <- eigen(sq_b, symmetric = TRUE)$values

# クロネッカー積の固有値(大きい順)
eigen_ab <- eigen(kronecker(sq_a, sq_b), symmetric = TRUE)$values

# 固有値どうしのすべての積(大きい順)
eigen_prod <- sort(c(outer(eigen_a, eigen_b)), decreasing = TRUE)

round(eigen_ab, 6)
[1] 15.970775 10.854102  6.100293  5.737429  4.145898  2.191503
round(eigen_prod, 6)
[1] 15.970775 10.854102  6.100293  5.737429  4.145898  2.191503
all.equal(eigen_ab, eigen_prod)
[1] TRUE

非可換 \(A \otimes B \neq B \otimes A\)

kron_ab_sq <- kronecker(sq_a, sq_b)
kron_ba_sq <- kronecker(sq_b, sq_a)

# 行列としては一致しない
isTRUE(all.equal(kron_ab_sq, kron_ba_sq))
[1] FALSE
# しかし、要素を並べ替えれば同じ集合になる
all.equal(sort(c(kron_ab_sq)), sort(c(kron_ba_sq)))
[1] TRUE

ベクトルのクロネッカー積と外積の関係

ベクトルを「1列の行列」とみなした場合、ベクトル同士のクロネッカー積は、外積を列方向に並べたものと一致します。

vec_a <- c(1, 2, 3)
vec_b <- c(1, 10)

# ベクトルのクロネッカー積
kronecker(vec_a, vec_b)
[1]  1 10  2 20  3 30
# 外積を作り、列ごとに縦に並べたもの
outer(vec_b, vec_a)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]   10   20   30
c(outer(vec_b, vec_a))
[1]  1 10  2 20  3 30

vec 作用素との関係

\(\mathrm{vec}(AXB) = (B^{\top} \otimes A)\,\mathrm{vec}(X)\) を確かめます。

Rの as.vector() は列ごとに縦に並べますので、\(\mathrm{vec}\) そのものになります。

# A (2行3列)、X (3行2列)、B (2行4列)
mat_l <- matrix(c(
  1, 2, 0,
  0, 1, 3
), nrow = 2, byrow = TRUE)
mat_x <- matrix(c(
  1, 2,
  3, 4,
  5, 6
), nrow = 3, byrow = TRUE)
mat_r <- matrix(c(
  1, 0, 2, 1,
  3, 1, 0, 2
), nrow = 2, byrow = TRUE)

# 左辺: vec(AXB)
lhs_vec <- as.vector(mat_l %*% mat_x %*% mat_r)

# 右辺: (B^T ⊗ A) vec(X)
rhs_vec <- as.vector(kronecker(t(mat_r), mat_l) %*% as.vector(mat_x))

lhs_vec
[1] 37 84 10 22 14 36 27 62
rhs_vec
[1] 37 84 10 22 14 36 27 62
all.equal(lhs_vec, rhs_vec)
[1] TRUE

こちらの関係を利用しますと、行列に対する線形な演算を「ベクトルへの行列の掛け算」に書き換えられます。

クロネッカー積を繰り返して図形を作る

クロネッカー積を同じ行列で繰り返しますと、全体と部分が相似になる模様が現れます。

\(P = \begin{pmatrix} 1 & 0 \\ 1 & 1 \end{pmatrix}\) を \(k\) 回繰り返したクロネッカー積 \(P^{\otimes k}\) を描画してみます。

ここで、\(P^{\otimes k}\) は \(2^k \times 2^k\) 行列です。

library(ggplot2)

# 基本となる行列P
mat_p <- matrix(c(
  1, 0,
  1, 1
), nrow = 2, byrow = TRUE)

# 同じ行列をk回クロネッカー積する関数
kron_power <- function(mat, k) {
  Reduce(kronecker, rep(list(mat), k))
}

# 行列をggplot2で描ける縦長のデータに変換する関数
to_long <- function(mat) {
  data.frame(
    row = rep(seq_len(nrow(mat)), times = ncol(mat)),
    col = rep(seq_len(ncol(mat)), each = nrow(mat)),
    value = factor(as.vector(mat), levels = c(0, 1))
  )
}

# k = 5 の場合(32行32列)
mat_p5 <- kron_power(mat_p, 5)
dim(mat_p5)
[1] 32 32
to_long(mat_p5) |>
  ggplot(aes(x = col, y = row, fill = value)) +
  geom_tile() +
  scale_y_reverse() +
  scale_fill_manual(
    values = c("0" = "white", "1" = "steelblue"),
    labels = c("0" = "0(空白)", "1" = "1(塗りつぶし)")
  ) +
  coord_equal() +
  labs(
    title = "行列Pを5回クロネッカー積した結果",
    x = "列番号",
    y = "行番号",
    fill = "要素の値"
  ) +
  theme_minimal()
Figure 1

1回繰り返すごとに、行列全体が「\(P\) の形に並んだ、1つ前の行列のブロック」になります。

クロネッカー積が、各要素をブロックに置き換える操作であることが、イメージ出来ます。

以上です。