JP6810003B2 - 行列単純化装置、プログラム、および行列単純化方法 - Google Patents

Description

本発明は、行列単純化装置、プログラム、および行列単純化方法に関する。

データや物理現象に潜在する低ランク性に基づくデータ解析手法（以降，低ランクモデリング）が、コンピュータビジョン、画像処理、ゲノムデータ解析などの多くの分野で近年活発に研究されている。ここで、低ランク性とは、データや物理現象から導かれる行列Ｘについて、その階数ｒａｎｋ（Ｘ）がＸのサイズ（ここでは、行列Ｘの行数または列数のうちの小さい方）と比較して小さい性質である。低ランクモデリングでは、着目した行列の階数(rank) を最小化する数理計画問題を立式し、この問題の解を求めて所望の解析を実現する。

関数rankは不連続で微分不能、非凸であり、これに基づく計画問題はＮＰ困難な組合せ最適化となる。このため、rank関数の代わりに核型ノルムを正則化する緩和アプローチが広く用いられる。核型ノルムはrank関数の凸包絡であるため、核型ノルムを最小化することで間接的に低ランク性を高められる。核型ノルムを含む最適化問題は、Alternating Direction Method of Multipliers（ＡＤＭＭ）等の１次法の反復計算で解を求められ、核型ノルムの近接写像である特異値閾値処理(Singular Value Thresholding，ＳＶＴ）が繰り返し実行される。

しかし、ＳＶＴ算出には計算量の大きい特異値分解（Singular Value Decomposition, ＳＶＤ) の算出が必要なため、解析結果を得るのに多くの計算時間を要する。ここで、各解析法で立式される核型ノルム正則化問題を、（１）少数の大型行列を正則化する問題と、（２）多数の小型行列を正則化する問題に分類する。前者（１）の用途は、例えば、ロバスト主成分分析や、欠損地推定・補間や、オプティカルフロー推定や、ダイナミックＭＲＩ解析や、ゲノム解析等である。また、後者（２）の用途は、グラフ単純化や、偽色除去等である。

前者（１）の問題の場合には、計算量を抑えてＳＶＴを高速化する手法がいくつか提案されている。J.F.Caiらは、行列を事前にComplete Orthogonal Decomposition（ＣＯＤ，完全直交分解）した後にニュートン法で反復更新し、ＳＶＤを行わずにＳＶＴを求める高速ＳＶＴ（Fast SVT，ＦＳＶＴ）を提案している。またT.H.Ohらは、大型行列を直交行列と小型のコア行列の積に近似することで、ＳＶＤの入力サイズを小さくして高速化する高速ランダム化ＳＶＴ（Fast Randomized SVT，ＦＲＳＶＴ）を提案している。ここに挙げるいずれの手法も、入力行列のサイズが大きい(行数および列数がそれぞれ５００〜２０００程度）ときに計算量を抑え、大幅な速度改善結果を示す。

一方で後者（２）の問題の場合、上記手法による高速化の効果は限定的であると考えられる。ＦＳＶＴを用いると、入力が小型の場合、直接のＳＶＤ計算と比較してＣＯＤとニュートン法の計算量が大きいという問題がある。また、ＦＲＳＶＴを用いると、入力が小型の場合、コア行列の縮退効果が小さいために高速化できず、また近似法であるため計算誤差が大きいという問題がある。加えて、これらの手法は多数の行列を同時に処理するデータ並列のアプローチを取れず、近年の並列アーキテクチャの計算資源を活用できないという問題もある。

非特許文献１では、計算量を抑えながらデータ並列にＳＶＴを算出する高速マルチプルＳＶＴ（Fast Multiple SVT，ＦＭＳＶＴ）が提案されている。ＦＭＳＶＴでは、特異値を用いずにＬ_∞，２混合ノルムで核型ノルムを表現することで、ＳＶＤが不要なＳＶＴ計算を実現し、計算量を削減している。加えて、データを並列に処理する並列化アルゴリズムを容易に実現でき、多数の行列について同時処理が可能である。

佐々木崇元，北原正樹，清水淳，「領域情報符号化における核型ノルム最適化の高速計算法」，第３１回画像符号化シンポジウム（ＰＣＳＪ２０１６），ｐｐ．１４０−１４１，２０１６年１１月．

しかしながら、非特許文献１の技術では、対象とする行列のサイズが２×２に限られていた。より大きなサイズの行列に関して、高速にＳＶＴを算出することのできる技術が求められている。

本発明は、上記の課題認識に基づいて行なわれたものであり、より大きなサイズの行列を対象として、少ない計算量で高速に、行列の低ランク化をすることのできる行列単純化装置、プログラム、および行列単純化方法を提供しようとするものである。

［１］本発明の一態様は、Ｍ行２列の入力行列から各列に対応するＭ次元の第１ベクトルおよび第２ベクトルを抽出し、または２行Ｍ列の入力行列から各行に対応するＭ次元の第１ベクトルおよび第２ベクトルを抽出する（ただし、Ｍは２以上の整数）ベクトル化部と、前記第１ベクトルおよび前記第２ベクトルを含む平面内で前記第２ベクトルを第１の方向に所定角回転させて得られる第３ベクトルと前記第１ベクトルとの距離である第１距離と、前記平面内で前記第２ベクトルを第２の方向に前記所定角回転させて得られる第４ベクトルと前記第１ベクトルとの距離である第２距離とを求め、前記第１距離および前記第２距離に基づいて前記入力行列をより低ランクで近似した低ランク近似行列を求める低ランク近似化部と、を具備する行列単純化装置である。

［２］本発明の一態様は、上記の行列単純化装置であって、前記入力行列ＹがＭ行２列であるときにはＹ＝［ｙ_１，ｙ_２］とし、前記入力行列Ｙが２行Ｍ列であるときにはＹ^Ｔ＝［ｙ_１，ｙ_２］とし、前記第１ベクトルをｙ_１とし、前記第２ベクトルをｙ_２としたとき、前記低ランク近似化部は、前記所定角回転させる回転行列をＲとして、［ｙ_２，−ｙ_１］＝Ｙ（Ｙの上にバー）として、内分比パラメータδ（０≦δ≦１）を用いて、前記入力行列Ｙと行列ＲＹ（Ｙの上にバー）との前記内分比パラメータδによる内分である内分行列と、振幅パラメータγと、に基づいて前記低ランク近似行列を求める。

［３］本発明の一態様は、上記の行列単純化装置であって、前記入力行列ＹがＭ行２列であるときにはＹ＝［ｙ_１，ｙ_２］とし、前記入力行列Ｙが２行Ｍ列であるときにはＹ^Ｔ＝［ｙ_１，ｙ_２］とし、前記第１ベクトルをｙ_１とし、前記第２ベクトルをｙ_２としたとき、前記低ランク近似化部は、前記所定角回転させる回転行列Ｒを式（１４）として前記第２ベクトルを回転させる。
なお、式（１４）自体は、後の実施形態において記載する。

［４］本発明の一態様は、上記の行列単純化装置であって、前記入力行列Ｙの特異値をσ_１，σ_２（σ_１≧σ_２≧０）として、閾値をμとしたとき、前記低ランク近似化部は、前記行列Ｙの階数が２であり且つσ_１≠σ_２の場合には、式（１２）および式（１３）により、式（１２）の右辺を求めて前記低ランク近似行列とするものであり、且つ、前記低ランク近似化部は、式（１５）の右辺の計算を行うことによって、式（１２）内のＲＹ（Ｙの上にバー）を求める。
なお、式（１２）、式（１３）、式（１５）自体は、後の実施形態において記載する。

［５］本発明の一態様は、上記の行列単純化装置であって、前記低ランク近似化部は、「入力行列Ｙの階数が２且つσ_１≠σ_２」以外の場合には、式（１９）の右辺を求めて前記低ランク近似行列とするものである。
なお、式（１９）自体は、後の実施形態において記載する。

［６］本発明の一態様は、上記の行列単純化装置であって、前記ベクトル化部は複数の前記入力行列を基にそれぞれの前記入力行列の前記第１ベクトルおよび前記第２ベクトルを抽出するものであり、前記低ランク近似化部は、各入力行列から抽出された前記第１ベクトルおよび前記第２ベクトルを用いて前記低ランク近似行列を求めるものであり、単一命令列を、前記入力行列にそれぞれ対応する複数のデータに適用して並列処理を行うものである。

［７］本発明の一態様は、コンピューターを、Ｍ行２列の入力行列から各列に対応するＭ次元の第１ベクトルおよび第２ベクトルを抽出し、または２行Ｍ列の入力行列から各行に対応するＭ次元の第１ベクトルおよび第２ベクトルを抽出する（ただし、Ｍは２以上の整数）ベクトル化部と、前記第１ベクトルおよび前記第２ベクトルを含む平面内で前記第２ベクトルを第１の方向に所定角回転させて得られる第３ベクトルと前記第１ベクトルとの距離である第１距離と、前記平面内で前記第２ベクトルを第２の方向に前記所定角回転させて得られる第４ベクトルと前記第１ベクトルとの距離である第２距離とを求め、前記第１距離および前記第２距離に基づいて前記入力行列をより低ランクで近似した低ランク近似行列を求める低ランク近似化部と、として機能させるためのプログラムである。

［８］本発明の一態様は、Ｍ行２列の入力行列から各列に対応するＭ次元の第１ベクトルおよび第２ベクトルを抽出し、または２行Ｍ列の入力行列から各行に対応するＭ次元の第１ベクトルおよび第２ベクトルを抽出する（ただし、Ｍは２以上の整数）ベクトル化過程と、前記第１ベクトルおよび前記第２ベクトルを含む平面内で前記第２ベクトルを第１の方向に所定角回転させて得られる第３ベクトルと前記第１ベクトルとの距離である第１距離と、前記平面内で前記第２ベクトルを第２の方向に前記所定角回転させて得られる第４ベクトルと前記第１ベクトルとの距離である第２距離とを求め、前記第１距離および前記第２距離に基づいて前記入力行列をより低ランクで近似した低ランク近似行列を求める低ランク近似化過程と、を含む行列単純化方法である。

本発明によれば、より大きなサイズの行列を対象として、少ない計算量で高速に、行列の低ランク化をすることが可能となる。

本発明の実施形態を説明する図であり、空間ＩｍＹにおけるベクトルｙ_１，ｙ_２，Ｒｙ_２，−Ｒｙ_２の関係を示すグラフである。 同実施形態における、閾値μに応じた係数γ（１−δ），γδの値を示すグラフである。 同実施形態における、閾値μの値に応じた値（ＳＶＴ，下の式の通り）の軌跡を示すグラフである。 アルゴリズム１（Ｍ≧３の場合のＳＶＴ算出法）を示す概略図である。 アルゴリズム２（Ｍ＝２の場合のＳＶＴ算出法）を示す概略図である。 アルゴリズム３（Ｌ個の行列Ｙ_１，Ｙ_２，・・・，Ｙ_ＬについてのＳＶＴ算出法）を示す概略図である。 第１実施形態による行列単純化装置の概略機能構成を示すブロック図である。 第１実施形態による行列単純化装置が処理する主要なデータの構成を示す概略図である。 第２実施形態による行列単純化装置の概略機能構成を示すブロック図である。 第３実施形態による行列単純化装置の概略機能構成を示すブロック図である。

次に、本発明の複数の実施形態について、図面を参照しながら説明する。以下において、全実施形態に共通する低ランクモデリングについて説明し、その後で各実施形態による装置等の具体的な構成を説明する。

［１．核型ノルム正則化に基づく低ランクモデリング］
本実施形態が対象とする核型ノルム正則化問題は、次に説明する通りの問題である。
行列Ｘ∈Ｒ^Ｍ×Ｎ（ここで、Ｒは実数の集合。Ｒ^Ｍ×Ｎは実数を要素とするＭ行Ｎ列の行列の集合。以下においても同様。）の特異値をσ_ｉ（Ｘ），ｉ＝１，２，・・・，Ｋとする。ただし、Ｋ＝ｍｉｎ（Ｍ，Ｎ）であり、ｍｉｎ（）は引数の最小値を返す関数である。なお、以後において、Ｍ行Ｎ列の行列のサイズを単に「Ｍ×Ｎ」と表す場合がある。
このとき、核型ノルムは特異値の和として定義される。即ち、下の式（１）の通りである。

核型ノルムは関数ｒａｎｋ（Ｘ）の凸包絡であるため、核型ノルムを正則化することで低ランク性を推進できる。

本実施形態は、多数の小型行列が核型ノルムで正則化される最適化問題を取り扱う。典型的にはヒルベルト空間Χ（カイ）の変数χ∈Χ（カイ）に関する最適化問題であり、この問題は下の式（２）の形式を有する。

式（２）で、目的関数の第２項は核型ノルムによる低ランク性の正則化関数であり、第１項はその他の正則化や忠実化関数である。関数Φ_ｉ：Χ（カイ）→Ｒ^Ｍ×Ｎは低ランク性を規約する行列の生成関数で、典型的には線形写像である。

例えば、グラフ単純化の問題では、Ｍ次元空間に埋め込まれたグラフＧ＝（Υ，Ε）の頂点座標χ∈Ｒ^{Ｍ×｜Υ｜}を変数とし、隣接辺のベクトルを２つ１組で行列に押し込め、Φ_ｉ（χ）∈Ｒ^Ｍ×２を作成する．この行列の階数はグラフの曲折回数に比例することに着眼し、核型ノルムの和を正則化することで、各頂点を局所線形に整列させている。なお、核型ノルムの和は、次の通り表される。

背景技術において述べた「（２）多数の小型行列を正則化する問題」を解くためには、関数ｆと線形写像Φ_ｉに応じたアルゴリズムを選択する。そのアルゴリズムとして、文献［佐々木崇元，谷田隆一，清水淳，グラフ信号の局所線形近似によるグラフ形状単純化，第１５回情報科学技術フォーラム(ＦＩＴ２０１６)，第３分冊，ｐｐ．１−４，Sept.2016．］や、文献［S. Ono and I. Yamada, Color-line regularization for color artifact removal，IEEE Transactions on Computational Imaging, vol.2, no.3, pp.204-217, 2016.］では、関数ｆが微分可能あるいは近接写像計算可能（proximable）であるため、ＡＤＭＭやPrimal Dual Splitting等の１次法を用いている。

なお、ここで「近接写像計算可能」の意味は次の通りである。即ち、関数ｇと正実数μ＞０に対し、点ｙ（太字）での近接写像ｐｒｏｘ_μｇは、下の式の様に定義される。

核型ノルムを含むいくつかの関数については近接写像の効率的な計算法が知られており、そのような関数を近接写像計算可能であると呼ぶ。

上記の１次法は、最適解への収束列を生成する反復計算手続きで、各反復では関数の勾配や近接写像などの１次情報に基づき解を更新し、目的関数を最小化する。

（２）の「多数の小型行列を正則化する問題」の核型ノルム正則化関数については、補助変数Ｙ_ｉ＝Φ_ｉ（χ），ｉ＝１，２，・・・，Ｌを導入し、下の式で表されるｇ（Ｙ_１，Ｙ_２，・・・，Ｙ_Ｌ）の近接写像を計算すれば良い（「Ｙ_ｉ」の「Ｙ」は太字。以下においても同様。）。

この右辺の和を構成する各要素は、独立に計算可能であり、下の式（３）が成立する。

右辺各要素は、核型ノルムの近接写像で、特異値閾値処理（ＳＶＴ）に等しい。核型ノルムは、上の式（３）において、次の様に表されている。

ＳＶＴは、入力行列をＹ、閾値をμ＞０として、以下の式（４）および式（５）で算出される。

ここで、関数ＳＶＤ（・）は、特異値分解（ＳＶＤ）である。また、（・）_＋は、入力の各要素を非負値にクリッピングするランプ関数である。
なお、本実施形態での特異値分解は、「thin SVD」である。即ち、Ｙ＝ＵΣＶ^Ｔと分解した時、行列Ｕ，ＶはＫ＝ｍｉｎ（Ｍ，Ｎ）個の正規直交なベクトルから成り、（Ｋ＋１）番目以降の特異ベクトルの算出は省略する（Ｕ∈Ｒ^Ｍ×Ｋ，Ｖ∈Ｒ^Ｎ×Ｋ）。また、行列ΣはＫ×Ｋ対角行列である。数値解析ソフトＭＡＴＬＡＢ（MathWorks社，米国）で記述する場合は、「[U,S,V]=svd(Y,'econ');」である。

即ち、関数ｇ（Ｙ_１，Ｙ_２，・・・，Ｙ_Ｌ）の近接写像（式（３））は、Ｌ回のＳＶＴ（式（４））により算出でき、それぞれのＳＶＴでは小型行列を入力とするＳＶＤ算出（式（５））が必要である。１次法で反復演算する際、計算時間の多くは計算量の大きいＳＶＤ算出に費やされる。

さて、式（３）では、Ｌ個の行列Ｙ_１，Ｙ_２，・・・，Ｙ_Ｌの各々の間における依存関係が無いため、タスク並列処理が可能である。つまり、行列Ｙ_１，Ｙ_２，・・・，Ｙ_Ｌの各々についてＬ並列での処理が可能である。なお、扱う問題において各行列のサイズが小さい場合には、処理割当やメモリロードのオーバーヘッドが相対的に大きく、並列化による改善の効果が限定的である場合もある。

そこで、オーバーヘッドが少ない，データ並列処理について説明する。アルゴリズムをデータ並列化できれば、Single Instruction Multiple Data（ＳＴＭＤ）等のデータ並列アーキテクチャを用いた実装により処理を高速化できる。しかし、ＳＶＤの算出は逐次的であり、且つ行列要素の参照位置や処理内容が入力Ｙ_ｉに依存するため、行列間で共通に
計算できる処理が少なく、データ並列化は本質的に困難である。

以上より，多数の小型行列を低ランクに正則化するモデルは最適解を高速に得るのが困難で、その原因は、ＳＶＤの多大な計算量と並列化の難しさにある。

［Ｆａｓｔ Ｍｕｌｔｉｐｌｅ ＳＶＴ］
ここでは、式（３）に示した多数のＳＶＴ計算を高速に算出するＦＭＳＶＴを導出する。この手法は核型ノルムがある部分空間上のベクトルの距離で特徴付けられるという、幾何的性質の発見に基づいて導かれる。この性質により、特異値を用いずに核型ノルムを表現できる（後で、命題１，系２に記載）。また、ＳＶＤを用いずにＳＶＴを表現できる（後で、命題３に記載）。そして、このＳＶＴは、ほとんどが線形変換で記述できる（後で、命題４に記載）ため、データ並列なアルゴリズムを導ける。また、定理５、定理６では、ＳＶＴ算出式を示す。そして、最後に並列化したアルゴリズムを示す。本実施形態で得られるＳＶＴ算出法では、式（４）、式（５）を用いた算出と比較して計算量が削減されており、かつデータ並列なアルゴリズムとして実行可能である。

従来技術（非特許文献１）が入力行列のサイズを２×２に限定していたのに対し、本実施形態は、入力行列のサイズをＭ×２および２×Ｎに拡張する。

なお、近接写像の計算において上の式の通りであるから、入力が２×Ｎのサイズの行列のＳＶＴの算出は、Ｎ×２のサイズの行列のＳＶＴ算出の前後に、転置処理を施すことで実現できる。従って以降では、縦長のＭ×２のサイズの行列のＳＶＴについて説明する。

［１．１ 核型ノルムの幾何的性質とSVT の新しい表現］
ここでは、ＦＭＳＶＴ導出の中核となる、核型ノルムの幾何的性質とＳＶＤ不要なＳＶＴの表現について説明する。本節で提示する表現についての証明は、後で記載（付録）する。以後において、単一の入力行列をＹ＝［ｙ_１，ｙ_２］∈Ｒ^Ｍ×２，ｙ_１，ｙ_２∈Ｒ^Ｍとし（つまり、ｙ_１，ｙ_２はそれぞれＭ次元の縦ベクトル）、その特異値をσ_１，σ_２（σ_１≧σ_２≧０）とする。

＜命題１＞：特異値和σ_１＋σ_２と特異値差σ_１−σ_２は、下の式（６）に表す通りです。なお、式（６）において、複号同順である。

式（６）におけるＲ（太字）は回転行列であり、この回転行列Ｒ（太字。以下においても同様。）は、像ＩｍＹ上のベクトルを、原点回りにＩｍＹに沿ってπ／２［ｒａｄ］回転させる。ただし回転方向の正負については、図１に示すように、ｙ_１からｙ_２に最短で辿り着く方向を回転の正方向（第１の方向）とし、その反対方向を負方向（第２の方向）とする。

なお、上記の像ＩｍＹについて、次の通りである。即ち、行列Ａ∈Ｒ^Ｍ×Ｎに対し、部分空間ＩｍＡ＝｛Ａｘ｜ｘ∈Ｒ^Ｎ｝⊂Ｒ^ＭをＡの像という．
また、上記の回転行列Ｒ（太字）に関して、任意のｙ∈ＩｍＹについて、下の式に示す条件（１）から条件（５）までが成立する。

なお、命題１の証明を、後で付録Ａにおいて説明する。命題１は、特異値の和や差が、ベクトルｙ_１と±Ｒｙ_２との間の、Ｌ_２距離であることを表しており、それは、図１にも示す通りである。また、式（６）の中間式は、トレースや行列式で構成されるものであり、Ｙ^ＴＹが２×２の行列であるため，容易に算出可能なものである。

図１は、空間ＩｍＹにおけるベクトルｙ_１，ｙ_２，Ｒｙ_２，−Ｒｙ_２の関係を示すグラフである。ベクトルｙ_１，ｙ_２は、行列Ｙをベクトル化したものである。θは、ベクトルｙ_１からベクトルｙ_２への回転角である。なお、空間ＩｍＹにおいて、ｙ_１からベクトルｙ_２に最短で辿りつく側の回転方向を正方向としている（図中で「＋」で示す方向）。Ｒは回転角π／２［ｒａｄ］の回転行列である。ベクトルＲｙ_２とベクトル−Ｒｙ_２とを破線で示している。また、ベクトルｙ_１とベクトルＲｙ_２との距離（Ｌ_２ノルム）は、特異値σ_１とσ_２の和である。また、ベクトルｙ_１とベクトル−Ｒｙ_２との距離（Ｌ_２ノルム）は、特異値σ_１とσ_２の差である。

命題１と式（１）（核型ノルムの定義式）とから、直ちに以下の系を求められる。
＜系２＞：核型ノルムは、下の式（７）で表現できる。なお、式（７）の左辺が核型ノルムの表記であり、右辺がその定義である。

ただし、行列Ｂ∈Ｒ^{２Ｍ×２Ｍ}（Ｂは太字。以下においても同様。）は、下の式（８）に示す通りである。

また、式（７）の右辺に現れる上の表記は、Ｌ_∞，２混合ノルムである。Ｌ_∞，２混合ノルムは、Ｌ_∞ノルムとＬ_２ノルムの合成関数である。この合成関数の入力を、下の式の通りとする。

上記の入力に対し、Ｌ_∞，２混合ノルムの値、即ち上記合成関数が返す値は、下の式の通りである。

また、式（７）における関数ｖｅｃは、行列をベクトル化する関数である。つまり、関数ｖｅｃは、入力行列Ｙ＝［ｙ_１，ｙ_２］∈Ｒ^Ｍ×２を並べ替えて、下の式で表される列ベクトルを出力する線形変換である。

なお、後の式（９）に現れるｖｅｃ^Ｔは、ｖｅｃの逆変換である。つまり、Ｙ＝ｖｅｃ^Ｔ（ｖｅｃ（Ｙ））である。

系２は，特異値を使わずに、線形変換Ｂｖｅｃ（・）とＬ_∞，２混合ノルムの合成で、核型ノルムを表現できることを表している。これは、ＳＶＤを用いることなく、核型ノルムの近接写像（即ちＳＶＴ）を表現できる可能性を示している。系２と、後の付録Ｂに記載の補題Ｂ１，Ｂ２とに基づいて、次の命題３が得られる。

＜命題３＞：行列ＹのＳＶＴについて、下の式（９）が成立する。

命題３の証明を、後の付録Ｂに記載する。
この命題３により、ＳＶＴを、線形変換Ｂｖｅｃ（・）とＬ_∞，２ノルムの近接写像の合成という新しい形式で表現することができた。なお、Ｌ_∞，２ノルムの近接写像は、下の式で表されるものである。

この近接写像は非線形であるが、入力依存の対角行列による線形変換として記述することができる。

＜命題４＞：下の式で表されるｘに対して、その下の式（１０），式（１１）が成立する。
なお、ｘ_ｉ＝０のとき、ｋ_ｉｘ_ｉ＝（０／０）０となるが、例外的にこのときは、ｋ_ｉｘ_ｉ＝０とする。

この命題４の証明は、Moreau直交分解とＬ_１，２混合ノルム球への射影に基づくものである。この証明を、後の付録Ｃに記載する。

以上により、係数ｋ_１，ｋ_２の算出を除いて、ＳＶＴは全て線形変換で表現できることを明らかにした。次のセクションでは、上記の式を展開し、ＳＶＴ算出式およびアルゴリズムについて説明する。

［１．２ ＳＶＴ算出法］
前セクションで説明した命題から、ＳＶＴを算出する下記の定理５が得られる。
＜定理５＞：rankＹ＝２、且つσ_１≠σ_２のとき、下の式（１２）が成立する。

上記の定理５は、式（９）のｐｒｏｘ計算を式（１０）で展開し、さらに式（６）と式（１１）とを適用して得られる。定理５は、単純な式変形と、場合分けとで確認できるものであり明らかであるので、証明を省略する。
定理５は、ＳＶＴを、入力行列Ｙと行列ＲＹ（Ｙの上にバー）との線形結合で算出できることを表している。
式（１２）における結合係数γ（１−δ）およびγδは、振幅パラメータγと内分比パラメータδから構成され、いずれもＳＶＴの閾値μの関数である．係数γ（１−δ），γδの関数形を、図２にプロットして示す。

図２は、閾値μに応じた係数γ（１−δ），γδの値を示すグラフである。γ（１−δ）のグラフを実線で示し、γδのグラフを破線で示している。
図示するように、μ＝０のとき、γ（１−δ）＝１である。０≦μ≦σ_２においてγ（１−δ）はリニアに変化し、μ＝σ_２のときγ（１−δ）＝σ_１／（σ_１＋σ_２）である。また、σ_２≦μ≦σ_１においてγ（１−δ）はリニアに変化し、μ＝σ_１のときγ（１−δ）＝０である。そして、σ_１≦μにおいてγ（１−δ）＝０である。
また、μ＝０のとき、γδ＝０である。０≦μ≦σ_２においてγδはリニアに変化し、μ＝σ_２のときγδ＝σ_２／（σ_１＋σ_２）である。また、σ_２≦μ≦σ_１においてγδはリニアに変化し、μ＝σ_１のときγδ＝０である。そして、σ_１≦μにおいてγδ＝０である。
図２に示した係数に基づき、μの増加に伴うＳＶＴの軌跡を図３に示す。図３は、μの値に応じた値（ＳＶＴ，下の式の通り）の軌跡を示すグラフである。

つまり、行列ＹのＳＶＴである行列Ｚ（Ｍ行２列）をベクトル化したものが、ベクトルｚ_１，ｚ_２である。
図３は、ベクトルｚ_１，ｚ_２（ｚは太字）を上の式のように定義したときの、ＩｍＹ上のベクトルｚ_１，ｚ_２の軌跡を表している。図示するように、０＜μ≦σ_２では、ＹとＲＹ（Ｙの上にバー）の内分となる。また、σ_２＜μ＜σ_１では振幅が線形に減少し、μ≧σ_１ではゼロ行列Ｏ（太字）となる。

定理５を適用するには、π／２［ｒａｄ］回転行列Ｒによる変換ＲＹ（Ｙの上にバー）の具体的な値を求める必要がある。この求め方を、Ｍ≧３の場合とＭ＝２の場合の２通りについて説明する。

［Ｍ≧３の場合］ 回転行列Ｒに関して先に説明した条件（１）から条件（５）までのすべての条件が満たされる必要がある。また、回転行列Ｒは、命題１で述べた回転方向への回転の作用を有するものである必要がある。式（１４）に示す行列Ｒは、これらの条件を満たすものである。なお、式（１４）において、「ｄｅｔ Ｘ」は、行列Ｘの行列式である。

ただし、式（１４）によって行列Ｒを算出してから行列積ＲＹ（Ｙの上にバー）を求める手順を撮った場合、Ｏ（Ｍ^２）のオーダーの計算量が必要である。そこで、計算順序を変えて、下の式（１５）による計算を行うようにする。

式（１５）による計算では、内積を先に展開する手順をとることができる。これにより、計算量のオーダーをＯ（Ｍ）に削減できる。さらに、式（１２）について、式（１５）を適用してＹの積で括る形とすると、下の式（１６）に変形できる。

この式（１６）によって計算することにより、さらに計算量を削減することができる。

［Ｍ＝２の場合］ 符号（正、負、０）に応じて＋１，−１，０のいずれかの値をとる符号関数ｓｇｎを用いると、下の式（１７）が成立する。

即ち、行列要素の入れ替えと、符号反転とによってＲＹ（Ｙの上にバー）を算出できる。つまり、下の式（１８）の通りＲＹ（Ｙの上にバー）を算出できる。

以上、定理５に基づいてＳＶＴを求める方法を説明した。なお、この定理は、ｒａｎｋＹ＝２且つσ_ｉ≠σ_２という条件を前提としている。この条件が満たされない場合には、μにより定まる振幅をＹに乗じて、下の定理６の通り算出できる。
＜定理６＞：「ｒａｎｋＹ＝２且つσ_ｉ≠σ_２」以外のとき、下の式（１９）の通りである。

なお、式（１９）に含まれる下の表現は、行列Ｙのフロベニウスノルムであり、即ち行列Ｙの要素の２乗和平方根である。

定理６の証明は、後の付録Ｄに記載する。以上説明したように、定理５および定理６を用いてＳＶＴを算出するためには、行列Ｙの階数ｒａｎｋＹと、特異値σ_ｉ，σ_２に基づいて適切に場合分けする必要がある。階数ｒａｎｋＹについては、Ｙがゼロ行列Ｏ（太字）であるか否か、またｄｅｔＹ^ＴＹが０であるか否かにより、下の表１の通り判別できる。また特異値については、式（６）の中間式よりσ_ｉ，σ_２を算出すれば良い。

以上の説明に基づき、次に述べるアルゴリズム１および２によって、ＳＶＴを算出することができる。アルゴリズム１はＭ≧３の場合に適用可能な手順である。アルゴリズム２はＭ＝２の場合に適用可能な手順である。これらのアルゴリズムにおいて、場合分けのために必要なｄｅｔＹ^ＴＹ，σ_ｉ−σ_２や、係数γ，δはすべて内積ｙ_ｉ ^Ｔｙ_ｊを用いて算出している。

図４は、アルゴリズム１を示す概略図である。このアルゴリズム１は、疑似的なコードによって記述されている。アルゴリズム１は、Ｍ≧３の場合のＳＶＴ算出法である。以下、この図に沿ってアルゴリズムを説明する。
本アルゴリズムにおいて、入力は、Ｍ×２（Ｍ行２列）の行列Ｙ、およびＳＶＴの閾値μである。Ｙの１列目、２列目の列ベクトル（Ｍ次元）をそれぞれｙ_１，ｙ_２（ｙは太字。以下においても同様。）と表す。また、μ＞０である。また、出力は、閾値μにより算出されたＳＶＴである。出力される行列をＺ（太字）と表す。

以下では、図の左側に付した行番号を参照しながら説明する。
第１行のｉｆ文の条件節において、行列Ｙが、Ｍ×２のサイズのゼロ行列であるか否かを判定する。これは、Ｙの階数が０であるか否かの判定と等価である。
第２行は、第１行の条件が真である場合に実行される処理であり、Ｍ×２のサイズのゼロ行列をＺに代入する。この場合、処理を終了してＺを出力する。

第３行は、第１行の条件が偽である場合に対応するｅｌｓｅ節であることを表す。このｅｌｓｅ節は第１７行まで続く。
第４行から第６行までは、第３行からのｅｌｓｅ節の一部である。
第４行において、変数ａ，ｂ，ｃへの代入が行われる。変数ａには、ｙ_１ ^Ｔｙ_１を代入する。変数ｂには、ｙ_１ ^Ｔｙ_２を代入する。変数ｃには、ｙ_２ ^Ｔｙ_２を代入する。言うまでもなく、変数ａ，ｂ，ｃにはそれぞれスカラー値が代入される。
第５行において、変数ｄに（ａｃ−ｂ^２）の値を代入する。また、変数ｅに、ｄの値の平方根を代入する。第５行の「％」以後はコメントであり、実行コードではない。以後においても同様である。第５行のコメントで示すように、変数ｄに代入された値は、行列Ｙ^ＴＹの行列式の値である。
第６行において、変数ｆに、ａ＋ｃの値を代入する。

第７行のｉｆ文の条件節において、ｄ＝０であるか否かを判定する。これは、Ｙの階数が１であるか否かの判定と等価である。
第８行は、第７行の条件が真である場合に実行される処理であり、（１−μ／ＳＱＲＴ（ｆ））_＋をＹに乗じた行列をＺに代入する。この場合、処理を終了してＺを出力する。なお、ここで、ＳＱＲＴ（）は、引数の平方根を返す関数である。

第９行は、第７行の条件が偽である場合に対応するｅｌｓｅ節であることを表す。このｅｌｓｅ節は第１７行まで続く。
第１０行において、変数ｗ_２に、ＳＱＲＴ（ｆ−２ｅ）の値を代入する。コメントに示すように、変数ｗ_２には、特異値σ_１とσ_２の差が代入される。

第１１行のｉｆ文の条件節において、ｗ_２の値が０であるか否かを判定する。これは、特異値σ_１とσ_２とが等しいか否かの判定と等価である。
第１２行は、第１１行の条件が真である場合に実行される処理であり、（１−（（ＳＱＲＴ（２）・μ）／ＳＱＲＴ（ｆ）））_＋をＹに乗じた行列をＺに代入する。この場合、処理を終了してＺを出力する。

第１３行は、第１１行の条件が偽である場合に対応するｅｌｓｅ節であることを表す。このｅｌｓｅ節は第１７行まで続く。
第１４行において、変数ｗ_１に、ＳＱＲＴ（ｆ＋２ｅ）の値を代入する。コメントに示すように、変数ｗ_１には、特異値σ_１とσ_２の和が代入される。
第１５行において、変数σ_２に、（ｗ_１−ｗ_２）／２の値を代入する。
第１６行において、変数γに、（１−（（μ−σ_２）_＋／ｗ_２））_＋の値を代入する。また、変数δに、ｍｉｎ（μ，σ_２）／ｗ_１の値を代入する。前述の通り、γは振幅パラメータであり、δは内分比パラメータである。
そして、第１７行において、下の式で示される行列をＺに代入する。そして、処理を終了してＺを出力する。

以上、説明したように、ＳＶＤを用いず、少ない計算量でＳＶＴを算出することができる。

図５は、アルゴリズム２を示す概略図である。このアルゴリズム２も、疑似的なコードによって記述されている。アルゴリズム２は、Ｍ＝２の場合のＳＶＴ算出法である。以下、この図に沿ってアルゴリズムを説明する。
本アルゴリズムにおいて、入力は、２×２（２行２列）の行列Ｙ、およびＳＶＴの閾値μである。Ｙの１列目、２列目の列ベクトル（２次元）をそれぞれｙ_１，ｙ_２と表す。また、行列Ｙの各要素を、行番号および列番号をこの順で並べたサフィックス（添え字）を用いて、ｙ_１１，ｙ_１２，ｙ_２１，ｙ_２２で表す。また、μ＞０である。また、出力は、閾値μにより算出されたＳＶＴである。出力される行列をＺ（太字）と表す。

以下では、図の左側に付した行番号を参照しながら説明する。
第１行のｉｆ文の条件節において、行列Ｙが、２×２のサイズのゼロ行列であるか否かを判定する。これは、Ｙの階数が０であるか否かの判定と等価である。
第２行は、第１行の条件が真である場合に実行される処理であり、２×２のサイズのゼロ行列をＺに代入する。この場合、処理を終了してＺを出力する。

第３行は、第１行の条件が偽である場合に対応するｅｌｓｅ節であることを表す。このｅｌｓｅ節は第１７行まで続く。
第４行から第６行までは、第３行からのｅｌｓｅ節の一部である。
第４行において、変数ａ，ｃへの代入が行われる。変数ａには、ｙ_１ ^Ｔｙ_１を代入する。変数ｃには、ｙ_２ ^Ｔｙ_２を代入する。
第５行において、変数ｄにＹの行列式の値を代入する。また、変数ｅに、ｄの値の絶対値を代入する。
第６行において、変数ｆに、ａ＋ｃの値を代入する。

第７行のｉｆ文の条件節において、ｄ＝０であるか否かを判定する。これは、Ｙの階数が１であるか否かの判定と等価である。
第８行は、第７行の条件が真である場合に実行される処理であり、（１−μ／ＳＱＲＴ（ｆ））_＋をＹに乗じた行列をＺに代入する。この場合、処理を終了してＺを出力する。

第９行は、第７行の条件が偽である場合に対応するｅｌｓｅ節であることを表す。このｅｌｓｅ節は第１７行まで続く。
第１０行において、変数ｗ_２に、ＳＱＲＴ（ｆ−２ｅ）の値を代入する。コメントに示すように、変数ｗ_２には、特異値σ_１とσ_２の差が代入される。

第１１行のｉｆ文の条件節において、ｗ_２の値が０であるか否かを判定する。これは、特異値σ_１とσ_２とが等しいか否かの判定と等価である。
第１２行は、第１１行の条件が真である場合に実行される処理であり、（１−（（ＳＱＲＴ（２）・μ）／ＳＱＲＴ（ｆ）））_＋をＹに乗じた行列をＺに代入する。この場合、処理を終了してＺを出力する。

第１３行は、第１１行の条件が偽である場合に対応するｅｌｓｅ節であることを表す。このｅｌｓｅ節は第１７行まで続く。
第１４行において、変数ｗ_１に、ＳＱＲＴ（ｆ＋２ｅ）の値を代入する。コメントに示すように、変数ｗ_１には、特異値σ_１とσ_２の和が代入される。
第１５行において、変数σ_２に、（ｗ_１−ｗ_２）／２の値を代入する。
第１６行において、変数γに、（１−（（μ−σ_２）_＋／ｗ_２））_＋の値を代入する。また、変数δに、ｍｉｎ（μ，σ_２）／ｗ_１の値を代入する。
そして、第１７行において、下の式で示される行列をＺに代入する。そして、処理を終了してＺを出力する。

ここで本実施形態のアルゴリズムによるＳＶＴ算出のための計算量について考察する。
式（４）および式（５）の通りＳＶＤを用いるＳＶＴ算出法（従来技術）では、ＳＶＤを求めるために２４Ｍ＋１６０回、閾値処理に２回、行列積を求めるために６Ｍ＋４回の浮動小数点演算が必要である。即ち合計で、３０Ｍ＋１６６回の浮動小数点演算が必要である。
本実施形態のアルゴリズム１による方法では１２Ｍ＋２６回、アルゴリズム２による方法では３６回の浮動小数点演算でＳＶＴを算出することができる。つまり、本実施形態による計算量は、従来技術の１６％〜４０％程度である。即ち、本実施形態により計算量を従来技術よりも６０％〜８４％削減することができる。

［１．３ 行列間並列化］
次に、行列間での処理を並列化して実施する方法について説明する。つまり、上記のアルゴリズムを用いながら、Ｌ個の行列入力Ｙ_１，Ｙ_２，・・・，Ｙ_ＬのＳＶＴを同時に求める方法について説明する。アルゴリズム１，２の処理の大半は，ベクトルや行列に関する基本演算（和、積、定数倍）より構成され，並列化の効果が高い。そこでＬ個の行列Ｙ_１，Ｙ_２，・・・，Ｙ_Ｌに対してデータ並列にＳＶＴを算出するため、アルゴリズム３を用いる。なお、このアルゴリズム３を、Fast Multiple SVT（ＦＭＳＶＴ）と呼ぶ。
ＦＭＳＶＴの計算量は，単純にアルゴリズム１，２の計算量のＬ倍である。

図６は、アルゴリズム３を示す概略図である。このアルゴリズム３も、疑似的なコードによって記述されている。アルゴリズム３は、ＦＭＳＶＴを用いている。以下、この図に沿ってアルゴリズムを説明する。
本アルゴリズムにおいて、入力は、Ｌ個の行列Ｙ_１，Ｙ_２，・・・，Ｙ_Ｌ、およびＳＶＴの閾値μである。また、出力は、Ｌ個の行列Ｚ_１，Ｚ_２，・・・，Ｚ_Ｌである。Ｚ_１，Ｚ_２，・・・，Ｚ_Ｌは、それぞれ、入力行列Ｙ_１，Ｙ_２，・・・，Ｙ_ＬのＳＶＴである。

以下では、図の左側に付した行番号を参照しながら説明する。
第１行において、各ｉ（ｉ＝１，２，・・・，Ｌ）について、アルゴリズム１あるいは２にしたがって、変数値ａ^（ｉ），ｂ^（ｉ），ｃ^（ｉ），ｄ^（ｉ），ｅ^（ｉ），ｆ^（ｉ），ｗ_１ ^（ｉ），ｗ_２ ^（ｉ），σ_２ ^（ｉ）を算出する。
第２行において、各ｉ（ｉ＝１，２，・・・，Ｌ）について、第３行から第１０行までに示す処理を実行する。つまり、入力される行列Ｙ_ｉに基づく条件により分岐し、分岐先において変数γ^（ｉ）およびδ^（ｉ）に値を代入する。

第３行のｉｆ文の条件節において、Ｙ_ｉがゼロ行列であるか否かを判定する。つまり、Ｙ_ｉの階数が０であるか否かを判定する。
第４行は、第３行の条件が真である場合に実行される処理であり、γ^（ｉ）およびδ^（ｉ）にそれぞれ０を代入する。そして、γ^（ｉ）およびδ^（ｉ）への代入後に、第１１行の処理に移る。

第５行のｅｌｓｅ ｉｆは、第３行の条件が偽である場合に、別の条件判定を行うためのものである。第５行のｉｆ文の条件節において、ｄ^（ｉ）が０であるか否かを判定する。つまり、行列Ｙ_ｉの階数が１であるか否かを判定する。
第６行は、第５行の条件が真である場合に実行される処理であり、γ^（ｉ）に、（１−（μ／ＳＱＲＴ（ｆ^（ｉ））））_＋の値を代入する。また、δ^（ｉ）に０を代入する。そして、γ^（ｉ）およびδ^（ｉ）への代入後に、第１１行の処理に移る。

第７行のｅｌｓｅ ｉｆは、第５行の条件が偽である場合に、別の条件判定を行うためのものである。第７行のｉｆ文の条件節において、ｗ_２ ^（ｉ）が０であるか否かを判定する。
第８行は、第７行の条件が真である場合に実行される処理であり、γ^（ｉ）に、（１−（（ＳＱＲＴ（２）・μ）／ＳＱＲＴ（ｆ^（ｉ））））_＋の値を代入する。また、δ^（ｉ）に０を代入する。そして、γ^（ｉ）およびδ^（ｉ）への代入後に、第１１行の処理に移る。

第９行のｅｌｓｅは、第７行の条件が偽である場合に対応する。この場合の処理は、第１０行に記述されている。
そして、第１０行において、て、γ^（ｉ）およびδ^（ｉ）にそれぞれ値を代入する。γ^（ｉ）には、（１−（μ−σ_２ ^（ｉ））_＋／ｗ_２ ^（ｉ））_＋の値を代入する。そして、δ^（ｉ）には、ｍｉｎ（μ，σ_２ ^（ｉ））／ｗ_１ ^（ｉ）の値を代入する。そして、第１１行の処理に移る。

条件分岐の結果に応じてγ^（ｉ）およびδ^（ｉ）の値が設定された状態で、第１１行の処理を実行する。
第１１行においては、各ｉ（ｉ＝１，２，・・・，Ｌ）について、アルゴリズム１（Ｍ≧３の場合）あるいはアルゴリズム２（Ｍ＝２の場合）の第１７行の処理を実行する。その処理により求められたＺを、Ｚ_iとする。
そして、処理を終了し、行列Ｚ_１，Ｚ_２，・・・，Ｚ_Ｌを出力する。

アルゴリズム３では、上述したように、係数γ^（ｉ），δ^（ｉ）の参照のみを条件分岐させ、その他の処理をインデクスｉ＝１，２，・・・，Ｌ間で共通化する。また、異なるｉの間でデータが相互に干渉しあわない。これにより、係数γ^（ｉ），δ^（ｉ）を算出するために費やす４Ｌ回の浮動小数点演算以外の処理を、すべて並列化できる。具体的には、Ｍ≧３では９３．５％以上、Ｍ＝２では８８．９％の浮動小数点演算をデータ並列に処理することができる。

また、アルゴリズム３はＳＩＭＤ型アーキテクチャとの親和性が特に高い。ＳＩＭＤは、「single instruction multiple data」の略であり、単一の命令列を複数のデータに適用する処理形態である。ＳＩＭＤでは条件分岐をマスク演算で実現でき、係数γ^（ｉ），δ^（ｉ）の算出もデータ並列化することが可能である。

なお、処理を並列化する手法自体には様々なものがあり、アルゴリズム３を実施するために適宜並列化手法を選択して用いればよい。

［２． 第１実施形態］
次に、第１実施形態による装置構成等について説明する。
図７は、本実施形態による行列単純化装置の概略機能構成を示すブロック図である。図示するように、行列単純化装置１は、入力部１１と、ベクトル化部１２と、低ランク近似化部１５と、出力部１６とを含んで構成される。行列単純化装置１が有する各部の機能は、例えば、電子回路を用いて構成される。また、各部は、必要に応じてデータを記憶するための記憶部を内部に備える。この記憶部は、半導体メモリーや磁気ハードディスク装置などといった記憶手段を用いて実現される。また、各部の機能は、コンピューターとプログラムとで実現されてもよい。

行列単純化装置１は、Ｍ行２列（ただし、Ｍ≧２）または２行Ｍ列の行列を入力し、その行列を低ランク近似化し、低ランク近似化した行列を出力する装置である。

入力部１１は、Ｍ行２列または２行Ｍ列の行列のデータを外部から取得する。なお、Ｍは、２以上の整数である。この行列の要素は、数値（スカラー値）である。入力部１１は、必要に応じて、入力した行列を転置する。つまり、後段のベクトル化部１２および低ランク近似化部１５がＭ行２列または２行Ｍ列のいずれか一方の形式の行列のみを処理するように構成されているときであって、入力された行列がその形式に合わないとき（つまり、縦と横が逆）に、入力部１１は、入力された行列を転置する。これにより、行列単純化装置１は、Ｍ行２列または２行Ｍ列のいずれの行列をも処理することができるようになる。

以下において、ベクトル化部１２と低ランク近似化部１５とは、Ｍ行２列の行列を処理するものとして説明する。但し、これが２行Ｍ列の行列を処理するものであってもよく、本質的な処理内容は変わらない。
ベクトル化部１２は、Ｍ行２列の行列（Ｙとする）をベクトル化する。ここでのベクトル化とは、Ｍ行２列の行列Ｙを、２個のＭ次元の列ベクトルｙ_１（第１ベクトル）とｙ_２（第２ベクトル）とに分割して出力する処理である。つまり、Ｙ＝［ｙ_１，ｙ_２］である。ベクトル化部１２は、これらのベクトルｙ_１，ｙ_２を、低ランク近似化部１５に渡す。
つまり、ベクトル化部１２は、Ｍ行２列の入力行列から各列に対応するＭ次元の第１ベクトルおよび第２ベクトルを抽出し、または２行Ｍ列の入力行列から各行に対応するＭ次元の第１ベクトルおよび第２ベクトルを抽出する。

低ランク近似化部１５は、ベクトル化部１２から渡されたベクトルｙ_１，ｙ_２に基づき、行列Ｙを低ランク近似化する。言い換えれば、低ランク近似化部１５は、行列Ｙを単純化する。具体的には、低ランク近似化部１５は、Ｍの値に応じて、前述のアルゴリズム１または２のいずれかを用いて、行列Ｙの低ランク近似化を行う。具体的には、Ｍ≧３の場合には、低ランク近似化部１５は、アルゴリズム１を用いる。また、Ｍ＝２の場合には、低ランク近似化部１５は、アルゴリズム２を用いる。これにより、低ランク近似化部１５は、低ランク化された行列Ｚを出力する。
つまり、低ランク近似化部１５は、ｙ_１（第１ベクトル）およびｙ_２（第２ベクトル）を含む平面（像ＩｍＹ）内でｙ_２を正方向に所定角（π／２）回転させて得られるＲｙ_２（第３ベクトル）と第１ベクトルとの距離である第１距離と、平面内で前記第２ベクトルを負方向に前記所定角回転させて得られる−Ｒｙ_２（第４ベクトル）と前記第１ベクトルとの距離である第２距離とを求め、第１距離および第２距離に基づいて前記入力行列をより低ランクで近似した低ランク近似行列Ｚを求める。ここで、角の正方向とは、ベクトルｙ_１からベクトルｙ_２まで最短で辿り着ける回転の方向である。
また、低ランク近似化部１５は、入力行列ＹがＭ行２列であるときにはＹ＝［ｙ_１，ｙ_２］とし、入力行列Ｙが２行Ｍ列であるときにはＹ^Ｔ＝［ｙ_１，ｙ_２］とし、次の方法により低ランク近似行列Ｚを求める。即ち、低ランク近似化部１５は、所定角回転させる回転行列をＲとして、［ｙ_２，−ｙ_１］＝Ｙ（Ｙの上にバー）として、内分比パラメータδ（０≦δ≦１）を用いて、前記入力行列Ｙと行列ＲＹ（Ｙの上にバー）との内分比パラメータδによる内分である内分行列に基づいて低ランク近似行列Ｚを求める。この「内分行列」を式で表すと、（１−δ）Ｙ＋δＲＹ（Ｙの上にバー）である。そして、さらに、振幅パラメータγを乗じて、Ｚ＝γ（１−δ）Ｙ＋γδＲＹ（Ｙの上にバー）が得られる。つまり、低ランク近似化部１５は、さらに、振幅パラメータγにも基づいて低ランク近似行列Ｚを求める。
また、低ランク近似化部１５は、所定角回転させる回転行列Ｒを前に示した式（１４）として第２ベクトルを回転させるものである。
また、低ランク近似化部１５は、入力行列Ｙの特異値をσ_１，σ_２（σ_１≧σ_２≧０）としたとき、行列Ｙの階数が２であり且つσ_１≠σ_２の場合には、式（１２）および式（１３）により、式（１２）の右辺を求めて低ランク近似行列Ｚとするものである。且つ、低ランク近似化部１５は、式（１５）の右辺の計算を行うことによって、式（１２）内のＲＹ（Ｙの上にバー）を求める。なお、式（１２）、式（１３）、式（１５）に関しては、前に示した通りである。
また、低ランク近似化部１５は、「入力行列Ｙの階数が２且つσ_１≠σ_２」以外の場合には、式（１９）の右辺を求めて低ランク近似行列とする。式（１９）自体は、前に示した通りである。

出力部１６は、低ランク近似化部１５によって求められた行列Ｚのデータを外部に出力する。
つまり、行列単純化装置１は、入力される行列Ｙを低ランク単純化し、その結果である行列Ｚを出力する。
なお、入力部１１が行列の転置を行った場合、出力部１６は、低ランク近似化部１５によって得られた行列Ｚを再び転置してから出力する。これにより、入力される行列のサイズ（行および列の数）と、出力される行列のサイズとを合わせることができる。

次に、行列単純化装置１が扱う行列およびベクトルのデータの構造について説明する。
図８は、行列単純化装置１が処理する主要なデータの構成を示す概略図である。同図（ａ）は、入力行列Ｙのデータ構成を示す。ここでは、Ｍ行２列の場合の行列を示している。この図では行番号と列番号を付して示しており、行列Ｙの要素であるｙ_ｉｊ（ｉ＝１，・・・，Ｍ、ｊ＝１，２）が各領域に格納される。なお、２行Ｍ列の行列の場合には、行と列の方向が入れ替わる。同図（ｂ）は、行列Ｙを基にベクトル化されたベクトルｙ_１およびｙ_２のデータ構成を示す。ベクトルｙ_１およびｙ_２は、それぞれ、Ｍ次元の列ベクトルである。行列Ｙ（同図（ａ））の第１列がベクトルｙ_１に対応し、第２列がベクトルｙ_２に対応する。同図（ｃ）は、出力行列Ｚのデータ構成を示す。行列Ｚについても、Ｍ行２列の場合の行列を示している。行列Ｚの要素であるz_ｉｊ（ｉ＝１，・・・，Ｍ、ｊ＝１，２）が各領域に格納される。

［３． 第２実施形態］
次に、第２実施形態による装置構成等について説明する。なお、前実施形態において既に説明した事項については説明を省略する場合がある。ここでは、本実施形態に特有の事項を中心に説明する。
図９は、本実施形態による行列単純化装置の概略機能構成を示すブロック図である。図示するように、行列単純化装置２は、入力部２１と、ベクトル化部２２−１，２２−２，・・・，２２−Ｌと、低ランク近似化部２５−１，２５−２，・・・，２５−Ｌと、出力部２６とを含んで構成される。

行列単純化装置２は、第１実施形態の行列単純化装置１と同様の処理によって、入力行列Ｙを低ランク化して出力するものである。ただし、行列単純化装置２は並列処理のための機構を有しており、Ｌ個（Ｌは自然数）の入力行列を同時に処理することができる。行列単純化装置２は、前述のアルゴリズム３を用いるものであり、例えばＳＩＭＤによる処理を行う。つまり、行列単純化装置２は、アルゴリズム３を実現する命令列を、Ｌ個のデータ列（入力行列Ｙ_１，・・・，Ｙ_Ｌ）に適用する形態をとる。言い換えれば、行列単純化装置２は、単一命令列を、入力行列にそれぞれ対応する複数のデータに適用して並列処理を行うものである。

入力部２１は、Ｍ行２列または２行Ｍ列の、Ｌ個の行列Ｙ_１，・・・，Ｙ_Ｌのデータを外部から取得する。また、第１実施形態における入力部１１と同様に、入力部２１は、必要に応じて入力された各行列の転置を行う。そして、入力部２１は、行列Ｙ_１，・・・，Ｙ_Ｌのデータを、それぞれ、ベクトル化部２２−１，２２−２，・・・，２２−Ｌに渡す。

ベクトル化部２２−１，２２−２，・・・，２２−Ｌは、入力部２１から、それぞれ、行列Ｙ_１，Ｙ_２，・・・，Ｙ_Ｌのデータを受け取る。そして、ベクトル化部２２−１，２２−２，・・・，２２−Ｌは、その行列をベクトル化する。つまり、ベクトル化部２２−ｉ（ｉ＝１，・・・，Ｌ）は、行列Ｙ_ｉのデータをベクトル化し、ベクトルｙ₁ ^（ｉ），ｙ_２ ^（ｉ）を出力する。ただし、ベクトルｙ_ｊ ^（ｉ）は、行列Ｙ_ｉの第ｊ番目の列の要素で成る列ベクトルである。そして、各々のベクトル化部２２−ｉは、ベクトルｙ₁ ^（ｉ），ｙ_２ ^（ｉ）のデータを、低ランク近似化部２５−ｉに渡す。
つまり、ベクトル化部２２−１，２２−２，・・・，２２−Ｌは、複数の入力行列を基にそれぞれの入力行列についてベクトルｙ_１およびベクトルｙ_２を抽出するものである。

低ランク近似化部２５−１，２５−２，・・・，２５−Ｌは、それぞれ、前述のアルゴリズム３（図６）を用いて、入力行列に対するＳＶＴの処理を実行する。なお、アルゴリズム３の第１１行においてアルゴリズム１，２を参照しているが、行列のサイズＭに応じて、アルゴリズム１（Ｍ≧３の場合）またはアルゴリズム２（Ｍ＝２の場合）のいずれを用いるかを決定する。つまり、低ランク近似化部２５−ｉ（ｉ＝１，・・・，Ｌ）は、行列Ｙ_ｉのＳＶＴを実行し、出力行列Ｚ_ｉを算出する。そして、低ランク近似化部２５−ｉは、得られた行列Ｚ_ｉを出力部２６に渡す。
つまり、低ランク近似化部２５−１，２５−２，・・・，２５−Ｌは、各入力行列から抽出されたｙ_１（第１ベクトル）およびｙ_２（第２ベクトル）を用いて低ランク近似行列Ｚを求めるものである。

出力部２６は、低ランク近似化部２５−１，２５−２，・・・，２５−Ｌから、それぞれ、出力行列Ｚ_１，Ｚ_２，・・・，Ｚ_Ｌを受け取り、外部に出力する。

本実施形態によれば、複数の行列について、同時並列的に高速に低ランク近似化を行うことが可能となる。

［４． 第３実施形態］
次に、第３実施形態による装置構成等について説明する。なお、前実施形態までにおいて既に説明した事項については説明を省略する場合がある。ここでは、本実施形態に特有の事項を中心に説明する。
図１０は、本実施形態による行列単純化装置の概略機能構成を示すブロック図である。図示するように、行列単純化装置３は、画像入力部３１と、ベクトル化部１２と、低ランク近似化部１５と、画像出力部３６とを含んで構成される。

行列単純化装置３は、第１実施形態の行列単純化装置１と同様の処理によって、入力行列Ｙを低ランク化して出力するものである。ただし、行列単純化装置３が入力する行列は画像のデータであり、行列単純化装置３はその画像のデータを単純化する処理を行う。
画像入力部３１は、画像データを外部から取得する。なお、画像入力部３１は、Ｍ行２列または２行Ｍ列の行列の形式の画像データを取得する。つまり、例えば、画像入力部３１は、Ｍ行２列または２行Ｍ列に配列された２Ｍ個の画素から成る画像を入力する。各画素値はスカラーである。画像入力部３１は、取得した画像を入力画像Ｙとしてベクトル化部１２に渡す。
ベクトル化部１２は、第１実施形態の場合と同様の機能を有する。
低ランク近似化部１５は、第１実施形態の場合と同様の機能を有する。低ランク近似化部１５は、低ランク化された行列Ｚを画像出力部３６に渡す。
画像出力部３６は、入力された画像に対応するサイズの画像を出力する。出力画像は、入力画像を単純化した（低ランク化した）ものである。

本実施形態によれば、行列単純化装置３は、Ｍ×２または２×Ｍのサイズを有する画像データを単純化することができる。

上述した各実施形態における行列単純化装置の機能の全部または一部を、コンピューターで実現するようにしても良い。その場合、この機能を実現するためのプログラムをコンピューター読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピューターシステムに読み込ませ、実行することによって実現しても良い。なお、ここでいう「コンピューターシステム」とは、ＯＳや周辺機器等のハードウェアを含むものとする。また、「コンピューター読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ−ＲＯＭ、ＤＶＤ−ＲＯＭ、ＵＳＢメモリー等の可搬媒体、コンピューターシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピューター読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間の間、動的にプログラムを保持するもの、その場合のサーバーやクライアントとなるコンピューターシステム内部の揮発性メモリーのように、一定時間プログラムを保持しているものも含んでも良い。また上記プログラムは、前述した機能の一部を実現するためのものであっても良く、さらに前述した機能をコンピューターシステムにすでに記録されているプログラムとの組み合わせで実現できるものであっても良い。

以上、この発明の実施形態について図面を参照して詳述してきたが、具体的な構成はこの実施形態に限られるものではなく、この発明の要旨を逸脱しない範囲の設計等も含まれる。

以上、説明した複数の実施形態のいずれかによる行列単純化装置は、特異値閾値処理（ＳＶＴ）の高速算出アルゴリズムを実行可能である。これにより、多数の小型行列を低ランク正則化する最適化問題の解を高速に求めることが可能となる。
具体的には、Ｍ行２列または２行Ｍ列（Ｍ≧２）の行列を対象として、低ランク化処理を、少ない計算量で高速に実行することが可能となる。また、複数の行列を対象として、低ランク化処理を並列して実行することができる。つまり、Ｍ行２列または２行Ｍ列の行列に対するＦＭＳＶＴアルゴリズムが実現される。
さらに具体的には、核型ノルムを部分空間上のベクトル距離で表現できる発見に基づき、特異値分解（ＳＶＤ）が不要でデータ並列なＳＶＴ算出法を実現した。また、実際のデータを用いた評価実験の結果、従来手法以上の計算精度を持ちつつ、最大９５．８２倍高速にＳＶＴを算出できることを確認した。

以下において、命題等の証明を記載する。
［付録Ａ 命題１の証明］
行列Ｙ^ＴＹの固有値をλ_１，λ_２（λ_１≧λ_２≧０）とする。

λ_１，λ_２に関する上の式より、下の式（２０）が成立する。

ここでＹを表現行列とする線形写像の像ＩｍＹを考えると、高々２次元の部分空間である（∵ｄｉｍ（ＩｍＹ）＝ｒａｎｋＹ＝２）。そこで図１のようにＩｍＹを示すと、
ベクトルｙ_１，ｙ_２∈Ｒ^Ｍの幾何的性質に着目する。ｙ_１，ｙ_２の成す角をθ∈［０，π］［ｒａｄ］とする。

すると、上の式の通りであるから、下の式（２１）が得られる。

また図１に示したように、ベクトルｙ_２をπ／２［ｒａｄ］回転させたＲｙ_２を考えると（Ｒは太字であり、前述の通り回転行列）、ｙ_１とＲｙ_２の成す角はθ＋π／２［ｒａｄ］である。

上の式の通りであるから、下の式（２２）が成立する。

これを式（２１）に代入することにより、下の式（２３）が得られる。つまり、定理１が証明された。

［付録Ｂ 命題３の証明］
命題３の証明に先立って、次の２つの補題Ｂ１およびＢ２を示す。

＜補題Ｂ１＞：式（８）の行列Ｂに関して、任意のｘ＝［ｘ_１ ^Ｔ，ｘ_２ ^Ｔ］^Ｔ，ｘ_１，ｘ_２∈ＩｍＹに対し、下の式（２４）が成立する。

補題Ｂ１の証明は、下記の通り。
前述の Ｒ^ＴＲｙ＝ＲＲ^Ｔｙ＝ｙ の性質より、下の式（２５）および式（２６）の通りである。つまり、補題Ｂ１が証明された。

＜補題Ｂ２＞： ベクトルｚ_１，ｚ_２∈Ｒ^Ｍを、下の式の通りとする。

このとき、ｚ_１，ｚ_２∈ＩｍＹである。

補題Ｂ１の証明は、下記の通り。
（Ｕ，Σ，Ｖ）＝ＳＶＤ（Ｙ）とし、また下の式の通りとする。

ＳＶＤの定義より、ｙ_ｉ＝σ_１ｖ_ｉ１ｕ_１＋σ_２ｖ_ｉ２ｕ_２なので、下の式（２７）の通りである。

同様にＳＶＴの定義（式（４））より、ｚ_ｉ＝（σ_１−μ）_＋・ｖ_ｉ１ｕ_１＋（σ_２−μ）_＋・ｖ_ｉ２ｕ_２ なので、式（２７）における場合分け毎に、ｚ_ｉ∈ＩｍＹが確認できる。つまり、補題Ｂ２が証明された。

以上の補題Ｂ１，Ｂ２を踏まえて、命題３を以下の通り証明する。
系２と、近接写像の定義より、下の式（２８）が成立する。

ここで、ｚ＝ｖｅｃ（Ｚ）と変数置換しているが、補題Ｂ２より探索範囲をｚ∈ＩｍＹに限定できる。そこで変数を再度置換して、下の式（２９）の通りとする。

両辺にＢ^Ｔを掛けると、補題Ｂ１より、Ｂ^Ｔｘ＝ｚとなる。また、再び補題Ｂ１より、ｖｅｃ（Ｙ）＝（１／２）・Ｂ^ＴＢｖｅｃ（Ｙ）が成立する。以上より、下の式（３０）の通りとなる。なお、式（３０）の最右辺への変形にも補題Ｂ１を使用している。

式（２８）に式（２９），式（３０）を適用すると、式（３１）の通りとなる。つまり、命題３が証明された。

［付録Ｃ 命題４の証明］
命題４の証明を下に記す。
Ｌ_∞，２ノルムの双対ノルムはＬ_１，２ノルムなので、Moreauの直交分解より、下の式（３２）の通りである。

なお、Moreauの直交分解についてここで補足する。ある凸関数ｆ（ｘ）とそのルジャンドル変換ｆ^＊（ｘ）について、ｘの分解ｘ＝ｐｒｏｘ_ｆ（ｘ）＋ｐｒｏｘ_ｆ＊（ｘ）をMoreauの直交分解と言う。特に、凸関数がノルムであるとき、双対ノルムの球Ｂ_ｄへの射影を用いて、ｘに関して、下記の通りである。
即ち、ノルムと双対ノルムと射影は、次の通り表される。

そして、ｘに関して、下の式が成立する。

式（３２）に戻る。ここで、式（３２）の右辺のｐｒｏｊ_{μＢ１，２}（ｘ）は、半径μのＬ_１，２球であるμＢ_１，２へのＬ_２距離射影で、下の式（３３）で表す最適化問題の解である。

Ewout van den Bergらによれば、式（３３）で表される問題はＬ_１球射影とＬ_２球射影の問題に分離でき、下の式（３４），式（３５）で表される２つの問題と等価である。

式（３４）のＬ_２ノルム射影は、下の式の通り計算できる。

一方で式（３５）のＬ_１ノルム射影は、下の式（３６）の通りとなる。

以上より、式（３７）、式（３８）が成立する。つまり、命題４が証明された。

［付録Ｄ 定理６の証明］
行列Ｚを行列ＹのＳＶＴとする。即ち、下の式の通りとする。

ｒａｎｋＹ＝０の場合、Ｙ＝Ｏ_Ｍ×２であるため、自明である。
ｒａｎｋＹ＝１の場合、σ_１＞０，σ_２＝０であるため、下の式（３９）の通りである。

ここで、下の式（４０）の通りであるので、階数１の場合についても証明された。

ｒａｎｋＹ＝２、且つσ_１＝σ_２の場合、下の式（４１）の通りである。

ここで、下の式（４２）の通りであるので、階数２の場合についても証明された。
つまり、定理６が証明された。

１，２，３ 行列単純化装置
１１ 入力部
１２ ベクトル化部
１５ 低ランク近似化部
１６ 出力部
２１ 入力部
２２−１，２２−２，・・・,２２−Ｌ ベクトル化部
２５−１，２５−２，・・・,２５−Ｌ 低ランク近似化部
２６ 出力部
３１ 画像入力部
３６ 画像出力部

Claims (8)

1. Ｍ行２列の入力行列から各列に対応するＭ次元の第１ベクトルおよび第２ベクトルを抽出し、または２行Ｍ列の入力行列から各行に対応するＭ次元の第１ベクトルおよび第２ベクトルを抽出する（ただし、Ｍは２以上の整数）ベクトル化部と、
   前記第１ベクトルおよび前記第２ベクトルを含む平面内で前記第２ベクトルを第１の方向に所定角回転させて得られる第３ベクトルと前記第１ベクトルとの距離である第１距離と、前記平面内で前記第２ベクトルを第２の方向に前記所定角回転させて得られる第４ベクトルと前記第１ベクトルとの距離である第２距離とを求め、前記第１距離および前記第２距離に基づいて前記入力行列をより低ランクで近似した低ランク近似行列を求める低ランク近似化部と、
   を具備する行列単純化装置。
2. 前記入力行列がＭ行２列の行列ＹであるときにはＹ＝［ｙ_１，ｙ_２］とし、前記入力行列が２行Ｍ列の行列ＹであるときにはＹ^Ｔ＝［ｙ_１，ｙ_２］とし、前記第１ベクトルをｙ_１とし、前記第２ベクトルをｙ_２としたとき、
   前記低ランク近似化部は、前記所定角回転させる回転行列をＲとして、［ｙ_２，−ｙ_１］＝Ｙ（Ｙの上にバー）として、内分比パラメータδ（０≦δ≦１）を用いて、前記入力行列Ｙと行列ＲＹ（Ｙの上にバー）との前記内分比パラメータδによる内分である内分行列と、振幅パラメータγと、に基づいて前記低ランク近似行列を求める、
   請求項１に記載の行列単純化装置。
3. 前記入力行列がＭ行２列の行列ＹであるときにはＹ＝［ｙ_１，ｙ_２］とし、前記入力行列が２行Ｍ列の行列ＹであるときにはＹ^Ｔ＝［ｙ_１，ｙ_２］とし、前記第１ベクトルをｙ_１とし、前記第２ベクトルをｙ_２としたとき、
   前記低ランク近似化部は、前記所定角回転させる回転行列Ｒを式（１４）として前記第２ベクトルを回転させる、
   

   

   （ただし、「ｄｅｔ Ｘ」は行列Ｘの行列式）
   請求項１に記載の行列単純化装置。
4. 前記行列Ｙの特異値をσ_１，σ_２（σ_１≧σ_２≧０）として、閾値をμとしたとき、
   前記低ランク近似化部は、前記行列Ｙの階数が２であり且つσ_１≠σ_２の場合には、式（１２）および式（１３）により、式（１２）の右辺を求めて前記低ランク近似行列とするものであり、
   

   

   且つ、前記低ランク近似化部は、式（１５）の右辺の計算を行うことによって、式（１２）内のＲＹ（Ｙの上にバー）を求める、
   

   

   （ただし、スカラー値の右下に付ける「＋」はランプ関数を表す）
   請求項３に記載の行列単純化装置。
5. 前記低ランク近似化部は、「入力行列Ｙの階数が２且つσ_１≠σ_２」以外の場合には、式（１９）の右辺を求めて前記低ランク近似行列とするものである、
   

   

   請求項４に記載の行列単純化装置。
6. 前記ベクトル化部は複数の前記入力行列を基にそれぞれの前記入力行列の前記第１ベクトルおよび前記第２ベクトルを抽出するものであり、
   前記低ランク近似化部は、各入力行列から抽出された前記第１ベクトルおよび前記第２ベクトルを用いて前記低ランク近似行列を求めるものであり、
   単一命令列を、前記入力行列にそれぞれ対応する複数のデータに適用して並列処理を行うものである、
   請求項１から５までのいずれか一項に記載の行列単純化装置。
7. コンピューターを、
   Ｍ行２列の入力行列から各列に対応するＭ次元の第１ベクトルおよび第２ベクトルを抽出し、または２行Ｍ列の入力行列から各行に対応するＭ次元の第１ベクトルおよび第２ベクトルを抽出する（ただし、Ｍは２以上の整数）ベクトル化部と、
   前記第１ベクトルおよび前記第２ベクトルを含む平面内で前記第２ベクトルを第１の方向に所定角回転させて得られる第３ベクトルと前記第１ベクトルとの距離である第１距離と、前記平面内で前記第２ベクトルを第２の方向に前記所定角回転させて得られる第４ベクトルと前記第１ベクトルとの距離である第２距離とを求め、前記第１距離および前記第２距離に基づいて前記入力行列をより低ランクで近似した低ランク近似行列を求める低ランク近似化部と、
   として機能させるためのプログラム。
8. コンピューターが、Ｍ行２列の入力行列から各列に対応するＭ次元の第１ベクトルおよび第２ベクトルを抽出し、または２行Ｍ列の入力行列から各行に対応するＭ次元の第１ベクトルおよび第２ベクトルを抽出する（ただし、Ｍは２以上の整数）ベクトル化過程と、
   コンピューターが、前記第１ベクトルおよび前記第２ベクトルを含む平面内で前記第２ベクトルを第１の方向に所定角回転させて得られる第３ベクトルと前記第１ベクトルとの距離である第１距離と、前記平面内で前記第２ベクトルを第２の方向に前記所定角回転させて得られる第４ベクトルと前記第１ベクトルとの距離である第２距離とを求め、前記第１距離および前記第２距離に基づいて前記入力行列をより低ランクで近似した低ランク近似行列を求める低ランク近似化過程と、
   を含む行列単純化方法。

Images