(57)【要約】【課題】特別な形をした行列たちのクロネッカー積を含む行列演算を高速化する。【解決手段】正方行列Ａ［ｉ］＝γ［ｉ］Ｕ（ｍ［ｉ］）＋ρ［ｉ］Ｅ（ｍ［ｉ］）とサイズがｍのベクトル（ｘ１，…，ｘｍ）との行列演算を行う。そのために、ｓ＝１，…，ｍ／ｍ［ｉ］についてｑ［ｉ，ｓ］＝γ［ｉ］Σｊ∈｛１，…，ｍ［ｉ］｝  ｘ（ｓ−１）ｍ［ｉ］+ｊを得る処理と、得られたｑ［ｉ，ｓ］を用いてｔ＝１，…，ｍ［ｉ］について値ｖ［ｉ，ｓ，ｔ］＝ｑ［ｉ，ｓ］＋ρ［ｉ］ｘ（ｓ−１）ｍ［ｉ］+ｔを得る処理とを行う。ただし、ｍ［ｉ］は正の整数であり、γ［ｉ］およびρ［ｉ］は定数であり、Ｕ（ｍ［ｉ］）は全要素１でサイズｍ［ｉ］×ｍ［ｉ］の正方行列であり、Ｅ（ｍ［ｉ］）はサイズｍ［ｉ］×ｍ［ｉ］の単位行列であり、ｍがｍ［ｉ］の倍数である。【選択図】図１
2015135536.tif 000002
【特許請求の範囲】
【請求項１】
正方行列Ａ［ｉ］＝γ［ｉ］Ｕ（ｍ［ｉ］）＋ρ［ｉ］Ｅ（ｍ［ｉ］）とサイズがｍのベクトル（ｘ１，…，ｘｍ）との行列演算を行う行列演算装置であって、
ｍ［ｉ］が正の整数であり、γ［ｉ］およびρ［ｉ］が定数であり、Ｕ（ｍ［ｉ］）が全要素１でサイズｍ［ｉ］×ｍ［ｉ］の正方行列であり、Ｅ（ｍ［ｉ］）がサイズｍ［ｉ］×ｍ［ｉ］の単位行列であり、ｍがｍ［ｉ］の倍数であり、
ｓ＝１，…，ｍ／ｍ［ｉ］についてｑ［ｉ，ｓ］＝γ［ｉ］Σｊ∈｛１，…，ｍ［ｉ］｝  ｘ（ｓ−１）ｍ［ｉ］+ｊを得る第１演算部と、
前記第１演算部で得られたｑ［ｉ，ｓ］を用いてｔ＝１，…，ｍ［ｉ］について値ｖ［ｉ，ｓ，ｔ］＝ｑ［ｉ，ｓ］＋ρ［ｉ］ｘ（ｓ−１）ｍ［ｉ］+ｔを得る第２演算部と、
を有する行列演算装置。
【請求項２】
請求項１の行列演算装置であって、
ｎが２以上の整数であり、ｍ＝Πｉ∈｛１，…，ｎ｝  ｍ［ｉ］であり、ｐ∈｛２，…，ｎ｝であり、
前記第１演算部の処理と前記第２演算部の処理とを行うループ処理がｎ回繰り返えされ、
１回目の前記ループ処理での前記サイズがｍのベクトル（ｘ１，…，ｘｍ）は入力ベクトルであり、
ｐ回目の前記ループ処理での前記サイズがｍのベクトル（ｘ１，…，ｘｍ）は、ｐ−１回目の前記ループ処理で得られた前記値ｖ［ｉ，ｓ，ｔ］を｛（ｔ−１）ｍ／ｍ［ｉ］＋ｓ｝番目の要素としたサイズがｍのベクトルである行列演算装置。
【請求項３】
請求項２の行列演算装置であって、
前記正方行列Ａ［１］，…，Ａ［ｎ］のクロネッカー積がＡであり、ベクトルν＝（ν１，…，νｍ）とベクトルβ＝（β１，…，βｍ）に対するν＊βが（ν１β１，…，νｍβｍ）であり、ν／βが（ν１／β１，…，νｍ／βｍ）であり、νｔがνの転置であり、撹乱後のテーブルから得られたクロス集計を表すサイズがｍのベクトルｙからクロス集計を再構築するための演算がχ（ｒ）＝χ（ｒ−１）＊（Ａ（ｙ／χ（ｒ−１）Ａ）ｔ）ｔを含み、χ（０）＝ｙであり、
少なくとも、前記入力ベクトルがχ（０）であって、かつ、前記ループ処理をｎ回繰り返して得られた前記ｖ［ｉ，ｓ，ｔ］を｛（ｔ−１）ｍ／ｍ［ｉ］＋ｓ｝番目の要素としたサイズがｍのベクトルがχｋＡであるか、または、前記入力ベクトル（ｘ１，…，ｘｍ）が（ｙ／χ（ｒ−１）Ａ）ｔであって、かつ、前記ループ処理をｎ回繰り返して得られた前記ｖ［ｉ，ｓ，ｔ］を｛（ｔ−１）ｍ／ｍ［ｉ］＋ｓ｝番目の要素としたサイズがｍのベクトルがＡ（ｙ／χ（ｒ−１）Ａ）ｔである、行列演算装置。
【請求項４】
請求項１から３の何れかの行列演算装置であって、
前記正方行列Ａ［ｉ］は遷移確率行列であり、前記ρ［ｉ］は０以上１以下の維持確率であり、γ［ｉ］＝（１−ρ［ｉ］）／ｍ［ｉ］である、行列演算装置。
【請求項５】
請求項１から４の何れかの行列演算装置であって、
或るｓについて前記第１演算部の処理を実行してから第２演算部の処理を実行する処理を、ｓ＝１，…，ｍ／ｍ［ｉ］についてそれぞれ実行する、行列演算装置。
【請求項６】
正方行列Ａ［ｉ］＝γ［ｉ］Ｕ（ｍ［ｉ］）＋ρ［ｉ］Ｅ（ｍ［ｉ］）とサイズがｍのベクトル（ｘ１，…，ｘｍ）との行列演算を行う行列演算方法であって、
ｍ［ｉ］が正の整数であり、γ［ｉ］およびρ［ｉ］が定数であり、Ｕ（ｍ［ｉ］）が全要素１でサイズｍ［ｉ］×ｍ［ｉ］の正方行列であり、Ｅ（ｍ［ｉ］）がサイズｍ［ｉ］×ｍ［ｉ］の単位行列であり、ｍがｍ［ｉ］の倍数であり、
第１演算部で、ｓ＝１，…，ｍ／ｍ［ｉ］についてｑ［ｉ，ｓ］＝γ［ｉ］Σｊ∈｛１，…，ｍ［ｉ］｝  ｘ（ｓ−１）ｍ［ｉ］+ｊを得る第１演算ステップと、
第２演算部で、前記第１演算ステップで得られたｑ［ｉ，ｓ］を用いてｔ＝１，…，ｍ［ｉ］について値ｖ［ｉ，ｓ，ｔ］＝ｑ［ｉ，ｓ］＋ρ［ｉ］ｘ（ｓ−１）ｍ［ｉ］+ｔを得る第２演算ステップと、
を有する行列演算方法。
【請求項７】
請求項１から５の何れかの行列演算装置としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、行列演算を行う技術に関し、特に、特別な形をした行列たちのクロネッカー積を含む行列演算を行うための技術に関する。
【背景技術】
【０００２】
行列のクロネッカー積とベクトルとの積を高速に計算する方法として非特許文献１に開示されたものがある。この方法では、クロネッカー積を直接計算せずに、個々の行列のベクトルへのある作用を定義して計算を行う。このような方法は、非特許文献２等に開示された撹乱再構築における再構築等に利用される。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】永井彰，五十嵐大，濱田浩気，松林達史，“クロネッカー積を含む行列積演算の最適化による効率的なプライバシー保護データ公開技術，”暗号と情報セキュリティシンポジウム2010 (SCIS2010).
【非特許文献２】五十嵐大，千田浩司，高橋克巳，“多値属性に適用可能な効率的プライバシー保護クロス集計，”コンピューターセキュリティシンポジウム2008 (CSS2008).
【発明の概要】
【発明が解決しようとする課題】
【０００４】
非特許文献２等に開示された再構築法では、特別な形をした行列たちのクロネッカー積を含む行列演算が必要になり、この演算が支配的な計算部分となる。
【０００５】
本発明の課題は、特別な形をした行列たちのクロネッカー積を含む行列演算を高速化することである。
【課題を解決するための手段】
【０００６】
正方行列Ａ［ｉ］＝γ［ｉ］Ｕ（ｍ［ｉ］）＋ρ［ｉ］Ｅ（ｍ［ｉ］）とサイズがｍのベクトル（ｘ１，…，ｘｍ）との行列演算を行う。そのために、ｓ＝１，…，ｍ／ｍ［ｉ］についてｑ［ｉ，ｓ］＝γ［ｉ］Σｊ∈｛１，…，ｍ［ｉ］｝  ｘ（ｓ−１）ｍ［ｉ］+ｊを得る処理と、得られたｑ［ｉ，ｓ］を用いてｔ＝１，…，ｍ［ｉ］について値ｖ［ｉ，ｓ，ｔ］＝ｑ［ｉ，ｓ］＋ρ［ｉ］ｘ（ｓ−１）ｍ［ｉ］+ｔを得る処理とを行う。ただし、ｍ［ｉ］は正の整数であり、γ［ｉ］およびρ［ｉ］は定数であり、Ｕ（ｍ［ｉ］）は全要素１でサイズｍ［ｉ］×ｍ［ｉ］の正方行列であり、Ｅ（ｍ［ｉ］）はサイズｍ［ｉ］×ｍ［ｉ］の単位行列であり、ｍがｍ［ｉ］の倍数である。
【発明の効果】
【０００７】
本発明では、特別な形をした行列たちのクロネッカー積を含む行列演算を高速化することができる。
【図面の簡単な説明】
【０００８】
【図１】図１Ａは、実施形態の再構築装置のブロック図であり、図１Ｂは、実施形態の行列演算装置のブロック図である。
【図２】図２は、実施形態の再構築演算を説明するためのフロー図である。
【図３】図３は、実施形態の｛ｘＡ演算｝を説明するためのフロー図である。
【図４】図４は、実施形態の｛Ａ［ｉ］．ｘ演算｝を説明するためのフロー図である。
【発明を実施するための形態】
【０００９】
以下、本発明の実施形態を説明する。
［撹乱再構築］
まず、プライバシー保護データ公開技術を実現させる技術の一つである、撹乱再構築の概要を述べる。情報提供者からのデータの集合であるテーブルに対して、各情報提供者からのデータをレコードと呼ぶ。各レコードはあらかじめ定められた項目に対する値から成り立っており、この項目を属性といい、属性の値（属性値）の取りうる範囲を属性の値域という。以下にテーブルの例を示す。以下の例では、「性別」「年代」が属性であり、それらの属性値からなる（女性，２０代）（男性，２０代）（女性，１０代）等がレコードである。以下の例では、属性「性別」の値域が「男性」「女性」からなり、属性「年代」の値域が「１０代」「２０代」からなる。
【表１】
2015135536.tif 000003
【００１０】
クロス集計とは、着目したすべての属性に関して属性値が等しいレコードの個数を集計したものである。クロス集計は大きさが各属性の値域の要素数の積からなる行ベクトルとして考えられる。表１の例の場合、クロス集計は、例えば、４個の要素からなる行ベクトル（｛男性，１０代｝，｛男性，２０代｝，｛女性，１０代｝，｛女性，２０代｝）＝（６１，４５，４３，５０）となる。
【００１１】
撹乱再構築は、情報提供者が情報を提供する際にプライバシー保護のために自らテーブルに対して行う撹乱と、マイニング者が撹乱されたテーブルからクロス集計のみの統計量を取り出す再構築とからなる。このようなことが可能なのは、撹乱がなされたテーブルで個々のデータの情報量が低下しても、統計量としては撹乱アルゴリズムごとに期待される値に収束するからである。
【００１２】
撹乱では各情報提供者がデータを確率的にランダムに変化させる。テーブルにおける各属性の属性値を一定の維持確率のもとで維持させ、維持しなければ属性値から一様にランダムに遷移させる。このように撹乱処理では属性ごとに維持確率を設定してレコード内の属性値を確率的に置換してプライバシー保護を図る。維持確率は、次のようにサイズが属性ｉの値域の要素数ｍ［ｉ］（ただし、ｍ［ｉ］は正の整数）の正方行列である遷移確率行列Ａ［ｉ］で表現できる。ここで、Ｕ（ｍ［ｉ］）とＥ（ｍ［ｉ］）は、それぞれ、サイズがｍ［ｉ］×ｍ［ｉ］の要素がすべて１の行列と単位行列であり、ρ［ｉ］（ただし、０≦ρ［ｉ］≦１）はその要素をどの程度維持するかを表す定数（維持確率）である。また、γ［ｉ］＝（１−ρ［ｉ］）／ｍ［ｉ］である。
【数１】
2015135536.tif 000004
【００１３】
先の表１の例の場合、性別の属性の維持確率を０．２とし、年代の属性の維持確率を０．６とすると、遷移確率行列Ａ［性別］およびＡ［年代］は以下のようになる。
【数２】
2015135536.tif 000005
例えば、Ａ［性別］の１行２列の要素（Ａ［性別］）１２＝０．４は男性から女性にデータが撹乱される確率が０．４であることを意味する。
【００１４】
複数の属性に対して遷移確率行列｛Ａ［ｉ］｝ｉ＝１，...，ｎ（ただし、ｎは２以上の整数）が与えられた際、レコードが撹乱後にあるレコードに遷移する確率を要素とする行列Ａは、以下のような遷移確率行列たちのクロネッカー積で表現できる。
【数３】
2015135536.tif 000006
【００１５】
＜クロネッカー積の定義＞
サイズがｍＢ×ｍＢである正方行列Ｂ＝（ｂφβ）とサイズがｍＣ×ｍＣである正方行列Ｃ＝（ｃφβ）に対して、クロネッカー積を次のように定義する。
【数４】
2015135536.tif 000007
ここで、ｂφβＣとは、行列Ｃのｂφβによるスカラー倍である。なお、このクロネッカー積のサイズはｍＢｍＣ×ｍＢｍＣである。なお、以降では、記載表記の制約上、正方行列ＢとＣとのクロネッカー積をＢ（×）Ｃと表記する場合がある。
【００１６】
先の表１の例の遷移確率行列Ａ［性別］とＡ［年代］とのクロネッカー積Ａ［性別］（×）Ａ［年代］は以下のようになる。
【数５】
2015135536.tif 000008
例えば、要素（Ａ［性別］（×）Ａ［年代］）３４＝０.１２は、ペア｛女性，１０代｝から｛女性，２０代｝にデータが撹乱される確率が０．１２であることを意味する。
【００１７】
再構築では、ベイズ推定法に基づき、撹乱で用いた遷移行列のクロネッカー積Ａ＝（×）ｉ∈｛１，...，ｎ｝  Ａ［ｉ］と撹乱後のテーブルから得られるクロス集計ｙ＝（ｙ１，…，ｙｍ）（サイズがｍのベクトルであり、ｍはｍ［ｉ］の倍数）とから、元のテーブルのクロス集計の推定を行う。推定のアルゴリズムは下記Ａｌｇｏｒｉｔｈｍ  １のようになる。ここで、ｒは反復回数に対応する整数であり、ベクトルν＝（ν１，…，νｍ）とベクトルβ＝（β１，…，βｍ）に対して、ν＊βが（ν１β１，…，νｎβｎ）であり、ν／βが（ν１／β１，…，νｎ／βｎ）である。ただし、νｉβｉはνｉとβｉとの積を表し、νｉ／βｉはβｉによるνｉの商を表す。また、νｔはνの転置であり、｜ν｜Ｌ１はベクトルνのＬ１距離、すなわちベクトルνの要素の絶対値の総和であり、εは予め定めた定数（収束半径）である。また、ν：＝βはνをβと定義すること（νにβを代入すること）を意味する。
＜Ａｌｇｏｒｉｔｈｍ  １＞
Ｉｎｐｕｔ：Ａ，ｙ＝（ｙ１，…，ｙｍ）
Ｏｕｔｐｕｔ：χ（ｒ）＝（χ１（ｒ），…，χｍ（ｒ））
χ（０）：＝（ｙ１，…，ｙｍ）
ｒ：＝０
ｗｈｉｌｅ  ｜χ（ｒ）−χ（ｒ−１）｜Ｌ１＞ε  ｄｏ
χ(ｒ)＝χ（ｒ−１）＊（Ａ（ｙ／χ（ｒ−１）Ａ）ｔ）ｔ
ｒ：＝ｒ＋１
ｅｎｄ  ｗｈｉｌｅ
ｒｅｔｕｒｎ  χ（ｒ）
【００１８】
ここでＡｌｇｏｒｉｔｈｍ  １を実装するにあたって支配的な計算は、演算χ(ｒ)＝χ（ｒ−１）＊（Ａ（ｙ／χ（ｒ−１）Ａ）ｔ）ｔに含まれるクロネッカー積Ａとベクトルｘ（ただし、ｘ＝χ（ｒ）およびｘ＝（ｙ／χ（ｒ−１）Ａ）ｔ）との積ｘＡ＝ｘ（（×）ｉ∈｛１，...，ｎ｝  Ａ［ｉ］）である。
【００１９】
［再構築アルゴリズムの高速化］
非特許文献１では、ｘＡ＝ｘ（（×）ｉ∈｛１，...，ｎ｝  Ａ［ｉ］）を計算するにあたり、以下のようなｘに対する各遷移行列Ａ［ｉ］の作用Ａ［ｉ］．ｘを与えている。まず、サイズがｍ＝Πｉ∈｛１，…，ｎ｝  ｍ［ｉ］のベクトルｘ＝（ｘ１，…，ｘｍ）に対して、ｍ／ｍ［ｉ］個のサイズがｍ［ｉ］のベクトル｛ｚｓ｝ｓ＝１，…，ｍ／ｍ［ｉ］を以下のように定める。
ｚｓ＝（ｘ（ｓ−１）ｍ[ｉ]＋１，…，ｘｓｍ[ｉ]）  …（１）
ただし、Πｉ∈｛１，…，ｎ｝  ｍ［ｉ］は
【数６】
2015135536.tif 000009
を表す。
【００２０】
また、行列Ａ［ｉ］に対して、ｍ［ｉ］個のサイズｍ［ｉ］のベクトル｛αｔ（ｉ）｝ｔ＝１，…，ｍ［ｉ］を以下のように定める。
αｔ（ｉ）＝（（Ａ［ｉ］）ｔ１，（Ａ［ｉ］）ｔ２，…，（Ａ［ｉ］）ｔｍ［ｉ］）  …（２）
ただし、（Ａ［ｉ］）ηωは行列Ａ［ｉ］のη行ω列の要素である。すなわち、αｔ（ｉ）は、行列Ａ［ｉ］のｔ行目の行ベクトルである。
【００２１】
これらの｛ｚｓ｝と｛αｔ（ｉ）｝とを用いて作用Ａ［ｉ］．ｘを以下のように定める。
Ａ［ｉ］．ｘ
＝（ｚ１・α１（ｉ），ｚ２・α１（ｉ），…，ｚｍ／ｍ［ｉ］・α１（ｉ），
ｚ１・α２（ｉ），ｚ２・α２（ｉ），…，ｚｍ／ｍ［ｉ］・α２（ｉ），
…，
ｚ１・αｍ［ｉ］（ｉ），ｚ２・αｍ［ｉ］（ｉ），…，ｚｍ／ｍ［ｉ］・αｍ［ｉ］（ｉ））  …（３）
ここで、ν・βはベクトルνとベクトルβとの内積を表す。
【００２２】
この作用Ａ［ｉ］．ｘにより、
ｘ（（×）ｉ∈｛１，...，ｎ｝  Ａ［ｉ］）＝Ａ［１］.（Ａ［２］.（...Ａ［ｎ］.ｘ）...））
となり、各作用を再帰的に計算することによって、クロネッカー積の計算を直接行うよりも計算量を軽減できる。実際、クロネッカー積の計算を直接計算を行う場合の計算量はＯ（ｍ２）であるが、各作用を再帰的に計算させる場合の計算量は、サイズがｍ［ｉ］のベクトルの内積をｍ回行うため、Ｏ（ｍΣｉ∈｛１，...，ｎ｝  ｍ［ｉ］）となる。ただし、Σｉ∈｛１，...，ｎ｝  ｍ［ｉ］は
【数７】
2015135536.tif 000010
を表す。
【００２３】
作用Ａ［ｉ］．ｘのアルゴリズムを形式的に表すと以下のＡｌｇｏｒｉｔｈｍ  ２のようになる。
＜Ａｌｇｏｒｉｔｈｍ  ２＞
Ｉｎｐｕｔ：Ａ［ｉ］，ｘ＝（ｘ１，…，ｘｍ）
Ｏｕｔｐｕｔ：Ａ［ｉ］．ｘ
｛ｚｓ：＝（ｘ（ｓ−１）ｍ[ｉ]＋１，…，ｘｓｍ[ｉ]）｝ｓ＝１，…，ｍ［ｉ］
｛αｔ（ｉ）＝（ａｔ１（ｉ），ａｔ２（ｉ），…，ａｔｍ［ｉ］（ｉ））｝ｔ＝１，…，ｍ／ｍ［ｉ］
ｖ：＝（ｖ１，…，ｖｍ）
ｆｏｒ  ｔ：＝１，…，ｍ／ｍ［ｉ］  ｄｏ
ｆｏｒ  ｓ：＝１，…，ｍ［ｉ］  ｄｏ
ｖ（ｔ−１）ｍ／ｍ［ｉ］＋ｓ：＝ｚｓ・αｔ（ｉ）
ｅｎｄ  ｆｏｒ
ｅｎｄ  ｆｏｒ
ｒｅｔｕｒｎ  ｖ
【００２４】
このＡｌｇｏｒｉｈｔｍ  ２を用いてクロネッカー積とベクトルとの積を計算させるアルゴリズムは下記のＡｌｇｏｒｉｔｈｍ３である。
＜Ａｌｇｏｒｉｔｈｍ  ３＞
Ｉｎｐｕｔ：｛Ａ［ｉ］｝ｉ＝１，...，ｎ，ｘ＝（ｘ１，…，ｘｍ）
Ｏｕｔｐｕｔ：ｘ（（×）ｉ∈｛１，...，ｎ｝  Ａ［ｉ］）
ｖ：＝ｘ
ｆｏｒ  i：＝ｎ，...，１  ｄｏ
ｖ：＝Ａ［ｉ］.ｖ
ｅｎｄ  ｆｏｒ
ｒｅｔｕｒｎ  ｖ
【００２５】
［作用Ａ［ｉ］．ｘの高速化］
本形態では、Ａ［ｉ］が全要素１でサイズｍ［ｉ］×ｍ［ｉ］の正方行列Ｕ（ｍ［ｉ］）と、サイズｍ［ｉ］×ｍ［ｉ］の単位行列Ｅ（ｍ［ｉ］）との一次結合（Ａ［ｉ］＝γ［ｉ］Ｕ（ｍ［ｉ］）＋ρ［ｉ］Ｅ（ｍ［ｉ］））であることを利用し、作用Ａ［ｉ］.ｘの計算高速化を図る。
【００２６】
ここでＡ［ｉ］＝γ［ｉ］Ｕ（ｍ［ｉ］）＋ρ［ｉ］Ｅ（ｍ［ｉ］）という形なので、前述の式（２）のベクトルαｔ（ｉ）は、ｔ番目の要素だけがγ［ｉ］＋ρ［ｉ］でその他の要素がγ［ｉ］のサイズｍ［ｉ］のベクトルとなる。
【数８】
2015135536.tif 000011
【００２７】
前述の式（３）の作用Ａ［ｉ］．ｘは、ベクトルｚｓとベクトルαｔ（ｉ）との内積ｚｓ・αｔ（ｉ）を１≦ｓ≦ｍ／ｍ［ｉ］と１≦ｔ≦ｍ［ｉ］に亙らせて計算すればよいが、本形態では、まずｓを固定し、ｔを亙らせてｚｓ・α１（ｉ），ｚｓ・α２（ｉ），…，ｚｓ・αｍ［ｉ］（ｉ）を計算させる。式（１）のようにｚｓ＝（ｘ（ｓ−１）ｍ[ｉ]＋１，…，ｘｓｍ[ｉ]）であり、ベクトルαｔ（ｉ）はｔ番目の要素だけがγ［ｉ］＋ρ［ｉ］でその他の要素がγ［ｉ］であるから、
【数９】
2015135536.tif 000012
となる。ここで、式（４）のγ［ｉ］Σｊ∈｛１，…，ｍ［ｉ］｝  ｘ（ｓ−１）ｍ［ｉ］+ｊの項はｔとは独立であり、ｓにしか依存しない。故に作用Ａ［ｉ］．ｘの計算量はＯ（ｍ＋ｍ［ｉ］）となり、ｘ（（×）ｉ∈｛１，...，ｎ｝  Ａ［ｉ］）＝Ａ［１］.（Ａ［２］.（...Ａ［ｎ］.ｘ）...））全体の計算量はＯ（ｎｍ＋Σｉ∈｛１，...，ｎ｝ｍ［ｉ］）＝Ｏ（ｎｍ）となる。
【００２８】
作用Ａ［ｉ］．ｘのアルゴリズムは以下のＡｌｇｏｒｉｔｈｍ  ４のようになる。
＜Ａｌｇｏｒｉｔｈｍ  ４＞
Ｉｎｐｕｔ：Ａ［ｉ］＝γ［ｉ］Ｕ（ｍ［ｉ］）＋ρ［ｉ］Ｅ（ｍ［ｉ］），ｘ＝（ｘ１，…，ｘｍ）
Ｏｕｔｐｕｔ：Ａ［ｉ］．ｘ
ｖ＝（ｖ１，…，ｖｍ）
ｆｏｒ  ｓ：＝１，...，ｍ／ｍ［ｉ］  ｄｏ
ｑ［ｉ，ｓ］＝γ［ｉ］Σｊ∈｛１，…，ｍ［ｉ］｝  ｘ（ｓ−１）ｍ［ｉ］+ｊ
ｆｏｒ  ｔ：＝１，...，ｍ［ｉ］  ｄｏ
ｖ（ｔ−１）ｍ／ｍ［ｉ］＋ｓ：＝ｑ［ｉ，ｓ］＋ρ［ｉ］ｘ（ｓ−１）ｍ［ｉ］+ｔ
ｅｎｄ  ｆｏｒ
ｅｎｄ  ｆｏｒ
ｒｅｔｕｒｎ  ｖ
【００２９】
次に、本発明の実施形態を図面を参照して詳細に説明する。
［第１実施形態］
＜構成＞
図１Ａに例示するように、本形態の再構築装置１は、行列演算装置１１、初期設定部１２、再構築部１３、制御部１４、およびメモリ１５を有する。図２Ｂに例示するように、本形態の行列演算装置１１は、第１演算部１１１、第２演算部１１２、制御部１１３、係数抽出部１１４、およびメモリ１１５を有する。再構築装置１は、ＣＰＵ（central processing unit）やＲＡＭ（random-access memory）等を有する汎用または専用のコンピュータに所定のプログラムが読み込まれて構成されたものである。再構築装置１は制御部１４の制御のもと、行列演算装置１１は制御部１１３の制御のもとで各処理を実行する。再構築装置１や行列演算装置１１に入力されたデータや各部で得られたデータは、逐一、メモリ１５やメモリ１１５に格納され、必要に応じて読み出される。
【００３０】
＜再構築＞
再構築装置１は、撹乱後のテーブルから得られるクロス集計ｙ＝（ｙ１，…，ｙｍ）と、その撹乱に用いられた遷移確率行列｛Ａ［ｉ］｝ｉ＝１，...，ｎ（ただし、ｎは２以上の整数）とを入力とし、ベイス推定法に基づいて元のクロス集計を再構築し、元のクロス集計の推定であるベクトルχ（ｒ）＝（χ１（ｒ），…，χｍ（ｒ））を出力する。以下、図２を用いて再構築装置１の処理を説明する。
【００３１】
初期設定部１２（図１Ａ）は、クロス集計ｙ＝（ｙ１，…，ｙｍ）を入力とし、χ（０）をｙ＝（ｙ１，…，ｙｍ）とし（χ（０）：＝（ｙ１，…，ｙｍ））、サイズがｍのベクトルχ（０）を出力する（ステップＳ１２）。制御部１４は、ｒを１にする（ｒ：＝１）（ステップｓ１４ａ）。再構築部１３は、遷移確率行列｛Ａ［ｉ］｝ｉ＝１，...，ｎおよびベクトルχ（０）を入力とし、行列演算装置１１を呼び出して図３の｛ｘＡ演算｝を実行させて、ベクトルｙからクロス集計を再構築するための演算χ（ｒ）＝χ（ｒ−１）＊（Ａ（ｙ／χ（ｒ−１）Ａ）ｔ）ｔを行う。すなわち再構築部１３は、χ（ｒ−１）Ａの演算に｛ｘＡ演算｝を用いるか、Ａ（ｙ／χ（ｒ−１）Ａ）ｔの演算に｛ｘＡ演算｝を用いるか、χ（ｒ−１）Ａの演算およびＡ（ｙ／χ（ｒ−１）Ａ）ｔの演算の両方に｛ｘＡ演算｝を用いる。この詳細は後述する（ステップＳ１３）。制御部１４は、｜χ（ｒ）−χ（ｒ−１）｜Ｌ１＞εであるかを判定し（ステップＳ１４ｂ）、｜χ（ｒ）−χ（ｒ−１）｜Ｌ１＞εであればｒ＋１を新たなｒとし（ｒ：＝ｒ＋１）（ステップＳ１４ｃ）、ステップＳ１３に戻す。一方、｜χ（ｒ）−χ（ｒ−１）｜Ｌ１＞εでなければ、再構築部１３が、ベクトルχ（ｒ）＝（χ１（ｒ），…，χｍ（ｒ））を出力する（ステップＳ１４ｄ）。
【００３２】
＜ｘＡ演算＞
図３を用いて、行列演算装置１１（図１Ｂ）の｛ｘＡ演算｝を説明する。｛ｘＡ演算｝では、行列演算装置１１が、ベクトルｘ＝（ｘ１，…，ｘｍ）（入力ベクトル）と遷移確率行列｛Ａ［ｉ］｝ｉ＝１，...，ｎとを入力とし、遷移確率行列｛Ａ［ｉ］｝ｉ＝１，...，ｎのクロネッカー積Ａとベクトルｘ＝（ｘ１，…，ｘｍ）との積ｘＡを計算し、その演算結果であるサイズがｍのベクトルｖ＝（ｖ１，…，ｖｍ）を出力する。なお、行列演算装置１１は、前述のχ（ｒ−１）Ａの演算を行う場合には、ｘ＝χ（ｒ−１）を入力ベクトルとしてｘＡ＝χ（ｒ−１）Ａを計算し、その計算結果を表すサイズがｍのベクトルｖ＝（ｖ１，…，ｖｍ）を出力する。前述のＡ（ｙ／χ（ｒ−１）Ａ）ｔの演算を行う場合には、行列演算装置１１は、ｘ＝（ｙ／χ（ｒ−１）Ａ）ｔを入力ベクトルとしてｘＡ＝Ａ（ｙ／χ（ｒ−１）Ａ）ｔを計算し、その計算結果を表すサイズがｍのベクトルｖ＝（ｖ１，…，ｖｍ）を出力する。以下では、ベクトルｘとクロネッカー積Ａとの積に一般化して｛ｘＡ演算｝の説明を行う。
【００３３】
まず、制御部１１３はｉをｎに設定する（ｉ：＝ｎ）（ステップＳ１１３ｈ）。次に、制御部１１３は、図４に例示する後述の｛Ａ［ｉ］．ｘ演算｝（ループ処理）を実行させる（ステップＳ１１３ｉ）。１回目（ｉ＝ｎ）の｛Ａ［ｉ］．ｘ演算｝でのベクトルｘ＝（ｘ１，…，ｘｍ）は、行列演算装置１１への入力ベクトル（例えば、ｘ＝χ（ｒ−１）やｘ＝（ｙ／χ（ｒ−１）Ａ）ｔ）であり、ｐ（ただし、ｐ∈｛２，…，ｎ｝）回目（ｉ＝ｎ−１，，…，１）の｛Ａ［ｉ］．ｘ演算｝でのベクトルｘ＝（ｘ１，…，ｘｍ）は、ｐ−１回目の｛Ａ［ｉ］．ｘ演算｝で得られた後述の値ｖ［ｉ，ｓ，ｔ］を｛（ｔ−１）ｍ／ｍ［ｉ］＋ｓ｝番目の要素ｘ（ｔ−１）ｍ／ｍ［ｉ］＋ｓとしたサイズがｍのベクトルである。次に、制御部１１３は、ｉ＝１であるかを判定し（ステップＳ１１３ｋ）、ｉ＝１でない場合には、ｉ＋１を新たなｉとし（ｉ：＝ｉ＋１）（ステップＳ１１３ｍ）、ステップＳ１１３ｉに戻す。ｉ＝１である場合には、このように｛Ａ［ｉ］．ｘ演算｝の処理をｎ回繰り返して得られたｖ［ｉ，ｓ，ｔ］を｛（ｔ−１）ｍ／ｍ［ｉ］＋ｓ｝番目の要素ｖ（ｔ−１）ｍ／ｍ［ｉ］＋ｓとしたサイズがｍのベクトルｖ＝（ｖ１，…，ｖｍ）を出力する（ステップＳ１１３ｎ）。
【００３４】
＜Ａ［ｉ］．ｘ演算＞
図４に例示するように、まず、係数抽出部１１４がＡ［ｉ］を入力とし、Ａ［ｉ］から前述のρ［ｉ］とγ［ｉ］とを抽出する（ステップＳ１１４）。制御部１１３はｓを１に設定する（ｓ：＝１）（ステップＳ１１３ａ）。次に、第１演算部１１１は、γ［ｉ］およびｘを入力としてｑ［ｉ，ｓ］＝γ［ｉ］Σｊ∈｛１，…，ｍ［ｉ］｝  ｘ（ｓ−１）ｍ［ｉ］+ｊを計算してｑ［ｉ，ｓ］を出力する（ステップＳ１１１）。次に、制御部１１３はｔを１に設定する（ｔ：＝１）（ステップＳ１１３ｂ）。第２演算部１１２は、ｑ［ｉ，ｓ］，ρ［ｉ］およびｘを入力とし、ｖ［ｉ，ｓ，ｔ］＝ｑ［ｉ，ｓ］＋ρ［ｉ］ｘ（ｓ−１）ｍ［ｉ］+ｔを計算する（ステップＳ１１２）。制御部１１３は、ｔ＝ｍ［ｉ］であるかを判定し（ステップＳ１１３ｃ）、ｔ＝ｍ［ｉ］でなければ、ｔ＋１を新たなｔとし（ｔ：＝ｔ＋１）（ステップＳ１１３ｄ）、ステップＳ１１２に戻す。ｔ＝ｍ［ｉ］であれば、制御部１１３は、ｓ＝ｍ／ｍ［ｉ］であるかを判定し（ステップＳ１１３ｅ）、ｓ＝ｍ／ｍ［ｉ］でなければ、ｓ＋１を新たなｓとし（ｓ：＝ｓ＋１）（ステップＳ１１３ｆ）、ステップＳ１１１に戻す。ｓ＝ｍ／ｍ［ｉ］であれば、第２演算部１１２は、ｖ［ｉ，ｓ，ｔ］を｛（ｔ−１）ｍ／ｍ［ｉ］＋ｓ｝番目（ただし、ｓ＝１，…，ｍ／ｍ［ｉ］，ｔ＝１，…，ｍ［ｉ］）の要素ｖ（ｔ−１）ｍ／ｍ［ｉ］＋ｓとしたサイズがｍのベクトルｖ＝（ｖ１，…，ｖｍ）を出力する（ステップＳ１１３ｇ）。
【００３５】
［特徴］
非特許文献１でのクロネッカー積を含む行列の演算の計算量はＯ（ｍΣｉ∈｛１，…，ｎ｝ｍ［ｉ］）であったが、本形態では計算量がＯ（ｎｍ）まで軽減され、値域の平均倍の高速化が行える。非特許文献２や参考文献（菊池亮，五十嵐大，千田浩司，濱田浩気，“データ分布依存処理によって高い有用性を実現する確率的ｋ−匿名化，”、暗号と情報セキュリティシンポジウム2013 (SICS2013).）に開示されている撹乱再構築法やデータ依存型撹乱再構築法等の再構築アルゴリズム、特に反復ベイズ法に対して本形態を適用しクロネッカー積を計算することで、高速に統計分析を行うことができる。例えば、行列のサイズが１００００×１００００の場合であれば、非特許文献１の計算方法に比べて行列のサイズの平均である１００００倍の高速化が図れる。実際に、非特許文献２を住所等の属性に適用する場合、行列のサイズはこの程度やそれ以上になるため、本形態の方式は再構築法において大きな効果を発揮する。また、高速に統計分析を行うことが可能となるため、統計データの範囲を大きくし、きめ細かな統計データを使用して撹乱再構築を行うことも可能となる。
【００３６】
［変形例等］
本発明は上述の実施の形態に限定されない。例えば、第１実施形態の｛Ａ［ｉ］．ｘ演算｝では、或るｓについて第１演算部１１１の処理（ステップＳ１１１）を実行してからｔ＝１，…，ｍ［ｉ］について第２演算部１１２の処理（ステップＳ１１２）を実行する処理を、ｓ＝１，…，ｍ／ｍ［ｉ］についてそれぞれ実行した。しかしながら、第１演算部１１１がｓ＝１，…，ｍ／ｍ［ｉ］についてｑ［ｉ，ｓ］＝γ［ｉ］Σｊ∈｛１，…，ｍ［ｉ］｝  ｘ（ｓ−１）ｍ［ｉ］+ｊを得、第２演算部１１２が第１演算部１１１で得られたｑ［ｉ，ｓ］を用いてｔ＝１，…，ｍ［ｉ］について値ｖ［ｉ，ｓ，ｔ］＝ｑ［ｉ，ｓ］＋ρ［ｉ］ｘ（ｓ−１）ｍ［ｉ］+ｔを得るのであれば、その他の順序で処理が成されてもよい。例えば、第１演算部１１１がすべてのｓ＝１，…，ｍ／ｍ［ｉ］についてｑ［ｉ，ｓ］＝γ［ｉ］Σｊ∈｛１，…，ｍ［ｉ］｝  ｘ（ｓ−１）ｍ［ｉ］+ｊを計算した後に、、第２演算部１１２が第１演算部１１１で得られたｑ［ｉ，ｓ］を用いてｔ＝１，…，ｍ［ｉ］について値ｖ［ｉ，ｓ，ｔ］＝ｑ［ｉ，ｓ］＋ρ［ｉ］ｘ（ｓ−１）ｍ［ｉ］+ｔを計算してもよい。その他、装置の処理能力あるいは必要に応じ、並列的にあるいは個別に各処理が実行されてもよい。
【００３７】
また、上記の実施形態では、再構築のために｛ｘＡ演算｝や｛Ａ［ｉ］．ｘ演算｝を行う例を示したが、その他のγ［ｉ］Ｕ（ｍ［ｉ］）＋ρ［ｉ］Ｅ（ｍ［ｉ］）の形をした
正方行列Ａ［ｉ］のクロネッカー積Ａ＝（×）ｉ∈｛１，...，ｎ｝  Ａ［ｉ］とサイズがｍのベクトルｘとの積のために｛ｘＡ演算｝や｛Ａ［ｉ］．ｘ演算｝が実行されてもよい。また、ｎ＝１であってもよい。
【００３８】
その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
【００３９】
上述の構成をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。このプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体の例は、非一時的な（non-transitory）記録媒体である。このような記録媒体の例は、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等である。
【００４０】
このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【００４１】
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。処理の実行時、このコンピュータは、自己の記録装置に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。
【００４２】
上記実施形態では、コンピュータ上で所定のプログラムを実行させて本装置の処理機能が実現されたが、これらの処理機能の少なくとも一部がハードウェアで実現されてもよい。
【符号の説明】
【００４３】
１  再構築装置
１１  行列演算装置
１１１  第１演算部
１１２  第２演算部
【図１】
2015135536.tif 000013
【図２】
2015135536.tif 000014
【図３】
2015135536.tif 000015
【図４】
2015135536.tif 000016