(57)【要約】【課題】精度よく、データを解析することができる。【解決手段】左因子行列Ｗ更新部３２は、左因子行列Ｗのｋ番目の列ベクトルｗｋの各々について、損失関数の値を最小化するように、近接点行列Ｐの重み付き平均ベクトルｑを用いて、列ベクトルｗｋを更新する。右因子行列Ｈ更新部３４は、右因子行列Ｈのｋ番目の列ベクトルｈｋの各々について、損失関数の値を最小化するように、近接点行列Ｐの重み付き平均ベクトルｑを用いて、列ベクトルｈｋを更新する。【選択図】図１
2017111754.tif 000002
【特許請求の範囲】
【請求項１】
Ｉ×Ｊの非負の要素からなる行列Ｘ、及び前記行列Ｘにおける各要素が欠損値であるか又は観測値であるかを表すＩ×Ｊの欠損値行列Ｍを受け付けるデータ入力部と、
前記データ入力部により受け付けた前記行列Ｘ、及び前記欠損値行列Ｍに基づいて、非負の要素からなるＩ×Ｋの左因子行列Ｗと、非負の要素からなるＪ×Ｋの右因子行列Ｈとの積の要素であって、かつ、前記観測値である要素と、前記行列Ｘの要素との距離を表す一般化カルバックライブラーダイバージェンス距離の和を用いて表される、非負値行列分解の損失関数の値を最小化するように、前記左因子行列Ｗ、及び前記右因子行列Ｈのパラメータを更新することを繰り返すパラメータ更新部と、
前記パラメータ更新部により更新された前記左因子行列Ｗ、及び前記右因子行列Ｈを出力するパラメータ出力部と、を含み、
前記パラメータ更新部は、前記行列Ｘ、前記欠損値行列Ｍ、前記左因子行列Ｗ、及び前記右因子行列Ｈに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記左因子行列Ｗ、及び前記右因子行列Ｈのパラメータを更新する解析装置。
【請求項２】
前記パラメータ更新部は、
前記左因子行列Ｗのｋ番目の列ベクトルｗｋの各々について、前記データ入力部により受け付けた前記行列Ｘ及び前記欠損値行列Ｍと、前記左因子行列Ｗと、前記右因子行列Ｈとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記列ベクトルｗｋを更新する左因子行列Ｗ更新部と、
前記右因子行列Ｈのｋ番目の列ベクトルｈｋの各々について、前記データ入力部により受け付けた前記行列Ｘ及び前記欠損値行列Ｍと、前記左因子行列Ｗと、前記右因子行列Ｈとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記列ベクトルｈｋを更新する右因子行列Ｈ更新部と、
予め定められた条件を満たすまで、前記左因子行列Ｗ更新部による前記列ベクトルｗｋの各々の更新、及び前記右因子行列Ｈ更新部による前記列ベクトルｈｋの各々の更新を繰り返させる繰り返し判定部と、
を含む請求項１に記載の解析装置。
【請求項３】
前記左因子行列Ｗ更新部は、前記列ベクトルｗｋを更新する際に、
前記近接点行列Ｐと前記重み付き平均ベクトルｑとに基づいて、並列近接法のための補助変数Ｙを演算し、
前記データ入力部により受け付けた前記行列Ｘ及び前記欠損値行列Ｍと、前記左因子行列Ｗと、前記右因子行列Ｈと、前記補助変数Ｙとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記列ベクトルｗｋを更新することを、所定回繰り返し、
前記右因子行列Ｈ更新部は、前記列ベクトルｈｋを更新する際に、
前記近接点行列Ｐと前記重み付き平均ベクトルｑとに基づいて、並列近接法のための補助変数Ｙを演算し、
前記データ入力部により受け付けた前記行列Ｘ及び前記欠損値行列Ｍと、前記左因子行列Ｗと、前記右因子行列Ｈと、前記補助変数Ｙとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記列ベクトルｈｋを更新することを、所定回繰り返す請求項２に記載の解析装置。
【請求項４】
データ入力部が、Ｉ×Ｊの非負の要素からなる行列Ｘ、及び前記行列Ｘにおける各要素が欠損値であるか又は観測値であるかを表すＩ×Ｊの欠損値行列Ｍを受け付けるステップと、
パラメータ更新部が、前記データ入力部により受け付けた前記行列Ｘ、及び前記欠損値行列Ｍに基づいて、非負の要素からなるＩ×Ｋの左因子行列Ｗと、非負の要素からなるＪ×Ｋの右因子行列Ｈとの積の要素であって、かつ、前記観測値である要素と、前記行列Ｘの要素との距離を表す一般化カルバックライブラーダイバージェンス距離の和を用いて表される、非負値行列分解の損失関数の値を最小化するように、前記左因子行列Ｗ、及び前記右因子行列Ｈのパラメータを更新することを繰り返すステップと、
パラメータ出力部が、前記パラメータ更新部により更新された前記左因子行列Ｗ、及び前記右因子行列Ｈを出力するステップと、を含み、
前記パラメータ更新部が更新するステップは、前記行列Ｘ、前記欠損値行列Ｍ、前記左因子行列Ｗ、及び前記右因子行列Ｈに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記左因子行列Ｗ、及び前記右因子行列Ｈのパラメータを更新する解析方法。
【請求項５】
前記パラメータ更新部が更新するステップは、
左因子行列Ｗ更新部が、前記左因子行列Ｗのｋ番目の列ベクトルｗｋの各々について、前記データ入力部により受け付けた前記行列Ｘ及び前記欠損値行列Ｍと、前記左因子行列Ｗと、前記右因子行列Ｈとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記列ベクトルｗｋを更新するステップと、
右因子行列Ｈ更新部が、前記右因子行列Ｈのｋ番目の列ベクトルｈｋの各々について、前記データ入力部により受け付けた前記行列Ｘ及び前記欠損値行列Ｍと、前記左因子行列Ｗと、前記右因子行列Ｈとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記列ベクトルｈｋを更新するステップと、
繰り返し判定部が、予め定められた条件を満たすまで、前記左因子行列Ｗ更新部による前記列ベクトルｗｋの各々の更新、及び前記右因子行列Ｈ更新部による前記列ベクトルｈｋの各々の更新を繰り返させるステップと、
を含む請求項４に記載の解析方法。
【請求項６】
前記左因子行列Ｗ更新部が更新するステップでは、前記列ベクトルｗｋを更新する際に、
前記近接点行列Ｐと前記重み付き平均ベクトルｑとに基づいて、並列近接法のための補助変数Ｙを演算し、
前記データ入力部により受け付けた前記行列Ｘ及び前記欠損値行列Ｍと、前記左因子行列Ｗと、前記右因子行列Ｈと、前記補助変数Ｙとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記列ベクトルｗｋを更新することを、所定回繰り返し、
前記右因子行列Ｈ更新部が更新するステップでは、前記列ベクトルｈｋを更新する際に、
前記近接点行列Ｐと前記重み付き平均ベクトルｑとに基づいて、並列近接法のための補助変数Ｙを演算し、
前記データ入力部により受け付けた前記行列Ｘ及び前記欠損値行列Ｍと、前記左因子行列Ｗと、前記右因子行列Ｈと、前記補助変数Ｙとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記列ベクトルｈｋを更新することを、所定回繰り返す請求項５に記載の解析方法。
【請求項７】
コンピュータを、請求項１〜請求項３のいずれか１項に記載の解析装置の各部として機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、解析装置、方法、及びプログラムに係り、特に、非負値行列分解を行うための解析装置、方法、及びプログラムに関する。
【背景技術】
【０００２】
非特許文献１に記載されている非負値行列因子分解法（ＮＭＦ）は、画像データ、音声データ、購買履歴データなどから属性の共起パターンを抽出する課題で頻繁に用いられている。特に非特許文献２に記載の一般化カルバックライブラー（ＫＬ）ダイバージェンスという距離規準をＮＭＦの損失関数に用いると、上記の課題では抽出された共起パターンの解釈性が高いと報告されており、この損失関数を最小化するための多数のアルゴリズムが提案されている（非特許文献３参照）。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】D. D. Lee and H. S. Seung. Algorithms for non-negative matrix factorization. In NIPS, 2000.
【非特許文献２】A. Cichocki, R. Zdunek, A. H. Phan, and S. Amari. Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation. Wiley, 2009.
【非特許文献３】N. Gillis. The why and how of nonnegative matrix factorization. arXiv preprint arXiv:1401.5226, 2014.
【非特許文献４】S. Sra, D. Kim, and B. Scholkopf. Non-monotonic poisson likelihood maximization.Technicalreport, 170, Max Planck Institute for Biological Cybernetics, 2008.
【非特許文献５】S. Sra and I. S. Dhillon. Generalized nonnegative matrix approximations with Bregman diver gences. In NIPS, 2005.
【非特許文献６】C.-J. Hsieh and I. S. Dhillon. Fast coordinate descent methods with variable selection for non-negative matrix factorization. In SIGKDD, 2011.
【非特許文献７】Dennis L. Sun and C´edric F´evotte. Alternating direction method of multipliers for non-negative matrix factorization with the beta-divergence. In Proc. ICASSP, pages 6201-6205. IEEE, 2014.
【発明の概要】
【発明が解決しようとする課題】
【０００４】
しかし、一般化ＫＬダイバージェンスは非リプシッツ連続という性質を持つため、既存のアルゴリズムでは計算が収束する保証がされていない（非特許文献４参照）。そのため、解析対象のデータに欠損値が多く含まれている場合に、アルゴリズムの挙動が悪化し、計算結果として得られる共起パターンの性能が悪化する場合がある。
【０００５】
本発明は、上記問題点を解決するために成されたものであり、精度よく、データを解析することができる解析装置、方法、及びプログラムを提供することを目的とする。
【課題を解決するための手段】
【０００６】
上記目的を達成するために、第１の発明に係る解析装置は、Ｉ×Ｊの非負の要素からなる行列Ｘ、及び前記行列Ｘにおける各要素が欠損値であるか又は観測値であるかを表すＩ×Ｊの欠損値行列Ｍを受け付けるデータ入力部と、前記データ入力部により受け付けた前記行列Ｘ、及び前記欠損値行列Ｍに基づいて、非負の要素からなるＩ×Ｋの左因子行列Ｗと、非負の要素からなるＪ×Ｋの右因子行列Ｈとの積の要素であって、かつ、前記観測値である要素と、前記行列Ｘの要素との距離を表す一般化カルバックライブラーダイバージェンス距離の和を用いて表される、非負値行列分解の損失関数の値を最小化するように、前記左因子行列Ｗ、及び前記右因子行列Ｈのパラメータを更新することを繰り返すパラメータ更新部と、前記パラメータ更新部により更新された前記左因子行列Ｗ、及び前記右因子行列Ｈを出力するパラメータ出力部と、を含んで構成され、前記パラメータ更新部は、前記行列Ｘ、前記欠損値行列Ｍ、前記左因子行列Ｗ、及び前記右因子行列Ｈに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記左因子行列Ｗ、及び前記右因子行列Ｈのパラメータを更新する。
【０００７】
また、第１の発明に係る解析装置において、前記パラメータ更新部は、前記左因子行列Ｗのｋ番目の列ベクトルｗｋの各々について、前記データ入力部により受け付けた前記行列Ｘ及び前記欠損値行列Ｍと、前記左因子行列Ｗと、前記右因子行列Ｈとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記列ベクトルｗｋを更新する左因子行列Ｗ更新部と、前記右因子行列Ｈのｋ番目の列ベクトルｈｋの各々について、前記データ入力部により受け付けた前記行列Ｘ及び前記欠損値行列Ｍと、前記左因子行列Ｗと、前記右因子行列Ｈとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記列ベクトルｈｋを更新する右因子行列Ｈ更新部と、予め定められた条件を満たすまで、前記左因子行列Ｗ更新部による前記列ベクトルｗｋの各々の更新、及び前記右因子行列Ｈ更新部による前記列ベクトルｈｋの各々の更新を繰り返させる繰り返し判定部と、を含むようにしてもよい。
【０００８】
また、第１の発明に係る解析装置において、前記左因子行列Ｗ更新部は、前記列ベクトルｗｋを更新する際に、前記近接点行列Ｐと前記重み付き平均ベクトルｑとに基づいて、並列近接法のための補助変数Ｙを演算し、前記データ入力部により受け付けた前記行列Ｘ及び前記欠損値行列Ｍと、前記左因子行列Ｗと、前記右因子行列Ｈと、前記補助変数Ｙとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記列ベクトルｗｋを更新することを、所定回繰り返し、前記右因子行列Ｈ更新部は、前記列ベクトルｈｋを更新する際に、前記近接点行列Ｐと前記重み付き平均ベクトルｑとに基づいて、並列近接法のための補助変数Ｙを演算し、前記データ入力部により受け付けた前記行列Ｘ及び前記欠損値行列Ｍと、前記左因子行列Ｗと、前記右因子行列Ｈと、前記補助変数Ｙとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記列ベクトルｈｋを更新することを、所定回繰り返すようにしてもよい。
【０００９】
第２の発明に係る解析方法は、データ入力部が、Ｉ×Ｊの非負の要素からなる行列Ｘ、及び前記行列Ｘにおける各要素が欠損値であるか又は観測値であるかを表すＩ×Ｊの欠損値行列Ｍを受け付けるステップと、パラメータ更新部が、前記データ入力部により受け付けた前記行列Ｘ、及び前記欠損値行列Ｍに基づいて、非負の要素からなるＩ×Ｋの左因子行列Ｗと、非負の要素からなるＪ×Ｋの右因子行列Ｈとの積の要素であって、かつ、前記観測値である要素と、前記行列Ｘの要素との距離を表す一般化カルバックライブラーダイバージェンス距離の和を用いて表される、非負値行列分解の損失関数の値を最小化するように、前記左因子行列Ｗ、及び前記右因子行列Ｈのパラメータを更新することを繰り返すステップと、パラメータ出力部が、前記パラメータ更新部により更新された前記左因子行列Ｗ、及び前記右因子行列Ｈを出力するステップと、を含んで実行することを特徴とし、前記パラメータ更新部が更新するステップは、前記行列Ｘ、前記欠損値行列Ｍ、前記左因子行列Ｗ、及び前記右因子行列Ｈに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記左因子行列Ｗ、及び前記右因子行列Ｈのパラメータを更新する。
【００１０】
また、第２の発明に係る解析方法において、前記パラメータ更新部が更新するステップは、左因子行列Ｗ更新部が、前記左因子行列Ｗのｋ番目の列ベクトルｗｋの各々について、前記データ入力部により受け付けた前記行列Ｘ及び前記欠損値行列Ｍと、前記左因子行列Ｗと、前記右因子行列Ｈとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記列ベクトルｗｋを更新するステップと、右因子行列Ｈ更新部が、前記右因子行列Ｈのｋ番目の列ベクトルｈｋの各々について、前記データ入力部により受け付けた前記行列Ｘ及び前記欠損値行列Ｍと、前記左因子行列Ｗと、前記右因子行列Ｈとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記列ベクトルｈｋを更新するステップと、繰り返し判定部が、予め定められた条件を満たすまで、前記左因子行列Ｗ更新部による前記列ベクトルｗｋの各々の更新、及び前記右因子行列Ｈ更新部による前記列ベクトルｈｋの各々の更新を繰り返させるステップと、を含むようにしてもよい。
【００１１】
また、第２の発明に係る解析方法において、前記左因子行列Ｗ更新部が更新するステップでは、前記列ベクトルｗｋを更新する際に、前記近接点行列Ｐと前記重み付き平均ベクトルｑとに基づいて、並列近接法のための補助変数Ｙを演算し、前記データ入力部により受け付けた前記行列Ｘ及び前記欠損値行列Ｍと、前記左因子行列Ｗと、前記右因子行列Ｈと、前記補助変数Ｙとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記列ベクトルｗｋを更新することを、所定回繰り返し、前記右因子行列Ｈ更新部が更新するステップでは、前記列ベクトルｈｋを更新する際に、前記近接点行列Ｐと前記重み付き平均ベクトルｑとに基づいて、並列近接法のための補助変数Ｙを演算し、前記データ入力部により受け付けた前記行列Ｘ及び前記欠損値行列Ｍと、前記左因子行列Ｗと、前記右因子行列Ｈと、前記補助変数Ｙとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、前記損失関数の値を最小化するように、前記列ベクトルｈｋを更新することを、所定回繰り返すようにしてもよい。
【００１２】
第３の発明に係るプログラムは、コンピュータを、上記第１の発明に係る解析装置の各部として機能させるためのプログラムである。
【発明の効果】
【００１３】
本発明の解析装置、方法、及びプログラムによれば、一般化カルバックライブラーダイバージェンス距離の和を用いて表される、非負値行列分解の損失関数の値を最小化するように、左因子行列Ｗ、及び右因子行列Ｈのパラメータを更新することを繰り返し、パラメータの更新において、行列Ｘ、欠損値行列Ｍ、左因子行列Ｗ、及び右因子行列Ｈに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、損失関数の値を最小化するように、左因子行列Ｗ、及び右因子行列Ｈのパラメータを更新することにより、精度よく、データを解析することができる、という効果が得られる。
【図面の簡単な説明】
【００１４】
【図１】本発明の実施の形態に係る解析装置の構成を示すブロック図である。
【図２】本発明の実施の形態に係る解析装置における解析処理ルーチンを示すフローチャートである。
【図３】本発明の実施の形態に係る解析装置における並列近接処理ルーチンを示すフローチャートである。
【発明を実施するための形態】
【００１５】
以下、図面を参照して本発明の実施の形態を詳細に説明する。
【００１６】
＜本発明の実施の形態に係る原理＞
【００１７】
まず、本発明の実施の形態における原理を説明する。
【００１８】
一般化ＫＬダイバージェンスを用いたＮＭＦは、入力された行列Ｘ、欠損値行列Ｍ、左因子行列Ｗ、及び右因子行列Ｈを用いて、以下（１）式の損失関数最小化問題を解くために左因子行列Ｗと右因子行列Ｈの更新を繰り返し行うものである。行列Ｘ、左因子行列Ｗ、及び右因子行列Ｈが非負の値のみを取ることを共用することが手法名の“非負値”の意味するところである。また属性の共起パターンは、因子数であるＫ個分、左因子行列Ｗ、右因子行列Ｈとして抽出される。左因子行列Ｗは各行属性の共起、右因子行列Ｈは各列属性の共起を表す。
【００１９】
【数１】
2017111754.tif 000003
【００２０】
このとき、ｗｋは左因子行列Ｗの第ｋ列ベクトル、ｈｋは右因子行列Ｈの第ｋ列ベクトル、Ｄは行列間の一般化ＫＬダイバージェンス距離、
【００２１】
【数２】
2017111754.tif 000004
【００２２】
は行列の要素毎の積である。また、要素ごとの一般化ＫＬダイバージェンス距離は以下（２）式のように定められる。
【００２３】
【数３】
2017111754.tif 000005
【００２４】
このときζは行列Ｘの要素値、ξはＷＨＴの要素値である。なお、行列間の一般化ＫＬダイバージェンス距離は要素ごとの一般化ＫＬダイバージェンス距離の和として定められる。
【００２５】
本発明の実施の形態に係る手法においては、参考文献１に記載の凸最適化技術である並列近接法を用いることで、上記の一般化ＫＬダイバージェンス距離の非リプシッツ連続性の問題を解決し、新たに収束が保証されたアルゴリズムを提案する。
【００２６】
［参考文献１］：P. L. Combettes and J.-C. Pesquet. A proximal decomposition method for solving convex variational inverse problems. Inverse problems, 24(6):065014, 2008.
【００２７】
＜本発明の実施の形態に係る解析装置の構成＞
【００２８】
次に、本発明の実施の形態に係る解析装置の構成について説明する。図１に示すように、本発明の実施の形態に係る解析装置１００は、ＣＰＵと、ＲＡＭと、後述する解析処理ルーチンを実行するためのプログラムや各種データを記憶したＲＯＭと、を含むコンピュータで構成することが出来る。この解析装置１００は、機能的には図１に示すようにデータ入力部１０と、演算部２０と、出力部５０とを備えている。
【００２９】
データ入力部１０は、非負の要素からなる行列Ｘ、及び行列Ｘにおける各要素が欠損値であるか又は観測値であるかを表す欠損値行列Ｍを受け付ける。行列Ｘの行と列の長さをＩ×Ｊとし、行列Ｘの要素は非負（０以上）の値のみを許容する。欠損値行列Ｍは行と列の長さが、行列Ｘと同じＩ×Ｊの場合のみを許容し、行列Ｍの要素は０あるいは１のみに制限する。値が０の要素は欠損値、１の要素は観測値として扱う。
【００３０】
また、データ入力部１０は、パラメータの初期値として、左因子行列Ｗと、右因子行列Ｈとの初期値の入力を受け付ける。左因子行列Ｗの行と列の長さはＩ×Ｋとし、右因子行列Ｈの行と列の長さはＪ×Ｋとする。左因子行列Ｗ、及び右因子行列Ｈの要素は非負（０以上）の値のみを許容する。また、データ入力部１０は、最大パラメータ更新繰り返し演算回数Ｎ、及び最大並列近接法繰り返し演算回数Ｔの入力を受け付ける。
【００３１】
演算部２０は、パラメータ更新部３０と、パラメータ出力部３８とを含んで構成されている。
【００３２】
パラメータ更新部３０は、左因子行列Ｗ更新部３２と、右因子行列Ｈ更新部３４と、繰り返し判定部３６とを含んで構成されている。
【００３３】
パラメータ更新部３０は、上記各部の処理によって、データ入力部１０により受け付けた行列Ｘ、及び欠損値行列Ｍに基づいて、非負の要素からなるＩ×Ｋの左因子行列Ｗと、非負の要素からなるＪ×Ｋの右因子行列Ｈとの積の要素であって、かつ、観測値である要素と、行列Ｘの要素との距離を表す一般化ＫＬダイバージェンス距離の和を用いて表される、非負値行列分解の損失関数の値を最小化するように、左因子行列Ｗ、及び右因子行列Ｈのパラメータを更新することを繰り返す。ここでは、行列Ｘ、欠損値行列Ｍ、左因子行列Ｗ、及び右因子行列Ｈに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、損失関数の値を最小化するように、左因子行列Ｗ、及び前記右因子行列Ｈのパラメータを更新する。
【００３４】
左因子行列Ｗ更新部３２は、左因子行列Ｗのｋ番目の列ベクトルｗｋの各々について、データ入力部１０により受け付けた行列Ｘ及び欠損値行列Ｍと、左因子行列Ｗと、右因子行列Ｈとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、損失関数の値を最小化するように、列ベクトルｗｋを更新する。また、左因子行列Ｗ更新部３２は、列ベクトルｗｋを更新する際に、近接点行列Ｐと重み付き平均ベクトルｑとに基づいて、並列近接法のための補助変数Ｙを演算し、データ入力部１０により受け付けた行列Ｘ及び欠損値行列Ｍと、左因子行列Ｗと、右因子行列Ｈと、補助変数Ｙとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、損失関数の値を最小化するように、列ベクトルｗｋを更新することを、所定回繰り返す。
【００３５】
左因子行列Ｗ更新部３２の処理の詳細を以下に説明する。左因子行列Ｗ更新部３２は、左因子行列Ｗの第ｋ列ベクトルを以下の計算法により、左因子行列Ｗの列ベクトルの各々（ｋ＝１からｋ＝Ｋまで）について更新することを、最大パラメータ更新繰り返し演算回数Ｎの回数分繰り返す。
【００３６】
左因子行列Ｗの第ｋ番目の列ベクトルｗｋに関して、ＮＭＦの損失関数は、行列Ｘ、欠損値行列Ｍ、左因子行列Ｗ、及び右因子行列Ｈに基づいて、以下（３）式のように定義される。
【００３７】
【数４】
2017111754.tif 000006
【００３８】
ここで、
【００３９】
【数５】
2017111754.tif 000007
【００４０】
はベクトルｗｋの少なくとも１つの要素が負であれば無限大を返す関数である。
【００４１】
上記（３）式の損失関数は、以下（４）式のように簡略化される。
【００４２】
【数６】
2017111754.tif 000008
【００４３】
このとき関数ｆｊ（ｗｋ）は、以下（５）式により表される。
【００４４】
【数７】
2017111754.tif 000009
【００４５】
本実施の形態では、上記（４）式において示した最小化問題を並列近接法によって解く。並列近接法はＭｏｒｅａｕ包という技術を使うことで一般化ＫＬダイバージェンスをリプシッツ連続にした複数の関数によって、それぞれの関数に対応した近接点Ｐｔが得られる。そして、近接点Ｐｔを得た上で最小化問題を解くことで非リプシッツ連続性の問題を回避することが出来る（上記参考文献１参照）。
【００４６】
具体的には、ｎ回目（ｎ≦Ｎ）の第ｋ番目の列ベクトルｗｋの更新において、並列近接法を適用し、以下の処理をＴ回繰り返す。
【００４７】
まず、並列近接法のためのＩ×Ｊの補助変数Ｙの各列をｗｋによって初期化する。ｔ回目（ｔ≦Ｔ）の繰り返しにおいて、並列近接法によって近接点Ｐｔを得る。
【００４８】
並列近接法は、行列Ｘ、欠損値行列Ｍ、左因子行列Ｗ、右因子行列Ｈ、補助変数Ｙ、及び重み行列Ｓに基づいて、以下（６）式の操作によって近接点行列Ｐｔを演算する。
【００４９】
【数８】
2017111754.tif 000010
【００５０】
このとき、
【００５１】
【数９】
2017111754.tif 000011
【００５２】
は要素毎の割り算、
【００５３】
【数１０】
2017111754.tif 000012
【００５４】
である。
【００５５】
上記（６）式の演算により、補助変数の各要素ｙｔに対する近接点行列Ｐｔのｊ番目の要素ｐj,tを得る。
【００５６】
【数１１】
2017111754.tif 000013
【００５７】
次に近接点行列Ｐｔの重み付き平均ベクトルｑｔを、以下（７）式に従って、予め定められた重み行列Ｓの各要素との積をとることにより求める。
【００５８】
【数１２】
2017111754.tif 000014
【００５９】
このときｅＪは長さＪの単位行列とする。
【００６０】
最後にｔ回目の第ｋ番目の列ベクトルｗｋの推定結果を以下（８）式に従って得る。
【００６１】
【数１３】
2017111754.tif 000015
【００６２】
ここで、λは正の値を取るステップサイズ変数である。
【００６３】
なお、ｔ＋１回目の演算では以下（９）式の補助変数Ｙｔ＋１を求める。
【００６４】
【数１４】
2017111754.tif 000016
【００６５】
このときＱｔはベクトルｑｔを並べた行列
【００６６】
【数１５】
2017111754.tif 000017
【００６７】
である。
【００６８】
そして、左因子行列Ｗ更新部３２は、ｎ回目のパラメータ更新において、ｋ＝１からｋ＝Ｋまで、左因子行列Ｗの第ｋ番目の列ベクトルｗｋのパラメータを、上記（８）式で得られた並列近接法により得られたｗｋに更新する。
【００６９】
右因子行列Ｈ更新部３４は、右因子行列Ｈのｋ番目の列ベクトルｈｋの各々について、データ入力部１０により受け付けた行列Ｘ及び欠損値行列Ｍと、左因子行列Ｗと、右因子行列Ｈとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、損失関数の値を最小化するように、列ベクトルｈｋを更新する。また、右因子行列Ｈ更新部３４は、列ベクトルｈｋを更新する際に、近接点行列Ｐと重み付き平均ベクトルｑとに基づいて、並列近接法のための補助変数Ｙを演算し、データ入力部１０により受け付けた行列Ｘ及び欠損値行列Ｍと、左因子行列Ｗと、右因子行列Ｈと、補助変数Ｙとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、損失関数の値を最小化するように、列ベクトルｈｋを更新することを、所定回繰り返す。
【００７０】
右因子行列Ｈ更新部３４では、上記（３）〜（９）式の演算においてＷとＨを入れ替えて実行する。つまり、左因子行列Ｗの第ｋ番目の列ベクトルｗｋを、右因子行列Ｈの第ｋ番目の列ベクトルｈｋに置き換え、右因子行列Ｈの第ｋ番目の列ベクトルｈｋを、左因子行列Ｗの第ｋ番目の列ベクトルｗｋに置き換えて実行する。
【００７１】
右因子行列Ｈ更新部３４の処理の詳細を以下に説明する。右因子行列Ｈ更新部３４は、右因子行列Ｈの第ｋ列ベクトルを以下の計算法により、右因子行列Ｈの列ベクトルの各々（ｋ＝１からｋ＝Ｋまで）について更新することを、最大パラメータ更新繰り返し演算回数Ｎの回数分繰り返す。
【００７２】
右因子行列Ｈの第ｋ番目の列ベクトルｈｋに関して、損失関数は、上記（３）式、（４）式において表される。
【００７３】
このとき、上記（５）式により表される関数ｆｊ（ｈｋ）におけるａｊ，ｋ、及びｂｊ，ｋは、それぞれ
【００７４】
【数１６】
2017111754.tif 000018
【００７５】
である。
【００７６】
本実施の形態では、上記（４）式において示した最小化問題を並列近接法によって解く。
【００７７】
具体的には、ｎ回目（ｎ≦Ｎ）の第ｋ番目の列ベクトルｈｋの更新において、並列近接法を適用し、以下の処理をＴ回繰り返す。
【００７８】
まず、並列近接法のためのＩ×Ｊの補助変数Ｙの各列をｈｋによって初期化する。ｔ回目の繰り返しにおいて、並列近接法によって近接点Ｐｔを得る。
【００７９】
並列近接法は、行列Ｘ、欠損値行列Ｍ、左因子行列Ｗ、右因子行列Ｈ、補助変数Ｙ、及び重み行列Ｓに基づいて、上記（６）式の操作によって演算する。
【００８０】
上記（６）式の演算により、補助変数の各要素ｙｔに対する近接点Ｐｔのｊ番目の要素を得る。
【００８１】
次に近接点Ｐｔの重み付き平均ベクトルｑｔを、以下（７）式に従って、重み行列Ｓの各要素との積をとることにより求める。
【００８２】
最後にｔ回目の第ｋ番目の列ベクトルｈｋの推定結果を以下（１０）式に従って得る。
【００８３】
【数１７】
2017111754.tif 000019
【００８４】
そして、右因子行列Ｈ更新部３４は、ｎ回目のパラメータ更新において、ｋ＝１からｋ＝Ｋまで、右因子行列Ｈの第ｋ番目の列ベクトルｈｋのパラメータを、上記（１０）式で得られたｈｋに更新する。
【００８５】
繰り返し判定部３６は、予め定められた条件を満たすまで、左因子行列Ｗ更新部３２による列ベクトルｗｋの各々の更新、及び右因子行列Ｈ更新部３４による列ベクトルｈｋの各々の更新を繰り返させる。本実施の形態では、予め定められた条件としては、最大パラメータ更新繰り返し演算回数Ｎの回数分繰り返すこととする。
【００８６】
パラメータ出力部３８は、パラメータ更新部３０によりＮ回更新された左因子行列Ｗ、及び右因子行列Ｈを出力する。
【００８７】
＜本発明の実施の形態に係る解析装置の作用＞
【００８８】
次に、本発明の実施の形態に係る解析装置１００の作用について説明する。データ入力部１０において行列Ｘ、欠損値行列Ｍ、左因子行列Ｗと、右因子行列Ｈとの初期値、最大パラメータ更新繰り返し演算回数Ｎ、及び最大並列近接法繰り返し演算回数Ｔを受け付けると、解析装置１００は、図２に示す解析処理ルーチンを実行する。
【００８９】
まず、ステップＳ１００では、データ入力部１０において受け付けた、左因子行列Ｗと、右因子行列Ｈとの初期値に基づいて、左因子行列Ｗ、及び右因子行列Ｈを初期化する。また、ｎをｎ＝１と初期化する。
【００９０】
次に、ステップＳ１０２では、データ入力部１０により受け付けた行列Ｘ及び欠損値行列Ｍと、左因子行列Ｗと、右因子行列Ｈと、並列近接法のための補助変数Ｙとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、上記（４）式の損失関数の値を最小化するように、列ベクトルｗｋの各々のパラメータを更新する。
【００９１】
ステップＳ１０４では、データ入力部１０により受け付けた行列Ｘ及び欠損値行列Ｍと、左因子行列Ｗと、右因子行列Ｈと、補助変数Ｙとに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、上記（４）式の損失関数の値を最小化するように、列ベクトルｈｋの各々のパラメータを更新する。
【００９２】
ステップＳ１０６では、ｎ＝Ｎであるかを判定し、ｎ＝ＮであればステップＳ１１０へ移行し、ｎ≠Ｎであれば、ステップＳ１０８へ移行してｎ＝ｎ＋１とカウントアップし、ステップＳ１０２に戻って処理を繰り返す。
【００９３】
ステップＳ１１０では、ステップＳ１０２、及びＳ１０４で更新された因子行列Ｗと、右因子行列Ｈとを出力し、処理を終了する。
【００９４】
上記ステップＳ１０２、及びＳ１０４は、図３に示す並列近接処理ルーチンによって実行される。以下では、ｋ番目の列ベクトルｗｋを更新する場合について説明する。
【００９５】
ステップＳ２００では、ｋ番目の列ベクトルｗｋを選択する。
ステップＳ２０２では、ｔ＝１と初期化する。また、補助変数Ｙを列ベクトルｗｋによって初期化する。
【００９６】
ステップＳ２０４では、ｔ＝２以降の場合、前回得られた近接点行列Ｐｔ、及び後述するステップＳ２０６で求められたｑｔに基づいて、上記（９）式により補助変数Ｙを求める。そして、行列Ｘ、欠損値行列Ｍ、左因子行列Ｗ、右因子行列Ｈ、補助変数Ｙ、及び重み行列Ｓに基づいて、上記（６）式に従って、ｔ回目の近接点行列Ｐｔを得る。
【００９７】
ステップＳ２０６では、ステップＳ２０４で得た近接点行列Ｐｔについて、上記（７）式に従って、重み行列Ｓの各要素との積をとることにより、近接点行列Ｐｔの重み付き平均ベクトルｑｔを求める。
【００９８】
ステップＳ２０８では、ステップＳ２０６で求められたｑｔを用いて、上記（８）式に従って、第ｋ番目の列ベクトルｗｋの推定結果を得る。
【００９９】
ステップＳ２１０では、ｔ＝Ｔであるかを判定し、ｔ＝ＴであればステップＳ２１４へ移行し、ｔ＝ＴでなければステップＳ２１２へ移行してｔ＝ｔ＋１とカウントアップし、ステップＳ２０４に戻って処理を繰り返す。
【０１００】
ステップＳ２１４では、ｋ＝Ｋであるかを判定し、ｋ＝Ｋであれば処理を終了し、ｋ＝ＫでなければステップＳ２１６へ移行してｋ＝ｋ＋１とカウントアップし、ステップＳ２００に戻って処理を繰り返す。
【０１０１】
上記ステップＳ１０４における、ｋ番目の列ベクトルｈｋを更新する処理についても、上記図３に示す並列近接処理ルーチンと同様の処理によって実行される。
【０１０２】
＜実験例における効果＞
【０１０３】
人工的に生成した欠損値を含む行列データで既存技術と本発明の実施の形態に係る手法との性能差を比較する。行列データの行と列の長さはＩ＝Ｊ＝１００、潜在因子数はＫ＝５とし、観測値の割合（ｄｅｎｓｉｔｙ）は１．００、０．７０、０．５０、０．３０、及び０．１０の場合を採用した。人工的なデータは、実データと似た傾向を持つようにガンマ分布から左因子行列Ｗと右因子行列Ｈを生成し、その内積ＷＨＴを平均値とするポアソン分布から行列Ｘを作成した。
【０１０４】
各手法の性能を比較するために欠損値の推定精度を利用する。実験では観測値の９割を学習データとして因子行列を学習し、因子行列の内積を推定値とした。既存法として、ＭＵ（上記非特許文献５）、ＣＯＤ（上記非特許文献６）、ＡＤＭＭ（上記非特許文献７）を利用した。なお、ＣＯＤとＡＤＭＭを欠損値に対応させるために非特許文献８の技術を用いている。
【０１０５】
［非特許文献８］：Y. Xu,W. Yin, Z.Wen, and Y. Zhang. An alternating direction algorithm for matrix completion with nonnegative factors. Frontiers of Mathematics in China, 7(2):365-384, 2012.
【０１０６】
残りの１割の観測値をテストデータとし、実際の値と、推定値の一般化ＫＬダイバージェンスの平均値を得た。各設定での推定精度を以下の表に示す。値が低いほど良い推定精度となる。本実験において、本実施の形態に係る手法は、観測値の割合が０．７０、０．５０、０．３０、及び０．１０の場合で最も良い推定精度を示した。また、観測値の割合が１．００の場合は既存のＣＯＤが良い性能を示した。以上から、本実施の形態に係る手法による推定精度の向上が確認できる。
【０１０７】
【表１】
2017111754.tif 000020
【０１０８】
以上説明したように、本発明の実施の形態に係る解析装置によれば、一般化カルバックライブラーダイバージェンス距離の和を用いて表される、非負値行列分解の損失関数の値を最小化するように、左因子行列Ｗ、及び右因子行列Ｈのパラメータを更新することを繰り返し、パラメータの更新において、行列Ｘ、欠損値行列Ｍ、左因子行列Ｗ、及び右因子行列Ｈに基づいて並列近接法により求めた近接点行列Ｐの重み付き平均ベクトルｑを用いて、損失関数の値を最小化するように、左因子行列Ｗ、及び右因子行列Ｈのパラメータを更新することにより、精度よく、データを解析することができる。
【０１０９】
なお、本発明は、上述した実施の形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。
【符号の説明】
【０１１０】
１０   データ入力部
２０   演算部
３０   パラメータ更新部
３２   左因子行列Ｗ更新部
３４   右因子行列Ｈ更新部
３６   繰り返し判定部
３８   パラメータ出力部
５０   出力部
１００ 解析装置
【図１】
2017111754.tif 000021
【図２】
2017111754.tif 000022
【図３】
2017111754.tif 000023