(57)【要約】【課題】より小さい観測率で（より少ない観測信号数で）より高い推定精度を実現できる求解技術を提供する。【解決手段】求めるべき変数を、前記変数の各要素の値を表現する状態変数と、前記変数の各要素の非ゼロ要素らしさを表現する構造変数との要素毎の積として表す。構造変数は、実数の値を取る制御変数の単調増加関数として０以上１以下の実数に制限されたものである。そして、コンピュータが、前記状態変数と前記構造変数又は前記状態変数との関数として定義される評価関数を用い、各反復において、前記評価関数の前記状態変数に関する情報に基づき、前記状態変数の各要素の値を更新する状態変数更新ステップと、前記評価関数の前記構造変数に関する情報に基づき、前記構造変数の各要素の値を更新する構造変数更新ステップと、を実行する。【選択図】図１
2015108994.tif 000002
【特許請求の範囲】
【請求項１】
劣決定系の連立線型方程式をコンピュータによる反復計算により解くための求解方法であって、
求めるべき変数を、前記変数の各要素の値を表現する状態変数と、前記変数の各要素の非ゼロ要素らしさを表現する構造変数との要素毎の積として表し、
前記構造変数は、実数の値を取る制御変数の単調増加関数として０以上１以下の実数に制限されたものであり、
前記状態変数と前記構造変数又は前記制御変数との関数として定義される評価関数を用い、
各反復において、
前記評価関数の前記状態変数に関する情報に基づき、前記状態変数の各要素の値を更新する状態変数更新ステップと、
前記評価関数の前記構造変数に関する情報に基づき、前記構造変数の各要素の値を更新する構造変数更新ステップと、をコンピュータが実行する
ことを特徴とする劣決定系線型方程式の求解方法。
【請求項２】
前記評価関数の前記状態変数に関する情報は、前記評価関数の前記状態変数に関する偏微分又は数値微分である
ことを特徴とする請求項１に記載の劣決定系線型方程式の求解方法。
【請求項３】
前記評価関数の前記構造変数に関する情報は、前記評価関数の前記制御変数に関する偏微分又は数値微分である
ことを特徴とする請求項１又は２に記載の劣決定系線型方程式の求解方法。
【請求項４】
前記状態変数更新ステップにおいて、前記状態変数の各要素の値を更新する処理が複数回実行される
ことを特徴とする請求項１〜３のうちいずれか１項に記載の劣決定系線型方程式の求解方法。
【請求項５】
前記構造変数更新ステップにおいて、前記構造変数の各要素の値を更新する処理が複数回実行される
ことを特徴とする請求項１〜４のうちいずれか１項に記載の劣決定系線型方程式の求解方法。
【請求項６】
各反復において、前記評価関数の値、前記連立線型方程式の残差ノルム、反復計算の回数、又は、これらのうちの二つ以上の組み合わせにより構成される収束条件を満足するか否かで、反復計算の終了を判断するステップ、をコンピュータが実行する
ことを特徴とする請求項１〜５のうちいずれか１項に記載の劣決定系線型方程式の求解方法。
【請求項７】
反復計算の終了後、各反復で得られた前記状態変数の値のうち、前記評価関数の値又は前記連立線型方程式の残差ノルムを最も小さくするものを、前記連立線型方程式の解に決定するステップ、をコンピュータが実行する
ことを特徴とする請求項１〜６のうちいずれか１項に記載の劣決定系線型方程式の求解方法。
【請求項８】
求めるべき変数をｘ、状態変数をｓ、構造変数をρ、制御変数をν、解くべき連立線型方程式をＡｘ＝ｂとしたとき、評価関数Ｌ２が下記式で定義される
ことを特徴とする請求項１〜７のうちいずれか１項に記載の劣決定系線型方程式の求解方
法。
【数１】
2015108994.tif 000018
【請求項９】
係数θの値を反復回数に応じて単調増加させる
ことを特徴とする請求項８に記載の劣決定系線型方程式の求解方法。
【請求項１０】
劣決定系の連立線型方程式を反復計算により解くための求解装置であって、
求めるべき変数を、前記変数の各要素の値を表現する状態変数と、前記変数の各要素の非ゼロ要素らしさを表現する構造変数との要素毎の積として表し、
前記構造変数は、実数の値を取る制御変数の単調増加関数として０以上１以下の実数に制限されたものであり、
前記状態変数と前記構造変数又は前記状態変数との関数として定義される評価関数を用いて、反復計算を行う計算手段を備え、
前記計算手段は、各反復において、
前記評価関数の前記状態変数に関する情報に基づき、前記状態変数の各要素の値を更新する状態変数更新ステップと、
前記評価関数の前記構造変数に関する情報に基づき、前記構造変数の各要素の値を更新する構造変数更新ステップと、を実行する
ことを特徴とする劣決定系線型方程式の求解装置。
【請求項１１】
観測系によって得られた観測データから被検体の内部の情報を表す画像データを再構成する被検体情報取得装置であって、
前記観測系の特性を表すデータ、及び、前記観測系によって得られた観測データを取得するデータ取得手段と、
請求項１０に記載の劣決定系線型方程式の求解装置と、を備え、
前記求解装置は、前記観測系の特性を表すデータを行列Ａ、前記観測系によって得られた観測データを定数ベクトルｂ、被検体の内部の情報を表す画像データを未知変数ｘとする、連立線型方程式Ａｘ＝ｂの解を計算する
ことを特徴とする被検体情報取得装置。
【請求項１２】
前記求解装置は、
予め与えられた変換行列Φを用いて連立線型方程式（ＡΦ−１）ｚ＝ｂにおける変数ｚ
を求解した後、
ｘ＝Φ−１ｚにより、変数ｘを求める
ことを特徴とする請求項１１に記載の被検体情報取得装置。
【請求項１３】
前記変換行列Φは、被検体の内部の情報を表す画像データのゼロ要素を増加させるフィルタである
ことを特徴とする請求項１２に記載の被検体情報取得装置。
【請求項１４】
請求項１〜９のうちいずれか１項に記載の劣決定系線型方程式の求解方法の各ステップをコンピュータに実行させることを特徴とするプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、連立線型方程式の求解方法及び求解装置に関し、特に、変数の数に対して方程式の数が少ない劣決定系線型方程式の求解技術に関する。
【背景技術】
【０００２】
Ｘ線ＣＴ（Computed Tomography）やＭＲＩ（Magnetic Resonance Imaging）等の医療
機器では、人体内部の３次元構造を、２次元平面への投影観測画像から再構成する。推定したい３次元構造を未知ベクトルｘ、Ｘ線投影のような観測過程を行列Ａ、観測信号をｂとおくと、この観測過程は次のように定式化できる：
【数１】
2015108994.tif 000003
Ａは係数行列、ｂは定数ベクトルあるいは右辺ベクトルとも呼ばれる。
【０００３】
式（１）の方程式を求解する場合、３次元構造の推定精度を向上させるために観測の回数を多くすることが一般的である。以後、推定すべき信号ｘの次元をＮ、観測信号ｂの次元をＭとし、Ｍ＜Ｎを仮定する。式（１）は、変数の次元Ｎに対して方程式の数Ｍが少ないので、劣決定系線型方程式と呼ばれる。
劣決定系線型方程式は解が無数に存在するので、一般に正則化と呼ばれる技術を適用することによって解を一意に決定する。
【０００４】
近年、推定すべき信号のスパース性を仮定する（つまり、ベクトルｘの多くの要素がゼロであると仮定する）ことによって、劣決定系線型方程式を精度よく求解する技術が提案された。この技術はスパース正則化と呼ばれ、例えば、画像のノイズ除去、医用画像の再構成等に有効な技術として注目されている。
スパース正則化のアルゴリズムは大きく２種類に分類することが出来る。「Greedy algorithm」と「Ｌ１正則化」である。Greedy algorithmで最も性能が良い方法が、非特許文献１に開示されたsubspace pursuit（以下ＳＰと略記する。）であり、Ｌ１正則化で最も性能が良い方法が非特許文献２に開示されたＦＩＳＴＡである。
【０００５】
ＦＩＳＴＡでは、以下の式（２）で定義される評価関数の最小値を探索することによって、式（１）を満足し、かつ非ゼロ要素数が少ない解ベクトルｘを得る：
【数２】
2015108994.tif 000004
評価関数の第１項は連立方程式の残差ベクトルのＬ２ノルムである。また、第２項は解ベクトルｘのＬ１ノルム（ベクトルｘの全要素の絶対値和）に定数λを乗じたものである。この評価関数Ｌ１は、解ベクトルｘが正解に近づくほど第１項が小さくなり、解ベクトルｘがスパースになるほど（ゼロ要素が増すほど）第２項が小さくなる、という性質をもつ。
【０００６】
ＦＩＳＴＡの処理の流れを、図３を用いて説明する。
ステップＳ３０１では、変数ｘ（０）とｙ（１）に初期値を設定し、反復回数ｋ＝１、ｔ（１）＝１とする。なお、カッコ書きの上付き添え字は反復回数（イテレーション回数）を表す。
ステップＳ３０２では、ステップ幅βを決定する。ステップ幅βには、∇Ｌ１：＝∂Ｌ
１／∂ｘのLipschitz定数の逆数を設定するが、次元が大きくなると計算時間が膨大にな
るので、近似値で代替してもよい。ただし、近似値はLipschitz定数の逆数を超えてはな
らない。
【０００７】
ステップＳ３０３では、次式で変数ｘ（ｋ）を計算する。
【数３】
2015108994.tif 000005
ここで、ｐｒｏｘβは下記式で定義される関数である。
【数４】
2015108994.tif 000006
【０００８】
ステップＳ３０４では、次式でｔ（ｋ＋１）を計算する。
【数５】
2015108994.tif 000007
ステップＳ３０５では次式でｙ（ｋ＋１）を計算する。
【数６】
2015108994.tif 000008
所定の収束条件を満足するまで、ステップＳ３０３〜Ｓ３０５を繰り返す。
【先行技術文献】
【非特許文献】
【０００９】
【非特許文献１】Wei Dai and Olgica Milenkovic : “Subspace Pursuit for Compressive Sensing Signal Reconstruction,” IEEE Transactions on Information Theory, VOL.55, NO.5, pp.2230-2249, MAY 2009.
【非特許文献２】Amir Beck and Marc Teboulle, “A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems,” SIAM Journal on Imaging Sciences, Vol.2, No.1, pp.183-202, 2009.
【発明の概要】
【発明が解決しようとする課題】
【００１０】
従来の２つの求解法ＳＰ及びＦＩＳＴＡの性能は、図４に示す完全再構成マップで見ることが出来る。推定する信号の数Ｎ、即ち変数の数Ｎは２５６である。図４は、正規分布に従って生成した行列Ａ及び正解１０ケースに対して各求解法を用い、９ケースが相対誤差５％以下であったときを完全再構成領域として作図した。図中横軸は疎性率（＝非ゼロ要素の数／推定信号の数）、縦軸は観測率（＝観測信号の数／推定信号の数）である。グラフ中の各線はそれぞれの解法に対するものであり、線の左上の領域が完全再構成領域である。線４０１は理論限界、線４０２はＳＰの限界、線４０３はＦＩＳＴＡの限界である。同図より、いずれの方法も理論限界に対して改善の余地があることがわかる。特に観測率が小さい領域では、疎性率の変化に対する観測率の変化割合が大きい。これは、信号中
の非ゼロ要素数が少し増加しただけでも、多くの追加観測信号を必要とすることを意味している。例えば、Ｘ線ＣＴやＭＲＩ等の医療機器の場合、被験者の身体的・心理的負担を減らすために、観測回数をより少なくし測定時間をより短くすることが望まれる。よって、必要な観測信号の数を削減することは、実用上、きわめて重要な課題である。
【００１１】
本発明は、上記課題に鑑みなされたものであって、より小さい観測率で（より少ない観測信号数で）より高い推定精度を実現できる求解技術を提供することを目的とする。
【課題を解決するための手段】
【００１２】
本発明の第一態様は、劣決定系の連立線型方程式をコンピュータによる反復計算により解くための求解方法であって、求めるべき変数を、前記変数の各要素の値を表現する状態変数と、前記変数の各要素の非ゼロ要素らしさを表現する構造変数との要素毎の積として表し、前記構造変数は、実数の値を取る制御変数の単調増加関数として０以上１以下の実数に制限されたものであり、前記状態変数と前記構造変数又は前記制御変数との関数として定義される評価関数を用い、各反復において、前記評価関数の前記状態変数に関する情報に基づき、前記状態変数の各要素の値を更新する状態変数更新ステップと、前記評価関数の前記構造変数に関する情報に基づき、前記構造変数の各要素の値を更新する構造変数更新ステップと、をコンピュータが実行することを特徴とする劣決定系線型方程式の求解方法である。
【００１３】
本発明の第二態様は、劣決定系の連立線型方程式を反復計算により解くための求解装置であって、求めるべき変数を、前記変数の各要素の値を表現する状態変数と、前記変数の各要素の非ゼロ要素らしさを表現する構造変数との要素毎の積として表し、前記構造変数は、実数の値を取る制御変数の単調増加関数として０以上１以下の実数に制限されたものであり、前記状態変数と前記構造変数又は前記状態変数との関数として定義される評価関数を用いて、反復計算を行う計算手段を備え、前記計算手段は、各反復において、前記評価関数の前記状態変数に関する情報に基づき、前記状態変数の各要素の値を更新する状態変数更新ステップと、前記評価関数の前記構造変数に関する情報に基づき、前記構造変数の各要素の値を更新する構造変数更新ステップと、を実行することを特徴とする劣決定系線型方程式の求解装置である。
【００１４】
本発明の第三態様は、観測系によって得られた観測データから被検体の内部の情報を表す画像データを再構成する被検体情報取得装置であって、前記観測系の特性を表すデータ、及び、前記観測系によって得られた観測データを取得するデータ取得手段と、本発明に係る劣決定系線型方程式の求解装置と、を備え、前記求解装置は、前記観測系の特性を表すデータを行列Ａ、前記観測系によって得られた観測データを定数ベクトルｂ、被検体の内部の情報を表す画像データを未知変数ｘとする、連立線型方程式Ａｘ＝ｂの解を計算することを特徴とする被検体情報取得装置である。
【００１５】
本発明の第四態様は、本発明に係る劣決定系線型方程式の求解方法の各ステップをコンピュータに実行させることを特徴とするプログラムである。
【発明の効果】
【００１６】
本発明によれば、より小さい観測率で（より少ない観測信号数で）より高い推定精度を実現できる求解技術を提供することができる。
【図面の簡単な説明】
【００１７】
【図１】本発明の実施形態に係る求解方法の処理を示すフローチャート。
【図２】本発明の実施形態に係る求解装置の構成を示すブロック図。
【図３】ＦＩＳＴＡの計算処理例を示すフローチャート。
【図４】従来技術と本発明の実施形態の求解方法の性能を示す完全再構成マップ。
【図５】本発明の求解方法を適用した被検体情報取得装置の構成例。
【発明を実施するための形態】
【００１８】
本発明は、劣決定系の連立線型方程式Ａｘ＝ｂをコンピュータによる反復計算により解くための求解技術に関し、前述したスパース正則化の改良アルゴリズムを提案するものである。
【００１９】
本方法の特徴の一つは、求めるべき変数ｘ（Ｎ個の要素をもつＮ次元ベクトルである。）のスパース性を評価するために、変数ｘを、ｘの各要素の値を表現する状態変数ｓと、ｘの各要素の非ゼロ要素らしさ（非ゼロ要素である確率）を表現する構造変数ρの要素毎の積として表現する点である。構造変数ρは状態変数ｓとは独立した変数である。構造変数ρも、状態変数ｓも、変数ｘと同じく、Ｎ個の要素をもつＮ次元ベクトルで表される。構造変数ρは確率的性質を持つので、０以上１以下の実数に制限される。そこで、微分等の取り扱いを容易にするために、実数の値を取る制御変数νを導入し、その単調増加関数として、例えば次式で構造変数ρを定義する。
【数７】
2015108994.tif 000009
式（８）中の係数θは構造変数ρの傾きを決定するパラメタであり、反復回数tの単調
増加関数とするとよい。係数θは例えば次式で与えられる（θ０は初期値である。）。
【数８】
2015108994.tif 000010
【００２０】
式（８）の変数ρにより、状態変数ｓの各要素の非ゼロ要素らしさを、０〜１に規格化された値（０がゼロ要素、１が非ゼロ要素）で表すことができる。そして、反復回数ｔの増加にしたがい係数θの値を次第に大きくしていくことで、ゼロ要素か非ゼロ要素かが不明な反復計算の初期段階では構造変数ρが０〜１の中間値を取りやすくし、反復計算が進むにつれ構造変数ρの０又は１への収束を促すことができる。
式（８）を用いることによって、構造変数ρは制御変数νで一意に決定することが出来るので、以下、構造変数ρの更新や微分と、制御変数νの更新や微分を、等価なものとして記載する場合がある。
【００２１】
本方法では、状態変数ｓと構造変数ρの関数として定義される評価関数Ｌ２（ｓ，ρ）を用いる（構造変数ρは制御変数νの関数であるから、この評価関数Ｌ２は、状態変数ｓと制御変数νの関数Ｌ２（ｓ，ν）で表すこともできる。）。具体的には、各反復において、（１）「Ｌ２のｓに関する情報」に基づき、状態変数ｓの各要素の値を更新する状態変数更新ステップと、（２）「Ｌ２のρに関する情報」に基づき、構造変数ρの各要素の値を更新する構造変数更新ステップを実行する。構造変数更新ステップでは、まず制御変数νの各要素の値を更新した後、式（８）によって構造変数ρの各要素の値を更新するとよい。
【００２２】
従来方法のＦＩＳＴＡでは、式（２）に示したように、評価関数Ｌ１の第２項に解ベクトルのＬ１ノルムを導入することで、解ベクトルのスパース性を評価している。しかし、この第２項のＬ１ノルムでは、解ベクトルの値そのものが小さい（例えば、信号レベルが小さい）のか、解ベクトルのゼロ要素が多いのかの区別ができない。それゆえ、評価関数Ｌ１の値が解ベクトル自体の値に影響を受け、推定精度や収束性の悪化を招くことがあっ
た。これに対し、本方法では、状態変数ｓとは独立した構造変数ρにより、状態変数ｓのスパース性を評価する。よって、状態変数ｓの値の大小が評価関数Ｌ２に与える影響をなくすことができ、従来方法に比べて、推定精度や収束性の向上を図ることが可能である。
【００２３】
「Ｌ２のｓに関する情報」としては、評価関数Ｌ２の状態変数ｓに関する偏微分（∂Ｌ２／∂ｓ）を用いるとよい。これにより、評価関数Ｌ２の値をより小さくするために、状態変数ｓを正／負のどちらの方向に変化させればよいか、を直接得ることができるからである。もし、評価関数Ｌ２の微分ができない場合又は難しい場合には、偏微分の代わりに数値微分を用いてもよい。数値微分とは、状態変数ｓの値を微小量Δｓだけずらしたときの評価関数の値Ｌ２（ｓ＋Δｓ，ρ）と、Ｌ２（ｓ，ρ）との差分から、近似的にＬ２のｓに関する微分を計算する操作である。数値微分以外にも、Ｌ２のｓに関する微分と等価な値や相関をもつ値であれば、「Ｌ２のｓに関する情報」として好ましく利用できる。
【００２４】
「Ｌ２のρに関する情報」としては、評価関数Ｌ２の制御変数νに関する偏微分（∂Ｌ２／∂ν）を用いるとよい。これにより、評価関数Ｌ２の値をより小さくするために、制御変数νを正／負のどちらの方向に変化させればよいか、を直接得ることができるからである。もし、評価関数Ｌ２の微分ができない場合又は難しい場合には、偏微分の代わりに数値微分を用いてもよい。数値微分とは、制御変数νの値を微小量Δνだけずらしたときの評価関数の値Ｌ２（ｓ，ν＋Δν）と、Ｌ２（ｓ，ν）との差分から、近似的にＬ２のνに関する微分を計算する操作である。数値微分以外にも、Ｌ２のνに関する微分と等価な値や相関をもつ値であれば、「Ｌ２のρに関する情報」として好ましく利用できる。
【００２５】
＜実施形態＞
以下、本発明に係る求解方法および求解装置の好ましい実施形態について、図面を参照して詳しく説明する。
【００２６】
（求解装置の構成）
図２は、本発明の実施形態に係る求解装置として機能するコンピュータの基本構成を示すブロック図である。このコンピュータは、大略、バス２０７を介して接続されたＣＰＵ２０１、表示装置２０２、入力装置２０３、ＲＡＭ２０４、ハードディスク装置２０５、ＮＩＣ２０６等から構成される。
【００２７】
ＣＰＵ（中央演算装置）２０１は、ＲＡＭ（主記憶装置）２０４に記憶されているプログラムやデータを用いて本装置全体の制御を行うと共に、後述する各処理を実行する。表示装置２０２は、液晶ディスプレイなどにより構成されており、ＣＰＵ２０１による処理結果（得られた解、再構成した画像）を文字や画像などでもって表示することができる。入力装置２０３は、キーボードやマウスなどの操作系の装置により構成されており、各種の指示をＣＰＵ２０１に対して入力することができる。ＲＡＭ２０４は，ハードディスク装置２０５からロードされたプログラムやデータを一時的に記憶する為のエリアを備えると共に、ＣＰＵ２０１が各種の処理を実行する際に必要なワークエリアを備える。ハードディスク装置２０５には、ＯＳ（オペレーティングシステム）や、コンピュータが行うべき後述の各処理をＣＰＵ２０１に実行させるためのプログラムやデータが保存されている。これらのプログラムがＲＡＭ２０４にロードされ、ＣＰＵ２０１によって実行されることで、本コンピュータが連立線型方程式を数値的に解くための求解装置として機能する。ＮＩＣ（ネットワークインターフェースコントローラ）２０６は、外部装置とのデータ通信を行う際のインターフェースとして機能するものである。本コンピュータはこのＮＩＣ２０６を介して外部装置とプログラムやデータの送受信を行う。
【００２８】
なお、ハードディスク装置２０５が保存しているプログラムやデータの代わりに、ＮＩＣ２０６を介して外部装置から受信したプログラムやデータをＣＰＵ２０１による処理対
象としても良い。すなわち、求解装置は、クライアント−サーバ、クラウドコンピューティング、グリッドコンピューティングなどのシステム構成を採ることもできる。あるいは、医療機器や画像分析機器などに搭載される組み込み用コンピュータに求解装置の機能を実装することもできる。あるいは、求解装置の機能の全部又は一部をＡＳＩＣやＦＰＧＡなどの回路で構成してもよい。
【００２９】
（求解方法）
次に、本実施形態に係る劣決定系線型方程式の求解方法について説明する。
【００３０】
本実施形態では、解くべき連立線型方程式Ａｘ＝ｂが与えられたときに、次式で定義される評価関数Ｌ２を用いる。
【数９】
2015108994.tif 000011
ただし、
【数１０】
2015108994.tif 000012
である。
【００３１】
式（１０）の評価関数Ｌ２の第１項は、連立方程式の残差ベクトルのＬ２ノルムの２乗である。この第１項は、解ベクトルｘの各要素（即ち状態変数ｓと構造変数ρとの要素毎の積）が正解に近づくほど小さい値をとる性質をもつ。また、評価関数Ｌ２の第２項は、変数ρの総和に対応し、解ベクトルｘがスパースになるほど（ゼロ要素が増すほど）小さい値をとる性質をもつ。また、第３項は、変数ρと（１−ρ）の積の総和に対応し、変数ρが０又は１に近付くほど（ゼロ要素か非ゼロ要素かの分離が進むほど）小さい値をとる性質をもつ。係数λ１、λ２は重みを調整するためのパラメタである。このような評価関数Ｌ２を用いることで、非ゼロ要素の抽出およびより良い解の探索を効率的に進めることができる。また、スパース性を評価する第２項、第３項の値は、解ベクトルｘの影響を受けないため、ＦＩＳＴＡのような問題は生じない。
【００３２】
図１は、求解装置によって実行される求解処理の流れを示すフローチャートである。図１に示す各処理の実行主体は、求解装置のＣＰＵ（計算手段）２０１である。
【００３３】
ステップＳ１０１において、ＣＰＵ２０１は、ユーザによって入力された問題、すなわち、係数行列Ａの値と定数ベクトルｂの値をそれぞれ設定する。また、ＣＰＵ２０１は、状態変数ｓのメモリ領域を確保し、初期値を設定する。
【００３４】
ステップＳ１０２において、ＣＰＵ２０１は、状態変数ｓのステップ幅βｓと制御変数νのステップ幅βνを設定する。本実施形態では、次式のようにステップ幅βｓ、βνを決める。
【数１１】
2015108994.tif 000013
ただし、Ｌｉｐ（∇Ｌ２）は、式（１０）の評価関数Ｌ２のgradientのLipschitz定数で
ある。
【００３５】
ステップＳ１０３では、ＣＰＵ２０１は、制御変数ν及び構造変数ρのメモリ領域を確保し、初期値を設定する。また、ＣＰＵ２０１は、反復回数ｋを１に設定する。
以上の初期化処理が完了したら、ＣＰＵ２０１はステップＳ１０４〜Ｓ１０６の反復計算を繰り返す。
【００３６】
ステップＳ１０４は、状態変数ｓの各要素ｓｊの値を更新する状態変数更新ステップである。本実施形態では、ＣＰＵ２０１は次式により状態変数ｓ（ｋ）からｓ（ｋ＋１）を計算する。
【数１２】
2015108994.tif 000014
ここで、Δｓ（ｋ＋１）は、評価関数Ｌ２の状態変数ｓに関する微分であり、次式で求まる値である。
【数１３】
2015108994.tif 000015
【００３７】
ステップＳ１０５は、構造変数ρの各要素ρｊの値を更新する構造変数更新ステップである。本実施形態では、まず、ＣＰＵ２０１は次式により制御変数ν（ｋ）からν（ｋ＋１）を計算する。
【数１４】
2015108994.tif 000016
ここで、Δν（ｋ＋１）は、評価関数Ｌ２の制御変数νに関する微分であり、次式で求まる値である。
【数１５】
2015108994.tif 000017
更に、ＣＰＵ２０１が式（８）を用いて制御変数ν（ｋ＋１）から構造変数ρ（ｋ＋１）を計算する。
【００３８】
ステップＳ１０７は、反復計算を終了するか否かを判断するステップである。ＣＰＵ２０１は、反復計算の状態が所定の収束条件を満足するか否かを評価し、収束条件を満たさない場合は反復回数ｋをインクリメントしてステップＳ１０４に戻り、収束条件を満たす場合は反復計算を抜けてステップＳ１０８に進む。収束条件としては、例えば、評価関数
Ｌ２の値が所定の閾値よりも小さくなること、連立線型方程式の残差ノルム（評価関数Ｌ２の第１項に相当）の値が所定の閾値よりも小さくなること、反復回数ｋが所定数に達すること、などの条件を用いることができる。あるいは、評価関数、残差ノルム、反復回数のうちの二つ以上を組み合わせた収束条件を設定してもよい。組み合わせとは、複数の条件（命題）の論理積又は論理和である。例えば、「反復回数が所定数を超え、且つ、評価関数の値が閾値よりも小さくなること」を収束条件としたり、「評価関数の値が第１の閾値よりも小さくなるか、又は、残差ノルムの値が第２の閾値よりも小さくなること」を収束条件としてもよい。あるいは、評価関数や残差ノルムの値が閾値を下回る状態が複数回続くこと、を収束条件としてもよい。
【００３９】
ステップＳ１０８は、反復計算終了後に、連立線型方程式の解を決定するステップである。簡単には、最後の反復計算で得られた状態変数ｓと構造変数ρとの要素毎の積の値を方程式の解に選べばよい。あるいは、ＣＰＵ２０１が、ステップＳ１０４〜Ｓ１０６の反復計算の中で、評価関数Ｌ２の値又は残差ノルムを最も小さくする変数ｓ，ρをメモリに記憶しておき、そのときの変数ρとｓとの要素毎の積の値を最適解として出力してもよい。
【００４０】
以上述べた処理によって、劣決定系の連立線型方程式の解を精度よく求めることができる。本実施形態の求解方法の性能を図４の完全再構成マップの線４０４に示す。従来方法のＳＰ（４０２）及びＦＩＳＴＡ（４０３）に比べて、本実施形態の方法の性能（４０４）は理論限界（４０１）に近づいている。特に疎性率が０．１の付近および０．４の付近において、従来方法では不可能だった問題の再構成が可能になっていることが分かる。
【００４１】
（変形例）
図３のフローでは、状態変数更新ステップＳ１０４のなかで状態変数ｓの更新処理（式（１３））を１回だけ実行したが、ステップＳ１０４のなかで同更新処理を複数回実行することも好ましい。つまり、構造変数ρの値を固定したまま状態変数ｓの値を複数回更新するのである。これにより、より良い性能が得られることがある。更新回数は５回や１０回のように予め設定してもよいし、動的に変えてもよい。動的に更新回数を変える方法としては、例えば、ステップＳ１０７での条件判定と同じように、評価関数Ｌ２や残差ノルムの値やその変化率などに基づき更新処理の繰り返しを止めればよい。
【００４２】
同じように、構造変数更新ステップＳ１０５のなかで構造変数ρの更新処理を複数回実行することも好ましい。つまり、状態変数ｓの値を固定したまま構造変数ρの値を複数回更新するのである。これにより、より良い性能が得られることがある。この更新回数も上記同様、予め設定してもよいし、評価関数Ｌ２や残差ノルムなどに基づき動的に変えてもよい。
【００４３】
（適用例）
本実施形態の求解方法を被検体情報取得装置における画像再構成処理に適用した例について説明する。被検体情報取得装置とは、観測系によって得られた観測データから被検体の内部の情報を表す画像データを再構成する装置であり、例えば、医用画像診断装置（Ｘ線ＣＴ、ＭＲＩ、光音響トモグラフィ、超音波診断装置など）や非破壊検査装置が該当する。
【００４４】
図５は、本実施形態の求解装置を搭載した医用画像診断装置の構成を模式的に示す図である。図に示すように、医用画像診断装置は、観測系（測定装置）５０１、投影データ生成部５０２、データ取得部５０３、再構成部５０４、表示装置５０５などを備える。
【００４５】
観測系５０１は、測定装置筐体５０６、信号発生源５０８、５０９、５１０、検出器５
１１、５１２、５１３を有している。信号発生源５０８、５０９、５１０は、筐体５０６の内部に置かれた被検体５０７に対し、測定用信号を照射する装置である。例えば、Ｘ線ＣＴであればＸ線源、ＭＲＩであれば磁界発生源、光音響トモグラフィであればレーザ光源、超音波診断装置であれば超音波送信機が該当する。被検体５０７を透過した信号、又は被検体５０７の内部で発生もしくは反射した信号が、検出器５１１、５１２、５１３で検出される。
【００４６】
投影データ生成部５０２は、複数の検出器５１１、５１２、５１３で検出された信号を所定のフォーマットに従って一つの観測データにまとめる。例えば２次元画像のラスタースキャンによるベクトル表現のようにである。これによって、それぞれの検出信号ｂ１、ｂ２、ｂ３が一つのベクトルｂとして統合される。
【００４７】
データ取得部５０３は、投影データ生成部５０２から取得した観測データ（ベクトルｂ）と、予め記憶されている観測系５０１の特性を表すデータ（係数行例Ａ）とから、連立方程式Ａｘ＝ｂを作成する。ｘは未知変数であり、Ｎ次元のベクトルである。係数行列Ａは、観測系５０１を構成する信号発生源や検出器の物理的特性などを元に理論的に作成してもよいし、内部構造が既知のサンプルを観測系５０１で実測した結果から作成してもよい。図５のＡ１、Ａ２、Ａ３はそれぞれ、信号発生源５０８と検出器５１１の特性、信号発生源５０９と検出器５１２の特性、信号発生源５１０と検出器５１３の特性を示している。
【００４８】
再構成部５０４は、図１に示した求解方法を用いて未知変数ｘを求める。そして、再構成部５０４は、得られた未知変数ｘの解に基づき、被検体５０７の内部構造を表す画像データを生成する。生成されたデータは表示装置５０５に表示される。このとき、被検体の断面（２次元画像）を表示してもよいし、被検体内部の３次元構造を表示してもよい。その際、ユーザーインターフェイスによって、画像の回転・並進などの操作を指示できるとよい。
【００４９】
ところで、本実施形態の求解方法は、未知変数ｘのうちの多くがゼロ要素とみなせるとの仮定をおいている。それゆえ、未知変数ｘの疎性が十分でない場合には、解の推定精度が低下するなど、期待する性能が得られない可能性がある。そこで、データ取得部５０３が、係数行列Ａの代わりに係数行列（ＡΦ−１）を用いて、連立方程式（ＡΦ−１）ｚ＝ｂを作成し、再構成部５０４が図１に示した求解方法を用いて未知変数ｚを求めた後、ｘ＝Φ−１ｚにより本来の未知変数ｘを計算してもよい。Φは、未知変数ｘの疎性を高める作用をもつ特徴量変換行列である。つまり、未知変数ｘに変換行列Φを作用させた変数であるｚ＝Φｘは、変数ｘに比べて高い疎性をもつので、本実施形態の方法による求解が容易となるのである。変換行列Φとしては、被検体の内部の情報を表す画像データ（変数ｘ）のゼロ要素を増加させるフィルタ、例えば、微分フィルタ、２次微分フィルタ、ハイパスフィルタ、カットオフフィルタなどを用いることができる。このような方法により、画像再構成の精度をより高めることが期待できる。
【符号の説明】
【００５０】
２０１…ＣＰＵ、２０２…表示装置、２０３…入力装置、２０４…ＲＡＭ、２０５…ハードディスク装置、２０６…ＮＩＣ、２０７…バス
【図１】
2015108994.tif 000019
【図２】
2015108994.tif 000020
【図３】
2015108994.tif 000021
【図４】
2015108994.tif 000022
【図５】
2015108994.tif 000023