特集「最適化技術に基づく統計的推論」［研究詳解］

統計数理（2013）
第 61 巻 第 1 号 135–146
c 2013 統計数理研究所
X 線自由電子レーザーによる分子の電子密度推定1池田 思朗 ・河野 秀俊

2（受付 2012 年 6 月 29 日；改訂 10 月 9 日；採択 10 月 10 日）

要旨

米国，欧州，日本及び韓国では X 線自由電子レーザーと呼ばれる装置が開発され，米国，日
本では運用段階にある．現在，この装置によって初めて可能となった X 線領域の波長を持つ高強度レーザー光を用いた，様々な実験が計画および実施されている．本稿では，そうした実験のひとつ，X 線回折を用いた生体単粒子の電子密度推定について説明し，2 次元回折画像から
2 次元電子密度を推定する位相回復問題に対する我々の研究成果について解説する．
キーワード： X 線回折，XFEL，位相回復，MAP 推定，疎性．
1. はじめに

近年開発された X 線自由電子レーザー（X-ray Free Electron Laser: XFEL）1) によって，これまでにない高強度 X 線レーザー光が出力できるようになった．現在，日本では XFEL の運用が始まり，このレーザーを用いる様々な実験が開始されている．本稿では，その実験のなかか
ら生体単粒子の 3 次元電子密度推定について説明し，その中で必要となる位相回復問題に関する我々の研究成果について解説する．
生体単粒子，その中でも特にタンパク質分子などの生体高分子の 3 次元構造を知ることは，医学を始め，生化学の様々な分野で重要である．生体高分子の詳細な立体構造は，主に X 線結晶構造解析によって決定されている．X 線結晶構造解析では，結晶中の分子の反復構造によって回折した X 線が作る膨大な数の回折点を取得する．この回折点を用いて分子の電子密度を推定し，立体構造を決定する．この方法によって，すでに 8 万以上のタンパク質分子の立体構造が決定されている．
しかし，全ての生体高分子の構造を X 線結晶構造解析によって決定するのは難しい．細胞中の約半数のタンパク質分子は結晶にならないからだ．そうした結晶にならない生体高分子に対し，単粒子の状態で X 線を当てて回折画像を計測し，電子密度を解析する実験が開始されている．単粒子であるがゆえに，回折データを得るためには位相の揃った非常に強い X 線を照射しなければならない．しかし，このような非常に強い X 線にさらされると，分子は数十フェムト秒で壊れてしまう．そのため，実験に用いる X 線レーザーは，高強度かつ短パルスでなければならない．XFEL の光は，位相の揃った高強度かつ短パルスの X 線であるため，分子が破壊される前に回折画像を撮影できると考えられている．様々な方向から生体高分子の回折画像を撮り，得られた複数の回折画像を処理することによって 3 次元電子密度を推定する方法が提案されている（Neutze et al., 2000; Huldt et al., 2003）．
1
2

統計数理研究所：〒190–8562 東京都立川市緑町 10–3
日本原子力研究開発機構：〒619–0215 京都府木津川市梅美台 8–1–7

統計数理

136

第 61 巻

第 1 号 2013

こうした方法を用い，実際に巨大ウィルスであるミミウィルスのコヒーレント X 線イメージングが行われている（Seibert et al., 2011）．しかし，生体内のタンパク質分子などは一般にミミウィ
ルスよりはるかに小さいため，XFEL であってもミミウィルスと同等のコントラストを持つ回折画像を撮るのは難しく，3 次元電子密度を推定するには様々な困難があると想像される．本稿ではそうした困難の中から 2 次元回折画像からの位相回復という問題をとりあげ，我々が提案するSparse Phase Retrieval（SPR）法（Ikeda and Kono, 2012）について説明する．SPR 法は LASSOなどと同様に情報源の疎性を仮定した推定法であり（Tibshirani, 1996; Zhao and Yu, 2006），ベ
イズ統計の枠組においては最大事後確率法（Maximum a posteriori: MAP）推定を行なうことに対応する．SPR 法は，従来の方法が適用できない回折画像面に到達する光子数が少ない状況下において，単粒子の電子密度推定を可能にする．
2. XFEL を用いた生体単粒子回折画像の計測
2.1 回折画像
まず，回折画像の基本的な性質について説明する．図 1 に実験の概略を示す．真空中に落下させた生体単粒子に X 線を当て，回折画像を得る．図中に示したように I（photons/pulse/mm2 ）は X 線の強度，L は分子の直径である．
生体単粒子の電子密度を X 線と垂直な面に射影したものを fxy ≥ 0，x, y = 1, . . . , M と表わす．x, y は 2 次元座標のインデックスである．X 線レーザー光によって得られる回折画像はその f のパワースペクトラム，すなわち 2 次元フーリエ変換したものの大きさに対応する量である．そこで，fxy のフーリエ変換を Fuv とする．定義は以下の通りである．

1 
2πi(ux + vy)
(2.1)
Fuv =
fxy exp.
M x,yM

通常，上式の右辺にある係数は 1/M ではなく 1/M 2 だが，fxy と Fuv の 2 乗ノルムを等しくするため，本稿では上式の定義を用いる．簡便のため，以下の記号を定義しておく．
f = {fxy } = {f11 , . . . , fM M },F = {Fuv } = {F11 , . . . , FM M }.

また，フーリエ変換を F = F(f ) と表わすことにする．
2 次元電子密度 f と理想的な回折画像 Suv の間には以下のような関係があることが知られ
ている．

2
λ
(2.2)
Suv = α|Fuv |2 cuv = Irc2
|Fuv |2 cuv .
σL

図 1.

生体単粒子を用いた XFEL による回折実験の概略図．

X 線自由電子レーザーによる分子の電子密度推定

図 2.

137

リゾチームの電子密度（左）と理想的な回折画像（右）．各画像は 308×308 の大きさの像である．

ここで rc は電子の古典半径（2.82×10−12 mm），λ は X 線の波長（0.1 nm），σ は各次元のオー
バーサンプリング比（2），cuv = c(u, v) = cos3 θ である（図 1 参照）．本来は回折画像は平面ではなく Ewald 球面上の画像として理解すべきだが，ここでは回折の角度が小さいとしてその影響は無視する．本研究の方法を Ewald 球面上のデータに適用することは可能である．
図 2 に理想的な状態での電子密度と回折画像をシミュレーションで作成したものを示す．結
晶化した生体単粒子を用いた X 線回折画像も同様の関係式に従う．結晶は非常に多くの粒子が規則的に並んだ構造を持つため，電子密度のフーリエ変換である回折画像はブラッグの回折条
件（Bragg’s Law）を満たした点のみで光子が観測される．一方，単分子の回折画像は図にあるような全体に連続的に拡がった像となる．
2.2 XFEL による回折画像
図 3 に示すのは XFEL による X 線レーザーをリゾチームというタンパク質分子にあてて得られる回折画像を計算機シミュレーションによって作った例である．図 2 に対して粗い画像であることがわかる．こうした画像となる理由は，（2.2）式に示したように，回折強度が古典半径の自乗に比例しているため小さい単粒子では回折強度が弱くなり，高強度の XFEL であっても観測面に到達する光子の数が少なくなってしまうからである．この図では 4630 個の光子のみが観測面に到達している．図 3 は約 90000 画素の図であるから，ほとんどの画素では光子が観測されていないことがわかる．また，各画素で観測されるのは光子の数，つまり理想的な場合のような実数値ではなく，非負の整数値である2) ．

図 3.

XFEL によって得られると予想される回折画像．強度を表す I は
5.0 × 1021 photons/pulse/mm2 であり，観測面に到達した光子の数は 4630 個である．

138

統計数理

第 61 巻

第 1 号 2013

画素 (u, v) で得られる光子数を Nuv と置く．Nuv の期待値は（2.2）式の Suv だが，実際に観測される光子数は全く同一の構造に同じ方向から X 線を照射したとしても同じになるとは限らず，実験を行う度に確率的に振舞う．全光子数 Nall = uv Nuv の期待値は強度 I に比例し，各画素で観測される光子数はポアソン分布に従うことが知られている．各画素で観測される光子数は独立で，各 Nuv が独立でポアソン分布に従うというのは光学で一般的に用いられている仮定である（Hau-Riege et al., 2005）．このとき Nuv の分布は以下のように書ける．
p(Nuv |Suv ) =Nuv
Suv
exp(−Suv ).
Nuv !

（2.2）式より Suv は α を正の定数として Suv = α|Fuv |2 cuv とかける．フーリエ変換によって得られる F は f の関数であることから，N = {Nuv } = {N11 , . . . , NM M } の分布は f の関数であり，以下のように書ける．
(2.3)

p(N |f ) = (|Fuv |2 cuv )Nuv exp(−|Fuv |2 cuv )
uv

Nuv !

.

なお，ここでは α は 1 とした．データだけから α の大きさと |Fuv |2 の大きさの両方を同時に決定できない．|Fuv |2 の大きさは全エネルギーなど別の情報をもとに決定する必要がある．
2.3 位相回復
本稿での位相回復とは，回折画像から電子密度を推定する問題である3) ．
まずは，光量は十分にあると仮定し，Suv から fxy を復元する際，なぜ位相回復問題が生じるかを説明する．
Suv が観測されるならば， |Fuv | が観測できるとして良い． Fuv は fxy のフーリエ変換にあたることから，逆フーリエ変換を用いて fxy が復元できそうである．しかし 観測されるのは，
Fuv の大きさのみであって位相情報は得られないため，逆フーリエ変換を適用できない．したがって fxy の復元には困難が伴う．
別の見方をしよう．
（2.1）式の Fuv の定義，そして fxy が実数をとることから，

1 2πi((M − u)x + (M − v)y)
(2.4)
FM −u M −v =fxy exp
M x,yM−2πi(ux + vy)
1 ∗
fxy exp, (∗ は複素共役を示す).
=
= Fuv
M x,y
M

すなわち |FM −u M −v | = |Fuv | となる．回折画像 Suv は M 2 画素の画像だが，そのうち半分は同じ値を持つことになり， Suv からは M 2 /2 個の画素の値しか得られない．これらの値を用いて fxy の M 2 の画素の値を復元しなければならない．したがって位相回復問題は不良設定問題である（図 4）．
2.4 既存法
ここでは，単粒子の 2 次元回折画像の解析における位相回復問題に対する既存の方法について説明する．
既存の方法では，2 次元空間の電子密度 f が正となりうる領域，すなわち正となりうる成分が限定されていると仮定する．言い換えると，ひとつの分子のみが真空状態に存在し，そのまわりには何もない状態を仮定する．その領域を γ と書くことにする．そして，次の関数の最適X 線自由電子レーザーによる分子の電子密度推定

図 4.

139

回折画像の対称性．同じ色の領域は同じ値をもっている．

図 5. （2.5）式 における制約．この図では正方形の領域を仮定したが，任意の形で構わない．

値を与える f を推定結果としている．
(2.5)

minimize

1/2 2(|Fuv |c1/2
uv − Nuv ) ,uvsubject tofxy ≥ 0, (x, y) ∈ γ and fxy = 0, (x, y) ∈/ γ.

上式の目的関数は電子密度を f としたときの回折画像 |Fuv |2 cuv と観測された回折画像の各画
素の値 Nuv との近さを表している．γ 上のみで正となりうる f によってこの目的関数を小さくすることが目的である（図 5）．
この問題を解くアルゴリズムとして広く用いられているのが hybrid input-output（HIO）法（Fienup, 1978; Fienup, 1982; Miao et al., 1998）である．XFEL による単粒子回折画像の解析でも HIO 法が用いられている（Neutze et al., 2000; Huldt et al., 2003）．HIO 法が
（2.5）
式の最適化
問題を解いていること，また，HIO 法のアルゴリズムが既存の最適化法である Douglas-Rachford
法の適用例とみなせることは既に報告されている（Fienup, 1982; Bauschke et al., 2002）．
単粒子の重さなど他の指標によって f の正となりうる領域 γ がわかっていて，理想的な回
折画像 |Fuv |2 cuv が正確に求まっている場合には，HIO 法によって電子密度が復元できる．しかし，理想的な回折画像を実験的に得るためには，十分な光量が必要である．ミミウィルスを対象とした XFEL による単粒子回折画像の解析に関する先行研究では HIO 法によって電子密度が復元されている（Seibert et al., 2011）．これは粒子が大きいために十分な光量が得られ，理想に近い回折画像が観測できたためだと考えられる．
3. 提案する方法
3.1 SPR 法
タンパク質分子などの単粒子に XFEL を用いて得られる回折画像は十分な光量が得られな
いため，前節で示した既存の方法は必ずしもうまく働かない．そこで我々は SPR 法を提案し

統計数理

140

第 61 巻

第 1 号 2013

た（Ikeda and Kono, 2012）．
（2.3）式にあるように，観測された回折画像は確率モデルによって記述できることから，電子密度 f を復元することは推定の問題である．しかし，2.3 節で説明したように，位相回復問題は
不良設定問題であり，最尤推定に関しても解は一意に求まらない．そこで，我々は電子密度に関する事前知識を事前分布として記述し，ベイズ統計の枠組における MAP 推定によってこの問題を扱った．位相回復に対してベイズ統計の方法を用いた方法は他にも提案されているが，我々とは問題設定，事前分布の選択が異なる（Irwan and Lane, 1998; Baskaran and Millane, 1999）．
事前分布は電子密度 f に関する知識を反映するものであることから，中心に高い確率で値を持ち，具体的にどの画素かは分らないが，多くの画素が 0 になると仮定するのが妥当だろう．
こうした仮定を満たす事前分布として，LASSO（Tibshirani, 1996; Zhao and Yu, 2006）で用いられたものと同様に指数分布を用いた事前分布を仮定する．具体的には以下のような事前分布を仮定する．
p(f ) =p(fxy ), p(fxy ) = ρxy exp(−ρxy fxy ),fxy ≥ 0, ρxy ≥ 0,xy

定義から分かるように，事前分布では各画素 fxy は独立であるとしている．また，ρxy はハイパーパラメータである．このような関数による事前分布を用いることによって，具体的に画素を指定することなく多くの画素が 0 となるが，ρxy の値が大きいほど fxy = 0 となり易い．我々は周辺のほうが 0 が多くなるように， ρxy を x, y の関数として以下のような関数 ρ(μ)xy と定
めた．

2 
2 
1+M
1+M
(3.1)
ρ(μ)xy = μwxy , wxy = a
x−
+ y−
+ b.
2
2
a と b は wxy が中心で 0，端で 1 となるように調整するパラメーターである（図 6 参照）．ここでは図 6 のように放物面の形状をした事前分布を仮定したが，電子密度に関してより詳細な事前知識があればそれを反映する関数を用いるべきである．なお，（3.1）式にある μ は事前分布
全体の大きさを調整するパラメータである．
回折画像 N を観測したときの電子密度 f の事後分布は，事前分布と尤度に対して以下の関係を持つ．
p(f |N ) ∝ p(N |f )p(f ).

（2.3）式の対数をとり p(f |N ) から f に関連する項を取り出す．そのMAP 推定を考えるために

図 6. 事前分布の重み wxy の 3 次元表示．周辺が 0 となり易い．

X 線自由電子レーザーによる分子の電子密度推定

141

結果，以下の関数が得られる．


(3.2)
wxy fxy .
Nuv ln |Fuv |2 − |Fuv |2 cuv − μµ (f |N ) =uvxy

第一項はポアソン分布の尤度に関する項である（同様の考察は他の研究でも見られる（Irwan and（3.1）式で導入した μ は 2
Lane, 1998））．第 2 項は指数分布による事前分布に関する項である．
つの項のバランスを調整するパラメータとみなすことができる．SPR 法では（3.2）式を最大にする f を求める．すなわち，以下の最適化問題を解く．
(3.3)

maximizeµ (f |N )subject tofxy ≥ 0(1 ≤ x, y ≤ M ).

一般に，目的関数 µ (f |N ) は f に関して上に凸な関数ではない．しかし，（3.3）式から，μ をひとつ決めればそれに応じた最適化問題が定まる．従って，どのように最適値を与える f を求めるかは別にして， μ をひとつ決めればそれに応じた f がひとつ決まる．我々は予め μ を決めることはせず，いくつかの μ に対してこの問題を解くことにする．その中からどの μ に対する電子密度を用いるかは後で議論する．
μ を非常に大きくとれば電子密度 f のほとんどが 0 と推定され，次節で示す勾配法が収束する f は比較的簡単に求まる．目的関数は凸でないが，この場合には最大値が求まる可能性が高いと考えられる．そこで最初は μ を大きく設定して問題を解くことで電子密度を推定し，次に μ を少し小さくし，前の推定結果を勾配法における初期値として用いて電子密度を推定す
る．これを繰り返し，複数の μ に対して電子密度を推定する．
3.2 アルゴリズム
前節で示した（3.3）式を解く方法として勾配法に基づく方法を提案する．まず，準備のために以下の関係式を示しておく．

22πi(ux + vy)∂|Fuv |2
Re Fuv exp −
=
.
∂fxy
M
M

逆フーリエ変換を以下のように定義する．
(F

すると

−1



1 
2πi(ux + vy)(F ))xy =Fuv exp −,M uv
Mµ (f |N ) の f による微分は以下のように書ける．

∂ µ (f |N )= 2Re F −1 (g(F ; N )) xy − μwxy ,∂fxy

ただし g(F ; N ) = {guv (Fuv ; Nuv )} は以下に定義する．

Nuv
guv (Fuv ; Nuv ) =−c
uv Fuv .
|Fuv |2µ (f |N ) の f による微分を計算するためには f のフーリエ変換が必要になり，さらに g(F ; N )に対して逆フーリエ変換を行なう必要がある．SPR 法は勾配法を用いるため，繰り返し演算が必要になるが，その 1 ステップの計算量は，同じく繰り返し演算を必要とする HIO 法の 1 ス
テップの計算量と同程度である（Fienup, 1982）．電子密度 fxy ，x, y = 1, . . . , M の勾配法に基づ
く更新則は以下のように書ける．
∂ µ (f t |N )
t+1
t
(3.4)
fxy = max 0, fxy + ηt
,
t
∂fxy

統計数理

142

第 61 巻

第 1 号 2013

ここで t は繰り返しの回数を示し，max 操作によって fxy は非負となる．正の定数 ηt がステップ幅を示すが，このステップ幅を直線探索によって変化させることによって収束速度が飛躍的に向上する．
一般に勾配法では初期値が重要である．最初 μ を非常に大きな値にとれば，多くの fxy が 0となり，初期値によらず解が比較的簡単に求まる．そうして求めた解を μ を少し小さくした問
題の初期値とする．このように，複数の μ に対して解を求める．
3.3 既存法のベイズ的解釈
既存法では（2.5）式の最適化問題を解いている．この最適化問題は形式的には SPR 法と同様に MAP 推定と考えることができる．
以下のように尤度と事前分布を定義する．

1/C, (x, y) ∈ γ
1/2
1/2 2
p(N |f ) ∝ exp −A(|Fuv |cuv − Nuv ) , p(fxy ) =0,(x, y) ∈/ γ,uv

上式で A, C は正の定数である．尤度は指数の肩に 2 乗誤差が乗った正規分布の形をしている
が |Fuv | が非負であるために規格化定数は正規分布のものとは異なる．また，事前分布は一様分布としているが，各 fxy は 0 ≤ fxy ≤ C に制限されると仮定している．ただし C は十分大きい正の実数とする．
こうして定義した尤度と事前分布を考えると，事後確率の対数をとったものは（2.5）式の目
的関数と一致する．したがって，既存法も形式的には MAP 推定としてとらえることができる．
既存法では γ を設定しなければならないが，これをあらかじめ知ることは一般に難しい．
SPR 法では γ を設定する必要はない．ただし，μ という新たなパラメータを調整する必要がある．
3.4 数値実験
本節では図 3 に示した回折画像に対して SPR 法を適用した結果を示す．図 7 は μ を変化させて復元したタンパク質分子の電子密度を示したものである．大きい μ に対しては多くの画素
で電子密度が 0 と推定されていることがわかる．ある程度小さい μ に対しては図 2 左にある
電子密度を十分に復元できる．
ここで問題となるのは μ の選択である．μ を小さくすると，一般に尤度は大きくなる．適
切な μ を求めるためには，尤度ではない何らかの規準によって最適値を与える μ を選ぶ必要がある．もし，同じ条件で測定したデータが複数あるのならば， cross validation などの方法でμ を選ぶことができるが，位相復元の問題では同じ条件で測定したデータはひとつしかないた

SPR 法で復元した電子密度 µ = 100 とした．
図 7.

SPR 法で復元した電子密度 µ = 0.1 とした．

SPR 法によって復元された 2 次元電子密度．

X 線自由電子レーザーによる分子の電子密度推定

143

め，他の方法を考える必要がある．
この点は LASSO を用いたデータ解析でも問題となる． Zou らは LASSO では，0 でないパラメータの個数が確率モデルの自由度（degree of freedom）の推定量となることを示している（Zou et al., 2007）．そして μ を選択するときに自由度を罰則項に含むモデル選択基準， Akaike
information criterion（AIC）や Bayesian information criterion（BIC）などを使うことを提案している．SPR 法においても，モデルの自由度を推定し，モデル選択基準によって μ を求めることが考えられる．しかし，ポアソン分布を尤度としている SPR 法では，0 でないパラメータの個数は自由度の推定値とはならない．自由度の推定法については今後の課題である．
光学の分野では伝統的に規格化された 2 乗誤差を規準にしている．SPR 法でも 2 乗誤差を基準として μ を選ぶことにする．具体的には，ある μ に対して SPR 法によって µ (f |N ) を最大化して求めた電子密度を f̂ (μ) とし，以下の関数によって誤差を定義する．
(|F̂ (μ)uv |cuv 1/2 − Nuv 1/2 )2(3.5)
ErrorF = uv.
 u v  Nu v
ここで F̂ (μ)uv は F̂ (μ)uv = F(f̂ (μ))uv とした．ErrorF = 0 となるのなら，|F̂ (μ)uv |cuv 1/2 =Nuv 1/2 であり，尤度は最大になる．しかし，ErrorF は SPR 法の尤度とは異なるため， μ を小さくしたとしても ErrorF は必ずしも減少しない．厳密な議論ではないが，μ を非常に大きな値から小さくしていくと，ErrorF は徐々に小さくなり，どこかで増加に転じる．そこで ErrorF
を最も小さくする μ を選ぶことにする．図 3 の回折画像に対して μ を変化させ，それぞれに対する ErrorF を表示したものが図 8 である．μ は 1×104 から 1×10−2 まで変化させた．ErrorF
を最も小さくする μ は μ = 0.1 であった．なお，このときの復元された電子密度画像は図 7 の
右，ErrorF の値は 0.209 であった．同じデータに対して HIO 法を用いて電子密度を復元した
とき（ただし，fxy が正となりうる領域 γ は 159 × 159 画素の正方形の領域とした．
）の ErrorF
は 0.220 であった．
最後に計算速度について述べる．SPR 法は Intel Xeon X5570（2.93 GHz 4 cores）を搭載した
デスクトップ機に C 言語で実装した．上記の実験では，全ての μ に対して 10 分以内で解が求まった．
4. まとめ

本稿では，我々が提案した位相回復の解法を解説した．提案する SPR 法はベイズ統計の MAP推定に基づく方法であり，今後 XFEL によって得られる生体単粒子などの回折画像からの位図 8.

図 3 の回折画像に対して µ を変化させて SPR 法を行った結果得られた ErrorF を示したもの．黒い三角形は最小の ErrorF に示した．横軸 µ は左から右に増加するが，実験を行なう際には大きい µ から始めて小さくしていった．

統計数理

144

第 61 巻

第 1 号 2013

相回復に対して有効な方法であると考えられる．
ベイズ統計の立場から考えると，SPR 法も既存法も MAP 推定とみなせる．ただし，二つの方法では尤度と事前分布，それぞれの分布が異なる．光子の数が十分に多い場合，SPR 法と既存法による電子密度の復元結果に大きな差は見られなかった（Ikeda and Kono, 2012）が，光子の数が少ない場合，復元結果には差が見られた．以下ではその差について簡単に考察する．
まず，尤度について考える．単粒子の回折画像では中心およびその周辺に大部分の光子が観測される．回折画像の中心部は電子密度の空間周波数を考えると低い周波に対応し，単粒子のおおまかな形を表現している．SPR 法では光子数 Nuv をポアソン分布として記述しているた
め，光子数の多い画素の方が尤度に大きく寄与し，回折画像の中心部分を重視することになる．
したがって電子密度の推定値は，単粒子のおおまかな形を反映したものになると考えられる．
一方，既存法は回折画像のどの画素も等価に扱っており，光子の観測されなかった周辺の画素の影響も中心と同じ重みとなる．このため，電子密度の空間周波数を考えると高周波のノイズの影響を受け易いと考えられる．
次に事前分布について考える．光子数の少ない場合，観測された回折画像は理想とは程遠いものとなる．既存法では γ 領域内での f で回折画像を表現しなければならないが，SPR 法では γ を設定しないため，正の値をもつ fxy は領域の全体に拡がる．μ を適切に設定すれば，観測された回折画像をうまく表現できる可能性がある．また，HIO 法を用いる場合，γ 上で定義される f によって回折画像が説明できなければ収束性が保証されない．一方，SPR 法は局所最適解しか得られない可能性はあるが，収束はする．計算量も現実的な範囲であり，既存の HIO
法よりも小さな誤差を持つ解を与えることが確認できた．
XFEL による生体単粒子の回折画像からの 3 次元構造解析には，この位相回復だけではなく，いくつもの問題を克服する必要がある．回折画像はある方向からの観測だけである．3 次元構造を知るためには数多くの画像を取り，それらを正しく組み合わせて 3 次元の回折像を構築して位相回復を行うか，2 次元回折像から 2 次元の実像を求めてそれらを 3 次元再構成する必要がある．輝度が低い画像を組み合わせる方法は未だ十分には検討されておらず，今後の研究が必要である．
注．
1)

独立法人理化学研究所播磨研究所 X 線自由電子レーザー，http://xfel.riken.jp を参照．
分子のサイズが大きい場合にはより鮮明な画像が得られるが，生体単粒子の回折画像はこの程度の光子数だと推定される．
3)「計測と制御」2011 年 5 月号にミニ特集「回折イメージング ∼ 位相回復の新展開 ∼ 」が
ある．
2)

謝

辞

本研究の一部は文部科学省の X 線自由電子レーザー利用推進研究課題の事業として行なわれました．

参 考 文 献Baskaran, S. and Millane, R. P.（1999）. Bayesian image reconstruction from partial image and aliased
spectral intensity data, IEEE Transaction on Image Processing, 8, 1420–1434.

X 線自由電子レーザーによる分子の電子密度推定

145

Bauschke, H. H., Combettes, P. L. and Luke, D. R.（2002）. Phase retrieval, error reduction algorithm
and Fienup variants: A view from convex optimization, Journal of the Optical Society ofAmerica A, 19, 1334–1345.
Fienup, J.（1978）. Reconstruction of an object from the modulus of its Fourier transform, Optics
Letters, 3, 27–29.
Fienup, J.（1982）. Phase retrieval algorithms: A comparison, Applied Optics, 21, 2758–2769.
Hau-Riege, S. P., London, R. A., Huldt, G. and Chapman, H. N.（2005）. Pulse requirements for x-ray
diﬀraction imaging of single biological molecules, Physical Review E, 71, 061919.
Huldt, G., Szőke, A. and Hajdu, J.（2003）. Diﬀraction imaging of single particles and biomolecules,
Journal of Structural Biology, 144, 219–227.
Ikeda, S. and Kono, H.（2012）. Phase retrieval from single biomolecule diﬀraction pattern, Optics
Express, 20, 3375–3387.
Irwan, R. and Lane, R. G.（1998）. Phase retrieval with prior information, Journal of the OpticalSociety of America A, 15, 2302–2311.
Miao, J., Sayre, D. and Chapman, H.（1998）. Phase retrieval from the magnitude of the Fouriertransforms of nonperiodic objects, Journal of the Optical Society of America A, 15, 1662–
1669.
Neutze, R., Wouts, R., van der Spoel, D., Weckert, E. and Hajdu, J.（2000）. Potential for biomolecularimaging with femtosecond X-ray pulses, Nature, 406, 752–757.
Seibert, M. M. et al.（2011）. Single mimivirus particles intercepted and imaged with an X-ray laser,
Nature, 470, 78–81.
Tibshirani, R.（1996）. Regression shrinkage and selection via the lasso, Journal of the Royal StatisticalSociety, Series B, 58, 267–288.
Zhao, P. and Yu, B.（2006）. On model selection consistency of lasso, Journal of Machine Learning
Research, 7, 2541–2563.
Zou, H., Hastie, T. and Tibshirani, R.（2007）. On the “degrees of freedom” of the lasso, Annals of
Statistics, 35, 2173–2192.

146

Proceedings of the Institute of Statistical Mathematics Vol. 61, No. 1, 135–146(2013)

Molecular Electron Density Estimation with X-ray Free Electron Laser
Shiro Ikeda1 and Hidetoshi Kono2
1 The Institute of Statistical Mathematics
2 Japan Atomic Energy Agency

New laser equipment x-ray free electron laser facilities are built in the United States,
Europe, Japan, and Korea. A variety of new experiments in physics, chemistry, and biology
are being carried out with this new X-ray laser, and new data will be collected. We explain
one such experiment, electron density estimation from a single biological particle, andshow our recent result on phase retrieval, which recovers the two-dimensional electron
density from a two-dimensional X-ray diﬀraction pattern.

Key words: Phase retrieval, X-ray diﬀraction imaging, MAP estimation, sparsity.