Mixtureモデルの周辺尤度をモンテカルロ積分で近似してみる
はじめに
この記事では、Mixtureモデルの周辺尤度をモンテカルロ積分で近似するアルゴリズムChib's methodを試します。 Mixtureモデルの周辺尤度は簡単には求まらない(数値的にも)という印象だったので、サンプリングを使って計算する手法がどこまでうまく行くのか検証しようと思いました。
経緯
Non-local prior というモデル選択のための事前分布の紹介qiitaを読んで、興味を持ちNon-local priorの論文も読みました。この論文はクラスタのダブりや空のクラスタが現れる事前確率を0にすることで、周辺尤度を使ったモデル選択(この場合はクラスタ数選択)の性能を上げるというものでした。この論文は面白かったのですが、周辺尤度の計算には疑問が残りました。見てみると周辺尤度の計算にはこの論文(Marin & Robert, 2008)の方法を使っており、この方法はこの論文内でChib's methodと呼ばれる方法をベースとしています。そこで、Chib's methodからMarin & Robertの方法からNon-local priorまですべて計算機実験でためそう!となったわけです。
そして、Chib's methodの検証までで飽きたのでそこまでを記事にすることにしました。
想定するモデル
想定するモデルは以下のようなものです。
- Mixtureモデル
- サンプルはi.i.d.
- パラメタの事前分布はconjugate
Chib's method
まず、Chib's methodを説明します。ベイズの定理より、
が成り立ちます。ここで、はデータ、はパラメタを指します。今回欲しいのは周辺尤度なので、事後分布さえ計算できればOKです。しかも左辺はに依存しないので、として任意の値を代入した確率密度が計算できればよいです。簡単簡単。と、思いきや、肝心なことを忘れています。そうですクラスタを示す潜在変数です。に含めて考えることもできますが、潜在変数とパラメタの同時分布の確率密度を得る必要が出てきて苦しいです。
そこで、潜在変数を周辺化して上式を得ることを考えます。とは簡単です。前者はの事前分布の割合で按分するだけで、後者はそもそも必要ないです。問題はです。簡単には求まりません。
解決法は単純です。MCMCによってのサンプルを得てそれを用いて近似して評価します。すなわち、
と近似します。ここで、はMCMCによる番目のの事後サンプルを示し、はサンプリング回数を示します。 これによって周辺尤度が近似できます。この方法がChib's methodです。 また、Marin & Robertの方法ではMCMCにおけるlabel switchingの問題に対処するためを各クラスタについてパーミュテーションした値全通りで計算して平均を取っています。 今回はとして各クラスタで同じ値を使っています。いい値を選ばないと精度が悪くなるとのことですが、今回は気にしていません。
検証実験
ここからは検証実験の内容になります。
モデル
用いるモデルは超簡単なGaussian mixture modelです。具体的な生成過程は、
です(がになってます)。厳密には記号の定義をするべきですが、しません。ハートで理解してください。
正解のための解析的な計算
ここからが問題です。Chib's methodを検証するためには正確な周辺尤度を計算しないといけません。しかし今回対象にできるモデルはmixtureモデルのみなのでそれができたら苦労しない...。なので、サンプルサイズを小さくし、による積分を厳密に計算できるデータで検証を行いました。が計算できるなら正確な周辺尤度を出せるはず...!つまりは、上の周辺事後分布の
の部分を計算します。 まずです。こちらは、単峰分布の場合の事後分布の計算とほぼ変わりません。すなわち、
のように計算できます。右辺はクラスタkに属するデータについてだけ(クラスタの分け方すなわちは与えられているので可能)単峰の場合の事後分布を計算することになるので、が共役事前分布であることを使って正規化し左辺を求められます。 今回の場合では、
となります。ただし、はクラスタkに属するデータ数です。 を計算するコードは以下になります。
def calc_ln_alpha(x, z, mu_star, D, K): ln_alpha = 0 for k in range(K): x_k = x[z==k] N_k = x_k.shape[0] if N_k == 0: ln_alpha += multivariate_normal.logpdf( mu_star, mean=np.zeros(D), cov=np.identity(D)) continue post_mu_ex = x_k.sum(0) / (N_k + 1) post_mu_cov = np.identity(D) / (N_k + 1) ln_alpha += multivariate_normal.logpdf( mu_star, mean=post_mu_ex, cov=post_mu_cov) return ln_alpha
次にです。ベイズの定理から、
が成り立ちます。ここで、
となり、kに所属するxについてだけ注目したときの、周辺尤度を計算すればよいことがわかります。これは単峰の場合の周辺尤度なので計算できます。具体的には、単峰を考えてを省略するとして、
def calc_ln_betas(x, N, D, K): ln_betas = [] for i in range(pow(K, N)): this_z = [] cluster_str = None if K == 1: cluster_str = ''.zfill(N) else: cluster_str = np.base_repr(i, K).zfill(N) for s in cluster_str: this_z.append(int(s)) this_z = np.array(this_z) ln_beta = 0 for k in range(K): x_k = x[this_z==k] N_k = x_k.shape[0] if N_k == 0: continue ln_betak = - D * N_k * np.log(2 * np.pi) / 2 ln_betak -= D * np.log(N_k + 1) / 2 x_term = 0 for i in range(N_k): x_term += x_k[i] @ x_k[i] x_term -= x_k.sum(0) @ x_k.sum(0) / (N_k+1) x_term /= -2 ln_betak += x_term ln_beta += ln_betak ln_betas.append(ln_beta) ln_betas -= logsumexp(ln_betas) return ln_betas最後に、掛け合わせてKN通りの和を取って、を得ます。 コードは以下です。としていくつかの値を計算して同じ値になったので実装に問題はない...と思ってます。
def calc_true_log_posterior(x, mu_star, N, D, K): ln_etazs = [] ln_betas = calc_ln_betas(x, N, D, K) for i in range(pow(K, N)): this_z = [] cluster_str = None if K == 1: cluster_str = ''.zfill(N) else: cluster_str = np.base_repr(i, K).zfill(N) for s in cluster_str: this_z.append(int(s)) this_z = np.array(this_z) ln_etaz = calc_ln_alpha(x, this_z, mu_star, D, K) + ln_betas[i] ln_etazs.append(ln_etaz) return logsumexp(ln_etazs)
Chib's methodのためのギブスサンプリング実装
こちらは普通にギブスサンプリングです。初期値から、を使った尤度と事前確率の値を正規化した分布からをサンプリング、をgivenとして共役事前分布を利用して計算した事後分布からをサンプリングします。まあ定式化はしなくていいでしょう。ギブスサンプリングとそれを使って対数事後確率を近似するコードは以下です。
def GMM_gibbs_sample(x, K, T): N, D = x.shape mu_sample = np.random.normal(0, 1, size=(K, D)) z_samples = [] mu_samples = [] mu_samples.append(mu_sample) for t in tqdm(range(T)): z_sample = np.zeros((N), dtype=int) for i in range(N): ln_posterior = np.zeros(K) for k in range(K): ll = multivariate_normal.logpdf(x[i], mean=mu_sample[k], cov=np.identity(D)) ln_prior = np.log(1 / K) ln_posterior[k] = ll + ln_prior ln_posterior -= logsumexp(ln_posterior) z_sample[i] = np.random.choice(range(K), p=np.exp(ln_posterior)) z_samples.append(z_sample) for k in range(K): x_k = x[z_sample==k] N_k = x_k.shape[0] mu_sample[k] = np.random.multivariate_normal( x_k.sum(0) / (N_k + 1), np.identity(D) / (N_k + 1)) mu_samples.append(mu_sample.copy()) return np.array(z_samples), np.array(mu_samples)def calc_GS_approx_log_posterior(x, mu_star, D, K): T = 1000 z_samples, _ = GMM_gibbs_sample(x, K, T) n_burn_in = 100 ln_alphas = [] for z_sample in z_samples[n_burn_in:]: ln_alphas.append(calc_ln_alpha(x, z_sample, mu_star, D, K)) lpos_approx = logsumexp(ln_alphas) - np.log(T-n_burn_in) return lpos_approx
カーネル密度推定
加えて、共役事前分布を利用しないケースも考えてみます。共役事前分布でないとを解析的に求められないため、MCMCののサンプルからカーネル密度推定で確率密度関数の値を近似することにします。 この場合のMCMCはギブスサンプリングを使ってしまうとズルなので(上記のギブスサンプリングは共役事前分布を利用している)、Stanでサンプリングすることにします。 Stanコード(を格納する文字列変数を定義するpythonコード)は以下です。
stan_code = """ data { int<lower=0> N; int<lower=0> D; int<lower=0> K; vector[D] X[N]; } parameters { vector[D] mu[K]; } transformed parameters{ vector[K] ps[N]; for (i in 1:N){ for (k in 1:K){ ps[i, k] = normal_lpdf(X[i] | mu[k], 1) + log(1.0/K); } } } model { for (k in 1:K){ mu[k] ~ normal(0, 1); } for (i in 1:N){ target += log_sum_exp(ps[i]); } }このStanコードを使って近似したを得るコードは以下です。
def calc_KDE_approx_log_posterior(x, mu_star, N, D, K): stan_data = {"N": N, "D": D, "K": K, "X": x, } random_seed = np.random.randint(0, 10e6) posterior = stan.build(stan_code, data=stan_data, random_seed=random_seed) fit = posterior.sample(num_chains=4, num_samples=10000) # クラスタを区別しないmuの事後サンプル mu_nuts = fit["mu"] mu_post_samples = mu_nuts.transpose((2, 0, 1)).reshape(-1, D) kde = KernelDensity(kernel='gaussian', bandwidth='silverman').fit(mu_post_samples) log_density = kde.score(mu_star.reshape(1, -1)) return log_density
対数周辺尤度の計算
これで、三者の方法による事後分布の確率密度関数の値の計算ができました。あとは対数周辺尤度の計算です。 事前分布と尤度は共通です。以下のコードで計算できます(として各クラスタで同じ値を用いる前提のコードです)。
def calc_log_prior(mu_star, D, K): lpri = 0 for k in range(K): lpri += multivariate_normal.logpdf(mu_star, np.zeros(D), np.identity(D)) return lpri対数周辺尤度は、事後分布を計算する関数を引数として、以下のように書きました。def calc_log_likelihood(x, mu_star, D): lls = multivariate_normal.logpdf(x, mu_star, np.identity(D)) return lls.sum()
def calc_lml(x, mu_star, N, D, K, calc_log_posterior): ll = calc_log_likelihood(x, mu_star, D) lpri = calc_log_prior(mu_star, D, K) lpos = calc_log_posterior(x, mu_star, N, D, K) return ll + lpri - lpos
実験結果
では、いよいよ三者の比較を行います。確認です。 True: 正解周辺尤度を厳密に計算。共役事前分布が必要で、Nは30が限度。 GS: Chib's methodを使った近似計算。共役事前分布が必要。 KDE: KDEを使ったChib's method。どんなモデルでも使える。KDEのためのハイパーパラメタは必要。
データは以下のコードで生成しました。まあ正解パラメタを当てたいわけではないのでサンプルサイズと次元さえ決められればどんな作り方でも良いです。
def generate_data_and_param(N, D, K): mu = np.random.multivariate_normal(np.zeros(D), np.identity(D), size=K) z = np.random.choice(range(K), size=N) x = [] for i in range(N): x.append(np.random.multivariate_normal(mu[z[i]], np.identity(D))) x = np.array(x) return x, z, muサンプルサイズは3, 5, 10、サンプルの次元は1, 3, 5、クラスタ数は1, 2の全ての組み合わせで作成したデータで実験を行いました。 比較実験は以下のコードで行いました。各条件で生成したデータ一つに対して、各近似手法で100回ずつ計算を行って近似のブレも調べました。
T = 100 settings = [(N, D, K) for N in [3, 5, 10] for D in [1, 3, 5] for K in [1, 2]] for N, D, K in settings: output_filename = f'N{N}D{D}K{K}.pickle' print(output_filename) if os.path.exists(output_filename): continue x, z, mu = generate_data_and_param(N, D, K) mu_star = x.mean(0) true_lml = calc_lml(x, mu_star, N, D, K, calc_true_log_posterior) GS_lmls = [] KDE_lmls = [] for t in tqdm(range(T)): GS_lml = calc_lml(x, mu_star, N, D, K, calc_GS_approx_log_posterior) KDE_lml = calc_lml(x, mu_star, N, D, K, calc_KDE_approx_log_posterior) GS_lmls.append(GS_lml) KDE_lmls.append(KDE_lml) result = [true_lml, GS_lmls, KDE_lmls] with open(output_filename, mode='wb') as f: pickle.dump(result, f)そして可視化はこちらのコードで行いました。
fig, axes = plt.subplots(3, 6, figsize=(30, 15)) fontsize = 20 settings = [(N, D, K) for N in [3, 5, 10] for D in [1, 3, 5] for K in [1, 2]] for t, (N, D, K) in enumerate(settings): i = t // 6 j = t % 6 ax = axes[i][j] condition = f'N{N}D{D}K{K}' output_filename = f'{condition}.pickle' with open(output_filename, mode='rb') as f: true_lml, GS_lmls, KDE_lmls = pickle.load(f) approx_lml_df = pd.DataFrame([GS_lmls, KDE_lmls], index=['GS', 'KDE']).T sns.violinplot(x='variable', y='value', data=pd.melt(approx_lml_df), scale='width', ax=ax) plt.sca(ax) plt.hlines(true_lml, -0.5, 1.5, color='blue') plt.title(condition, fontsize=fontsize) plt.xticks(fontsize=fontsize) plt.yticks(fontsize=fontsize) plt.xlabel('Method', fontsize=fontsize) plt.ylabel('Log marginal likelihood', fontsize=fontsize) plt.plot() plt.tight_layout()結果は以下になります。
まず、肝心のGS(Chib's method)ではそこそこ当たっていそうな感じです。あと、サンプルサイズや次元を大きくすると急激に精度が悪化するようなこともなさそうです(K=2までしか調べてないですし、各条件でデータセットも一つしか試していないので雰囲気です)。KDEだとバイアスが大きそうですね。やり方をちゃんと吟味すればうまく行くかもしれません。 周辺尤度の近似精度としてはそんな感じです。
やり残したこと
近似精度の解釈
ただ、周辺尤度を計算したいときというのはモデル選択を行いたいときが多いはずで、他のモデルでの周辺尤度との差も考慮して誤差を評価したいところです。例えば、近似値のほとんどが真の値の±0.1くらいに収まったとしても、比較したいモデルの真の値に0.001くらいの差しかなければ近似手法の有用性は低そうです(ノイズに差が埋もれてしまう)。
渡辺ベイズ理論における周辺尤度の近似
渡辺ベイズ(ベイズ統計の理論と方法、5.1.4 自由エネルギーの近似)では周辺尤度は逆温度を少しずつ変えた複数回のMCMCにより計算できる一方で、ナイーブなモンテカルロ積分ではうまくいかないと書かれています。 WBICも含めて、渡辺ベイズ理論の手法とも比較したいところでした。畜生、根気さえあれば...!
の吟味
元論文では最尤推定などで事後分布の大きいところを占めるの値をとして設定することが推奨されていました。今回は各クラスタで同じ値を使ってしまったので、理論の上では同じことだとしても、性能に差が出る可能性があります。
高次元データでの性能
今回は低次元データのみでの検証でした。高次元データではどうなるのかわかりません...ので、軽くやってみたところ... ん...?あれ...?若干分布が怪しいですね。中央値を見ればそこそこ近そうですが... さらなる調査が望まれる。
まとめ
いかがでしたか? タイトルの内容を試して、理論値(と言っていいのか?正解の値?)と比較しました。 やり残したことはありますが、個人的には満足です。 では。
ガウス過程の導関数を求める
はじめに
ガウス過程回帰による予測値の導関数をGPyのpredictive_gradients
メソッドによって簡単に得られることがわかったので記しておこうと思います.
ガウス過程の導関数
まずガウス過程によって推定された関数の導関数が求められることを確認します.こちらのページを参考にしました.
まず,ガウス過程回帰における予測値(事後期待値)は
で計算できます.ここで,
これにより,ガウス過程回帰の事後期待値の導関数はカーネル値のベクトルの導関数を用いて求められることが分かります.
GPyによる導関数値の取得
まずはデータを作って可視化します.サインカーブにでもしてみましょう.
import GPy import numpy as np import matplotlib.pyplot as plt N=100 np.random.seed(0) x = np.random.normal(0, 1, size=N) y = np.sin(x) + np.random.normal(0, 0.1, size=N) kernel = GPy.kern.RBF(1) model = GPy.models.GPRegression(x.reshape(-1, 1), y.reshape(-1, 1), kernel=kernel) model.optimize() model.plot() plt.show()
GPyではガウス過程回帰の導関数がpredictive_gradients
メソッドで簡単に求められます.
このデータの導関数値を数値解析的に求めます.
x_new = np.array([[0.0]]) epsilon = 1.0e-5 x0_minus_0 = model.predict_quantiles(x_new - epsilon, quantiles=(5, 50, 95))[1][0][0] x0_plus_0 = model.predict_quantiles(x_new + epsilon, quantiles=(5, 50, 95))[1][0][0] print((x0_plus_0 - x0_minus_0) / (2 * epsilon)) # 0.9633798891961475
ここでは小さい値だけ前後にずらした2点の傾きを求めています.
今度はpredictive_gradients
メソッドを使ってみます.
print(model.predictive_gradients(x_new)[0][0][0][0]) # 0.9633798901678498
ちゃんと求められていますね(数値解析的に求めたものを答えとするのってどうなんだろう…).GPyのgithubからも分かるように予測値の分散の導関数(導関数の分散ではない)も一つ目の添字を1
にすることで得られます.
まとめ
pystanでglassoを実装してみる
はじめに
変数間のinteractionを発見するのに使われるglassoの拡張モデルを作る準備としてstanでglassoを実装してみます.ただglassoを使うだけなら
sklearn.covariance.graphical_lasso — scikit-learn 0.23.1 documentation
を使うのが一番手っ取り早いです.
stanファイルの実装
stanでパラメタ推定を行う場合はモデルを記述するstanファイルを用意する必要があります. 詳細は省きますが,glassoに置けるL1正則化項は正規分布の精度行列パラメタの事前分布と見なせます. この事前分布は精度行列の各要素を対角成分と非対角成分に分けて,それぞれ独立同分布とみなすものです. ただし,この事前分布は精度行列としての条件(逆行列=分散共分散行列が半正定値行列)を満たすことが保証されないのでstanでの実装においてはcov_matrix型を定義することで対処します. このcov_matrix型で宣言されたパラメタは分散共分散行列としての条件を満たすことが保証されます. そして実装したのが以下になります.
// glasso.stan data { int<lower=0> n; int<lower=0> p; vector[p] x[n]; real<lower=0> lambda; // real<lower=0> prior_parameter_for_lambda; } transformed data { } parameters { cov_matrix[p] invsigma; vector[p] mu; // real<lower=0> lambda; } transformed parameters { real<lower=0> invlambda; real<lower=0> lambda_half; invlambda = 1/lambda; lambda_half = lambda/2; } model { for (i in 1:n) { x[i] ~ multi_normal_prec(mu, invsigma); } for (j1 in 1:p) { for (j2 in 1:j1) { if (j1==j2) { invsigma[j1,j2] ~ exponential(lambda_half); } else { invsigma[j1,j2] ~ double_exponential(0,invlambda); } } } //Prior for lambda // lambda ~ exponential(prior_parameter_for_lambda); } generated quantities { }
こちらの論文&コードを参考にしました(なんでバイオインフォやねんというのはたまたま思い浮かんだのがこの論文でして...)
人工データによる実験
生成過程から人工データを作成して真のパラメタを当てられるか検証します. まず,importとデータの次元を決めます.
import pandas as pd import numpy as np import sys import seaborn as sns sns.set() import matplotlib.pyplot as plt import pystan N = 1000 p = 5
精度行列を作ります.ただし,精度行列をglassoの生成過程通りに作るとスパースにならないので,categorical分布で0か否かを生成してから値を生成します.
prec_matrix = np.zeros([p,p]) for i in range(p): for j in range(i+1): if i == j: prec_matrix[i, j] = np.random.gamma(100,1) else: prec_matrix[i, j] = np.random.choice(2, p=[0.8, 0.2]) prec_matrix[i, j] = prec_matrix[i, j] * np.random.normal(0,100) prec_matrix[j, i] = prec_matrix[i, j]
固有値の最小値を出力し,負でないことを確認します(負であれば再生成します).
np.linalg.eigvals(np.linalg.inv(prec_matrix)).min()
# 0.005406751950232075
精度行列を可視化します.
sns.heatmap(prec_matrix, cmap='coolwarm', center=0)
データを生成します.平均は0にします.
x = [] for i in range(N): x.append(np.random.multivariate_normal(np.zeros(5), np.linalg.inv(prec_matrix))) x = pd.DataFrame(x)
罰則項の係数は1にしてみます.
stan_data = {'n': N, 'p': p, 'x': x, 'lambda': 1} fit = pystan.stan(file='./glasso.stan', data=stan_data, iter=1000, chains=1) print('Sampling finished.')
一応MAP推定をしたいので,サンプル列からposteriorが最大のindexを取ってきて,該当するサンプルを可視化します.
la = fit.extract(permuted=True) names = fit.model_pars map_index = pd.Series(fit.get_logposterior()[0]).idxmax() la = fit.extract(permuted=True) invsigma_sample = pd.DataFrame(la[names[0]].reshape(-1, la[names[0]][0].size)) # burn_in期間の500を引いたindexを使用 estimated_prec_matrix = pd.DataFrame(np.array(invsigma_sample.iloc[map_index - 500]).reshape(p, p)) sns.heatmap(estimated_prec_matrix, cmap='coolwarm', center=0, annot=True)
うまく推定できていますね. ただ,posteriorが最大になるサンプルを取ってきているだけなので0に見える部分は0ではなかったりします.
おわりに
glassoはposteriorを最大化した場合にスパースになりやすいという手法なのでMCMCで行うのは邪道かもしれませんね.精度行列の非対角成分にラプラス分布の代わりに正規分布をおいてもたいして変わらないと思います.しかし,MCMCで近似された事後分布の形はおおよそ正しい訳ですから,考案したモデルをスクラッチから実装する前の確認としては十分有用なのでは無いかと思います.
特異値分解(singular value decomposition)でベクトルを直交化する
はじめに
特異値分解を使った直交化をする機会があったので書き留めておきます.数式で説明するか文章で説明するかはフィーリングです.
モチベーション
行列が与えられたとします.の列ベクトルを,変換後の空間のスケール(厳密な用語が分かりません...)を変えないように直交したベクトルに変換することを行いたいという話です.なぜかについては説明しません(面倒ですので・・・).
説明
間違ってたらご一報頂けると幸いです.
特異値分解
行列分解の一種.行列をのように分解します.ここで,は直交行列,は対角行列です.
直交行列の性質
行列が直交行列である,とは以下のことと同値です.
今回の直交化においてこの性質が大事になってきます.
特異値分解(SVD)による直交化
SVDによる直交化は を利用して のように直交化します.ここで,性質2より,Vの逆行列をかける代わりに転置行列をかけています.
特異値分解(SVD)による直交化の性質
SVDによる直交化では,被変換行列による変換後の単位円を変えずに各ベクトルを直交にすることができます.以下,説明です.全てのノルムが1のベクトルを考えたとき,座標とみなしたベクトルの集合によって円ができます.これを単位円と言います.によってその各座標を変換すると単位円の形が変わります(楕円もしくは円に).SVDによる直交化では,この楕円の形を変えないベクトルに直交化できるのです.もちろん,各ベクトルについて見てみれば直交化する前と後では値が変わってしまいます.ただ,ノルムが一定のベクトルによってできる,集合レベルでは変わらないようになっているということですね. この話の直感的な図がWikipediaにあります. en.wikipedia.org
この直交化を行列に対して行うと,
となり,SVDによって得られるとで変換後の行列を得られることになります(もちろん元の行列にをかけてもよい). 更にを見てみると直交行列と対角行列の行列積になっているので,性質1より,変換後の行列は列ベクトルが直交していることが分かります(行ベクトルは直交するとは限らない).
要するに,による変換をしてあげると,変換後の空間のスケール(厳密な用語がry)を変えずに列ベクトルを直交化させることができる,というわけですね.
実験
実際にSVDを使って得られるベクトルを眺めてみます.
実験条件
pythonで動かしていきます. まずimportします.numpyと可視化関連です(個人的な好みでseabornを使ってますがimportしなくてもこの後のコードは実行できます,figureの見た目が良くなるだけです).
import numpy as np import matplotlib.pyplot as plt import seaborn as sns sns.set()
適当な行列を作ります.
M = np.array([[1,4],[3,4]])
SVDしてを計算することでMを直交化します(numpyの関数だとは対角成分のみのベクトルで返ってくるためnp.diagで行列にしている).
U, sigma, V = np.linalg.svd(M, full_matrices=True)
M_prime = np.dot(U, np.diag(sigma))
一応直交しているか確認.
np.dot(M_prime.T, M_prime) # array([[ 4.04164878e+01, -1.21399145e-17], # [ -1.21399145e-17, 1.58351216e+00]])
大丈夫そう.最後に可視化します.
fig, ax = plt.subplots(1, 1, figsize=(5, 5)) ax.set_xlim([-6,5]) ax.set_ylim([-6,5]) ax.plot([0,M.transpose()[0][0]], [0,M.transpose()[0][1]], color='b') ax.plot([0,M.transpose()[1][0]], [0,M.transpose()[1][1]], color='b') ax.plot([0,M_prime.transpose()[0][0]], [0,M_prime.transpose()[0][1]], color='r') ax.plot([0,M_prime.transpose()[1][0]], [0,M_prime.transpose()[1][1]], color='r') plt.show()
直交化後のの列ベクトルが描く単位円が変わってなさそうなことが分かります(適当).