2つの部分空間の共通部分の基底を求める

はじめに

しばらくブログを更新してなかったので線形代数のネタを1つ。やりたいのは線形独立な $m$ 個のベクトル ${\bf u}_{1}, \cdots, {\bf u}_{m} \in R^{d}$ が張る部分空間と線形独立な $n$ 個のベクトル ${\bf v}_{1}, \cdots, {\bf v}_{n} \in R^{d}$ が張る部分空間の共通部分 ${\rm span} \left\{ {\bf u}_{1}, \cdots, {\bf u}_{m} \right\} \cap {\rm span} \left\{ {\bf v}_{1}, \cdots, {\bf v}_{n} \right\}$ の基底を求めること。

以降、行列 ${\bf U}$ を ${\bf U} = \left[ {\bf u}_{1}, \cdots, {\bf u}_{m} \right]$ 、行列 ${\bf V}$ を ${\bf V} = \left[ {\bf v}_{1}, \cdots, {\bf v}_{n} \right]$ とする。

方法1：零空間を使う方法

まず、 ${\rm span} \left\{ {\bf u}_{1}, \cdots, {\bf u}_{m} \right\}$ の任意の元は適当なベクトル ${\bf a} \in R^{m}$ を使って ${\bf U} {\bf a}$ と表せる。同様に ${\rm span} \left\{ {\bf v}_{1}, \cdots, {\bf v}_{n} \right\}$ の任意の元は適当なベクトル ${\bf b} \in R^{n}$ を使って ${\bf V} {\bf b}$ と表せる。共通部分の元は、

$\displaystyle {\bf U} {\bf a} = {\bf V} {\bf b}$

を満たす。これを変形すると、

$\displaystyle \begin{aligned} {\bf U} {\bf a} - {\bf V} {\bf b} &= &{\bf 0} \\ \begin{bmatrix} {\bf U} & {\bf V} \end{bmatrix} \begin{bmatrix} {\bf a} \\ -{\bf b} \end{bmatrix} & = & {\bf 0} \\ {\bf W} {\bf c} &=& {\bf 0} \end{aligned}$

最後の式は ${\bf W} = \left[{\bf U} \ \ {\bf V}\right]$ 、 ${\bf c} = \left[ {\bf a}^T \ \ -{\bf b}^T \right]^T$ とおいた。最後の式より、この条件を満たす ${\bf c}$ は ${\bf W}$ の零空間 $\ker ({\bf W})$ の元である。従って、この $\ker ({\bf W})$ の基底を求め、それぞれの基底について上から $m$ 次元分取り出し（これを ${\hat {\bf a}}_1, \cdots, {\hat {\bf a}}_k$ とする）、 ${\bf U}{\hat {\bf a}}_1,, \cdots, {\bf U}{\hat {\bf a}}_k$ とすれば共通部分の基底が求まる。

これをpythonで実装する場合は、scipy.linalg.null_spaceを使えば一発で零空間の基底を求められるのでこれを活用すればよい。ちなみに零空間の基底は、行列を特異値分解してその特異値ゼロに対応する右特異ベクトルとして得られる。実際、scipy.linalg.null_spaceも中では特異値分解を使ってるとのこと。これを試してみたコードは以下の通り（コメントにある通りこのコードの実験だと共通部分は3次元になるはず*1）。

import numpy as np
from scipy import linalg

def make_test_matrix(d, m, rand_seed = 1):
    np.random.seed(rand_seed)

    #Uはランダムに設定
    U = np.random.randn(d, m)
    #Vの2つのベクトルはランダム
    v1 = np.random.randn(d, 1)
    v2 = np.random.randn(d, 1)
    #他の3つは既存のUとVから構成(これで共通部分は3次元になるはず)
    v3 = ((7 * U[:, 0]) + (5 * U[:, 1]))[:, np.newaxis]
    v4 = ((1 * U[:, 1]) + (2 * U[:, 3]))[:, np.newaxis]
    v5 = ((3 * U[:, 3]) + (-2 * v1[:, 0]))[:, np.newaxis]
    V = np.concatenate([v1, v2, v3, v4, v5] , axis = 1)
    return U, V

if __name__ == '__main__':
    d = 15
    m = 5

    #テストデータの生成
    U, V = make_test_matrix(d, m, 1)

    #UとVを連結した行列の零空間の基底を計算
    W = np.concatenate([U, V], axis=1)
    C = linalg.null_space(W)
    #Uの係数のみ取り出して基底を計算
    A_hat = C[:m, :]
    basis = U @ A_hat

    #結果の表示
    print("A_hat")
    print(A_hat)
    print("basis of intersection")
    print(basis)

結果

A_hat
[[-7.16152592e-01  4.38912368e-01  1.34809890e-01]
 [-6.53812424e-01 -2.10948597e-01 -7.03781271e-02]
 [-1.63931368e-16 -2.84223600e-16 -1.73499453e-16]
 [-1.05903826e-01 -6.16391057e-01  6.45398134e-01]
 [-3.74700271e-16  5.09575021e-17 -3.66568755e-16]]
basis of intersection
[[-0.64967372  1.50336269 -0.43045986]
 [ 0.47368583 -1.57489224 -0.22715953]
 [ 0.34052603  1.31304872  0.09422749]
 [ 0.89595509 -0.47240257 -0.10889637]
 [-0.01343938 -1.03428625  0.09537105]
 [ 0.59837101 -0.10904928 -0.25641896]
 [ 0.84424754  0.30109334 -0.61081408]
 [ 0.56380177 -0.79295082  1.14816012]
 [ 0.53848895 -0.94016879  1.12891538]
 [ 0.31863787 -0.39392353 -0.02175922]
 [ 0.0523343   0.42138681 -0.16020853]
 [-0.99889429 -0.09553967  0.20435405]
 [-0.24730755 -0.41166361 -0.38226338]
 [-0.91721327 -1.61905413  1.32073652]
 [ 1.27124225 -1.06748573  0.40630033]]

期待通り3つの基底ベクトルが得られた。

方法2：正準角を使う方法

たぶん上で書いた方法がスタンダードなやり方だと思うが、次のような方法でも求められる。簡単のため各列ベクトルはそれぞれが張る部分空間の正規直交基底とする。つまり、 ${\bf U}^T {\bf U} = {\bf I}, \ \ {\bf V}^T {\bf V} = {\bf I}$ 。もし与えられているベクトルが正規直交じゃない場合でもQR分解を使って直交化すればよい*2。

今、両部分空間の元 ${\bf U} {\bf a}$ と ${\bf V} {\bf b}$ のなす角 $\theta$ を考えると

$\displaystyle \cos \theta = \frac{{\bf a}^T {\bf U}^T {\bf V} {\bf b}}{\parallel {\bf U} {\bf a} \parallel_2^{2} \parallel {\bf V} {\bf b} \parallel_2^{2}} = \frac{{\bf a}^T {\bf U}^T {\bf V} {\bf b}}{ ({\bf a}^T {\bf a}) ( {\bf b}^T {\bf b} )}$

ここで、 ${\bf U}$ と ${\bf V}$ に共通部分があるなら、上式の $\cos \theta$ を $1$ （つまりなす角を $0$ ）にできる ${\bf a}, {\bf b}$ が存在するはずであり、それを ${\hat {\bf a}}_1, {\hat {\bf b}}_1$ とすると ${\bf U}{\hat {\bf a}}_1$ と ${\bf V}{\hat {\bf b}}_1$ は共通部分の元である。さらに、 ${\hat {\bf a}}_1, {\hat {\bf b}}_1$ と直交して、かつ、上式 $\cos \theta$ を $1$ にできる ${\hat {\bf a}}_2, {\hat {\bf b}}_2$ が存在するなら、 ${\bf U}{\hat {\bf a}}_2, {\bf V}{\hat {\bf b}}_2$ も共通部分の元である*3。同様にして、互いに直交し上記 $\cos \theta$ を $1$ にできる ${\bf a}, {\bf b}$ が $k$ 組み見つかったとすると、 ${\bf U}{\hat {\bf a}}_1, \cdots, {\bf U}{\hat {\bf a}}_k$ が共通部分の基底となる。このなす角 $\theta$ は正準角（principal angle）と呼ばれ、このブログでも過去に少しだけ触れている（9年越しでタイトルの疑問に回答・・・）。

yamagensakam.hatenablog.com

この記事の内容より、結局のところ ${\bf U}^T {\bf V}$ を特異値分解して、その特異値が $1$ に対応する左特異ベクトルが上記の ${\hat {\bf a}}_1, \cdots, {\hat {\bf a}}_k$ に当たるため、後はそれぞれ左から ${\bf U}$ をかければ良い。

こちらも以下のコードで実験してみた（テストデータの生成関数make_test_matrixは方法1のコードと同じ）。

if __name__ == '__main__':
    d = 15
    m = 5
    #テストデータの生成
    U, V = make_test_matrix(d, m)

    #直交化
    UQ , _ = np.linalg.qr(U)
    VQ , _ = np.linalg.qr(V)

    #特異値分解でcos(theta)と係数a, bを計算
    A_hat, D, B_hat = np.linalg.svd(UQ.T @ VQ)

    #特異値が1の左特異ベクトルを使って共通部分の基底ベクトルを計算
    A_hat = A_hat[:, np.abs(D - 1.0) < 1e-15]
    basis = UQ @ A_hat

    #結果の表示
    print("singular values")
    print(D)
    print("A_hat")
    print(A_hat)
    print("basis of intersection")
    print(basis)

結果

singular values
[1.         1.         1.         0.90457795 0.29262578]
A_hat
[[ 6.44393662e-01 -1.84065423e-01 -7.42210703e-01]
 [ 6.08863932e-01  7.10720557e-01  3.52364871e-01]
 [-2.07087630e-01  3.03916180e-01 -2.55165300e-01]
 [ 4.13710326e-01 -6.07150036e-01  5.09757726e-01]
 [ 2.49800181e-16 -9.54097912e-16  8.53483950e-16]]
basis of intersection
[[-0.37542124  0.14849663  0.16837461]
 [ 0.22018893 -0.18075024 -0.51658548]
 [ 0.00558178  0.36511633  0.39239376]
 [ 0.26960343  0.1627379  -0.16594012]
 [ 0.0955871  -0.22341439 -0.22149631]
 [ 0.13222487  0.16072238 -0.14863213]
 [ 0.10764232  0.333674   -0.21546023]
 [ 0.41641461 -0.06119148  0.37256877]
 [ 0.41790505 -0.09802718  0.32441074]
 [ 0.11692351  0.01010546 -0.1092835 ]
 [-0.04649923  0.11043067  0.03039536]
 [-0.23646815 -0.31480919  0.06458789]
 [-0.10098034 -0.13695922 -0.30038741]
 [ 0.09864014 -0.66361076  0.22550194]
 [ 0.50983159  0.12044854 -0.0594418 ]]

こちらも3つの特異値が1になり、期待通り3次元分の基底が得られた。

*1:ランダムでたまたま線形従属なベクトルが生成されなければ

*2:直交化しなくても、そのまま共通部分の基底を求めることも可能だが、少し計算が煩雑になるのでここでは書かない。このあたりについては、例えば以下の書籍の6章なんかに書いてあるので興味あればそちらを参照。

入門機械学習による異常検知―Rによる実践ガイド

作者:剛, 井手
発売日: 2015/02/19
メディア: 単行本

*3: ${\bf U}$ が直交行列なので ${\hat {\bf a}}_1\perp {\hat {\bf a}}_2$ ならば ${\bf U}{\hat {\bf a}}_1 \perp {\bf U}{\hat {\bf a}}_2$ 。同様に ${\bf V}$ が直交行列なので ${\hat {\bf b}}_1\perp {\hat {\bf b}}_2$ ならば ${\bf V}{\hat {\bf b}}_1 \perp {\bf V}{\hat {\bf b}}_2$ 。