特集／ LATEX 2ε クラスファイル

大域的最適化を用いた先進結晶構造解析
特集記事用

富 安

亮 子

1. はじめにの比である．各 p j は，ある物体の座標 z j における密度の大小に対応し，物体の形の詳細が分かる．

大域的最適化手法である半正定値計画法 (SDP)
は，緩和 (SDR) と呼ばれる手法により，2 次計画問題を中心とする非線形最適化を短時間で実施できる．さらに，間違ったモデルの確度 100%での棄却，状況によっては，大域的最適解の一意性の保証や，複数ある場合にその全てを求めることができる，といった性質が知られている．
結晶構造解析では，位相問題という様々な計測の問題に表れる問題のため，出力（= 実験データ）
が入力 (= 求めたい構造) の 2 次式となる状況で，出力から入力を求める逆問題を解く．フーリエ変換など調和解析に関係する問題だが，近年は圧縮
センシングの適用も議論された1) ．
結晶学・数理計画法の背景に加え，著者と実験家の共同研究で得られた SDR の応用を紹介する2) ．

図1

1 次元トーラス (= 円周, 左) と 2 次元トーラス (右)

考えている状況では，各差ベクトル zi − z j (i, j =1, . . . , m) に対し，pi p j の値を測定できるが，解像度のため，差ベクトルの座標が近いと，それぞれの pi p j を足した値になる．
状況を正確に表すため，差ベクトル zi − z j がほぼ同じかそうでないかで，ペア集合 {(i, j) :i, j = 1, . . . , m} を部分集合 P0 , P1 , . . . , Pt に分割
する．各ペアはただ１つの P s に属し，二つのペア(i1 , j1 ), (i2 , j2 ) が同じ P s に属すならばそのときに
2. 本稿で扱う問題
背景の説明に入る前に，一般の位相回復の問題
を含む差ベクトル情報からの密度割り当て問題を説明する．節 4 で説明される位相回復問題と見た目が大きく異なるが，フーリエ変換すればほぼ同じ問題になることを節 5 で紹介する．

限り，zi1 − z j1 , zi2 − z j2 は観測の解像度でほとんど一致するようグループ分けできたとする．このとき以下の I s (s = 0, . . . , t) の比が計測されるので，比 p1 : · · · : pm を推定せよ，という問題になる．

Is =

∑

pi p j .

(1)(i, j)∈P s

Rd もしくは d 次元トーラス上 (図 1) の有限個のほとんどの P s が一つの (i, j) しか含まない場合

代表点を z1 , . . . , zm とする．知りたいのは既知のを考えると，この問題は簡単に解けそうだが，解各座標点 z j に割り当てられた密度 p1 , . . . pm ∈ C

像度の問題がなくても，そうとは限らないことは，

数理科学

NO. 677, NOVEMBER 20191

位相回復のような問題も含まれることから分かる．
興味深いことに，密度割り当て問題の中には，
上記の自明なケースと，位相回復の非常に困難なケースの間に，SDR により解の一意性の判定や複数解の取得が個別に実施できる場合が含まれる．
サイエンス分野で，結晶構造解析における原子種・占有率の解析，および，磁気構造解析と呼ばれる解析はこの中間層に属し，SDR の優れた機構
が実際に役に立つことが確かめられた2) ．ただし磁気構造解析は， p j が 3 次ベクトルで，各 I s は p jの成分の 2 次式となる特別な場合である．

基本的な方針になる．

)2
N ( obs∑
I (Ki ) − I cal (Ki ; p)
Err[I obs (Ki )]

i=1

.

(2)

I cal (Ki ; p) は p の線形式で表せないため，非線形
最小二乗法を実施する．
求めた構造の信頼性が高いと認められる基準に，R 因子が 0.05–0.1 以下，というものがあり4, 5) ，そのため解析用ソフトウェアを用いて I obs の I cal のグラフをきっちり合わせることが求められる．それでも構造が間違って報告されることがあることに加え，解析の自動化の観点から，結晶学には大域的最適化の需要が存在する．

3. 結晶学分野における大域的最適化の需要

しかし正解にかなり近い初期パラメータが得ら

結晶学では，結晶構造を中心とする物質のミクれている場合を除き，局所最適化である非線形最ロ構造を扱う．求める構造は図 2 のようなものだ

小二乗法で結晶構造を決めるのは難しい．主な原が，各原子が見えるのは Å(= 10−10 m) の解像度な因として局所最小値の問題が知られているが，最ので，結晶構造解析においては，1Å 前後の波長を小化する目的関数が p の凸関数なら自然に回避さ持つ X 線が用いられる．

れる（図 3）．ここで f が凸関数とは，(p1 , f (p1 )),(p2 , f (p2 )) を結ぶ線分が f のグラフ (x, f (x)) の下を通らないことである．
y = f (x)図2結晶構造の例 (球は原子，結晶構造は図の平行六面体を格子のマスとして 3 方向に無限に並べたものを考える．描画は Vesta3) 使用)

•

局所最小値

•

真の最小値

(p1 , f (p1 ))•• (p2 , f (p2 ))結晶学における大域的最適化の需要に著者が出会ったのは，結晶学を知って間もない 2006 年ご

図3

非凸関数の局所最小値 (左図) と凸関数 (右図)

ろのことで，関連する国際会議において simulated

annealing や遺伝的アルゴリズム, モンテカルロ法

結晶構造モデルとしては，独立原子モデルと呼の結果が紹介されていた．その理由をまず紹介する．

ばれるものがよく用いられる; 結晶構造は 3 方向

X 線回折で得られる実験データは，強度と呼ばに周期的なので，ある 1 次独立な b1 , b2 , b3 ∈ R3

れるヒストグラムデータ I obs (K1 ), . . . , I obs (KN ) とが存在して，+b だけ平行移動しても元の構造に一

その見積もり誤差 Err[I obs (Ki )] からなる．図 2 の致する b ∈ R3 全体は，m1 b1 + m2 b2 + m3 b3 (mi :結晶構造を表す密度関数を考え，m 次元ベクトル整数) と書ける点集合，すなわち 3 次元格子となp の各成分を密度関数のパラメータとする．このとる．結晶構造の原子座標は，この格子 L の平行移

き強度の理論値 I cal (Ki ; p) は，Ki , p の関数となる．
最小二乗法の考え方に基づいて以下の値を最小化する最適解 p を求めることは，結晶構造解析の

2

動 z j + L ( j = 1, . . . , m) の合併∪m

j=1 (z j + L) の各点になる．
結晶構造自体は以下のモデル式で表せる．

f (x) =

m
∑

pj

j=1

∑

f j (x − z j − l),

(3)

l∈L

ただし， p j ∈ C は座標 z j + L に存在する原子の種−x /2σ2j類に固有な値で，f j は，正規分布 e2√/ 2πσ2j ，もしくは x が原点から離れるにしたがって単調減少する他の確率分布（積分面積が 1）になる．異なる原子座標の間には距離があり，ピークオーバー

2.

f (x) 7→ f (−x),

3.

f (x) 7→ f (x + x0 ) (x0 だけ平行移動).
特に 1 と 3 から生じる複数解は，構造情報として区別する必要がない．そこで 1–3 の変換で写り合う解は 1 つの同じ解と数え，さらに観測誤差はないと仮定しても，１次元 (d = 1) の位相回復では，解の一意性に深刻な問題が生じることが分かっている6, 7) ．しかし 2 次元以上では，複数解はずっ

ラップの影響はそれほどではない．
通常の結晶構造解析では p j , z j ともに未知変数
だが，節 2 の冒頭に述べたように，以下では p j が未知，z j が既知の場合を一般的に取り扱う．

と生じにくいと考えられる8, 9) ，観測データの解析においては，C1 , C2 は無限次元の関数空間より，有限次元の線形空間と考えた方が，実際の状況に近い．例えば，関数 f ∈ C1 , g ∈ C2 を

Rd の座標点集合 {x1 , . . . , x M }, {k1 , . . . , kN } における
4. 位相問題と結晶構造解析

値で表すと，( f (x1 ), . . . , f (x M )), (g(k1 ), . . . , g(kN ))位相問題が生じることも結晶構造解析の特徴のとなるので，C1 , C2 はそれぞれ C M , CN と同一視一つだが，この問題は X 線回折顕微鏡などのイメージング研究，宇宙望遠鏡の解析でも生じる．
物体（の密度など）を表す関数を f とし，単色波を当てて生じる回折波を遠方から観測する．これは Fraunhofer 回折と呼ばれ，回折波は，以下のフーリエ変換 F f の，ある定数倍に等しい．

∫

F f (k) =2πix·kRdf (x)edx.

(4)

ここで，x · k は，Rd 上のユークリッド内積で，通常は d = 2, 3 である．F f が観測できれば，逆フー
リエ変換

∫

Rd

F f (k)e−2πix·kdk により， f を復元できる．しかし F f の位相（= 複素数の偏角）を観測することは難しく，絶対値のみが観測できる状
況で，位相回復の必要が発生する．
様々な計測の状況を含むよう，問題を定式化すできる．観測値の (|T f (k1 )|2 , . . . , |T f (kN )|2 ) から，( f (x1 ), . . . , f (x M )) の各成分を何らかの方法で推定することが位相回復にあたる．
イメージングの分野で代表的な位相回復手法は，

Gerchberg-Saxton アルゴリズム10) , HIO11) といった
フーリエ反復法で，非晶質の構造解析のため提案されたコヒーレント回折イメージング12, 13) の基盤
技術にもなっている．
初期パラメータを用いない結晶構造解析は未知
構造解析と呼ばれる．反復フーリエ法も近年の手法として存在するが14) ，直接法15, 16) が代表的である．
これらの既存手法においては収束性や解の一意性の証明といった理論上の問題が残る一方，解像度の向上，高分子の解析といった新しい課題に向けた改善も求められている．

ると以下のようになる．Rd 上の関数空間 C1 , C2 の間の線形写像 T : C1 → C2 を考える．入力 f ∈ C1に対し，出力として T f ∈ C2 の絶対値を取った
関数 |T f | が観測される．|T f | から f を決められるか？

C1 = C2 = L2 (Rd ) で，T = F の場合，式 (4) の形から， f を以下のように変換しても，|F f | は変
わらない（証明は読者に任せる）．

1.

5. 最適化問題の導出
本節では結晶構造解析における原子種・占有率の解析が，節 2 の密度割り当て問題と同じ問題であることを示し，さらに 2 次制約つき 2 次計画問題（QCQP) と呼ばれる以下の最適化問題に帰着できることを説明する．

f 7→ ζ f (ζ は絶対値 1 の定数),
数理科学

NO. 677, NOVEMBER 2019

3






最小化: f (x)制約条件: gi (x) ≥ 0. (i = 1, . . . , N)

用いるなどで回避し，デルタ関数の仮定を用いる
ことは実際によく行われる．
式 (5) の f j を デ ル タ 関 数 で 置 き 換 え れ ば ，∑m

22πiz j ·K
dx である．そこで，
j=1 p j e

ここで，目的関数 f , 制約条件 gi は最適化する変I (K) =

数 x ∈ Rm の 2 次以下の多項式とする．磁気構造I cal ≈ I obs となる p j ∈ C を求めるには以下の 2

解析も位相回復も，各座標の密度が未知という状

次計画問題をとけばよい:況で，数学上の定式がほぼ同じため，同様に 2 次計画問題に帰着できる．
結晶構造は式 (3) で表現するが，z1 , . . . , zm ∈R3 /L は原子座標の候補でもよい．なぜなら，求めた密度 p j がほぼ 0 に等しければ，座標 z j + L に

cal


N
∑

I obs (K s ) − I cal (K s )



最小化:



Err[I obs (K s )]


s=1

2 (6)


m

∑


2πiz
·K

cal
j
s

dx .
p je
制約条件: I (K s ) =



j=1

は原子が存在しない，もしくは X 線をほとんど散式 (2) のような最小二乗法ではなく，絶対値の乱せず解析の対象から除外される種類の原子しか和を最小化する効用は色々知られているが，ここないことが分かるからである．

では，SDR を実施した後の最適化問題の変数や制観測値として，強度データ I obs (K1 ), . . . , I obs (KN )約条件の数を抑えることが主な理由になる．

とその見積り誤差 Err[I obs (K s )] が与えられたとす

5.2 設定その 2 ( f j の広がりを考慮する場合)

る．各 I obs (K s ) の理論式 I cal (K s ) は，式 (3) をフー

式 (5) の I cal (K) の以下のフーリエ級数を考える:∫

リエ級数展開したときの各係数 R3 /L f (x)e2πix·Ks dxの絶対値の 2 乗に，未知定数 c > 0 をかけたもの

I (K) = c

m
∑

∫

2

f j (x − z j )e

2πix·K

pj

R3

j=1

√

I (K) =

m
∑
j=1

dx .

の逆格子）である．一般に F cal (x) は，以下で定義される f の自己相関関数に等しい．

∫

f (z) f (z − x)dz

∫
pj

2

f j (x − z j )e

2πix·K

R3

dx . (5)

I obs から求められるのは新しくおいた方の p j の値なので，一般に強度からは p1 , · · · , pm の比しか求められない．しかし，求めた比に加えて，化学式を用いれば，各 p j の値は得ることができる．
以下ではまず， f j をデルタ関数とみなせる状況で最適化問題を定式化する (節 5.1)．デルタ関数の場合も含む，より一般の場合，分布関数 f j は既知だが，ピークの広がりを表す σ は未知である．このケースを節 5.2 で取り扱う．

=

R3 /L
m ∑
m
∑

各 f j がデルタ関数とみなせないのは，熱振動の影響などで原子の存在する領域が広がっている場
合で，試料の温度を下げる，規格化構造因子17) を

p j1 p j2 ( f j1 ∗ f j2 )(x − z j1 + z j2 ).

j1 =1 j2 =1

ただし， f j1 ∗ f j2 は実数値関数 f j1 と f j2 の自己相関関数 f j1 ∗ f j2 (x) =

∫

R3 /L

f j1 (z) f j2 (z − x)dz である．

f j1 , f j2 が原点で最大となるピーク形状の確率分布であることに注意すれば， f j1 ∗ f j2 (x) も式の形から，同様の確率分布になる．このことから F cal (x)は，差ベクトル z j1 − z j2 の各座標に積分面積がp j1 p j2 となるピークの極大点が存在する多峰ピー
ク関数になる．

F cal は，観測値 I obs を I cal の代わりに使っても近い関数を得ることができる．

5.1 設定その 1 ( f j はデルタ関数とみなせる場合)

4

I cal (K)e−2πix·K .

K∈L∗

cp j を p j と置きなおせば結局，
cal∑

ただし，L∗ = {k ∈ R3 : l · k ∈ Z for any l ∈ L} （L

に等しい:
cal

F cal (x) =F obs (x) =N∑

I obs (K s )e−2πix·Ks .

s=1

F obs に表れるピークは，節 2 で座標ごとにグルー
最小化: C • X







制約条件: Ak • X = bk (k = 1, . . . , N), (10)







X ⪰ 0.

プ分けした差ベクトルの集合 P1 , . . . , Pt のどれか１つに対応し，P s の座標に観測されるピークの積
分面積 Iobs
s が，式 (1) の I s の観測値を与えるので， p j を強度データから求める問題が，節 (2) の密度割り当て問題と同じ状況だと，これで分かる．


∑N


最大化: k=1
bk yk






∑N

制約条件: Y + k=1
yk Ak = C,







Y ⪰ 0.

特に， p j ∈ C を変数とする以下の 2 次計画問題に帰着される:


t
cal
∑

Iobs

s − Is


最小化:



Err[Iobs

s ]
s=1

∑




制約条件: Ical
pi p j .

s =



(7)

ただし，Ak , C, X, Y は同じサイズの対称行列で，bk ∈ R である．X ⪰ 0 は, X が半正定値行列という
∑n ∑n
i=1
j=1 Xi j Yi j で(i, j)∈P s

意味で，X • Y = Trace(XY) =6. 半正定値計画緩和 (SDR) の適用ある．式 (10), (11) は Ak , C, bk が共通なら，SDP

問題 (6), (7) はほぼ同じ形なので，以下のようの主問題・双対問題と呼ばれるものになる．
半正定値計画緩和 (SDR) とは最適化問題の制に１つにまとめておく．

t∑
Ikobs − Ikcal




最小化
:


Err[Ikobs ]

k=1





制約条件: Ikcal = pT S k p.

(11)

約条件を緩めて，式 (10), (11) の問題に帰着することである．目的関数・制約条件がともに変(8)数 x = (x1 , . . . , xn ) の多項式となる多項式最適化 (POP) も含め，一般的な方法で実施できるが，

ただし，p ∈ C の各成分は実部と虚部に分け，

QCQP の場合，緩和後の問題の変数と制約条件の数

p∈Rとなるように取り直しておく．また p を

を最も小さくするには，以下のようにすればよい;変数とする 2 次多項式を表すのに，対称行列 S k を

X = ppT とおくと，今，p は m 次の列ベクトルなの用いている．

で，X は階数 1 の m 次対称行列で，pT S k p = S k •Xm2m

問題 (8) が 2 次計画問題になることは一見分かである．逆に X が階数 1 ならば，ある p ∈ Rm にりにくいが，新たな変数 sk , ϵk を， sk = (Ikobs −対し，X = ppT と書けるので，以下の問題は問題

Ikcal )/Err[Ikobs ]，−ϵk ≤ sk ≤ ϵk を満たすとして加え(9) と同値である．

t
∑





ϵk
最小化
:





k=1




制約条件: sk = (I obs − I cal )/Err[I obs ],
(12)
k
k
k





cal


Ik = S k • X, X は階数 1,








−ϵk ≤ sk ≤ ϵk .

た以下の問題は，問題 (8) と同値で，しかも線形式・2 次式しか含まないことから分かる．


t
∑





最小化
:
ϵk





k=1

(9)



制約条件: sk = (Ikobs − Ikcal )/Err[Ikobs ],








Ikcal = pT S k p, −ϵk ≤ sk ≤ ϵk .

次に，制約条件 “X は階数 1” を “X ⪰ 0” に置き換

半正定値計画問題は，典型的には以下の式 (10),えて，緩和する．あとは機械的な手順で問題 (10)(11) のいずれかの形の最適化問題である．最適化の形に書き直せる．緩和後の最適化変数は，対称

変数 X, Y が動く範囲（許容領域）は半正定値対称

行列 X = (Xi j ) の各成分 Xi j (1 ≤ i ≤ j ≤ m) とな

行列からなる凸錐とアフィン線形空間の交わりで，り，元の m 個から m(m + 1)/2 個に増える．

目的関数も線形なので，凸計画問題の一種である．

弱双対定理より，主問題 (10) と双対問題 (11)の最適値は，最適解 X opt , Y opt , yopt において常に

∑N

opt
opt
となる．これらの最適値のk=1 bk yk ≤ C • X

数理科学

NO. 677, NOVEMBER 20195

差は双対ギャップと呼ばれる．
考えている問題では，以下の強双対定理の仮定

7. 応用例

が成立するため，理論上の双対ギャップは 0 とな

以下に登場する結晶学の解析は，原子種・占有率，り，出力される双対ギャップが 0 に近いかどうか

磁気構造は未知だが，原子座標の候補 z1 , . . . , zm はで，大域的収束の確認ができる．

ある程度判明している状況で実施されるものにな

定理 1（強双対定理） 主問題 (10) の目的関数が下に有界で，内点許容解を持つ (すなわち全ての等式
制約条件を満たす正定値な X が存在する) ならば，
最適解 X opt , Y opt , yopt が存在して，

C•X

opt

∑N

opt
k=1 bk yk=となる．双対問題 (11) において同様の仮定が成立する場合も，同じ等式が成り立つ．
一般に緩和を行えば，変数 X が動ける範囲は広がり，最適値は減少する（目的関数を最小化する

る．使用した SDP ソルバーは SDPA18) である．

7.1 原子種・占有率の大域的最適化と解の一意性
この解析で原子座標 z j に割り当てられる密度 p jの値は，以下の c j , o j の積である．化学式から各原子種の c j は既知となるため，求めた p j = c j o jの値から，c j , o j を導出できる．

• z j に存在する原子種に固有の値 c j（原子散乱因子から求まる），• 占有率 0 ≤ o j ≤ 1 (zi に実際に原子が存在する場合）ため，緩和問題の最適値が大きいと，元の単位胞の割合）．

最適値は必ずそれ以上の値になる．
「確度 100%で無機結晶構造データベースのうち，原子座標のモデル棄却ができる」とは，問題 (8) のような目

数が 3 ≤ m ≤ 49 であった結晶構造 3572 個について，IS の値をシミュレーションし，2 次計画問題

的関数を取った場合に成立する話である．
最適解も緩和後に変わる可能性はあるが，少（7）に SDR を適用し，p j の割り当てを実施した．

なくとも最適解 X opt の階数が 1 なら，X opt =差ベクトルを同一/異なると判定する際の解像度は) と分解して得られる ±pは，元の問0.01Å と小さく設定したので，複数解が見つかっ

題の大域的最適解になる．元の問題が他にも大域た場合も，観測精度が理由というわけではない．

的最適解を持つなら，内点法の機構より，X の階

最適値および双対ギャップはどのケースでも 0

数 > 1 となることから，解の一意性が証明される．

に近い値に収束し，SDP で得られた最適解 X の

opt

p

opt T

(p

opt

また階数 > 1 でも，代数的に，X

opt

=

∑

opt
i Xi

階数が 1, 2–4, 5 以上となった結晶構造はそれoptと分解して，Xi がそれぞれ階数 1 の緩和問題

ぞれ 3019, 377, 176 個であった． したがって，の大域的最適解にできることがある．このとき，3019/3572 ≈ 84.5% の結晶構造で，密度 p j は一

opt T
opt
Xiopt = popti (pi ) と分解して得られる ±pi は，元意に決まった．それ以外のケースでは密度を割りの問題の大域的最適解を全て与えると考えられる．

当てる方法が複数あったことになる．ただし解の節 (7) で紹介するが，これまで実施した密度割り当ての具体例では，９割以上が X

opt

が階数 1 か，対称性が理由のものも含まれている．
顕著にみられる傾向として，複数解が（たくさ

opt
もしくは階数 1 の大域的最適解 Xi に分解されるん）見つかる結晶構造の多くは，差ベクトル集合

ケースであった．

が R3 /L のある部分群を含むことが分かった．

この高い成功率の理由は，まず目的関数が ≥ 0

7.2 磁気構造解析の大域的最適化と下から抑えられていること，さらに結晶構造解

磁気構造解析で最適化する p j は 3 次ベクトル析の状況では，解の存在証明は必要なく，目的関（磁気モーメント）で，I cal (K) も以下の式になる:

数を 0–1 程度に小さくする解は必ず存在する，ということがあるのでは，と考えている．

I cal (K) =

m (
∑
j=1

pj · K
pj −
K
K·K

)∫
R3

2

f j (x − z j )e2πix·K dx .

ただし |∗|2 の中身は複素数成分の 3 次ベクトルな6ので，各成分の絶対値の 2 乗和が I cal (K) になる．
磁気構造解析では座標数 m の 3 倍の最適化変数があるが，利用できる I obs の数は少ないため（以8. まとめ著者の最近の研究結果から，SDR の機構を用いは 8 個)，既約表現解析て結晶構造解析における構造保証がどのように実が用いられる．すなわち結晶構造の対称性にラン施できるかを紹介した．また，位相回復との関連，

ダウ理論を適用し，磁気構造の対称性を複数個に理論的な解の一意性の調査に使用できることも説絞る．各対称性が要請する p j の線形関係式を用い明した．本稿で触れられなかったフーリエ反復法・てパラメータ数を減らし，それぞれの最適化結果

圧縮センシングに関わる近年の理論研究をまとめを総合して，最も良い対称性・モデルを求める．

たものに，19), 20) がある．

下の例では m = 8 で I

obs

以下に示すのはオークリッジ国立研究所の熱中

SDR は数学の理論に裏打ちされた強力な最適化

性子 3 軸分光計で計測された Nd2 Ir2 O7 の磁気散乱手法であるが，他分野で具体的な応用事例が得ら

強度データを解析した結果である．詳細は省くが，れたことで，使用範囲が拡大することを期待する．

既約表現解析により，磁気構造は 4 つの既約表現

Γ3 , Γ5 , Γ7 , Γ9 に絞られた．試料温度は 1.4K(ケル

参考文献ビン) で，デルタ分布の仮定を用いることができ

1) E. J. Candès, Y. C. Eldar, T. Strohmer and V. Voroninski,
SIAM J. Imaging Sci. 6(1), 199–225 (2013).
2) K. Tomiyasu, R. Oishi-Tomiyasu, M. Matsuda and K.
Matsuhira, Scientific Reports 8:16228 (2018).
3) K. Momma and F. Izumi, J. Appl. Crystallogr. 44, 1272–
1276 (2011).
4) Online dictionary of crystallography, http://reference.
iucr.org/dictionary/R_factor.
5) 大橋 裕二, 有機合成化学協会誌 58(4) 334–338 (2000)
6) E. J. Akutowicz, Trans. Amer. Math. Soc. 83 179–192(1956).
7) E. J. Akutowicz, Proc. Amer. Math. Soc. 8(2) 234–238(1957).
8) M. H. Hayes, Chapter 6 in Image Recovery: Theory andApplication, Academic Press (1987)
9) R. Barakat and G. Newsam, J. Math. Phys. 25(11) 3190–
3193 (1984).
10) R. W. Gerchberg and W. O. Saxton, Optik 35 237–246(1972)
11) J. R. Fienup, Opt. Lett. 3(1) 27–29 (1978)
12) J. Miao, P. Charalambous, J. Kirz and D. Sayre, Nature
400 342–344 (1999).
13) J. Miao, T. Ishikawa, I. K. Robinson, M. M. Murnane,
Science 348(6234) 530–535 (2015).
14) L. Palatinus, Acta Cryst. B69(1) 1–16 (2013).
15) Sheldrick, Acta Cryst. A64(1) 112–122 (2008).
16) 桜井 敏雄, X 線結晶解析 (物理化学選書), 裳華房 (1995).
17) C. Giacovazzo, Phasing in Crystallography: A Modern
Perspective, Oxford University Press (2013).
18) M. Yamashita, K. Fujisawa, K. Nakata, M. Nakata, M.
Fukuda, K. Kobayashi, and K. Goto, Research Report B460, Tokyo Institute of Technology (2010).
19) D. R. Luke, SIAG/OPT Views and News 25(1) 1–5(2017).
20) Y. C. Eldar, N. Hammen and D. G. Mixon, IEEE Signal
Processing Magazine 33(5) 158–162, (2016).

たため，2 次計画問題（6）を解いて，密度の割り当てを実施した．
双対ギャップは 0 に収束し，得られた最適値を自由度（= データ数-独立なパラメータ数) で割った値は，以下のようになった:

Γ3 : 0.75,

Γ5 : 4.8,

Γ7 : 16,

Γ9 : 14.

モデルが観測値を十分に表現できている場合の期待値（カイ二乗値の 1 にあたる）は√2/π ≈ 0.798

だが，Γ3 のみがこれに近い値となった．また，Γ3
optのモデルでは最適解 XΓ3 が以下の階数 1 の行列に収束した．


 (−1.20)2


XΓopt
= 
3−1.20 · 0.19


−1.20 · 0.19 .
(0.19)2 

他の既約表現では目的関数がこれ以上小さくならないこともこの時点で分かるので，Γ3 の大域的最適解が目的関数を最小化する唯一の大域的最適解と結論される．しらみつぶし探索では Γ3 , Γ5 , Γ7 ,

Γ9 の各既約表現に対し 2, 4, 4, 6 次元空間の探索が必要で，何時間かかっても得られなかった結論に，SDPA は同じパソコンを用いた計算で，1 秒もかからずたどりついた．

（とみやす・りょうこ，九州大学）

数理科学

NO. 677, NOVEMBER 20197