対称クロネッカー積におけるvecからsvecに変換する行列Qの実装について

0. はじめに

「21世紀の線形計画」といわれる半正定値計画(SDP)を実装し理解を深めようと思います。多くの教科書では半正定値計画の理論的なところを触れているので、実装例を備忘録の連載という形で書いていきます。21世紀の皆様に少しでもお役に立てれば幸いです。

1. vecからsvecに変換する行列Q

$vec$から$svec$に変換する行列Qの実装でつまづいたので備忘録として残す。

変換行列Qの用途をわかりやすいように、2x2の実対称行列$A$を実例として説明する。
実対称行列なので$a_{21}=a_{12}$である。

A = \left(
    \begin{array}{cc}
      a_{11} & a_{12}  \\
      a_{21} & a_{22} 
    \end{array}
  \right)

ここで、実対称行列$A$を$vec$作用素でベクトル化すると$vec(A)$をえる。

vec(A) := (a_{11},a_{21},a_{12},a_{22})^T

実対称行列の対称的な構造に着目し、変数を減らすことができる$svec$を導入する。例えば、2x2実対称行列において、$vec$の場合変数は4であるが、svecの場合変数が3と減らせることがわかる。

svec(A) := (a_{11},\sqrt{2} a_{21},a_{22})^T

$vec$と$svec$の内積が等しいことを確認するため、vecの内積を下記に示す。

trace(AA^T)=vec(A)^Tvec(A)=a_{11}^2 + a_{21}^2 + a_{12}^2 + a_{22}^2 = a_{11}^2 + 2a_{21}^2 + a_{22}^2 

$svec$の場合でも下記のように$vec$と同じ内積値を得られることがわかる。

svec(A)^Tsvec(A)= a_{11}^2 + 2a_{21}^2 + a_{22}^2 

このように$svec$が次元削減において便利であることがわかったと思う。大規模な実対称行列のサイズであれば、より演算回数が減らすうえで有用であることが予想できる。

便利さがわかれば$vec$から$svec$に変換する変換行列$Q$を生成する関数が欲しくなると思う。
2x2の実対称行列の場合の変換行列$Q$は下記で与えられる。

Qvec(A) = svec(A)

Q^Tsvec(A) = vec(A)

Q = \left(
    \begin{array}{cccc}
      1 & 0 & 0 & 0 \\
      0 & \sqrt{2} & \sqrt{2} & 0 \\
      0 & 0 & 0 & 1
    \end{array}
  \right)

2. 実装

任意の実対称行列における変換行列$Q$の実装例を下記に示す。

import numpy as np

def set_vec2svec_matrix(n):
  Q = np.zeros((np.int(n*(n+1)/2) , np.int(n**2)))

  for k in range(n):
      for l in range(n):
          if k == l:
              i = np.int(n*l -l*(l+1)/2 + k)
              j = np.int(n*l + k)
              Q[i,j] = 1
          elif k > l:
              i = np.int(n*l -l*(l+1)/2 + k )
              j = np.int(n*l + k)
              Q[i,j] = 1/np.sqrt(2)
          elif k < l:
              i = np.int(n*k -k*(k+1)/2 + l)
              j = np.int(n*l + k)
              Q[i,j] = 1/np.sqrt(2) 
  return Q

2x2の実対称行列における変換行列Qにおける実行方法を下記に示す。

$ n = 2
$ Q = set_vec2svec_matrix(n)

#[[1.         0.         0.         0.        ]
# [0.         0.70710678 0.70710678 0.        ]
# [0.         0.         0.         1.        ]]

参照

• On the Kronecker