2017-12-26

NMFでMovieLensのレコメンドを実装する

前回の投稿( SVDでMovieLensのレコメンドを実装する - け日記 )に引き続き、今回はNMFで映画のレコメンドを実装します。

NMF

非負値行列因子分解(Non-negative Matrix Factorization: NMF)は、ある行列XをW・Hに近似的に分解する方法の一つです。その名の通り、各行列の要素は全て0以上となります。

SVDなどの行列分解では、負の値が現れることがあります。
しかし、実世界で使われる値(画素、頻度、評価値など)の多くは非負値で、分解した結果に解釈を与えるなど非負の性質を保ちたいケースもあります(負の画素値などは存在しませんし)。そうした場合にNMFが使われるようです。

NMFでは、M×Nの行列Xを、N×kの行列Wとk×Mの行列Hに分解します。kはM,Nよりも小さな値が使われます。

$X \approx W \cdot H$

未知のWとHを求めるために、Xとの距離を定義します。
ここでは一番直感的なフロベニウスノルムを使います。

$D= || X-WH ||^2=\sum^{M}_{i=1}\sum^{N}_{j=1}(X_{ij}-(WH)_{ij})^2$

学習では、この距離(つまり誤差)が小さくなるように、HとWを更新します。

$H_{n+1} \leftarrow H_n \frac{W^T_n V_n}{W^T_n W_n H_n}, W_{n+1} \leftarrow W_n \frac{V_n H^T_{n+1}}{W_n H_{n+1} H^T_{n+1}}$

上述以外にもいくつか距離の定義と更新式が提案されていますが、割愛します。こちらの方の投稿がよくまとまっています。

非負値行列因子分解（NMF : Non-negative Matrix Factorization）を改めてやり直してみた - Qiita

NMFを推薦システムに適用する

以下のように、5本の映画に対して、4人のユーザが5段階(1〜5)で評価した、というケースを想定します。
恣意的に選択しているのでそうなっているのですが、これら5本の映画にはロマンス要素とアニメーション要素の2つがあり、ユーザもこの2つのグループに分けられそうです(AとBはロマンス要素が強い映画、CとDはアニメーション要素が強い映画をそれぞれ好んでいる)。

f:id:ohke:20171223190013p:plain

この表は5×4の行列(ここではX)で表現できます。なお、未評価となっている要素は0で埋めています。

$X=\left( \begin{matrix} 5 & 5 & 0 & 3 \\ 4 & 4 & 0 & 0 \\ 0 & 3 & 4 & 0 \\ 0 & 1 & 4 & 5 \\ 1 & 0 & 5 & 4 \end{matrix} \right)$

XをNMFでWとHに分解します。ここではk=2としています。

$X \approx W \cdot H=\left( \begin{matrix} 2.7 & 0.3 \\ 2.1 & 0.0 \\ 0.6 & 0.8 \\ 0.2 & 1.6 \\ 0.1 & 1.6 \end{matrix} \right) \left( \begin{matrix} 1.8 & 1.9 & 0.0 & 0.4 \\ 0.0 & 0.4 & 3.1 & 2.5 \end{matrix} \right) =\left( \begin{matrix} 4.9 & 5.3 & 0.9 & 1.9 \\ 3.8 & 4.0 & 0.0 & 0.9 \\ 1.1 & 1.5 & 2.4 & 2.2 \\ 0.3 & 0.9 & 4.8 & 4.1 \\ 0.2 & 0.8 & 4.9 & 4.1 \end{matrix} \right)$

分解されたWとHを見ると、映画の持つ要素やユーザの嗜好が反映された行列になっていることがわかります。

Wは映画数×kの行列
- 上2行(ローマの休日、タイタニック)は1列目、下2行(サマーウォーズ、AKIRA)は2列目の値がそれぞれ大きくなっており、映画の持つ要素が各列に反映されているように見えます
Hはk×ユーザ数の行列
- 左2列は1行目、右2列は2行目の値がそれぞれ大きくなっており、ユーザごとの嗜好がグループ分けされているように見えます
kによって集約された特徴と捉えることができます

レコメンドの実装

それではNMFを使ったレコメンドを実装します。データセットの準備とランキングの作成については前回の投稿( SVDでMovieLensのレコメンドを実装する - け日記 )とほぼ同じです。

データセットの準備

今回もMovieLens 100K Datasetを使います。

以下からダウンロードして、今回作成するPythonファイルと同じディレクトリに解凍しておきます。以前の投稿( word2vecでitem2vecを実装して映画を推薦する - け日記 )も参考にしてみてください。

MovieLens | GroupLens

評価値をnumpy行列にして持ちます。欠測値は0で補完しています。あわせてアイテムIDとタイトルの一覧もDataFrameにロードしておきます。

import numpy as np
import pandas as pd
import codecs


# 評価値を持つファイルはu.data。楽に行列化するために、まずはDataFrameでロードする
ratings = pd.read_csv('ml-100k/u.data', delimiter='\t', header=None).iloc[:, :3]
ratings.rename(columns={0: 'user_id', 1: 'item_id', 2: 'rating'}, inplace=True)
print('ratings.shape: {}'.format(ratings.shape)) # ratings.shape: (100000, 3)

# item_idを行、user_idを列とする1682×943の行列(欠測値は0埋め)
rating_matrix = ratings.pivot(index='item_id', columns='user_id', values='rating').fillna(0).as_matrix()
print('rating_matrix.shape: {}'.format(rating_matrix.shape)) # rating_matrix.shape: (1682, 943)

# pd.read_csvで読み込もうとすると文字コードエラーになるため、
# codecs.openでエラーを無視してファイルを読み込んだ後にDataFrameにする
with codecs.open('ml-100k/u.item', 'r', 'utf-8', errors='ignore') as f:
    item_df = pd.read_table(f, delimiter='|', header=None).ix[:, 0:1]
    item_df.rename(columns={0: 'item_id', 1: 'item_title'}, inplace=True)
    item_df.set_index('item_id', inplace=True)

items = pd.Series(item_df['item_title'])
print('items.shape: {}'.format(items.shape)) # items.shape: (1682, 1)

評価ランキングの予測と表示

指定したユーザに対して、未評価のアイテムの評価値を予測し、ランキングを返す関数を作ります。

未評価の映画については0で補完しており、また映画やユーザごとの評価値の平均なども考慮していないため、絶対値の予測精度は高くないと予想されます。そのため、ランキングにしています。

def predict_ranking(user_index, reducted_matrix, original_matrix, n):
    # 対象ユーザのSVD評価値
    reducted_vector = reducted_matrix[:, user_index]
    
    # 評価済みのアイテムの値は0にする
    filter_vector = original_matrix[:, user_index] == 0
    predicted_vector = reducted_vector * filter_vector

    # 上位n個のアイテムのインデックスと予測した評価値を返す
    return [(i, predicted_vector[i]) for i in np.argsort(predicted_vector)[::-1][:n]]


def print_ranking(user_ids, items, reducted_matrix):
    for user_id in user_ids:
        predicted_ranking = predict_ranking(user_id - 1, reducted_matrix, rating_matrix, 10)
        print('User: {}:'.format(user_id))
        for item_id, rating in predicted_ranking:
            # アイテムID, 映画タイトル, 予測した評価値を表示
            print('{}: {} [{}]'.format(item_id, items[item_id + 1], rating))

NMFによる行列分解

NMFによる行列分解ですが、sklearn.decomposition.NMFでお手軽にできますので、これを使います。

n_componentsにk=30を指定しています(前回のSVDと同じ)
initで初期値の設定方法を指定できます(今回はランダム)
MAE(絶対平均誤差)は約0.26
- XとWHの各要素は平均0.26の差がある、ということです

from sklearn.decomposition import NMF


def reduct_matrix_nmf(matrix, k):
    # NMFで分解
    model = NMF(n_components=k, init='random', random_state=1)
    w = model.fit_transform(matrix)
    h = model.components_
    
    return np.dot(w, h)


reducted_matrix_nmf = reduct_matrix_nmf(rating_matrix, k=30)

print('MAE: {}'.format(np.average(np.absolute(rating_matrix - reducted_matrix_nmf))))
# MAE: 0.25617962729902044

評価

それでは実際に推薦をさせてみます。

前回のSVDの結果と比較するために、前回と同じ2人のユーザ(ユーザID 267と868)を対象にランキングを生成しました。
概ねSVDの時と同じようなタイトルが得られましたが、負値を含まない分、予測評価値は全体的にNMFの方が高いことがわかります。

この2人は"Akira(1988)"と"Ghost in the Shell(Kokaku kidotai)(1995)"を5で評価しています(つまりSF+アニメーション好きです)

User: 267:
95: Terminator 2: Judgment Day (1991) [4.355731842820824]
10: Seven (Se7en) (1995) [3.5465670297733514]
231: Young Guns (1988) [3.3581447476132875]
78: Fugitive, The (1993) [3.3281046778126173]
153: Monty Python's Life of Brian (1979) [3.3028868626155874]
41: Clerks (1994) [3.1889273065852977]
172: Princess Bride, The (1987) [3.158916000181793]
356: One Flew Over the Cuckoo's Nest (1975) [3.1348566836324423]
0: Toy Story (1995) [3.128039827666597]
3: Get Shorty (1995) [2.8694586381163]

User: 868:
174: Brazil (1985) [4.270610564574103]
257: Contact (1997) [2.965456387373195]
356: One Flew Over the Cuckoo's Nest (1975) [2.914023120253882]
143: Die Hard (1988) [2.5656501189688483]
430: Highlander (1986) [2.5194608160565846]
195: Dead Poets Society (1989) [2.319622366933584]
21: Braveheart (1995) [2.2844329680205355]
41: Clerks (1994) [2.212755642335326]
379: Star Trek: Generations (1994) [2.212123031118453]
27: Apollo 13 (1995) [2.1517770562746565]

まとめ

今回はNMFを使った映画のレコメンドを実装しました。

SVDが含んでいた負値を取り除くことができましたが、欠測値や、ユーザごとのバイアス(甘め・辛め)が考慮されたモデルには、まだなっていません。

次は、欠測値やバイアスを考慮した推薦システムに取り組みたいと思います。

参考にした文献・記事

岩波データサイエンス Vol.5

発売日: 2017/02/16
メディア: 単行本（ソフトカバー）

qiita.com

http://lab.adn-mobasia.net/?p=1220lab.adn-mobasia.net

非負値行列因子分解(NMF)でブログ記事のレコメンドをしてみる | かものはしの分析ブログ

2017-12-20

SVDでMovieLensのレコメンドを実装する

Python Recommender system

前回の投稿( Pythonで特異値分解(SVD)を理解する - け日記 )で特異値分解(SVD)についてPythonで(車輪の再発明的に)実装してみました。
今回はSVDを使って、映画のレコメンドシステムを作ります。データセットはMovieLens 100Kを用います。

データセット

今回データセットに用いるMovieLens 100K Datasetを準備します。

以下からダウンロードして、今回作成するPythonファイルと同じディレクトリに解凍します。以前の投稿( word2vecでitem2vecを実装して映画を推薦する - け日記 )も参考にしてみてください。

grouplens.org

評価値をnumpy行列にして持ちます。欠測値は0で補完しています

あわせてアイテムIDとタイトルの一覧もDataFrameにロードしておきます。

import numpy as np
import pandas as pd
import codecs


# 評価値を持つファイルはu.data。楽に行列化するために、まずはDataFrameでロードする
ratings = pd.read_csv('ml-100k/u.data', delimiter='\t', header=None).iloc[:, :3]
ratings.rename(columns={0: 'user_id', 1: 'item_id', 2: 'rating'}, inplace=True)
print('ratings.shape: {}'.format(ratings.shape)) # ratings.shape: (100000, 3)

# item_idを行、user_idを列とする1682*943の行列(欠測値は0埋め)
rating_matrix = ratings.pivot(index='item_id', columns='user_id', values='rating').fillna(0).as_matrix()
print('rating_matrix.shape: {}'.format(rating_matrix.shape)) # rating_matrix.shape: (1682, 943)

# pd.read_csvで読み込もうとすると文字コードエラーになるため、
# codecs.openでエラーを無視してファイルを読み込んだ後にDataFrameにする
with codecs.open('ml-100k/u.item', 'r', 'utf-8', errors='ignore') as f:
    item_df = pd.read_table(f, delimiter='|', header=None).ix[:, 0:1]
    item_df.rename(columns={0: 'item_id', 1: 'item_title'}, inplace=True)
    item_df.set_index('item_id', inplace=True)

items = pd.Series(item_df['item_title'])
print('items.shape: {}'.format(items.shape)) # items.shape: (1682, 1)

SVDで行列分解

まずは行列分解するロジックを実装します。

numpy.linarg.svdで分解して、指定のk次元まで近似した行列を返します。
k=2(全ての特異値を使う)では同じ行列、k=1では元の行列と少し値が異なる行列がそれぞれ得られていることがわかります。

今回は評価値行列(rating_matrix)をk=30で次元圧縮します。この値は決め打ちです。

def reduct_matrix(matrix, k):
    # SVDで分解
    # full_metrices=Falseとすると、uがM×M(正方行列)ではなく、M×Nとなる
    u, s, v = np.linalg.svd(matrix, full_matrices=False)
    
    # k次元まで圧縮
    s_k = np.array(s)
    s_k[k:] = 0
    sigma_k = np.diag(s_k)
    
    return np.dot(u, np.dot(sigma_k, v))


reduct_matrix(np.matrix([[1, 2],[3, 4],[5, 6]]), k=2)
# matrix([[ 1.,  2.],
#         [ 3.,  4.],
#         [ 5.,  6.]])

reduct_matrix(np.matrix([[1, 2],[3, 4],[5, 6]]), k=1)
# matrix([[ 1.357,  1.718],
#         [ 3.097,  3.923],
#         [ 4.838,  6.128]])

# k=30で次元圧縮
reducted_matrix = reduct_matrix(rating_matrix, k=30)

評価ランキングの予測

対象ユーザが未評価のアイテムに対して、評価値を予測してランキングを作ります。

SVDでは負の値をとることもあります。また、今回の単純なモデルでは、アイテムやユーザごとの評価値の平均を考慮しておらず、評価していないアイテムについては0で補完しています。
そのため、絶対値の予測(例えば「ユーザAの映画Xの評価は3.7」という予測)の精度は期待できず、ランキングにしています。

def predict_ranking(user_index, reducted_matrix, original_matrix, n):
    # 対象ユーザのSVD評価値
    reducted_vector = reducted_matrix[:, user_index]
    
    # 評価済みのアイテムの値は0にする
    filter_vector = original_matrix[:, user_index] == 0
    predicted_vector = reducted_vector * filter_vector

    # 上位n個のアイテムのインデックスを返す
    return [(i, predicted_vector[i]) for i in np.argsort(predicted_vector)[::-1][:n]]


def print_ranking(user_ids, items):
    for user_id in user_ids:
        predicted_ranking = predict_ranking(user_id - 1, reducted_matrix, rating_matrix, 10)
        print('User: {}:'.format(user_id))
        for item_id, rating in predicted_ranking:
            # アイテムID, 映画タイトル, 予測した評価値を表示
            print('{}: {} [{}]'.format(item_id, items[item_id + 1], rating))

ここまで実装して、具体的に2人のユーザの推薦ランキング(TOP10)を作ります。
この2人(ユーザIDで267と868)は、"Akira(1988)"と"Ghost in the Shell(Kokaku kidotai)(1995)"をいずれも5で評価してます(つまり、SF+アニメーション好きです)。
結果は以下の通りです。

267の人は、1位が"Terminator 2: Judgment Day (1991)"でSFだったり、6位が"Toy Story (1995)"でアニメだったりです
- 全体的にはサスペンス物が多い印象です
868の人は、1位の"Brazil (1985)"や9位の"Apollo 13 (1995)"などのSFが含まれます
- 全体的には、ドラマのような、カルトのような、、、という感じです

User: 267
95: Terminator 2: Judgment Day (1991) [4.682164278048727]
133: Citizen Kane (1941) [3.3170890216340365]
10: Seven (Se7en) (1995) [3.2720405176668845]
78: Fugitive, The (1993) [3.1124326890407046]
3: Get Shorty (1995) [2.96920110165377]
172: Princess Bride, The (1987) [2.754536418480666]
0: Toy Story (1995) [2.6663561243251217]
231: Young Guns (1988) [2.6232331423202835]
635: Escape from New York (1981) [2.618586732939522]
153: Monty Python's Life of Brian (1979) [2.4655795942476715]

User: 868
174: Brazil (1985) [3.9280576476220705]
195: Dead Poets Society (1989) [2.7738253899665537]
430: Highlander (1986) [2.764929365245002]
356: One Flew Over the Cuckoo's Nest (1975) [2.573415252362662]
143: Die Hard (1988) [2.4954484877473875]
41: Clerks (1994) [2.29187936816679]
181: GoodFellas (1990) [2.132322532767669]
287: Scream (1996) [2.089763076093557]
27: Apollo 13 (1995) [2.0590921066379777]
82: Much Ado About Nothing (1993) [1.920902451839764]

まとめ

SVDを使ってMovieLens 100Kのデータでレコメンドを行いました。

素のSVDを推薦システムへ適用するといくつか問題があります。

負の値を含む
- MovieLensでも1〜5で評価されるため、予測値がマイナスとなるのは適合していない
欠測値も含めて計算される
- 今回は0で補完しているため、「値0で評価している」という誤った状態で計算されている
各ユーザのバイアス(甘め辛め)が考慮されていない

1番目の問題を解決するためには、Matrix Factorizationという別の行列分解手法が必要となります。
また、2番目と3番目の問題には、モデルを拡張して対応する必要があります。以下の記事などが参考になりそうです。

ameblo.jp

参考文献

岩波データサイエンス Vol.5

作者: 岩波データサイエンス刊行委員会
出版社/メーカー: 岩波書店
発売日: 2017/02/16
メディア: 単行本（ソフトカバー）
この商品を含むブログ (3件) を見る

2017-12-14

Pythonで特異値分解(SVD)を理解する

数学 Python

以前の投稿( 論文メモ: Item2Vec: Neural Item Embedding for Collaborative Filtering - け日記 )で比較対象になっていた特異値分解(SVD)についてまとめ、Pythonで実装してみます。

SVDとは

特異値分解(singular value decomposition: SVD)は、数学的にはM×N(M行N列)の行列を分解する方法の一つです。
コンピュータサイエンスでは、行列で表現される特徴(情報検索の分野では、例えば文書毎の単語の出現頻度の行列など)を可能な限り損ねること無く次元を圧縮するために利用され、多次元の特徴を扱う画像処理、自然言語処理、検索や推薦に応用されています。
他の次元圧縮手法としては、主成分分析(PCA)や非負値行列因子分解(NMF)などがあります。

SVDで何ができるのか

任意のM×Nの行列Cは、以下の形に分解することができます。

$C=U \Sigma V^T$

左特異行列Uは、M×Mの正方行列で、各列は ${CC^T}$ と直交する固有ベクトル
右特異行列Vは、N×Nの正方行列で、各列は ${C^TC}$ と直交する固有ベクトル
ΣはN×M行列で、対角行列成分が (ただし、 )で、それ以外が0
- ${\sigma_i}$ は、行列Cの特異値とよばれる
- rは行列Cのランク(一次独立な列または行ベクトルの数)

このとき、特異値σは ${CC^T}$ の固有値の平方根となります。

${CC^T=U \Sigma V^T V \Sigma^T U^T=U \Sigma\Sigma^T U^T}$ となる
${CC^T}$ は正方対称行列で、対称対角化定理(正方対称行列Sは ${Q \Lambda Q^T}$ に分解でき、このときのΛはSの固有値になる)により、 ${\Sigma\Sigma^T}$ の対角要素は ${CC^T}$ の固有値 ${\lambda_1 \dots \lambda_r}$ にあたる
Σは対角行列のため ${\Sigma\Sigma^T=\Sigma^2}$ であるため、 ${\sigma_i=\sqrt{\lambda_i}}$ となる

具体的な例を挙げます。

以下の3つの文章から単語(ここでは名詞と動詞に限定します)を抽出して、出現頻度行列Cを作ります。

I drive my car.
I read this document.
I read this paper.

行列Cは、行は単語(上から"I", "drive", "car", "read", "document", "paper")、列は文書で、値は単語の出現頻度にそれぞれ対応しています。

$$ C=\left( \begin{matrix} 1 & 1 & 1 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 1 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix} \right) $$

これを特異値分解( ${C=U \Sigma V^T}$ )するとこうなります。

UのN列目以降、ΣがN行目以降は削除してます(ΣのN行目以降は全て0のため)

$$ C=\left( \begin{matrix} 0.717 & -0.158 & 0 \\ 0.192 & -0.59 & 0 \\ 0.192 & -0.59 & 0 \\ 0.525 & 0.432 & 0 \\ 0.262 & 0.216 & -0.707 \\ 0.262 & 0.216 & 0.707 \end{matrix} \right) \left( \begin{matrix} 2.394 & 0 & 0 \\ 0 & 1.506 & 0 \\ 0 & 0 & 1 \end{matrix} \right) \left( \begin{matrix} 0.46 & 0.628 & 0.628 \\ -0.888 & 0.325 & 0.325 \\ 0 & -0.707 & 0.707 \end{matrix} \right) $$

実践では、データ量が減らせることもあって、 ${U\Sigma}$ と ${V^T}$ の2つの行列で持つことが多いとのことです。

$$ C=\left( \begin{matrix} 1.716 & -0.238 & 0 \\ 0.46 & -0.888 & 0 \\ 0.46 & -0.888 & 0 \\ 1.256 & 0.65 & 0 \\ 0.628 & 0.325 & -0.707 \\ 0.628 & 0.325 & 0.707 \end{matrix} \right) \left( \begin{matrix} 0.46 & 0.628 & 0.628 \\ -0.888 & 0.325 & 0.325 \\ 0 & -0.707 & 0.707 \end{matrix} \right) $$

numpyでSVDを実装するとこんな感じになりますが、numpy.linarg.svdがありますので、実践ではそちらを使いましょう。

import numpy as np
import pandas as pd


np.set_printoptions(precision=3, suppress=True)

# 分解する行列
c = np.array([[1,1,1],[1,0,0,],[1,0,0],[0,1,1],[0,1,0],[0,0,1]])

# 特異値行列Σ(sigma) → 右特異行列V(v) → 左特異行列U(u) の順に求める

# C^TCの固有値と固有ベクトルの計算
ctc = np.dot(c.T, c)
eigen_values, eigen_vectors = np.linalg.eig(ctc)

# 特異値の計算
singular_values = np.sqrt(eigen_values)
singular_index = np.argsort(singular_values)[::-1]

# 特異値行列の計算
sigma = np.diag(singular_values[singular_index])

# 右特異行列の計算
v = eigen_vectors[:,singular_index]

# 左特異行列の計算
u = np.array([np.dot(c, v[:,i]) / sigma.diagonal()[i] for i in range(len(sigma.diagonal()))]).T

SVDを使うとどうして次元圧縮できるのか

${C=U \Sigma V^T}$ においてΣは対角行列のため、対角要素σが大きければ、Cに与える影響も大きくなります(逆に言いますと、σが小さくなれば影響も小さくなります)。
そこで、上位k(<M)個のσを残し、それ以外は0とした ${C_k=U_k \Sigma_k V_k^T}$ で近似することで、次元を削減します。低ランク近似と呼ばれており、小さな特異値から削減することで、 ${C-C_k}$ の各要素の2乗誤差が最小になることが知られています。

上の例でk=2として ${C_k=U_k \Sigma_k V_k^T}$ を求めると、以下の形になり、概ね行列Cと ${C_k}$ の値が近いことがわかります。
またこの近似によって、 ${C_k}$ の5行目("document")と6行目("paper")の値はいずれも(0, 0.5, 0.5)となっています。同一の文書に現れるわけではない"document"と"paper"が、SVDによって類似した単語と判定できるようになった、と言うことができます。
SVDを使ってレコメンドなど欠測値が多いデータの補完に用いられているケースもあります。

$$ C_k=\left( \begin{matrix} 0.717 & -0.158 \\ 0.192 & -0.59 \\ 0.192 & -0.59 \\ 0.525 & 0.432 \\ 0.262 & 0.216 \\ 0.262 & 0.216 \end{matrix} \right) \left( \begin{matrix} 2.394 & 0 \\ 0 & 1.506 \end{matrix} \right) \left( \begin{matrix} 0.46 & 0.628 & 0.628 \\ -0.888 & 0.325 & 0.325 \end{matrix} \right)=\left( \begin{matrix} 1 & 1 & 1 \\ 1 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 1 \\ 0 & 0.5 & 0.5 \\ 0 & 0.5 & 0.5 \end{matrix} \right) $$

Pythonの実装は以下のとおりです。

# U_kの計算(k列目以降を削除)
u_k = np.delete(u, range(k, u.shape[1]), axis=1)

# Σ_kの計算(k×k行列にする)
sigma_k = np.delete(sigma, range(k, sigma.shape[0]), axis=0)
sigma_k = np.delete(sigma_k, range(k, sigma.shape[1]), axis=1)

# V_kの計算(k列目以降を削除
v_k = np.delete(v, range(k, v.shape[1]), axis=1)

c_k = u_k.dot(sigma_k.dot(v_k.T))

参考文献

理解・実装にあたり、こちら2冊が参考になりました。

岩波データサイエンス Vol.5

作者: 岩波データサイエンス刊行委員会
出版社/メーカー: 岩波書店
発売日: 2017/02/16
メディア: 単行本（ソフトカバー）
この商品を含むブログ (3件) を見る

情報検索の基礎

作者: Christopher D.Manning,Prabhakar Raghavan,Hinrich Schutze,岩野和生,黒川利明,濱田誠司,村上明子
出版社/メーカー: 共立出版
発売日: 2012/06/23
メディア: 単行本
購入: 2人クリック: 69回
この商品を含むブログ (5件) を見る