(1)3階直交テンソル積展開の性質とその計算法の改良に関する研究

全文

(1)3階直交テンソル積展開の性質とその計算法の改良に関する研究. 2011年. 3月. 博士(理学). 大. 隈. 千春. 熊本大学大学院自然科学研究科.

(2) 3 階直交テンソル積展開の性質とその計算法の改良に関する研究論文要旨多次元データを低次元のデータで近似することは信号処理や，画像処理，音声認識やパターン認識といった多くの分野で利用され，画像データの圧縮，信号処理でのディジタルフィルタ設計などに応用されている．多次元データを低次元分解するための手法として，2 次元データの場合には行列の固有値分解 (Eigen-Value Decomposition – EVD) や特異値分解 (Singular Value Decomposition – SVD)を利用する手法がよく知られている． 3 次元以上のデータに対しては， SVD を多次元に拡張して定義した多次元外積展開が提案されている．その数値計算には非線形最適化法が用いられるが，非線形最適化法は行列のサイズが大きくなると多大な計算時間が必要である．そのため，べき乗法を利用した 3 階テンソル積展開 (3 rd -Order Tensor Product Expansion – 3-TPE)や，さらに展開項間の直交性を満たすようにそれを修正した 3 階直交テンソル積展開 (3 rd -Order Orthogonal Tensor Product Expansion – 3-OTPE)が考案されている． 3 階テンソル積展開は，信号処理の分野では 3 次元ディジタルフィルタの設計問題，画像処理の分野では動画像データ圧縮やカラー画像からの文字データ抽出など幅広い応用が可能であるが，より良い計算精度や処理時間が要求されている．そこで本研究では，3 階テンソル積展開によって得られる展開項の性質を調べることと，それを利用してアルゴリズムの計算精度と計算時間の改良を行うことを目的とした．まず，村上によって提案されている 3-TPE と 3-OTPE の計算アルゴリズムを検証して，計算精度と計算時間の改良を図ることにした．またべき乗法による計算の過程において，反復ベクトルに条件を付加して修正できることを利用して，展開項を構成する展開ベクトルの成分を非負とする分解アルゴリズムを開発し，3 次元ディジタルフィルタの振幅応答特性の例について数値. (i).

(3) 実験を試みた．次に，近年パターン認識などの分野において注目されている Lathauwer らによって提案された高次特異値分解 (Higher-Order Singular Value Decomposition – HOSVD)に着目し，HOSVD と 3-TPE について展開計算法やそれぞれの計算によって得られる展開項の違い，計算精度や計算時間について調べることにした． HOSVD との比較により，これまでの 3-OTPE の計算アルゴリズムには多次元データの各次元のサイズが全て同じではない場合，展開項の打ち切り和が必ずしも元の多次元データとの残差を最小にしないことが分かった．そこでこの問題を解決する新たな計算アルゴリズムを考案し，計算精度の向上を図った．本論文の第 1 章では，研究の背景と目的および概要，さらに論文の構成を述べる．第 2 章では，3 階テンソル積展開の定義，およびべき乗法を用いた 3 階テンソル積展開と 3 階直交テンソル積展開の計算アルゴリズムについて説明した後，3 階直交テンソル積展開の最終項の計算方法の改良について述べる．次いで，3 階テンソル積展開によって得られる各展開項を構成する展開ベクトルの成分を非負とする非負分解のアルゴリズムを述べ，3 次元ディジタルフィルタの振幅応答特性を例に取って行った数値実験の結果を示す．第 3 章では，HOSVD についてその計算法を説明する． HOSVD には最良ランク(R 1, R 2 , … , R N )近似，最良ランク 1 近似と呼ばれる分解があるが，3 階テンソル積展開の手法および計算結果と比較してその相違点を述べる．第 4 章では，初めに 3-OTPE の計算アルゴリズムの問題点として，テンソルの各次元のサイズが全く同じではない場合には残項の計算が必要になってアルゴリズムが複雑になることと，残項の計算に Gram-Schmidt の直交化法を用いているため，初期値の選び方に依存して展開項には一意性がないことを指摘する．これらの問題を解決するために考案した改良法を提示し，従来の 3-OTPE と改良した 3-OTPE のアルゴリズムを用いて数値実験を通して改良法の有効性を示す．第 5 章では，本研究の総括を行い，得られた結論をまとめるとともに，今後に残された課題について述べる． (ii).

(4) 目. 次. 第 1 章序論 ................................................................................................. - 1 1.1. 研究の背景と目的 ............................................................................. - 1 -. 1.2. 本研究の概要 .................................................................................... - 2 -. 1.3. 各章の構成 ........................................................................................ - 3 -. 第 1 章の参考文献 ...................................................................................... - 5 第 2 章 3 階テンソル積展開 ........................................................................ - 7 2.1. 緒言 ................................................................................................... - 7 -. 2.2. テンソルの表記と演算 ..................................................................... - 7 -. 2.3. べき乗法を用いた 3 階テンソル積展開 ......................................... - 10 -. 2.3.1. べき乗法を用いた固有値分解 .................................................. - 10 -. 2.3.2. べき乗法を用いた特異値分解 .................................................. - 12 -. 2.3.3. 3 階テンソル積展開の定義 ....................................................... - 15 -. 2.3.4. べき乗法を用いた 3 階テンソル積展開の計算 ........................ - 16 -. 2.3.5. 3 階テンソル積展開の計算結果 ............................................... - 19 -. 2.4. 3 階直交テンソル積展開 ................................................................. - 25 -. 2.4.1. 3 階直交テンソル積展開の定義 ............................................... - 25 -. 2.4.2. べき乗法を用いた 3 階直交テンソル積展開の計算 ................ - 25 -. 2.4.3. 3 階直交テンソル積展開の計算結果 ........................................ - 28 -. 2.4.4. 3 階直交テンソル積展開の計算アルゴリズムの改良 ............. - 31 -. 2.5. 3 階テンソル積展開の非負分解 ..................................................... - 32 -. 2.5.1. 非負分解の計算 ........................................................................ - 33 -. 2.5.2. 数値実験 .................................................................................... - 34 -. 2.6. 結言 ................................................................................................. - 37 -. 第 2 章の参考文献 .................................................................................... - 39 第 3 章高次特異値分解と 3 階テンソル積展開 ....................................... - 41 3.1. 緒言 ................................................................................................. - 41 -. (iii).

(5) 3.2. 高次特異値分解 ............................................................................... - 42 -. 3.2.2. n -モード展開行列 ................................................................... - 42 n -モード積 ............................................................................... - 43 -. 3.2.3. N 階テンソルの高次特異値分解 ............................................. - 44 -. 3.2.1. 3.3. 最良ランク(R 1 , R 2 , … , R N )近似 ....................................................... - 45 -. 3.4. 最良ランク 1 近似 ........................................................................... - 47 -. 3.5. 高次特異値分解と 3 階テンソル積展開の比較 .............................. - 49 -. 3.5.1. 展開項の比較 ............................................................................ - 49 -. 3.5.2. 計算時間 .................................................................................... - 52 -. 3.6. 高次特異値分解と 3 階直交テンソル積分解の比較 ...................... - 52 -. 3.6.1. 展開項の比較 ............................................................................ - 53 -. 3.6.2. 計算時間 .................................................................................... - 59 -. 3.7. 結言 ................................................................................................. - 60 -. 第 3 章の参考文献 .................................................................................... - 61 第 4 章 3 階直交テンソル積展開の計算法の改良 .................................... - 63 4.1. 緒言 ................................................................................................. - 63 -. 4.2. 従来の計算アルゴリズムの問題点 ................................................. - 63 -. 4.3. 3 階直交テンソル積展開の計算アルゴリズムの改良 .................... - 64 -. 4.4. 数値実験 .......................................................................................... - 66 -. 4.4.1. 立方体テンソルの数値実験 ...................................................... - 66 -. 4.4.2. 非立方体テンソルの数値実験 .................................................. - 71 -. 4.5. 計算時間の比較 ............................................................................... - 77 -. 4.6. 画像データへの適用 ....................................................................... - 78 -. 4.7. 結言 ................................................................................................. - 87 -. 第 4 章の参考文献 .................................................................................... - 89 第 5 章結論 ............................................................................................... - 91 5.1. 本研究の総括 .................................................................................. - 91 -. 5.2. 今後の課題 ...................................................................................... - 92 -. 謝辞. (iv).

(6) 図目次図 2-1. 3 階テンソルのイメージ図 ....................................................... - 9 -. 図 2-2. 3 次元ディジタルフィルタの振幅応答特性のイメージ図 .... - 19 -. 図 2-3. 3-TPE による相対残差 ( L  M  N  12) ................................. - 24 -. 図 2-4. 3-OTPE による相対残差の対数グラフ ( L  M  N  3) .......... - 31 -. 図 2-5. 3-OTPE の計算時間の比較 ( L  M  N ) .................................. - 32 -. 図 2-6. 非負分解の収束性 ................................................................... - 37 -. 図 3-1. 3 階テンソル A  R I1I 2 I3 の展開行列への変換 ........................ - 42 -. 図 3-2. 計算時間の比較(3-TPE，最良ランク 1 近似) ........................ - 52 -. 図 3-3. 展開係数の対数グラフ(HOSVD，3-OTPE) ............................ - 56 -. 図 3-4. 相対残差(HOSVD，3-OTPE) .................................................. - 58 -. 図 3-5. 相対残差の対数グラフ(HOSVD，3-OTPE) ............................ - 58 -. 図 3-6. 計算時間の比較(HOSVD，3-OTPE) ....................................... - 59 -. 図 4-1. 展開係数( 3  3  3 の立方体テンソルの例の場合) ................... - 69 -. 図 4-2. 相対残差( 3  3  3 の立方体テンソルの例の場合) ................... - 70 -. 図 4-3. 相対残差の対数グラフ( 3  3  3 の立方体テンソルの例の場合) ...... ....................................................................................................... - 70 図 4-4. 残差最小値との差( 3  3  3 の立方体テンソルの例の場合 ) .... - 71 -. 図 4-5. 展開係数( 6  5  3 の非立方体テンソルの例の場合 ) ............... - 75 -. 図 4-6. 相対残差( 6  5  3 の非立方体テンソルの例の場合 ) ............... - 75 -. 図 4-7. 相対残差の対数グラフ( 6  5  3 の非立方体テンソルの例の場合 ) .. ....................................................................................................... - 75 図 4-8. 残差最小値との差( 6  5  3 の非立方体テンソルの例の場合 ) - 76 -. 図 4-9. 計算時間( m  m  m 立方体テンソルの場合) ............................ - 78 -. 図 4-10. 標準画像(Lenna) .................................................................... - 79 -. 図 4-11. 標準画像(Lenna)の RGB 分割画像 ........................................ - 79 -. 図 4-12. 相対残差( 16  16 の画像 Lenna の場合) .................................. - 80 (v).

(7) 図 4-13. 相対残差の対数グラフ ( 16  16 の画像 Lenna の場合) ........... - 80 -. 図 4-14. 残差最小値との差( 16  16 の画像 Lenna の場合) .................. - 81 -. 図 4-15. 実験画像(アルファベット順) ............................................... - 84 -. 図 4-16. 計算時間 ................................................................................ - 87 -. (vi).

(8) 表目次表 2-I. 3-TPE による展開ベクトル ( L  M  N  12) ........................... - 21 -. 表 2-II. 展開項間の内積 ( L  M  N  12) ............................................ - 22 -. 表 2-III. 3-TPE による展開係数 ( L  M  N  12) ................................ - 23 -. 表 2-IV. 3-OTPE による展開ベクトル ( L  M  N  3) ....................... - 29 -. 表 2-V. 展開項間の内積 ( L  M  N  3) ............................................. - 29 -. 表 2-VI. 3-OTPE による展開係数と相対残差 ( L  M  N  3) ............ - 30 -. 表 2-VII. 非負分解による展開ベクトル ( L  M  N  3) .................... - 35 -. 表 2-VIII. 非負分解の展開係数と相対残差 ........................................ - 36 -. 表 3-I. 最良ランク 1 近似と 3-TPE の展開ベクトル .......................... - 50 -. 表 3-II. 最良ランク 1 近似と 3-TPE の展開係数と相対残差 ............. - 51 -. 表 3-III. HOSVD と 3-OTPE による展開ベクトル.............................. - 53 -. 表 3-IV. HOSVD の展開項間の内積 ................................................... - 54 -. 表 3-V. HOSVD のコアテンソル成分値および 3-OTPE の展開係数 . - 55 -. 表 3-VI. 相対残差(HOSVD，3-OTPE) ................................................ - 57 -. 表 4-I. 展開ベクトル( 3  3  3 の立方体テンソルの例の場合 ) ............ - 67 -. 表 4-II. 展開係数( 3  3  3 の立方体テンソルの例の場合) .................. - 68 -. 表 4-III. 展開ベクトル( 6  5  3 の非立方体テンソルの例の場合 ) ..... - 73 -. 表 4-IV. 計算時間( m  m  m 立方体テンソルの場合 ) .......................... - 77 -. 表 4-V. 残差 5%以下となる項数と従来法に対する比率 (Lenna) ....... - 81 -. 表 4-VI. 展開項数の圧縮率(Lenna 5%近似の場合) ............................ - 82 -. 表 4-VII. 展開による復元画像 (その 1) ............................................... - 83 -. 表 4-VIII. 残差 5%以下となる項数と従来法に対する比率 ............... - 84 -. 表 4-IX. 展開による復元画像(その 2) ................................................ - 85 -. 表 4-X. 計算時間 ................................................................................. - 87 -. (vii).

(9)

(10) 第1章序論. 第1章. 序論. 1.1 研究の背景と目的多次元データを低次元のデータで近似することは，信号処理や画像処理，音声認識，パターン認識といった多くの分野で利用されている．信号処理の分野では波形処理に用いられる多次元ディジタルフィルタの設計法として，多次元の設計仕様を低次元化し，1次元ディジタルフィルタの組み合わせとして実現することが可能である[1, 2, 3]．画像処理では画像データ圧縮に応用されており圧縮率の向上が期待できる[4]．多次元データを低次元分解するための手法として，2次元データの場合には固有値分解(Eigen-Value Decomposition – EVD)や特異値分解(Singular Value Decomposition – SVD) を利用する手法がよく知られている[5, 6]．3次元以上のデータに対しては，SVDを多次元に拡張して定義した多次元外積展開(Multi-Dimensional Outer Product Expansion)が斎藤らによって提案されている[4]．その数値計算には非線形最適化法が用いられているが，非線形最適化法はデータのサイズが大きくなると多大の計算時間が必要である．そのため，村上らはべき乗法を利用して計算するアルゴリズムを提案し，3次元データに拡張した3階テンソル積展開(3rd-Order Tensor Product Expansion – 3-TPE)のアルゴリズムを提案している[7]．さらに村上は，斎藤らによる多次元外積展開の計算アルゴリズムによる結果の展開項間に直交性が成り立たないことを指摘して，直交性を満たすように修正した3階直交テンソル積展開(3rd-Order Orthogonal Tensor Product Expansion – 3-OTPE) のアルゴリズムも考案している．本研究の目的は，3階テンソル積展開によって得られる展開項の性質を調べることと，それを利用して計算アルゴリズムの計算精度と計算時間を向上させることである．そこで，村上によって提案された3-TPEと3-OTPEの計算アルゴリズムを検証し，計算精度と計算時間の改良を図ることにした． -1-.

(11) 第1章序論 3-OTPEによって得られる展開項は，3-TPEの計算の過程において反復ベクトルに展開項間の直交性を満たす条件を課して修正して計算されたものである．このように，べき乗法の過程で次々と計算される反復ベクトルに対して種々の拘束条件を課すことができる．村上はこのことに着目して2次元データの場合にべき乗法を用いた非負分解を計算するアルゴリズムを開発し，応用例として2次元ディジタルフィルタの設計問題への適用例を挙げている[8]．ディジタルフィルタの設計問題においては，フィルタの構造を簡単にするために，展開式の項数はなるべく少ないほうが望ましく，各外積項を構成するベクトルは1次元の振幅仕様となるので，ベクトルの各成分は非負でなければならないからである[3]．今回，3次元データを展開する3階テンソル積展開においても非負分解を行う計算アルゴリズム(3rd-Order Non-Negative Tensor Product Expansion – 3-NTPE)を開発し，3次元ディジタルフィルタの振幅応答特性の例について数値実験を試みた[9]．近年， 3階テンソル積展開の応用分野の一つでもあるパターン認識などの分野において， Lathauwerらによって提案された高次特異値分解 (Higher-Order Singular Value Decomposition – HOSVD)[10, 11]が注目されている[12, 13, 14]．そこで，HOSVDと3-TPE について展開計算法やそれぞれの計算によって得られる展開項の違い，計算精度や計算時間について比較を行うことにした[15]． HOSVDとの比較を通して，これまでの3-OTPEの計算アルゴリズムには多次元データの各次元のサイズが異なる場合に，展開項の和が元の多次元データの最小2乗近似になるという条件に合わないことが分かった．この問題を解決する新たな計算アルゴリズムを考案し，計算精度の向上を図ることにした[16]．. 1.2 本研究の概要本研究では，3階テンソル積展開の性質を調べることと，それを利用して計算アルゴリズムの計算精度と計算時間を改良することの2つを目的とし，特に展開項間の直交性を満足させた3-OTPEの計算アルゴリズムの改良を行った．本研究で行ったことをまとめると次のようになる． 3-OTPEの計算においては，各次元の次元数のうち最も小さい次元数までの展開項を計算するが，最後の項についてはそれまでに計算した展開項に直交するベクトルを求 -2-.

(12) 第1章序論めればよいから，べき乗法で計算する必要がないため最後の項を求める過程を Gram-Schmidtの直交化法[6, 17]に置き換えることで計算時間を改良した[9]．また， 3-OTPEのべき乗法の過程の中で次々と計算される反復ベクトルに対して直交条件を課して修正することを応用して，反復ベクトルを非負に制約する条件を課した3-NTPE の計算アルゴリズムを提案した[9]．次に，近年パターン認識などの分野で注目されているLathauwerらのHOSVDについて展開計算法を調べることにした．HOSVDには，多次元データの各次元のランクに注目した展開である，最良ランク(R1, R2, …, RN)近似(Best Rank (R1, R2, …, RN) Approximation) や最良ランク1近似(Best Rank 1 Approximation)と呼ばれる展開が提案されており[11]， 3-TPEと3-OTPEと展開項の性質や計算精度の違いを検証した[15]． 3階テンソル積展開の定義では元のテンソルを展開項の和で求める過程で残差を最小にする展開項を求めなくてはならないが，従来の3-OTPEの計算アルゴリズムでは，多次元データの各次元のサイズが異なる場合に，残差を最小とする展開項とならないことが分かり，この問題を解決するために新たな計算アルゴリズムを開発した[16]．. 1.3 各章の構成本論文は，全体で5章から構成されている．1.1節に示した研究の背景に基づき，前節の概要に記述した内容についてそれぞれ1つの章で述べる．第1章では，研究の背景と目的及び概要，さらに論文の構成を述べる．第2章では，3階テンソル積展開の定義，およびべき乗法を用いた3階テンソル積展開と3階直交テンソル積展開の計算アルゴリズムについて説明した後，3階直交テンソル積展開の最終項の計算方法の改良について述べる．また3階テンソル積展開によって得られる各展開項を構成する展開ベクトルに非負拘束条件を課した非負分解のアルゴリズムを述べ，3次元ガウス型ディジタルフィルタの振幅応答の周波数特性[18]を例にとって数値実験例として示す．第3章では HOSVDについてその計算法を説明する． HOSVDには最良ランク(R1, R2, …, RN)近似，最良ランク1近似と呼ばれる分解があり，3階テンソル積展開の手法および計算結果を比較してその類似性と相違点を列挙する．. -3-.

(13) 第1章序論第4章では，初めに，従来の3-OTPEの計算アルゴリズムの問題点として，テンソルの各次元のサイズが異なる場合に残項の計算が必要になりアルゴリズムが複雑になること，残項の計算にGram-Schmidtの直交化法を用いているため，初期値の選び方に依存して展開項には一意性がないことを指摘している．アルゴリズムの複雑さの軽減と残項の算出における初期値による問題を解決するために新たに考案した改良法を示し，従来の3-OTPEと改良した3-OTPEのアルゴリズムを用いた数値実験を通して改良の効果を示す．第5章では，本研究の総括を行い，得られた結論をまとめると共に，今後に残された課題について述べる．. -4-.

(14) 第1章序論. 第1章の参考文献 [1]. 川又政征, 樋口龍雄: 多次元ディジタル信号処理の基礎[Ⅲ], 電子情報通信学会誌, Vol.74, No.7, pp.757-762, 1991.. [2]. Tian-bo Deng, Masayuki Kawamata: Design of Two-Dimensional Recursive Digital Filters Based on the Iterative Singular Value Decomposition, Transaction IEICE, Vol.E73, No.6, pp.882-892, 1990.. [3]. 大木真, 川又政征, 樋口龍雄: インパルス応答仕様の外積展開に基づく 3 次元ディジタルフィルタの設計, 電子情報通信学会論文誌, Vol.J72-A, No.4, pp.648-655, 1989.. [4]. 斎藤弘, 小松隆, 原島博, 宮川洋: 多次元外積展開による静止画像の符号化, 電子情報通信学会論文誌, Vol.J68-B, No.4, pp.547-548, 1985.. [5]. James H. Wilkinson: The Algebraic Eigenvalue Problem, Oxford University Press, 1965.. [6]. 村田健郎: 線形代数と線形計算法序説, 第 4 章, サイエンス社, 1986.. [7]. 村上純, 山本直樹, 田所嘉昭: べき乗法を用いた 3 階テンソル積展開の高速計算法,電子情報通信学会論文誌, Vol.J82-A, No.8, pp.1351-1359, 1999.. [8]. 村上純, 山本直樹, 田所嘉昭: べき乗法による 2 次元ディジタルフィルタの設計仕様行列の非負分解, 情報処理学会論文誌, Vol.33, No.5, pp717-721, 1992.. [9]. Chiharu Okuma, Jun Murakami, Naoki Yamamoto: Calculation of 3-D Nonnegative Outer. Product Expansion by the Power Method and its Application to Digital Signal Processing, Proceedings of the Twelfth International Symposium on Artificial Life and Robotics, GS14-4, 2007. [10] Lieven De Lathauwer, Bart De Moor, Joos Vandewalle: A Multilinear Singular Value Decomposition, SIAM Journal on Matrix Analysis and Applications., Vol.21, No.4, pp.12531278, 2000. [11] Lieven De Lathauwer, Bart De Moor, Joos Vandewalle: On the Best Rank-1 and Rank-(R1, R2, …, RN) Approximation of Higher-Order Tensors, SIAM Journal on Matrix Analysis and Applications, Vol.21, No.4, pp.1324-1342, 2000. -5-.

(15) 第1章序論 [12] Berkant Savas, Lars Eldén: Handwritten Digit Classification using Higher Order Singular Value Decomposition, Pattern Recognition, Vol.40, pp.993-1003, 2007. [13] 森垣潤一, 片山薫: 高次特異値分解の画像分類への応用, DEWS2009 E10-2, 第 19 回データ工学ワークショップ, 2008. [14] 井上光平, 浦浜喜一: 行列データの主成分分析 MPCA の近似解法, 電子情報通信学会技術報告, PRMU2004-145, pp.67-70, 2004. [15] Chiharu Okuma, Jun Murakami, Naoki Yamamoto: Comparison between Higher-order SVD and Third-Order Orthogonal Tensor Product Expansion, International Journal Electronics, Communications and Computer Engineering, Vol.1, No.2, pp.131-137, 2009. [16] Chiharu Okuma, Naoki Yamamoto, Jun Murakami: An Improved Algorithm for Calculation of the Third-Order Orthogonal Tensor Product Expansion by Using Singular Value Decomposition International Journal of Electronics, Communications and Computer Engineering, Vol.2, No.1, pp.11-20, 2010. [17] Gilbert Strang 著, 山口昌哉監訳, 井上昭訳: 線形代数とその応用, 第 3 章, 産業図書, 1976. [18] 大木真, 川又政征: 周波数領域使用に対する並列完全分離型 3 次元ディジタルフィルタの設計, 電子情報通信学会論文誌, Vol.J72-A, No. 8, pp.1247-1252, 1989.. -6-.

(16) 第2章. 第2章. 3階テンソル積展開. 3階テンソル積展開. 2.1 緒言多次元データを低次元分解する手法として，2次元データの場合には固有値分解 (Eigen-Value Decomposition – EVD)や特異値分解(Singular Value Decomposition – SVD)がある[1, 2]．多次元データについては多次元外積展開[3]が知られているが，その計算には非線形最適化法が用いられており，データのサイズが大きい場合に計算時間も増大する[4, 5]．そこで本研究で扱う3階テンソル積展開の計算にはべき乗法を用いて高速化する方法が提案されている[6]．行列の固有値計算はべき乗法を用いると容易に求めることができるが，一般に行列は実対称行列であるとは限らず，任意の実行列について低次元分解を行うには特異値分解の計算が必要となる[7]．この章では，初めにべき乗法を用いた固有値分解および特異値分解の計算法を説明する．そして，それを応用して開発されたべき乗法を用いた3階テンソル積展開，展開項間に直交性を有する3階直交テンソル積展開について述べる．次に，べき乗法で次々に計算される反復ベクトルに非負拘束条件を課した3 階テンソルの非負分解の計算アルゴリズムを述べ，数値実験の結果を示す．. 2.2 テンソルの表記と演算多次元データは，その次元が N 次元のとき N 次元データと表現される．0次元データはスカラー，1次元データはベクトル，2次元データは行列を意味する．また多次元データはテンソル(tensor)とも呼ばれており，本論文では，主にテンソルの表現を用いている[8]．多次元データの次数はテンソルの表現では階数と呼ばれ， N 階テンソル -7-.

(17) 第2章. 3階テンソル積展開. と表現する．多次元データと同様に，0階テンソルはスカラー，1階テンソルはベクトル，2階テンソルは行列を意味する．高階テンソル(higher order tensor)と表現されるテンソルは3階以上のテンソルである． I I 2. 1階テンソルを A  R I1 ，2階テンソルを A R 1. ，3階テンソルを A R I1 I 2 I 3 のよ. I  I 2  I N. と表わされ，その (i1 , i2 , , iN ) 成分. うに表現すると， N 階テンソルは A  R 1. は各次元のインデックスを用いて， ai1i2 in ， (1  i1  I1 , , 1  in  I N ) と表す．本論文で扱うテンソルは3階までのテンソルであるから，表記を明確にするために，ベクトルである1階テンソルと，行列である2階テンソルは A, B のように太字の書体で表し，3階のテンソルのみ A, B のように手書き様の書体を用いて示すものとする．以下にテンソルの表記について示しておく． 1階テンソル(ベクトル) A  R I1 において， i1 を成分の配置を示すインデックス記号とするとき， A  R I1 の (i1 ) 成分を ai1 と表す．1階テンソル A は I 1 個の成分 ai1 を持つ．. A  [a1 a2  aI1 ]  [ai1 ]. (i1  1， 2, , I1 ).. (2.1). 2階テンソル(行列) A  R I1 I 2 において， i1 をテンソルの行方向の配置， i 2 を列方向の配置を示すインデックス記号とするとき， A  R I1 I 2 の( i1 , i2 )成分を ai1i2 と表す．2 階テンソル A は I 1  I 2 個の成分 ai1i2 を持つ．.  a11 a12 a a22 21 A     aI11 aI1 2.  a1I 2   a2 I 2   [ai1i2 ]      aI1I 2 . (i1  1， 2, , I1 , i2  1， 2, , I 2 ).. (2.2). 行列 A の i1 行は. [ai11 ai1 2  ai1I 2 ]. (i1  1， 2, , I1 ). (2.3). と表され，行列 A の i2 列は.  a1i2  a   2i2       a I1i2 . (i2  1， 2, , I 2 ). (2.4). と表せる．従って，行列 A は行方向に I 2 個の成分を持つ1階テンソルが I 1 個並んだもの，あるいは列方向に I 1 個の成分を持つ1階テンソルが I 2 個並んだものと考えることができ， I 1  I 2 行列と呼ぶ． -8-.

(18) 第2章. 3階テンソル積展開. 同様に3階テンソル A R I1 I 2 I 3 は成分に2階テンソル A R I1  I 2 を I 3 個並べたもの，あるいは成分に2階テンソル A  R I 2 I3 を I 1 個並べたもの，または成分に2階テンソル. A  R I1I3 を I 2 個並べたものとも考えることができる．図2-1に3階テンソルのイメージ図を示す．. I3 i3. ai1i2i3. I1 i1. i2. A R I I I 1. I2. 2. 3. 図2-1 3階テンソルのイメージ図図2-1において，3階テンソル A  R I1I 2 I3 は成分として2階テンソル A R I1  I 2 を I 3 個並べたものとみなせるから A  R ( I1I 2 )I3 とも表せる．この2階テンソル A R I1  I 2 をインデックス i3 により A i3 と表すことにすると，3階テンソル A は次のように表せる．. A  [A1 | A2 |  | A I 3 ]  [Ai3 ]. (i3  1, 2, , I3 ).. (2.5). ここで， | は単にデータの区切りを表す記号として用いている． A i3 の( i1 ,i2 )成分を ai1i2i3 と表すと，式(2.5)は，. A  [ai1i2i3 ]. (i1  1, 2, , I1，i2  1, 2, , I 2，i3  1, 2, , I3 ). (2.6). と表せる．ただし， ai1i2i3 は3階テンソル A  R I1I 2 I3 の (i1 , i2 , i3 ) 成分を意味し，テンソル A は I1  I 2  I 3 個の成分を持つ．次に，テンソルの演算の定義を述べる．和とスカラー倍に関しては，2階テンソルは行列であるから通常の行列の演算とし，3階テンソルの場合は文献[8]に基づいて次の通りに定める．. A  [ai1i2i3 ] と B  [bi1i2i3 ] が共に I1  I 2  I 3 テンソルであるとき， A と B の和を (i1 , i2 , i3 ) 成分が ai i i  bi i i である I1  I 2  I 3 テンソルと定義し A  B と書く．さらに c 1 2 3. 1 2 3. -9-.

(19) 第2章. 3階テンソル積展開. をスカラーとするとき， A と c のスカラー倍を (i1 , i2 , i3 ) 成分が cai1i2i3 である I1  I 2  I 3 テンソルと定義し cA と書く．また， A  (1)B を A と B の差と呼び A  B と書く． 2つのテンソル A と B があるとき，それぞれの成分の全ての組み合わせの積を成分とするテンソルを A と B のテンソル積と呼び， A  B と書く．. 2.3 べき乗法を用いた3階テンソル積展開多次元データの分解として2次元の場合にはEVDやSVDが知られている．多次元データの場合は，多次元外積展開が知られているが，その計算には非線形最適化法が用いられている．ただし，この手法は計算量が多く，大規模サイズの3次元データに対しては適用が困難である．そこで村上らはこの問題を解決するために，べき乗法による特異値分解の3次元データへの拡張として，べき乗法を用いた3階テンソル積展開のアルゴリズムを考案している．この章では初めにべき乗法を用いた固有値分解と特異値分解について示し，その拡張として提案されているべき乗法を用いた3階テンソル積展開を示す．. 2.3.1. べき乗法を用いた固有値分解. 成分が実数からなる対称行列(実対称行列)の固有値分解について述べる．サイズが N  N の実対称行列 A のEVDが次のように表されるとするとき，. Ax i  i xi. i  1,2,, N. (2.7). i は A の固有値， x i は i に対応する固有ベクトルと呼ばれる．固有値は， | 1 |  | 2 |    | N |. (2.8). となるように番号付けしておく．べき乗法による固有値分解は，まず，反復に使用する初期ベクトルを任意に選び. x ( 0) とする．このベクトルに行列 A を左から乗じて x (1) を求め，その x (1) にまた A を乗じて x( 2) を求める．この操作を繰り返して，絶対値最大の固有値とそれに対応する固有ベクトルを求める．さらに行列 A から求まった項を減じた行列を対象にして，次の固有値および固有値ベクトルを同様に求めることを繰り返す手法である．. - 10 -.

(20) 第2章. 3階テンソル積展開. 行列 A は半単純，即ち全ての固有値の重複度は1であると仮定する．任意の初期ベクトル x( 0) は固有ベクトル x1 , x 2 , , x N の線形結合 N. x( 0)  c1x1  c2 x 2    cN x N   ci xi. (2.9). i 1. として表される． c1 , c2 , , cN は定数である．.  (0)  || x(0) || と置き x(0)  x(0) /  (0) とする． x ( 0) に行列 A を左から乗じて， x (1)  Ax ( 0)  . c11. . ( 0). N.  i 1. c1. . x1 . ci i. (0). Ax1 . c22. . (0). c2. . (0). x2   . Ax 2    c N N.  (0). cN.  (0). Ax N. xN. (2.10). xi.  (0). となる．さらに  (1)  || x(1) || と置いて x(1)  x(1) /  (1) とし， x (1) に行列 A を乗じると， N. x( 2)  Ax (1)   i 1. ci i. N.  (0) (1). Ax i   i 1. ci i2.  (0) (1). xi. (2.11). となり，次々と行列 A を乗じることを繰り返す． x( k 1) が求まったとして，.  ( k 1)  || x ( k 1) || と置き x( k 1)  x( k 1) /  ( k 1) とすると， x(k ) は， N c k x ( k )   ( 0 ) (1)i i ( k 1) x i   i 1  (2.12) k k N   1   i   ci c1   x1   ( 0) (1)   xi     ( 0 ) (1) ( k 1)    ( k 1)  1   i 2      1  となる．反復回数 k が十分大きくなると，固有値の大きさの関係から (i / 1 ) k は無視 k 1. できるほど小さくなるから， x(k ) は c11k (2.13) x  ( 0) (1)  ( k 1) 1 (k ) (k ) となり， x は固有ベクトル x1 の定数倍に収束する．従って，固有ベクトル x1 は x x( k ) . をノルムが1となるように正規化すればよいから， ( k )  || x( k ) || とすると， x1 . x( k ).  (k ). (2.14). となる．実際の計算では，. x( k )  x( k 1)  . (2.15). となるまで繰り返す．ここで，  は設定した微小値とする．通常は，計算機イプシロン[9]を採用する．固有値 1 はRayleigh(レイリー)商[10] - 11 -.

(21) 第2章. 3階テンソル積展開. (x( k ) , Ax ( k ) ) (x( k ) , x( k 1) ) 1  ( k ) ( k )  ( k ) ( k ) (x , x ) (x , x ). (2.16). から計算でき， 1 は絶対値最大の固有値の近似値である． 2 以降を求めるには，元の行列 A を変形して階数を下げる．代表的な方法として次のHotelling(ホテリング)の方法がある[10]．. A2  A  1x1x1T. (2.17). と変形されるとき， x1T x1  1 となるように正規化されているから，行列 A 2 の固有値.  , N， 0 であり，それに対応する固有ベクトルは， x2 , x3 , , x N , x1 である．は， 2 , 3，. 3 以降の固有値についても同様に， A i 1  A i  i x i xTi. (2.18). を新しい行列と考えてべき乗法を繰り返すことで，全ての固有値および固有ベクトルを計算することができる．. 2.3.2. べき乗法を用いた特異値分解. 任意の実行列についてSVDは対称行列に限らず実数の範囲で分解が可能である[7]．サイズが M  N , M  N の実行列 A のSVDは次のように表される． r. A    i xi yTi .. (2.19). i 1. ここで，  i は特異値，ベクトル x i と y i は  i に対応する左，右特異ベクトルと呼ばれ， r は行列 A のランクを表す．. 特異ベクトル xi , y i は， Ay i   i xi , AT xi   i y i. (2.20). を満たすので， AT Ay i   i y i , AA T xi   i xi 2. 2. (2.21). の関係があり，特異値  i の2乗は AT A, AA T の固有値となり  i は実数となる．特異値は， | 1 |  |  2 |    |  r |. となるように番号付けしておく．式(2.19)を2.2節で定義したテンソル積で表現すると - 12 -. (2.22).

(22) 第2章. 3階テンソル積展開. r. A   i ( x i  y i ). (2.23). i 1. と表せる．べき乗法によるSVDは，EVDと同様に最大から数個の特異値と特異ベクトルを求めるのに有効である．いま，行列 A の特異値が全て異なると仮定する．任意の M 次元初期ベクトル x( 0) と N 次元初期ベクトル y ( 0) は特異ベクトル x1 , x2 , , xr および. y1, y 2 , , y r の線形結合で， r. x ( 0)  a1x1  a2 x 2    ar x r   ai x i , i 1 r. y. (0). (2.24).  b1y1  b2 y 2    br y r   bi y i i 1. と表せる．ここで， ai , bi は定数である．.  (0)  || x(0) ||,  (0)  || y (0) || と置き， x (0)  x (0) /  (0) , y (0)  y (0) /  (0) とする． x ( 0) , y ( 0) に行列 A と AT をそれぞれ左から乗じて， b b b x (1)  Ay ( 0)  (10 ) Ay 1  (20 ) Ay 2    (r0) Ay r.    b b b  1 ( 01) x1  2 ( 0 )2 x 2    r ( 0)r x r    r b   i ( 0i) xi , i 1  a a a y (1)  AT x ( 0 )  (10 ) AT x1  (20 ) AT x 2    (r0 ) AT x r    a1 1 a2 2 ar r  ( 0 ) y1  ( 0 ) y 2    ( 0 ) y r    r a   i ( 0 )i y i i 1 . (2.25). となる．さらに  (1)  || x(1) ||,  (1)  || y (1) || と置いて， x(1)  x(1) /  (1) , y (1)  y (1) /  (1) とし，. x(1) , y (1) に行列 A と AT を乗じると， r. x( 2)  Ay (1)  . r. Ay i  . ai i2. x,. ( 0 ) (1) i  ( 0)  (1)  i 1  r b  b 2  AT x (1)   ( 0i) i (1) AT xi   ( 0i ) i (1) y i   i 1  i 1  i 1 r. y ( 2). ai i. (2.26). となり，次々とこれを繰り返す．. x( 2 k 1) , y ( 2k 1) が求まったとして，  ( 2k 1)  || x( 2k 1) ||,  ( 2k 1)  || y ( 2k 1) || と置き， x( 2k 1)  x( 2k 1) /  ( 2k 1) , y ( 2k 1)  y ( 2k 1) /  ( 2k 1) とすると， x( 2 k ) , y ( 2 k ) は. - 13 -.

(23) 第2章. 3階テンソル積展開 ai i2 k   ( 0) (1) xi   ( 2 k 2)  ( 2 k 1) i 1  r. x. (2k ). . r   1 2k    12 k a ( ) x  ai ( i ) 2 k x i ,   1 1 ( 0 ) (1) ( 2 k  2 ) ( 2 k 1) 1     i 2  1 . b  2k   ( 0) (1) i i( 2 k 2 ) ( 2 k 1) y i    i 1  r. y. (2k ). (2.27). r   1 2k    12 k  ( 0 ) (1) b ( ) y  bi ( i ) 2 k y i    1 1 ( 2 k  2 ) ( 2 k 1) 1     i 2  1 . となる．反復回数 k が十分大きくなると，特異値の大きさの関係から ( i /  1 ) 2 k は無視できるほど小さくなるから， x( 2 k ) , y ( 2 k ) は，. a1 12 k x  ( 0) (1) x,    ( 2 k  2)  ( 2 k 1) 1 b1 12 k (2k ) y  ( 0) (1) y     ( 2 k  2) ( 2 k 1) 1 (2k ). (2.28). となり，特異ベクトル x1 , y1 の定数倍に収束する．従って，特異ベクトル x1 , y1 は. x( 2 k ) , y ( 2 k ) をノルムが 1 となるように正規化すればよいから，  ( 2k )  || x( 2k ) ||,.  ( 2k )  || y ( 2k ) || とすると x( 2k ). x1  x ( 2 k ) .  (2k ). y1  y ( 2 k ) . ,. y (2k ). (2.29).  (2k ). となる．実際の計算では，. x ( 2 k )  x ( 2 k 2)   , y ( 2 k )  y ( 2 k  2)  . (2.30). となるまで繰り返す．ここで，  は設定した微小値とする．特異値  1 は固有値の場合と同様に，Rayleigh商 (x ( 2 k ) , Ax ( 2 k ) ) (x ( 2 k ) , x( 2 k 1) )  , (x( 2 k ) , x ( 2 k ) ) (x( 2 k ) , x ( 2 k ) ) (y ( 2 k ) , Ay ( 2 k ) ) (y ( 2 k ) , y ( 2 k 1) ) 1   (y ( 2 k ) , y ( 2 k ) ) (y ( 2 k ) , y ( 2 k ) ). 1 . (2.31). から計算でき，  1 は絶対値最大の特異値の近似値である．  2 を求めるには，固有値の計算の場合と同様に元の行列 A を次のように書き換えてHotellingの方法を適用する．. A 2  A  1x1y 1T ..  3 以降も順次同様に求めることができる． - 14 -. (2.32).

(24) 第2章. 2.3.3. 3階テンソル積展開. 3階テンソル積展開の定義. 2.2節では本論文で扱うテンソルの定義と表記について述べた．ここでは斎藤らの多次元外積展開[3]に倣って3階テンソル積展開の定義を示す．サイズが L  M  N の3階テンソル A は，式(2.23)の行列 A のSVDの表現を拡張して，テンソル積によって次のように表せる． r. A    i (ui  v i  w i ) ．. (2.33). i 1. 3階テンソルのこの形の分解を3階テンソル積展開と呼ぶ．ここで，ベクトル. ui , vi , wi はそれぞれ， ui  vi  wi . L.  u ( j) j 1. 2.  1,. 2.  1,. i. M.  v ( j) j 1. i. N.  w ( j) j 1. 2. i. (2.34). 1. を満たす．ただし， u i ( j ) はベクトル u i の第 j 成分を表す．ベクトル v i , w i についても同様である．本論文ではこれらのベクトルのことを展開ベクトル(expansion vector) と呼び，また  i はSVDの特異値に相当するスカラー値で展開係数(expansion coefficient) と呼ぶ．行列のランクに相当する r は3階テンソルのランクと表現し，展開係数  i は絶対値の大きさの順に番号付けられている． | 1 |  |  2 |    |  r | .. (2.35). さらに，式(2.33)の展開式には以下の性質がある． (1) 展開式の各項は直交し，異なる展開項間の内積は， (u j  v j  w j )  (u k  v k  w k )  0,. jk. (2.36). となる．ただし，サイズが L  M  N の2つの3階テンソル A と B の内積は，2つのテンソルの対応する成分値の積の総和により次式で定義する．. A  B   A (i, j, k )B (i, j, k ). i , j ,k. - 15 -. (2.37).

(25) 第2章. 3階テンソル積展開. ~ (2) 展開式をある項数 p で打ち切ったとき，その項までの展開項の和 A p を便宜的. に p 項による近似テンソル，あるいは単に近似テンソルと呼ぶことにすると，近似テンソルは p. ~ A p   i (u i  v i  w i ). (2.38). i 1. で得られ，元のテンソルに対する最小2乗近似を与える．元のテンソル A から ~ 近似テンソル A p を引いたものを残差テンソルと呼び，次式で定義する． ~ E p  A  Ap. (2.39) ~ また，元のテンソル A と近似テンソル A p との誤差を残差テンソル E p のノルム || E p || と元のテンソルのノルム || A || を用いて，. ep . || E p ||. (2.40). || A ||. と定義する．これを，元のテンソルと近似テンソルとの相対残差と呼ぶ．ただし，3階テンソル A のノルムは次式で定義する． || A ||.  A (i, j, k ) . 2. (2.41). i,j,k. 2.3.4. べき乗法を用いた3階テンソル積展開の計算. べき乗法を用いたSVDの計算は，まず2個の初期ベクトルを選び，一方に対して行列自身を，他方に対して行列の転置行列を乗じ，得られた2個のベクトルに対して前者に転置行列を，後者に行列自身を乗じる計算を繰り返すことで，最終的に特異値と特異ベクトルが求められた．3階テンソル積展開においては，まずテンソルを行列に変換し，得られた行列についてSVDと同様な反復計算を行う．まず，3階テンソル A  R I1I2I3 とベクトル z  R I3 の積 Z R I1I2 を次のように定義する． I3. Z  Az   A i3 z i3 .. (2.42). i3 1. この積を A と z の縮約[8]と呼ぶ．縮約によって得られる行列 Z の (i, j ) 成分は次式で表せる． N. Z(i, j )   A (i, j, k )z (k ). k 1. - 16 -. (2.43).

(26) 第2章. 3階テンソル積展開. Algorithm 2.0に，べき乗法を用いた3階テンソル積展開によるテンソル A の絶対値最大の展開係数  1 および展開ベクトル u1, v1, w1 を求めるアルゴリズムを示す． Algorithm 2.0: 絶対値最大の展開係数と展開ベクトルの計算 (1) p  0 として初期ベクトル u1( p ) , v1( p ) , w1( p ) を任意に選び，正規化して. u1( p ) , v1( p ) , w1( p ) とする． (2) A とベクトル w1( p ) の縮約により， L  M 行列 F を次式により求める．. F  Aw1( p ) .. (2.44). 次に， F とベクトル v1( p ) の積からベクトル u1'( p ) ，さらに F の転置行列 F T と先に得たベクトル u1'( p ) の積からベクトル v1'( p ) を次のように求める．. u1'( p )  Fv1( p ) , v1'( p )  FT u1'( p ) .. (2.45). 同様にして， M  N 行列 G と N  L 行列 H からベクトル w1'( p ) , u1"( p ) , v1"( p ) , w1"( p ) を得る．. G  Av1'( p ) , w1'( p )  Gu1'( p ) , u1"( p )  G Tw1'( p ) ,. (2.46). H  Au1"( p ) , v1"( p )  Hw1'( p ) , w1"( p )  HT v1"( p ) .. (2.47). u1"( p ) , v1"( p ) , w1"( p ) を正規化して u1"( p ) , v1"( p ) , w1"( p ) と置き， u1"( p )  u1( p )   ,. v1"( p )  v1( p )   ,. w1"( p )  w1( p )  . (2.48). となるまで u1( p 1)  u1"( p ) , v1( p 1)  v1"( p ) , w1( p 1)  w1"( p ) として，さらに p の値を1増やして(2)を繰り返す．ここで，  は微小値とする． (3) (2) の条件が成立したときのベクトル u1"( p ) , v1"( p ) , w1"( p ) により展開ベクトルを近似できるものとして u1, v1, w1 と表す．展開係数  1 は， u1, v1, w1 のテンソル積と元のテンソル A との内積で次式により求める．.  1  A  (u1  v1  w1 ).. (2.49). (4) 展開係数 1 は絶対値最大の展開係数となり，展開ベクトル u1, v1, w1 のテンソル積から得られる展開項，. - 17 -.

(27) 第2章. 3階テンソル積展開. 1 (u1  v1  w1 ). (2.50). は，元のテンソル A の最小2乗近似を与える． (Algorithm 2.0 終わり). べき乗法を用いた3階テンソル積展開の第 n 項を求めるアルゴリズムは，元のテンソル A から既に得られている n  1 項の和を引いた残差テンソル E n を n 項目の3階テンソル積展開の対象として，絶対値最大の展開係数と展開ベクトルを求めるものである．残差テンソル E n は n 1. E n  A   i (ui  v i  w i ). (2.51). i 1. で表される．そこで，Algorithm 2.0を拡張して第 n 項を求めるアルゴリズムとして Algorithm 2.1に示す．これをべき乗法を用いた3階テンソル積展開の計算アルゴリズムとする．. Algorithm 2.1: べき乗法を用いた3階テンソル積展開. (1) 展開の対象とするテンソルを n 1. E n  A   i (ui  v i  w i ). (2.52). i 1. とする．ただし， E1  A とする． (2) p  0 として初期ベクトル u(np ) , v(np ) , w (np ) を任意に選び，正規化して un( p ) , v(np ) , wn( p ) とする．. (3) E n とベクトル w n( p ) の縮約により， L  M 行列 F を求める． F  E n w n( p ) .. (2.53). つぎに， F とベクトル v (np ) の積からベクトル u 'n( p ) ，さらに F の転置行列 F T と先に得たベクトル u 'n( p ) の積からベクトル v 'n( p ) を次式で求める． u'n( p )  Fv (np ) , v 'n( p )  FT u'n( p ) .. (2.54). 同様にして， M  N 行列 G と N  L 行列 H からベクトル w 'n( p ) , u "(n p ) , v "(n p ) , w "(n p ) を次式により得る． - 18 -.

(28) 第2章. 3階テンソル積展開. G  En v'n( p ) , w'n( p )  Gu 'n( p ) , u"n( p )  G Tw'n( p ) ,. (2.55). H  Enu"n( p ) , v"n( p )  Hw 'n( p ) , w"n( p )  HT v"n( p ) .. (2.56). u"n( p ) , v"n( p ) , w"n( p ) を正規化して un"( p ) , v"n( p ) , wn"( p ) と置き，. un"( p )  un( p )   ,. v "n( p )  v (np )   ,. w n"( p )  w n( p )  . (2.57). となるまで un( p1)  un"( p ) , v(np1)  v"n( p ) , wn( p1)  wn"( p ) として，さらに p の値を1増やして(3)を繰り返す．  は微小値とする． (4) (3)の条件が成立したときのベクトル un"( p ) , v"n( p ) , wn"( p ) を第 n 展開ベクトルの近似とし u n , v n , w n と表す．第 n 展開係数  n は， u n , v n , w n のテンソル積と残差テンソル E n との内積で次のようにして求める．.  n  E n  (u n  v n  w n ).. (2.58). 以上を指定された項数となるまで n の値を1増やして繰り返す． (Algorithm 2.1 終わり). 2.3.5. 3階テンソル積展開の計算結果. 3階テンソルの例として，図2-2に示す3次元ディジタルフィルタの振幅応答特性を持つ設計仕様行列[12]を扱うことにする．. 図2-2 3次元ディジタルフィルタの振幅応答特性のイメージ図. - 19 -.

(29) 第2章. 3階テンソル積展開. この例は，3階テンソルの例として信号処理の分野においてよく用いられるものである[4, 5]．このフィルタの振幅応答は，中心部が1で中心から離れるに従い，徐々に0 に近づく．3階テンソルの成分として次のように表せる． A (i, j , k )  h d ( xi , y j , z k ),  1.0  (0.6  r ) h d ( xi , y j , z k )    0.2  0.0. (2.59). (0.0  r  0.4), (0.4  r  0.6) , (0.6  r  1.0).. ここで， r は， r. 1. xi2  y 2j  z k2.  であり， xi , y j , zk は次式の値を取る．. (2.60). i  (0  i  L'1),  xi  L'1  j  (0  j  M '1), (2.61)  yj  M '  1   z k  k (0  k  N '1).  N '1  振幅仕様 hd ( xi , y j , zk ) は r  0.6 で0となることから， r  0.6 の部分を無視して考え. て構わないため，テンソルのサイズ L  M  N は，. L  L'0.6, M  M '0.6, N  N '0.6. (2.62). と縮小して考えることができる． . 展開ベクトル. この例についてテンソルのサイズを L  M  N  12 としてべき乗法を用いた 3 階テンソル積展開(3-TPE)で 30 項まで展開計算を行った．得られた 30 項の展開ベクトルのうち最初の 3 項までを表 2-I に示す．このテンソルは図 2-2 に示したように球状の振幅応答の 1/8 部分のデータを表すため，球体が持つ対称性がある．そのため展開項を構成するベクトルの組 ui , vi , wi のそれぞれのベクトルは符号を除いて同じ成分値を持っていることが確かめられる．. - 20 -.

(30) 第2章. 3階テンソル積展開. 表2-I 3-TPEによる展開ベクトル ( L  M  N  12) i. ui. vi. wi. 1.   3.798122 E  01   3.781495 E  01     3.730864 E  01     3.640552 E  01   3.499341 E  01     3.286772 E  01   2.977177 E  01     2.542092 E  01   1.937813 E  01     1.238805 E  01     6.338168 E  02   1.898468 E  02.   3.798122 E  01   3.781495 E  01     3.730864 E  01     3.640552 E  01   3.499341 E  01     3.286772 E  01   2.977177 E  01     2.542092 E  01   1.937813 E  01     1.238805 E  01     6.338168 E  02   1.898468 E  02.   3.798122 E  01   3.781495 E  01     3.730864 E  01     3.640552 E  01   3.499341 E  01     3.286772 E  01   2.977177 E  01     2.542092 E  01   1.937813 E  01     1.238805 E  01     6.338168 E  02   1.898468 E  02. 2.   4.376303 E  01   4.253504 E  01     3.887286 E  01     3.263883 E  01  2.318387 E  01    9.662223 E  02   6.783292 E  02     2.292411 E  01   3.222034 E  01     3.024091 E  01     2.166542 E  01   8.645486 E  02.   4.376303 E  01   4.253504 E  01     3.887286 E  01     3.263883 E  01   2.318387 E  01     9.662223 E  02   6.783292 E  02     2.292411 E  01   3.222034 E  01     3.024091 E  01     2.166542 E  01   8.645486 E  02.   4.376303 E  01   4.253504 E  01     3.887286 E  01     3.263883 E  01   2.318387 E  01     9.662223 E  02   6.783292 E  02     2.292411 E  01   3.222034 E  01     3.024091 E  01     2.166542 E  01   8.645486 E  02. 3.   8.276701E  02   6.575784 E  02     1.406363 E  02     7.435672 E  02   1.906600 E  01     3.174659 E  01   4.348816 E  01     5.038617 E  01   4.772634 E  01     3.582788 E  01     2.063133 E  01   6.523137 E  02.   8.276701 E  02   6.575784 E  02     1.406363 E  02     7.435672 E  02   1.906600 E  01     3.174659 E  01   4.348816 E  01     5.038617 E  01   4.772634 E  01     3.582788 E  01     2.063133 E  01   6.523137 E  02.   8.276701 E  02   6.575784 E  02     1.406363 E  02     7.435672 E  02   1.906600 E  01     3.174659 E  01   4.348816 E  01     5.038617 E  01   4.772634 E  01     3.582788 E  01     2.063133 E  01  6.523137 E  02. - 21 -.

(31) 第2章. 3階テンソル積展開. 展開項間の直交性を調べるために，展開項間の内積を計算した．得られた展開項のうち初めの5項までについて表2-IIに示したが，これらの展開項間には直交性がないことが分かる．残りの項についても同様であった．表2-II 展開項間の内積 ( L  M  N  12) 展開項の組. . 内積値. (u1 , v1 , w1 ) (u1 , v1 , w1 ). 1.000000E+00. (u1 , v1 , w1 ) (u 2 , v 2 , w 2 ). 1.700196E-01. (u1 , v1 , w1 ) (u3 , v 3 , w 3 ). 1.817749E-01. (u1 , v1 , w1 ) (u 4 , v 4 , w 4 ). 7.933135E-01. (u1 , v1 , w1 ) (u5 , v 5 , w 5 ). 6.047984E-01. (u 2 , v 2 , w 2 ) (u 2 , v 2 , w 2 ). 1.000000E+00. (u 2 , v 2 , w 2 ) (u3 , v 3 , w 3 ). 1.426595E-01. (u 2 , v 2 , w 2 ) (u 4 , v 4 , w 4 ). 1.891421E-01. (u 2 , v 2 , w 2 ) (u5 , v 5 , w 5 ). 1.632663E-01. (u3 , v 3 , w 3 ) (u3 , v 3 , w 3 ). 1.000000E+00. (u3 , v 3 , w 3 ) (u 4 , v 4 , w 4 ). 1.250297E-01. (u3 , v 3 , w 3 ) (u5 , v 5 , w 5 ). 9.676634E-03. (u 4 , v 4 , w 4 ) (u 4 , v 4 , w 4 ). 1.000000E+00. (u 4 , v 4 , w 4 ) (u5 , v 5 , w 5 ). 8.014975E-02. (u5 , v 5 , w 5 ) (u5 , v 5 , w 5 ). 1.000000E+00. 展開係数と相対残差. 式(2.58)によって展開係数を計算した．3-TPE による計算では，初期値の取り方によって展開係数は必ずしも絶対値の大きい順に求まらないことがあるため，絶対値の大きさの順に並べ替えた 20 項までの展開係数 p を表 2-III に示した． p が実際に求まった順番は No.の列に示してある．同表には文献[12]に示されている非線形最適化法による展開係数も載せたが，これら 20 項の展開係数は有効桁数７桁の精度で一致している．展開の一意性から 3-TPE と，非線形最適化法による展開は同じ展開を計算して. - 22 -.

(32) 第2章. 3階テンソル積展開. いると考えられ，非線形最適化法による展開ベクトルも展開項間の直交性を満たさないことが分かる．表2-III 3-TPEによる展開係数 ( L  M  N  12) べき乗法による 3階テンソル積展開 p. p. No.. 非線形最適化法による多次元外積展開 p. 1. 2.275862E+01. 1. 2.275862E+01. 2. 4.283573E+00. 2. 4.283573E+00. 3. 3.025678E+00. 3. 3.025678E+00. 4. 1.400982E+00. 4. 1.400982E+00. 5. 1.129514E+00. 5. 1.129514E+00. 6. 6.526013E-01. 6. 6.526013E-01. 7. 3.698252E-01. 11. 3.698252E-01. 8. 3.454422E-01. 7. 3.454422E-01. 9. 3.403220E-01. 8. 3.403220E-01. 10. 2.840046E-01. 10. 2.840046E-01. 11. 2.598434E-01. 12. 2.598434E-01. 12. 2.108770E-01. 9. 2.108770E-01. 13. 2.008905E-01. 14. 2.008905E-01. 14. 1.767691E-01. 13. 1.767691E-01. 15. 1.535205E-01. 18. 1.535205E-01. 16. 1.494289E-01. 15. 1.494289E-01. 17. 1.064934E-01. 22. 1.064934E-01. 18. 1.006265E-01. 19. 1.006265E-01. 19. 1.002772E-01. 23. 1.002772E-01. 20. 9.568394E-02. 26. 9.568394E-02. 3-TPEによって得られた展開ベクトルと展開係数を用いて p 項， p  1, 2,  , 30 による近似テンソルを求めた．図2-3は近似テンソルと元のテンソルとの相対残差をグラフで示している．最初の5項程度で元のテンソルを約5%まで近似できていることが分かり，実用的には十分少ない項数で元のテンソルを近似することが可能である．. - 23 -.

(33) 3階テンソル積展開. 3.00E-01 2.50E-01 2.00E-01 相対残差. 第2章. 1.50E-01. 1.00E-01 5.00E-02 0.00E+00 0. 3. 6. 9. 12 15 18 21 展開項数 p. 24. 図2-3 3-TPEによる相対残差 ( L  M  N  12). - 24 -. 27. 30.

(34) 第2章. 3階テンソル積展開. 2.4 3階直交テンソル積展開前節に示したように，3階テンソル積展開の計算をべき乗法や非線形最適化法で求めると，展開項間の直交性が全ての展開項の間では成立しない．村上は全ての展開項間の直交性が成立する分解法を提案している[13]．. 2.4.1. 3階直交テンソル積展開の定義. サイズが L  M  N の3階テンソル A は，サイズがそれぞれ L, M , N の正規直交基底の組{ u1, u2 ,  , u L }，{ v 1 , v 2 , , v M }，{ w 1 , w 2 , , w N }から1個ずつ取り出したベクトルのテンソル積の1次結合により，次のように展開できる．. A. . ijk. (ui  v j  w k ). (i  1,2,, L, j  1,2,, M , k  1,2, N ).. (2.63). i, j,k. この展開を3階直交テンソル積展開(3rd-Order Orthogonal Tensor Product Expansion – 3OTPE)と呼ぶ．展開ベクトルは正規直交ベクトルであるから，各展開項間には次のような関係があり直交する． (u j  v j  w j )  (u k  v k  w k )  {uTj u k }{vTj v k }{wTj w k }  0,. 2.4.2. j  k.. (2.64). べき乗法を用いた3階直交テンソル積展開の計算. 2.3.4で述べた3-TPEの計算アルゴリズムに，各展開ベクトルが正規直交基底となるような付加条件を課すことにより，3階直交テンソル積展開の計算を行うことができる．直交化にはGram-Schmidtの直交化法[2, 14]を用いる．いま，テンソル A のサイズを L  M  N とする． L, M , N の最小のサイズを m ，. m  min( L, M , N ) とするとき，まず m 個の展開ベクトルの組 { u i , v i , w i } ，. i  1, 2, , m を求め，さらにこれらに直交する残りの展開ベクトルを計算する．第 n 項を求める計算アルゴリズムをAlgorithm 2.2に示す．新たなベクトルが得られたとき，既に得られている展開ベクトルと直交するように常にGram-Schmidtの直交化法を用いて修正する．. - 25 -.

(35) 第2章. 3階テンソル積展開. Algorithm 2.2: べき乗法を用いた3階直交テンソル積展開.  m 個の展開ベクトルの組の計算 (1) これまでに得られている展開ベクトルのテンソル積と展開係数の1次結合により，元のテンソル A との残差を成分に持つ残差テンソル E n を求める． n 1. E n  A   i (u i  v i  w i ).. (2.65). i 1. ただし， E1  A とする． (2) p  0 として，初期ベクトル u (np ) , v (np ) , w (np ) を任意に選び，正規化して u n( p ) , v (np ) , w n( p ) とする．Gram-Schmidtの直交化法に従って，この初期ベクトル. から既に得られている項の n  1 個の展開ベクトルの成分を除いて，ベクトル ˆ (np ) を計算する．既に得られているベクトルは正規直交となってい uˆ (np ) , vˆ (np ) , w ˆ (np ) は，るから， uˆ (np ) , vˆ (np ) , w T T T uˆ (np )  un( p )  (u1 un( p ) )u1  (u 2 un( p ) )u 2    (u n1 un( p ) )u n1 ,. (2.66). T T T vˆ (np )  v (np )  ( v1 v (np ) ) v1  ( v 2 v (np ) ) v 2    ( v n1 v (np ) ) v n1 ,. (2.67). ˆ (np )  w n( p )  (w1T w n( p ) )w1  (w 2T w n( p ) )w 2    (w n1T w n( p ) )w n1 w. (2.68). ~ ( p) , ~ ~ ( p ) と置で得られる．これらを正規化して，正規直交初期ベクトル u v n( p ) , w n n. き直す．. ˆ (np ) uˆ (np ) vˆ (np ) w ( p) ( p) ( p) ~ ~ ~ un  ( p) , v n  ( p) , w n  . ˆ (np ) || || uˆ n || || vˆ n || || w. (2.69). (3) (1)で計算した残差テンソル E n に対して，以下を繰り返す． 3-TPEの計算アルゴリズムと同様に，残差テンソル E n と(2)で直交化した初期ベ ~ ( p) , ~ ~ ( p ) の縮約により， L  M 行列 F M  N 行列 G と N  L 行列クトル u v n( p ) , w ， n n H を求め，新たなベクトル u'n( p ) , v 'n( p ) , w 'n( p ) , u"n( p ) , v"n( p ) , w"n( p ) を次式で求める． ~ ( p ) , u '( p )  F~ F  En w v n( p ) , n n. v 'n( p )  FT uˆ 'n( p ) ,. (2.70). ˆ 'n( p ) , G  E n vˆ 'n( p ) , w 'n( p )  Guˆ 'n( p ) , u"n( p )  G T w. (2.71). H  E nuˆ "n( p ) ,. (2.72). ˆ 'n( p ) , w"n( p )  HT vˆ "n( p ) . v"n( p )  Hw. - 26 -.

(36) 第2章. 3階テンソル積展開. 式(2.70)，(2.71)，(2.72)で得られるこれらのベクトルは，(2)と同様にGramˆ 'n( p ) , uˆ "n( p ) , vˆ "n( p ) , w ˆ "n( p ) として常に既に得ら Schmidtの直交化法によって uˆ 'n( p ) , vˆ 'n( p ) , w. れている展開ベクトルと直交性を持つように修正する．これらのベクトル ~ "( p ) , ~ ~ "( p ) と置き， ˆ "( p ) を正規化して u uˆ "( p ) , vˆ "( p ) , w v "( p ) , w n. n. n. ~ ( p1) となるまで u n. n. n. n. ~" ( p)  u ~ ( p)   , u n n   ~" ( p) ~ ( p) (2.73)  vn  vn   ,  ~ " ( p) ~ ( p)  wn  wn    "( p ) ~ ~ ( p1)  w ~ "( p ) として，さらに p の値を1  un , ~ v n( p1)  ~ v "n( p ) , w n n. 増やして(3)を繰り返す．  は微小値とする． ~ "( p ) , ~ ~ "( p ) を第 n 展開ベクトルと (4) (3)の条件式が成立したときのベクトル u v "n( p ) , w n n. し u n , v n , w n と表す． n 項目の展開係数  n は次式で求める．.  n  E n  (u n  v n  w n ).. . (2.74). 残りの1次独立な展開ベクトルの組の計算. テンソルの各次元のサイズが全て同じ場合は，上記のプロセスで全ての直交ベクトルが求まるが，そうでない場合は1個あるいは2個以上の次元のベクトルの数は. m より多く，まだ得られていないベクトルが存在する．そこで，既に求めた展開ベクトルからGram-Schmidtの直交化法を用いて，残りの展開ベクトルを求める．この計算アルゴリズムは以下の通りである． L  m の場合についてのみ記述するが. M  m, N  m の場合についても同様に求めることができる． (1) n  m  1とする． (2) この時点までに得られている n  1 個の展開ベクトルに重複しない正規化された初期ベクトル u n を任意に定める． (3) (2)で定めた初期ベクトル u n からこの時点までに得られている n  1 個の展開ベクトル方向の成分を引いて， - 27 -.

(37) 第2章. 3階テンソル積展開 T T T uˆ n  un  (u1 un )u1  (u 2 un )u 2    (u n1 un )u n1. (2.75). ~  uˆ / || uˆ || として正規直交ベクトルとする．残りのベクトルが全とする． u n n n て求まるまで(1)から(3)を繰り返す． (Algorithm 2.2 終わり). . 展開係数の計算. Algorithm 2.2によって得られた展開ベクトルを用いて，展開係数を次の手順で計算する． Algorithm 2.3: 3-OTPEの展開係数の計算 (1) i  1, j  1, k  1 とする． (2) 展開ベクトル u i , v j , w k の3階テンソル積 u i  v j  w k を構成する． (3) 元の3階テンソル A と(2)のテンソル積との内積を取り，.  ijk  A  (u i  v j  w k ). (2.76). として  ijk を展開係数とする． i  1, 2, , L ， j  1, 2, , M ， k  1, 2, , N の全ての組み合わせについて展開係数  ijk が得られるまで，(2)，(3)を繰り返す． (Algorithm 2.3 終わり). 2.4.3. 3階直交テンソル積展開の計算結果. 2.3.5で示した3次元ディジタルフィルタの振幅応答特性の数値例について，. L  M  N  3 としたテンソルを3-OTPEで展開した．このサイズのテンソルは， L  M  N  27 項の和で元のテンソルを正確に表現できる． . 展開ベクトル. 表2-IVに，3-OTPEにより得られた展開ベクトルを示す．. - 28 -.

(38) 第2章. 3階テンソル積展開. 表2-IV 3-OTPEによる展開ベクトル ( L  M  N  3) i. ui. vi. wi. 1.  6.869175 E  01  6.222948 E  01    3.753580 E  01.  6.869175 E  01  6.222948 E  01    3.753580 E  01.  6.869175 E  01  6.222948 E  01    3.753580 E  01.  3.221110 E  01  2.022857 E  01    9.248379 E  01  6.514516 E  01  7.561938 E  01     6.149406 E  02.  3.221110 E  01  2.022857 E  01    9.248379 E  01  6.514516 E  01  7.561938 E  01     6.149406 E  02.  3.221110 E  01  2.022857 E  01    9.248379 E  01  6.514516 E  01  7.561938 E  01    6.149406 E  02. 2. 3. この展開ベクトルの全ての組み合わせについて内積を計算した結果を表2-Vに示す．互いに異なる展開項間の内積値は1016 のオーダーでほぼ0となっており，展開項間の直交性を満たしている．表2-V 展開項間の内積 ( L  M  N  3) 展開項の組. . 内積値. (u1 , v1 , w1 ) (u1 , v1 , w1 ). 1.000000E+00. (u1 , v1 , w1 ) (u 2 , v 2 , w 2 ). 1.480297E-16. (u1 , v1 , w1 ) (u3 , v 3 , w 3 ). -9.147775E-16. (u 2 , v 2 , w 2 ) (u 2 , v 2 , w 2 ). 1.000000E+00. (u 2 , v 2 , w 2 ) (u3 , v 3 , w 3 ). -8.211024E-16. (u3 , v 3 , w 3 ) (u3 , v 3 , w 3 ). 1.000000E+00. 展開係数と相対残差. 表 2-VI に Algorithm 2.3 によって求めた展開係数  ijk を絶対値の大きい順に示し，併せて元のテンソルと展開式を途中で打ち切って計算した近似テンソルとの相対残差を示す．展開係数  ijk に対応する展開ベクトル (u i , v j , w k ) の組み合わせはインデックスの欄に  (i, j, k ) として示した．図 2-4 に相対残差の対数グラフを示しているが，相対残差は 18 項で10 10 のオーダーに急速に減少し，さらに 24 項では10 16 のオーダーと - 29 -.

(39) 第2章. 3階テンソル積展開. なっている．3-OTPE では L  M  N 項の和で元のテンソルを正確に表現できることが分かる．表2-VI 3-OTPEによる展開係数と相対残差 ( L  M  N  3) No.. 展開係数  ijk. インデックス.  (i, j, k ). 相対残差. 1. 3.80071E+00. σ(1,1,1). 2.60734E-01. 2. 5.09803E-01. σ(2,1,2). 2.26304E-01. 3. 5.09803E-01. σ(2,2,1). 1.85593E-01. 4. -5.09803E-01. σ(1,2,2). 1.32952E-01. 5. 3.40569E-01. σ(2,2,2). 1.00959E-01. 6. 1.60926E-01. σ(3,2,1). 9.23133E-02. 7. 1.60926E-01. σ(3,1,2). 8.27700E-02. 8. -1.60926E-01. σ(1,2,3). 7.19722E-02. 9. 1.60926E-01. σ(2,1,3). 5.92378E-02. 10. -1.60926E-01. σ(1,3,2). 4.28747E-02. 11. 1.60926E-01. σ(2,3,1). 1.29368E-02. 12. -2.73712E-02. σ(3,1,3). 1.09098E-02. 13. 2.73712E-02. σ(1,3,3). 8.40750E-03. 14. -2.73712E-02. σ(3,3,1). 4.72746E-03. 15. -1.06581E-02. σ(2,3,3). 3.87553E-03. 16. -1.06581E-02. σ(3,3,2). 2.77320E-03. 17. -1.06581E-02. σ(3,2,3). 6.01255E-04. 18. 2.36707E-03. σ(3,3,3). 2.14418E-10. 19. -6.72771E-10. σ(1,2,1). 1.29506E-10. 20. -4.57480E-10. σ(2,1,1). 5.71710E-11. 21. -1.86121E-10. σ(1,3,1). 3.21477E-11. 22. -1.26562E-10. σ(3,1,1). 3.79868E-14. 23. 1.47633E-13. σ(3,2,2). 6.22127E-15. 24. -2.43817E-14. σ(2,3,2). 6.61145E-16. 25. 3.49377E-16. σ(1,1,3). 6.37790E-16. 26. -2.76688E-16. σ(1,1,2). 6.21132E-16. 27. -1.07912E-18. σ(2,2,3). 6.21214E-16. - 30 -.

(1)3階直交テンソル積展開の性質と その計算法の改良に関する研究

(1)3階直交テンソル積展開の性質とその計算法の改良に関する研究