Reduced Order Model とシミュレーション計算
名古屋大学 山本章夫
1.プロローグ .......................................................................................................................................2
2.本資料の全体構成 .........................................................................................................................3
3.シミュレーション計算におけるデータの圧縮(次元削減)の概念 ...................................................4
3.1 データ間の相関の利用 ..........................................................................................................4
3.2 Active subspace の概念 ...........................................................................................................7
4.シミュレーションの入出力と Reduced Order Model .......................................................................8
4.1 入力データの次元削減の概要 ..............................................................................................8
4.2 シミュレーションの簡略化(surrogate model)の概要 ..............................................................8
4.3 出力データの次元削減の概要 ..............................................................................................95.固有値分解、特異値分解及び低ランク近似 ...............................................................................10
5.1 固有値分解 ...........................................................................................................................10
5.1.1 固有値問題と対角化 .....................................................................................................10
5.1.2 固有値と固有ベクトルの意味 ....................................................................................... 11
5.2 特異値分解 ...........................................................................................................................13
5.3 低ランク近似 .........................................................................................................................19
6.出力データの次元削減とシミュレーションの簡略化 ...................................................................21
6.1 例 1：水の加熱 ......................................................................................................................21
6.2 例 2：一般的な時系列データ................................................................................................25
7.入力データの次元削減 ................................................................................................................28
7.1 入力データ間の相関を用いた次元削減 .............................................................................28
7.1.1 理論 ................................................................................................................................28
7.1.2 具体例 ............................................................................................................................32
7.2 入力データの感度係数を用いた次元削減 .........................................................................33
7.2.1 概要 ................................................................................................................................33
7.2.2 理論 ................................................................................................................................36
7.2.3 具体例 ............................................................................................................................37
8.エピローグ .....................................................................................................................................39

1

1.プロローグ
20xx 年、あなたは某メーカーのエンジニアとして勤務している。持ち前の能力を生かして、工学システムの性能評価についてバリバリ業務をこなしていたある日、上司に呼ばれた。
「【シミュレーション計算が得意】君、わが社の重要顧客から、例のシステムの性能評価に関する依頼が来た。この性能評価では、システムを使用する環境条件が変化したときに、どのよ
うにシステムの振る舞いが変化するか、徹底的に評価することが求められている。この仕事を頼めるのは、【シミュレーション計算が得意】君しかいない。日本政府からの認可取得の関係で、顧客は急いでいるようだ。1 か月後に検討結果が欲しいと言っている。
」「わかりました。【いろいろな仕事を頼んでくる上司】さん、やってみましょう。
」
あなたは早速抱えている仕事の再スケジューリングを行い、依頼された重要業務を開始した。
着手してすぐにわかったことは、この工学システムは、様々な異なった種類、しかも非常にミクロなスケールから、マクロなスケールまでの物理現象に基づいているものであり、そのシ
ミュレーションは、一筋縄ではいかない、ということであった。また、シミュレーションに入力するパラメータの数および、評価結果として考慮すべき出力パラメータの数も非常に多い。
顧客から求められているのは、システムの環境条件、すなわち入力パラメータを変化させた
ときのシステムの変化(応答)である。このシステムのシミュレーションコードは存在するが、非常に複雑なシステムを扱うことから、計算時間も長い。大型のクラスター計算機を利用しても、一回の計算結果を得るまでに日単位の時間がかかる。すなわち、愚直に入力パラメータを変更し、シミュレーションを実行し、その結果からシステムの応答をまとめる形では、１年かけても終わりそうにない。
そこで、シミュレーションに用いる計算モデルを大幅に簡略化し、試計算を行ってみた。その結果、分かったことは、大幅な簡略化をした計算モデルでは、システムの応答が適切に評価できず、使い物にならない、ということであった。
あなたはハタと困ってしまい、上司と相談することとした。
「【無茶な仕事を振ってくる上司】さん、ということで、どうすればよいか、悩んでいるんです。お知恵を拝借いただけないでしょうか？」「うむ、そうか。コードでメッシュなどを粗くするのはダメだったんだな？」「やってみましたが、ダメでした」「最近、計算時間を抜本的に減少させる手法として、Reduced Order Model という手法が流行っているらしい。そういえば、【この話は大学で聞いているはずだが】君、君の大学の指導教員も手掛けていたのでは？」「そ、そうでしたかね(汗)。
【意外に鋭い上司】さん、一度話を聞いてみることにします」ということで・・・。

2

2.本資料の全体構成
本資料では、Reduced Order Model (ROM)の基礎を学ぶため、以下の章構成としている。
第 3 章では、シミュレーション計算におけるデータ圧縮の概念について、概要を説明する。
また、ROM で重要となる active subspace の概念を説明する。
第 4 章では、シミュレーション計算の入出力データと ROM の関係の概要を解説する。
第 5 章では、ROM で重要となる線形代数の知識、特に固有値分解と特異値分解について説明を行う。また、これらの結果を用いて行列を近似する低ランク近似について解説する。
第 6 章では、シミュレーションの出力データに対して ROM を適用する例について説明する。
ここでは、入力データの変化に対して得られた出力データに相関があることに着目し、入力データから出力データを予測するモデルの作成方法を示す。
第 7 章では、シミュレーションの入力データに対して ROM を適用する例を二つ説明する。
一つは入力データ間の相関を利用し、独立な変動が少ない入力データを削減する方法、もう一つは、計算結果に大きな影響を与えない(感度が小さい)入力データを削減する方法である。

3

3.シミュレーション計算におけるデータの圧縮(次元削減)の概念
3.1 データ間の相関の利用
以下に示す二次元(x,y)のデータがある。
(x,y) = (0,0), (1,2), (2, 4), (3, 6), (4, 8), (5, 10), (6, 12), (7, 14), (8, 16), (9, 18)
このデータは、20 個の数字からなっているが、このデータを、もっと簡略に伝える方法はないだろうか？
データをプロットするまでもなく、これらのデータは、y=2x の関係を満たしている。したがって、x=(0, 1, 2, 3, 4, 5, 6, 7, 8, 9), y=2x という表記でも同じ内容のデータを伝えることができる。つまり、もともとは 2 次元のデータを、１次元のデータに置き換えられることを示している。これは、データ間に存在する関係性、すなわち「相関」や「規則性」を用いて、余分な(冗長な)データを削減できることを示している。
それでは、表 3.1～表 3.3 の例はどうであろうか？

4

表 3.1 データ 1X
0.778943
0.632191
0.6694180.911412
0.729727
0.122423
0.170492
0.000968348
0.0726456
0.789154
0.351568

Y
1.557891.26438
1.33884
1.82282
1.459450.244845
0.340984
0.0019367
0.145291
1.578310.703136

表 3.2 データ 2
XY
0.693185
0.508612
0.440180.350574
0.442068
1.55027
0.428305
0.611313
0.519177
0.0112239
0.116691
0.504184
0.965998
0.433146
0.135558
0.522415
0.473424
1.126510.474922
0.524639
0.403969
1.43150.415812
0.113456
0.0332364
1.726710.682517
1.479120.977801
1.08385Z
2.33683
1.89657
2.00826
2.73424
2.18918
0.367268
0.511476
0.00290504
0.217937
2.36746
1.0547Z
3.09678
2.02169
4.42675
2.50754
1.57998
1.358443.76429
1.45153.673292.47405
4.07491
1.47435
3.553135.00578
5.1011(以下省略)
表 3.3 データ 3X

0.319907
0.654223
0.0619923
0.151772
0.445456
0.478713
0.558450.612535
3.90043 10 6
0.399286
0.140148
0.806074
0.0114799

Y
0.236917
0.60694
0.843833
0.415480.850281
0.412266
0.578679
0.179420.000108227
0.015582
0.301477
0.0836212
0.521571(以下省略)

5Z
0.917349
0.451237
0.533013
0.896851
0.28034
0.775159
0.594360.769812
1.
0.916694
0.943117
0.585878
0.853131

2

1.5
1
0
0.0

1.0
0.2

0.5

0.4
0.6
0.8

0.0

図 3.1 データ 1 をプロットしたもの
0.0

0.5
6
4

0.0

1.0
6

2

0.5

4

0
2.0

1.0

2
0
2.0

1.5
1.0

1.5

0.5

1.0

0.0

0.5
0.0

図 3.2 データ 2 をプロットしたもの
0.0

0.5

1.0

0.5

0.5

1.0

0.0

1.0

0.5

1.0

1.0

1.0

0.5
0.5

0.0

1.0
0.5

0.0

1.0

0.0
0.5

0.51.0
0.0
0.5

1.0

1.0
0.5
10.0
1.0

0.5
0.5

0.0

0.0
0.5

0.5
1.0

1.0

図 3.3 データ 3 をプロットしたもの
これらの例は、もともとは三次元のデータであったものを、一次元あるいは二次元で示すことが可能であることを示している。すなわち、もともと三次元のデータを一次元あるいは二次
6

元に削減できることを示している。また、もともとのデータには、不要な次元が含まれていた、と言うこともできる。
「複雑なデータを、データ間の相関を活用して、簡略化する」ことが Reduced Order Model の
一つの本質である。

3.2 Active subspace の概念
図 3.1 の個々のデータは三次元の座標であらわされているものの、図 3.1 より明らかなように、一次元の直線状に並んでいる。すなわち、図 3.1 のデータを表現するためには、直線に沿った一次元を考えれば十分であり、この直線に直交する方向は、このデータを表現するためには全く不要である。
図 3.2 については、データがある平面上に並んでおり、これらのデータを取り扱うためには、
この二次元平面を考えれば十分である。この平面に直交する方向は、考えなくても全く差し支えない。
図 3.3 については、データが球面状に並んでおり、この球面上、すなわち二次元を考えればデータを取り扱うことが出来る。
これらの観察結果は、扱うデータに相関がある場合には、考慮しなくてもよい方向(次元)が存在することを示している。逆に、特定の方向は、データを表すために不可欠であるといえる。
図 3.1 の場合、直線の方向はデータを表すのに不可欠であり、図 3.2 の場合には、平面を構成
する 2 方向が不可欠である。図 3.3 の場合は、球面上のデータを表す極角と方位角の 2 つが不可欠となる。
このように、データを表すために不可欠(あるいは重要)な方向(次元)を active subspace と呼ぶ。Subspace は、部分空間と訳されるが、例えばある三次元空間を考えた場合、その中にある点(0 次元)、線(1 次元)、平面(2 次元)は、それぞれ、元の三次元空間の部分空間を構成することになる。このような部分空間は、無数に考えることができるが、そのうち、データを取り扱うために必要な部分空間を抜き出したものを active(有効な) subspace(部分空間)と呼んでいる。
シミュレーション計算を簡略化するためには、この active subspace をどのように見出すかがキーポイントとなる。

7

4.シミュレーションの入出力と Reduced Order Model
シミュレーション解析は、おおざっぱに言うと、入力データを計算モデルに与え、出力データを得る作業であるといえる。シミュレーションにおいては、図 4.1 に示す形で Reduced Order
Model を使用することができる。本章では、それぞれの項目について、概要を述べる。なお、
出力データの次元削減及びシミュレーションの簡略化については第 6 章で、入力データの次元削減については第 7 章で、より詳細な説明を行う。
入力データの次元削減

シミュレー
ション

シミュレーションの簡略化(Surrogate model)

出力データの次元削減

図 4.1 シミュレーション計算の概要と Reduced Order Model

4.1 入力データの次元削減の概要
4.2 シミュレーションの簡略化(surrogate model)の概要
一般にシミュレーションモデルは、実際の物理現象を模擬したモデルを設定し、数値計算を行う。例えば、原子炉の核特性解析では、中性子の振る舞いを Boltzmann 輸送方程式で記述し、この Boltzmann 輸送方程式の解を各種の数値解法により求める。あるいは、炉心損傷を含むよ
うな過酷事故の進展を評価する場合には、燃料溶融、燃料と被覆管の共晶反応、原子炉容器破損、コアコンクリート反応など、様々な物理現象を含むモデルを構築し、数値計算を行う。
このような数値計算は、対象とするシステムが複雑であるほど、一般的に長い計算時間を要する。このような場合、入力データを変えながら多数回シミュレーション計算を行うことは、
現実問題として難しい場合が多い。
しかし、入力データの変化と、計算結果である出力データの変化には、高い相関がみられる
ことが多い。例えば、燃料集合体の濃縮度に対する燃料集合体の無限増倍率の依存性を予測することを考える。計算精度の高い連続エネルギーモンテカルロコードを用い、統計誤差を非常に小さくなるようにして計算を実施すると、長い計算時間を要する。そのため、濃縮度を変えながら、数多くの計算をすることはかなりの計算時間を要するかもしれない。しかしながら、濃縮度と無限増倍率には、非常に良い相関があるため、代表的な濃縮度について連続エネルギ
ーモンテカルロコードの計算を実施し、計算していない濃縮度については、代表点の濃縮度における無限増倍率を内挿することで、実用上、十分な精度で無限増倍率を予測できると考えられる。
これは、シミュレーションモデルの簡略化(surrogate model)の一例である。Surrogate とは、「代8

理」を意味する言葉であり、surrogate model は、「代理モデル」すなわち、詳細なシミュレーションコードの「代わり」になるモデルであるといえる。
Surrogate model は、おおざっぱに言うと、入力データと詳細なシミュレーションの計算結果の相関(関係)を用い、入力データを用いて、出力データの内挿を行うものである。Surrogate model
は、メッシュ幅を広くしたり、考慮する物理モデルの数を減らしたりするなど、計算モデルの簡略化を行わないところに特徴がある。従って、入力データと計算結果の相関が高く、入力データから、高精度に計算結果が内挿できる場合、surrogate model は、低い計算コストで、高精度な結果を得ることが可能である。逆に、入力データと出力データの相関の発見が困難、両者の関係が極めて複雑で内挿モデルを作成することが困難な場合には、surrogate model の適用は制限される。

4.3 出力データの次元削減の概要
一般的に、シミュレーション計算においては、多数のデータが出力される。例えば、原子炉内の燃料棒ごとの三次元の出力分布を考える。燃料棒本数は数万本(104)のオーダーであることから、軸方向のメッシュ数を 100 すると、原子炉の一つの状態あたり、104×102=106 オーダーの燃料棒出力の計算結果が得られる。原子炉の安全評価では、数百～数千点の原子炉の状態を評価するため、106 オーダーの個数の燃料棒出力のデータが、数百～数千個出力されることになる。
しかし、これらのデータは、まったくの無相関ではない。例えば、制御棒配置、燃焼度や原子炉出力が似通っている場合、燃料棒出力のデータは、かなり似通ったものになることが容易に予想できる。これは、言い換えると、データ間に高い相関があることを示しており、この相関を用いて、得られたデータをより少数のデータであらわすことが可能となる。

9

5.固有値分解、特異値分解及び低ランク近似 1
本章では、データの次元削減を行うために必要となる数学的テクニックである特異値分解及び低ランク近似について解説を行う。

5.1 固有値分解
5.1.1 固有値問題と対角化
�⃗でないベクトル𝑣𝑣⃗が(5.1)式の関係を満ある正方行列𝐀𝐀による変換に対し、ある複素数𝜆𝜆および0たすとき、𝜆𝜆を固有値、ベクトル𝑣𝑣⃗を（縦）固有ベクトルとよぶ。

𝐀𝐀𝑣𝑣⃗ = 𝜆𝜆𝑣𝑣⃗

�𝑣𝑣⃗ ≠ �0⃗�

(5.1)(5.1)式は、行列𝐀𝐀を固有ベクトル𝑣𝑣⃗に作用させると、その固有ベクトルの方向は変化せず、大きさが𝜆𝜆倍（固有値倍）されることを示している。つまり固有値𝜆𝜆は、固有ベクトル𝑣𝑣⃗の行列𝐀𝐀による伸縮の度合を表す値である。

����⃗、
𝑣𝑣2 が存在したとする。ただし、一例として、2 × 2の正方行列𝐀𝐀に対して縦固有ベクトル𝑣𝑣
1 ����⃗この固有ベクトルに対応する固有値は𝜆𝜆1、𝜆𝜆2（𝜆𝜆1 ≠ 𝜆𝜆2）とする。

𝐀𝐀���⃗
𝑣𝑣1 = 𝜆𝜆1 ���⃗𝑣𝑣1
𝐀𝐀���⃗
𝑣𝑣2 = 𝜆𝜆2 ���⃗
𝑣𝑣2

����⃗
𝑣𝑣1 ≠ �0⃗�����⃗
𝑣𝑣2 ≠ �0⃗�

(5.2)
(5.3)

ここで、固有ベクトルからなる行列𝐏𝐏、固有値からなる行列𝚲𝚲を(5.4), (5.5)式で定義する。
(5.4)

𝐏𝐏 = (𝑣𝑣����⃗,
𝑣𝑣2
1 ����⃗)𝜆𝜆𝚲𝚲 = � 1
0

(5.5)

0
�
𝜆𝜆2

行列𝐏𝐏に行列𝐀𝐀を作用させると次式になる。
(5.6)式に、行列𝐏𝐏の逆行列𝐏𝐏

−𝟏𝟏

(5.6)

𝐀𝐀𝐏𝐏 = 𝐏𝐏𝐏𝐏

を両辺の左から作用させる。

𝐏𝐏 −1 𝐀𝐀𝐏𝐏 = 𝐏𝐏 −1 𝐏𝐏𝐏𝐏
𝐏𝐏 −1 𝐀𝐀𝐏𝐏 = 𝐄𝐄𝐄𝐄

(5.7)

𝐏𝐏 −1 𝐀𝐀𝐏𝐏 = 𝚲𝚲����⃗、𝑣𝑣
行列𝐀𝐀に対応する固有ベクトル𝑣𝑣
1 ����⃗2 を使えば、行列𝐀𝐀から(5.7)式のように対角成分だけの行列を取り出すことができる。なお、上記において、𝐄𝐄は、2 × 2の単位行列とする。
同様に、(5.6)式の両辺の右から逆行列𝐏𝐏−1を作用させる。
𝐀𝐀𝐏𝐏𝐏𝐏 −1 = 𝐏𝐏𝐏𝐏𝐏𝐏 −1

1

第 5 章及び第 6 章は、名古屋大学山本研究室の松下柾輝氏の研究会資料をベースに作成され

ている。
10

𝐀𝐀𝐀𝐀 = 𝐏𝐏𝐏𝐏𝐏𝐏 −1(5.8)

𝐀𝐀 = 𝐏𝐏𝐏𝐏𝐏𝐏 −1(5.8)式は、行列𝐀𝐀 を固有ベクトルに関する行列𝐏𝐏、逆行列𝐏𝐏分解できることを示している。

−𝟏𝟏と固有値に関する行列𝚲𝚲の 3 つに5.1.2 固有値と固有ベクトルの意味
固有値と固有ベクトルの本質的な意味を考えるために、具体的な行列を扱って説明していく。
2 次元の XY 座標系上の点を移動させる表現行列𝐀𝐀が次式で与えられるとする。
3 1
𝐀𝐀 = �
�
2 4

(5.9)

1
𝑥𝑥⃗ = � �
0

(5.10)XY 座標系の軸方向のベクトルを以下のように定義する。

0𝑦𝑦⃗ = � �
1

(5.11)1�⃗ = � � = 1𝑥𝑥⃗q
0(5.12)

3 1 1
3𝐀𝐀q�⃗ = �� � � = � � = 3𝑥𝑥⃗ + 2𝑦𝑦⃗
2 4 0
2

(5.13)

3 1 3
11𝐀𝐀(𝐀𝐀q�⃗) = ��⃗� � � = � � = 11x�⃗ + 14y
2 4 2
14

(5.14)

例えば、次の直交座標系上のデータ点q�⃗について考える。
このデータ点に対して、行列𝐀𝐀を作用させる。

ここで、もう一度行列𝐀𝐀を作用させる。

(5.13), (5.14)式からわかるように、XY 座標系では、行列𝐀𝐀による変換を単純に扱うことができない。もし行列𝐀𝐀による変換をn回繰り返す場合、変換後のベクトルを求めるためには行列とベクトルの演算をｎ回行う必要がある。
そこで、行列𝐀𝐀による変換を簡単に扱えるような新たな軸を考える。行列𝐀𝐀の固有値と固有ベ
クトルは以下のよう求められる。
𝜆𝜆𝚲𝚲 = � 1
0

0
2 0
�=�
�
𝜆𝜆2
0 5

1 1
𝐏𝐏 = (𝑣𝑣����⃗,
����⃗)
�
1 𝑣𝑣
2 =�
−1 2

(5.15)(5.16)

ここで、 q�⃗と 𝐀𝐀q
�⃗と 𝐀𝐀(𝐀𝐀q�⃗)の成分を直交座標系ではなく固有ベクトルを使って表現してみる。
2 1
1
1 1
2
1
�⃗ = 𝑥𝑥⃗ = � � = � � + � � = ����⃗q
𝑣𝑣1 + ����⃗
𝑣𝑣
0
3 −1
3 2
3 2
3
11

(5.17)

4 1
1
5 1
2
3
𝐀𝐀q
�⃗ = � � = � � + � � = 2 × ����⃗𝑣𝑣1 + 5 × ����⃗
𝑣𝑣
2
3 −1
3 2
3 2
3

8
2
1
25
����⃗ = �11� = � 1 � + �1� = 22 × ����⃗𝐀𝐀(𝐀𝐀q)𝑣𝑣1 + 52 × ����⃗
𝑣𝑣
14
3 −1
3
3 2
3 2

(5.18)(5.19)

この結果を見ると、あるベクトルに行列𝐀𝐀を作用させることは、固有ベクトル𝑣𝑣
����⃗1 方向に 2����⃗2 方向に 5 倍することに等しい。また、各倍率は(5.15)式の各固有ベクトルに対応する倍、𝑣𝑣固有値の値に等しい。以上のことを図 4.1 に示す。つまり、行列 A の変換は、あるベクトルを各固有ベクトルに分解し、その固有ベクトル上で固有値倍することに等しい。

」
図 5.1 行列 A による変換の固有値と固有ベクトルの解釈
行列𝐀𝐀による変換をn回繰り返すとき、(5.8)式を使うと以下のようになる。
2𝑛𝑛
𝐀𝐀𝑛𝑛 = 𝐏𝐏𝚲𝚲𝐏𝐏 −1 𝐏𝐏𝚲𝚲𝐏𝐏−1 ⋯ 𝐏𝐏𝚲𝚲𝐏𝐏 −1 = 𝐏𝐏𝚲𝚲𝒏𝒏 𝐏𝐏−1 = 𝐏𝐏 �
02𝑛𝑛
𝐀𝐀𝑛𝑛 q�⃗ = 𝐏𝐏 �00� 𝐏𝐏 −1
5𝑛𝑛

2𝑛𝑛 × 25𝑛𝑛 × 1
0
−1
�⃗ =
𝑣𝑣1 +
����⃗
𝑣𝑣2
����⃗
𝑛𝑛 � 𝐏𝐏 q
53
3(5.20)(5.21)

行列𝐀𝐀による変換は固有ベクトルを用いれば固有値で簡単に扱うことができる。言い換えると、行列𝐀𝐀による変換は、対応する各固有ベクトル上で固有値倍する操作である。以下に簡単
な行列 A の変換の解釈を示す。

12

図 5.2 行列 A の変換の解釈

5.2 特異値分解
特異値分解とは、任意の行列を適当な直交行列により対角行列に変換する操作である。固有値分解と同じ操作と考えて良いが、固有値分解が正方行列に対してのみ使用できるのに対し、特異値分解は任意の行数・列数の行列に対して適用可能となる。また、対象とする行列が対称行列(すなわち正方行列)である場合、固有値分解の結果と特異値分解の結果は一致する。特異値分解は、固有値分解の一般的なケースであると理解しても良い。ROM を作成する場合、正方でない行列を扱うことが一般的であり、従って、特異値分解が用いられる。
具体的に、𝑚𝑚 行 𝑡𝑡 列の実数行列 𝐀𝐀を正規直交基底 𝐔𝐔＝(𝑢𝑢�⃗1 , 𝑢𝑢�⃗2 , … 𝑢𝑢�⃗𝑚𝑚 ) と 𝐕𝐕＝(𝑣𝑣⃗1 , ����⃗,𝑣𝑣2 … 𝑣𝑣⃗𝑡𝑡 ) を用いて𝑚𝑚 行 𝑡𝑡 列の対角行列 𝚺𝚺 に変換する。
T(5.22)𝐔𝐔𝐀𝐀𝐕𝐕 = 𝚺𝚺
この式は、前節の(5.7)式に対応している。

ここで、行列𝚺𝚺 は対角行列として以下のように表される。
（m < t の時）
𝜎𝜎10𝚺𝚺 = �
⋮
0

0
𝜎𝜎2
⋱
⋯

⋯ ⋯
0 ⋱
⋱ 0
0 𝜎𝜎𝑚𝑚0⋮�⋮0(5.23)

対角成分は、𝜎𝜎1 > 𝜎𝜎2 > ⋯ > 𝜎𝜎𝑖𝑖 > ⋯ > 𝜎𝜎𝑚𝑚 ≥ 0 (1 ≤ 𝑖𝑖 ≤ m)で表される単調減少する正の値である。𝑚𝑚 は行列A のランクに一致する。

正規直交基底𝐔𝐔、𝐕𝐕 は共にユニタリ行列（転置行列と逆行列が等しい：T𝐔𝐔 = 𝐔𝐔 −1、T𝐕𝐕 = 𝐕𝐕 −1 ）であるので、両辺の左から𝐔𝐔、右から T𝐕𝐕をかけると、次式に変換できる。
𝐔𝐔 T𝐔𝐔𝐀𝐀𝐕𝐕 T𝐕𝐕 = 𝐔𝐔𝐔𝐔 T𝐕𝐕

𝐔𝐔𝐔𝐔 −1 𝐀𝐀𝐕𝐕𝐕𝐕 −1 = 𝐔𝐔𝐔𝐔 T𝐕𝐕
𝐄𝐄𝐀𝐀𝐄𝐄 = 𝐔𝐔𝐔𝐔 T𝐕𝐕
𝐀𝐀 = 𝐔𝐔𝐔𝐔 T𝐕𝐕(5.24)

つまり、行列 𝐀𝐀 は 3 つの行列に分解できる。ここで、行列 𝐔𝐔 を左特異ベクトル、行列 𝚺𝚺 を

特異値、行列 𝐕𝐕 を右特異ベクトルという。(5.24)式は、固有値問題の(5.8)式に対応している。
13

図 5.2 と同様にこの変換についての概念を図 5.3 に示す。内容としては、次元の縮約などあるが固有値分解とほぼ同じである。

図 5.3 行列 A の変換の解釈
特異値分解の具体例を以下の行列𝐁𝐁を用いて示す。
50
50
50
50
11√2
9√2
9√2
11√2
⎛
⎞
−
−
𝐁𝐁 = ⎜ 444
4 ⎟9√2 11√2
9√2
11√2
−
⎝− 4
4
4
4 ⎠

1 1
1
1
2
2 ⎞
1
0
0
⎛2 2
1
1 1
1
1
1
⎜
−
− ⎟
⎛0
⎞ 100 0 0 0
2 2
2⎟
𝐁𝐁 = ⎜√2 √2⎟ � 0
10 0 0� ⎜2
1 1
1
1
⎜
1
1
0
0 1 0
−
− ⎟
0 −
⎜2 2
2
2⎟
⎝
√2 √2⎠
11
1 1⎝2 − 2 − 2 2 ⎠(5.25)(5.26)(5.26)式の 3 つの行列を(5.24)式に対応させて考える。左特異ベクトルの行成分及び右特異ベクトルの列成分は正規化されていて、異なった行成分、列成分で内積をとると 0 となり直交していることが容易に分かる。また、行列𝐁𝐁と特異値 𝚺𝚺のランクは必ず一致しこの計算例では 3 になる。

この分解において、右特異ベクトルを観察してみる。右特異ベクトルを行方向に見ると 4 つのベクトルがあることが分かる。この 4 つのベクトルそれぞれを 4 つの時間ステップからなる変化としてみると、以下のようになる。

14

1 1
1
1
2
2 ⎞
𝑇𝑇
⎛2 2
𝑣𝑣1
����⃗
1
1 1
1
⎜
−
− ⎟
⎛ 𝑇𝑇����⃗
𝑣𝑣2 ⎞
2 2
2⎟
= ⎜2
⎜ 𝑇𝑇����⃗1
1
𝑣𝑣3 ⎟ ⎜1 1
−
− ⎟
𝑇𝑇
2
2⎟
𝑣𝑣4 ⎠ ⎜2 2⎝ ���⃗
11
1 1⎝2 − 2 − 2 2 ⎠(5.27)

����⃗𝟏𝟏 と 𝑻𝑻𝒗𝒗
����⃗𝟐𝟐 の概形
図 5.4 右特異ベクトル 𝑻𝑻𝒗𝒗

����⃗𝟑𝟑 と 𝐓𝐓𝒗𝒗
����⃗𝟒𝟒 の概形
図 5.5 右特異ベクトル 𝐓𝐓𝒗𝒗

ここで、行列 B の特徴を探してみると、1 行目は一定になっている。また、2 行目と 3 行目
を見ると正負が逆になっていて、右に行くほど全体としての値が小さくなっていることが分か
𝑣𝑣1 ,る。まさに、この特徴が特異ベクトルの T����⃗T𝑣𝑣
����⃗2,T𝑣𝑣
����⃗3 にそれぞれ対応している。特異ベクト

𝑣𝑣4 の性質は、行列 B の行方向のデータには存在しない。(5.26)式を見ると対応する特異値がル ���⃗T

𝑣𝑣4 は無視してもよい。
存在しないことが分かる。すなわち、この特異ベクトル T���⃗

ここで、ROM 作成における特異値分解の役割について説明する。簡単のために、状態が時間ステップごとに変化する 2 つのサンプルについて考える。

15

表 5.1 各サンプルの時間ステップごとの変化

𝑡𝑡1

15√2
4

1 サンプル ���⃗
𝑐𝑐1

5√2
4

2 サンプル ���⃗
𝑐𝑐2

𝑡𝑡2

𝑡𝑡3

5√2
4

15√2
4

15√2
4

5√2
4

このデータを行列𝐂𝐂に変換すると以下のようになる。
15√2 5√2 15√2 5√2⎛ 4444 ⎞ = �𝑐𝑐���⃗1 �𝐂𝐂 = ⎜𝑐𝑐2
���⃗5√2 15√2 5√2 15√2⎟
4
4
4 ⎠
⎝ 4

𝑡𝑡4

5√2
4

15√2
4

(5.28)

図 5.6 サンプル 1,2 の時間ステップ変化（行列 C）
ここで、行ベクトルが行列にどの程度寄与しているのかを表す寄与率の大きさは、次式のように定義される。

𝑐𝑐⃗2（行ベクトル��⃗
𝑐𝑐𝑖𝑖 の行列𝐂𝐂に対する寄与率）= � 𝑚𝑚 𝑖𝑖 2

(5.29)∑𝑖𝑖=1 𝑐𝑐⃗𝑖𝑖

ベクトル𝑐𝑐���⃗,𝑐𝑐
1 ���⃗の大きさは等しいので行列𝐂𝐂に対する寄与率の大きさは等しい。
2
表 5.2 ベクトルの寄与率

ベクトル

寄与率 [%]

𝑐𝑐1
���⃗

50
50

𝑐𝑐2
���⃗

ROM の作成の考え方は、この行列𝐂𝐂の傾向を表すベクトルの数をできるだけ減らすことであC の 50％の情報を削ることにる。現在の形では、ベクトル𝑐𝑐���⃗,𝑐𝑐
1 ���⃗の片方を無視することは行列2なる。そこで、寄与率が大きいベクトルと小さいベクトルに変換できれば、ベクトルを減らし
16

た際の行列 C の情報損失は小さくなることが考えられる。そのようなベクトルを探すために特異値分解を使用することができる。
行列 C を特異値分解すると、次式で書き表される。
1 1
1
1
2
2 ⎞
⎛2 2
1
1
1
1 1
1
⎜
⎟
⎛√2 √2 ⎞ 10 0 0 0 2 − 2 2 − 2⎜⎟
𝐂𝐂 = ⎜
��
1
1 ⎟ 0 5 0 0 ⎜1 1
1
1⎟
−
−
−
⎜2 2
2
2⎟
√2⎠
⎝√2
1
1
1 1
⎝2 − 2 − 2 2 ⎠(5.30)(5.30)式における特異値 10 に対応するベクトル T����⃗𝑣𝑣1 と特異値 5 に対応するベクトル 𝑇𝑇����⃗𝑣𝑣2 が特異値分解で得られた新たなベクトルである。
1 1𝑣𝑣
����⃗1
� 𝑇𝑇 � = �2 2
1
1
𝑣𝑣2
����⃗
−
2
2
𝑇𝑇

1 1
2 2 �
1
1
−
2
2

(5.31)

図 5.7 サンプル 1,2 の時間ステップ変化を表す特異ベクトル
ベクトル T����⃗𝑣𝑣1 は、ベクトル𝑐𝑐���⃗,𝑐𝑐
1 ���⃗の全体の傾向を表すベクトルだと考えることができる。ベク
2トル𝑐𝑐���⃗1 ,𝑐𝑐���⃗の平均値が一定なので全体の傾向として時間ステップ変化に対し一定ということを2

𝑣𝑣2 は、ベクトル𝑐𝑐���⃗,𝑐𝑐示しているベクトルである。ベクトル 𝑇𝑇����⃗1 ���⃗の時間ステップ変化ごとの状態変2化率を表すベクトルである。

𝑣𝑣1 に特異特異値は、行列 C による右特異ベクトルの拡大率である。これにより、ベクトル 𝑇𝑇����⃗𝑣𝑣2 に特異値 5 をかけると、軸変換した行列𝐂𝐂′を得ることができる。
値 10、ベクトル T����⃗5 5 5 5����⃗
10 × T����⃗
𝑣𝑣1
𝑐𝑐′
5
5 5
5� = � 1 �
𝐂𝐂′ = ��=
�(5.32)
T−−�����⃗
5 × ����⃗
𝑣𝑣2
𝑐𝑐′
2
2
2 2
2

17

図 5.8 サンプル 1,2 の時間ステップ変化（行列𝐂𝐂′）����⃗1 と𝑐𝑐′�����⃗2 の行列𝐂𝐂′に対する寄与率が与えられる。
(5.29)式より、ベクトル𝑐𝑐′
表 5.3 右特異ベクトルの寄与率
特異ベクトル

寄与率 [%]

����⃗𝑐𝑐′1�����⃗
𝑐𝑐′2

8020

上記の流れをまとめると、行列 C を表現する寄与が等しく与えられているベクトルを、特異値分解により変換し、一方のベクトルの寄与率をできるだけ大きくするようにすることができる。

図 5.9 特異値分解による軸変換の行列 C に対する寄与率の変化
右特異ベクトルは正規直交基底なので、2 乗ノルムは規格化されている。よって、新しく軸を取ったベクトルの寄与率は特異値のみで計算され次式のように与えることができる。なお、
特異値の 2 乗(𝜎𝜎𝑘𝑘 2)は、k 番目の特異値に対応する特異ベクトル方向の分散を表している。従っ2て、∑𝑚𝑚
𝑖𝑖=1 𝜎𝜎𝑖𝑖 はデータの全分散を表しており、寄与率は、データの全分散に対する k 番目の特異値に対応する特異ベクトル方向の分散の割合を表していることになる。

18

（各特異ベクトルの第k番目の直交基底の寄与率）=

𝜎𝜎𝑘𝑘 2𝑚𝑚∑𝑖𝑖=1
𝜎𝜎𝑖𝑖 2（各特異ベクトルの第k番目までの直交基底の累積寄与率）=

5.3 低ランク近似2∑𝑘𝑘𝑖𝑖=1 𝜎𝜎𝑖𝑖

𝑚𝑚∑𝑖𝑖=1
𝜎𝜎𝑖𝑖 2(5.33)(5.34)

前節の特異値分解で説明した通り、右特異ベクトルと左特異ベクトルの直交基底は行列 𝐀𝐀の寄与率が大きい順に並べられている。つまり、寄与率の小さい直交成分は無視しても行列 𝐀𝐀のデータ損失は小さいと考えられる。このように、寄与率の小さい直交成分を無視する技術は、
画像データ圧縮にも使われるものである。

特異値分解後の式を詳細にみると以下のようになっている。rank𝐀𝐀 = mとする。 (m < t の

時)

𝜎𝜎1
0
𝐀𝐀 = (𝑢𝑢�⃗1 , 𝑢𝑢�⃗2 , … 𝑢𝑢�⃗𝑚𝑚 ) �
⋮
0

0
𝜎𝜎2
0
…

… …
0 ⋱
⋱ 0
0 𝜎𝜎𝑚𝑚0⋮ T
� (𝑣𝑣⃗1 , 𝑣𝑣⃗2 , … 𝑣𝑣⃗𝑡𝑡 )⋮0(5.35)

ここで、行列 𝐀𝐀 のデータを特異値𝑘𝑘個（𝑘𝑘 ≤ m）まででおおよそ再現できると仮定すると行列 𝐀𝐀 は以下のように近似できる。特異値は、𝑘𝑘 × 𝑘𝑘の正方行列とする。
𝜎𝜎1 ⋯ 0𝐀𝐀′ = (𝑢𝑢�⃗1 , 𝑢𝑢�⃗2 , … 𝑢𝑢�⃗𝑘𝑘 ) � ⋮ ⋱ ⋮ � T(𝑣𝑣⃗1 , 𝑣𝑣2 , … 𝑣𝑣⃗𝑘𝑘 )0 ⋯ 𝜎𝜎𝑘𝑘m ≥ t の時も同様に近似できる。

(5.36)

具体例として(5.26)式の行列 B に対して低ランク近似を実施する。特異値の成分 100 に対して 1 は十分小さいと考えて低ランク近似を実施し、行列𝐁𝐁 ′を計算する。
1 1
1
1
2
2 ⎞
1
0
0
⎛2 2
1
1 1
1
1
1
100 0 0 0 ⎜
−
− ⎟
0
⎛
⎞
2 2
2⎟
𝐁𝐁 ′ = ⎜√2 √2⎟ � 0
10 0 0� ⎜2
1 1
1
1
⎜
1
1
0
0 0 0
−
− ⎟
0 −
⎜2 2
2
2⎟
⎝
√2 √2⎠
11
1 1⎝2 − 2 − 2 2 ⎠(5.37)(5.37)式の左特異ベクトル 3 列目及び右特異ベクトル 3、4 行目のベクトルは無視できる。
1

0
1

1 1
0
⎛
⎞
100
0
𝐁𝐁 ′ = ⎜� �2 2√2 ⎟ �
1
0
10 1
1
−
0 −
2
2
⎝
√2⎠

19

1 1
2 2 �
1
1
−
2
2

(5.38)

50
50
50
50
10√2 10√2
10√2
10√2
⎛
⎞
−
−
𝐁𝐁 ′ = ⎜ 444
4 ⎟10√2 10√2
10√2 10√2
−
⎝− 4
4
4
4 ⎠
50
50
50
50
9√2
11√2
9√2
11√2
⎛
⎞
−
−
≈⎜ 4
4
4
4 ⎟ = 𝐁𝐁9√2 11√2
9√2
11√2
−
⎝− 4
4
4
4 ⎠

(5.39)

データ損失は多少あるが(5.26)式の行列𝐁𝐁を概ね再現することができる。
さらに、特異値の成分 10 も落としてみる。

1 1
1
1
2
2 ⎞
1
0
0
⎛2 2
1
1 1
1
1
1
100 0 0 0 ⎜
−
− ⎟
0
⎛
⎞
2 2
2⎟
𝐁𝐁 ′′ = ⎜√2 √2⎟ � 0
0 0 0� ⎜ 2
1 1
1
1
⎜
1
1
0
0 0 0
−
− ⎟
0 −
⎜2 2
2
2⎟
⎝
√2 √2⎠
11
1 1⎝2 − 2 − 2 2 ⎠(5.40)(5.40)式の左特異ベクトル 2、3 列目及び右特異ベクトル 2、3、4 行目のベクトルは無視できる。
1
1 1 1 1𝐁𝐁 ′′ = �0� (100) ��(5.41)
2 2 2 2050 50 50 50
𝐁𝐁 ′′ = � 0(5.42)
0
0
0�
0
0
0
0
これより、最も大きな成分は保存されているが、それより小さな成分が失われている。この結果から分かる通り、低ランク近似の度合いが大きいほど、行列𝐁𝐁のデータ損失は大きくなる。
(5.33)式を使うと直交基底の寄与率を求めることができる。
表 5.4 各特異値に対応する直交基底の寄与率
特異値

特異値に対応する直交基底の寄与率 [%]100

99.0001100.99001

0.009920

6.出力データの次元削減とシミュレーションの簡略化
本章では、シミュレーションで得られた出力データの次元削減方法と、次元削減したデータを用いたシミュレーションの簡略化(surrogate model の作成)について説明を行う。一般的なシ
ミュレーションの流れは、入力データ→シミュレーション→出力データであるが、説明のしやすさの観点から、まず出力データの次元削減とシミュレーションの簡略化について説明する。
入力データの次元削減は、次章にて説明する。

6.1 例 1：水の加熱
本節では、低ランク近似を用い、簡単なシミュレーション例について、計算結果の次元削減を行う。また、次元削減した計算結果と入力データ(この例では水の重さ)との相関を用いることで、surrogate model の作成を行う。
やかんに水（初期温度 T0 ＝20 ℃）を注いで加熱させる場合を考える。加熱時間を320秒として、温度変化を80秒ごと観測する。注いだ水の重さを m =1kg、1.5kg、2kg とする。また、水の比熱をC = 4 k J/kg/℃、ガスコンロの仕事量をQ = 1 kW として仕事量はすべて水の温度上昇に消費されると仮定する。

温度を計算する理論式は以下のようになる。

t : 時間 [s]T : t 秒後の水の温度 [℃]

T = T0 +𝑄𝑄 × 𝑡𝑡
𝐶𝐶 × 𝑚𝑚(6.1)

この時の水の温度変化の時系列データは以下の表にまとめることができる。

21

表 6.1 各水の質量における水の温度変化
水の質量[kg]

時間[s]

1
1.5

水の温度変化[℃]

2

0

80

160

240

320

20

40

60

80

100

20

33.3

46.6

60

73.3

20

30

40

50

60

1.0 kg

1.5 kg

2.0 kg

80

160

240

水の温度 [℃]

100
80
60
40
20
0

0

320

時間 [s]

図 6.1 水の温度変化
例えば、水1.2kg の時の温度変化を計算する方法を考える。一般的には、比熱と質量と熱量を使って理論式(6.1)から求める。ここでは、温度変化を計算する別のアプローチとして ROMの理論を用いてみる。
表 6.1 を行列Cに変換すると以下のようになる 2。
20 40
60 80 100
𝐂𝐂 = �20 33.3 46.6 60 73.3�
(6.2)
20 30
40 50 60
この行列𝐂𝐂に対して特異値分解𝐂𝐂 = 𝐔𝐔𝐔𝐔 T𝐕𝐕を施すと特異値は以下のようになる。
208.763000 0𝚺𝚺 = � 0(6.3)
8.4441700 0�
000.0439405 0 0
3 行 3 列目の特異値の成分が十分小さいので無視し、特異値行列を 2×2 行列𝚺𝚺′で近似する。この場合、小さい特異値(0.0439405)に対応する左特異ベクトル及び右特異ベクトルを無視することになるため、左特異ベクトル𝐔𝐔及び右特異ベクトル T𝐕𝐕は以下のように近似される。
−0.710001 0.651234
𝐔𝐔′ = (����⃗𝑢𝑢2 ) = �−0.538726 −0.257268�
(6.4)
𝑢𝑢1 ����⃗−0.453512 −0.713938
2

この例では、3 系列のデータ列(それぞれ 5 点のデータを含む)を行方向に並べているが、列方向に並べることも出来る。その場合、特異値分解で得られる右特異ベクトルと左特異ベク
トルの役割が入れ替わることに留意。ROM 関係の文献を読む際には、データ列がどのように並べられているか、最初に確認しておくと良い。
22

𝑻𝑻
𝑣𝑣1
����⃗
𝐕𝐕′ = � 𝑻𝑻 �
𝑣𝑣2
����⃗

𝐓𝐓

−0.163079 −0.287143 −0.411208 −0.535532 −0.659597=��
−0.757855−0.4661−0.174345 0.114363
0.406118(6.5)

近似された特異値行列𝚺𝚺′、左特異ベクトル𝐔𝐔′及び右特異ベクトル T𝐕𝐕′を用いると、𝐂𝐂 ≈ 𝐂𝐂 ′ =𝐔𝐔′𝚺𝚺′ T𝐕𝐕′となる。

(6.5)式に示された右特異ベクトルをプロットすると、図 6.2 となる。図 6.1 に示された水温の上昇の時間変化の形は、図 6.2 に示された右特異ベクトルの重ね合わせ(重みをつけた後、足����⃗1 )の重みを大きくして、v2（ 𝑻𝑻����⃗𝑣𝑣2 ）し合わせる)で表現できることを示している。例えば、v1( 𝑻𝑻𝑣𝑣と足し合わせると、水温の上昇が緩やかになり、これは水の質量が大きくなることに相当する。

右特異ベクトル

ｖ1
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1

0

ｖ2

80

160

240

320

時間 [s]

図 6.2 (6.5)式に示す右特異ベクトルの形
ここで、𝐂𝐂 ′ = 𝐔𝐔′𝚺𝚺′ T𝐕𝐕′であることから、𝐔𝐔′𝚺𝚺′は、𝐂𝐂 ′、すなわち水温の時間変化を右特異ベクトル
T𝐕𝐕′で展開した場合の展開係数であるとみることが出来る。この展開係数は左特異ベクトルと特異値を用いて以下のように表される。
−0.710001 0.651234
208.7630����⃗1 ����⃗𝛃𝛃 = �𝛽𝛽�
𝛽𝛽2 � = �−0.538726 −0.257268� �
08.44417−0.453512 −0.713938−148.222 5.49912= �−112.466 −2.17241�
−94.6766 −6.02861
入力パラメータである水の質量と、この展開係数の相関を観察してみる。

23

(6.6)

8

-90

6

-100

4

展開係数β2

展開係数β1

-80

-110
-120
-130
-140
-150
-160

2
0
-2
-4
-6

1

1.2

1.4

1.6

1.8

-8

2

1

1.2

1.4

水の質量[kg]

展開係数𝛽𝛽11.61.8

水の質量[kg]

展開係数𝛽𝛽2

図 6.3 展開係数の水の質量依存性

𝑥𝑥を水の質量とし、図 5.2 の関係を近似曲線で表すと以下のようになる。
𝛽𝛽1 = −35.933𝑥𝑥 2 + 161.345𝑥𝑥 − 273.634

(6.7)
(6.8)

𝛽𝛽2 = 7.63066𝑥𝑥 2 − 34.41971𝑥𝑥 + 32.28817(6.5)式、(6.7)式、(6.8)式を使えば、水の質量 1 ≤ 𝑥𝑥 ≤ 2 の範囲であれば水の温度変化を瞬時に算出することができる。つまり、サンプルデータの数1kg、1.5kg、2kg の 3 つから様々な質
量の水の温度変化を再現することができる。
求めた展開係数の近似式に水の質量 1.2kg を代入して、対応する右特異ベクトル(6.9)式、β𝒗𝒗𝟏𝟏 + 𝜷𝜷𝟐𝟐 𝑻𝑻����⃗𝒗𝒗𝟐𝟐を計算すると入力データの関係を示す(6.7)式、(6.8)式を用いて(6.9)式のように𝜷𝜷𝟏𝟏 𝑻𝑻����⃗と、結果は表 6.2 となる。

(−35.933 × 1.22 + 161.345 × 1.2 − 273.634) ×(−0.163079 −0.287143 −0.411208 −0.535532 −0.659597)+(7.63066 × 1.22 − 34.41971 × 1.2 + 32.28817) ×(−0.757855 −0.4661 −0.174345 0.114363 0.406118)

表 6.2 水の質量 1.2kg の時の理論式と ROM の温度変化の比較
水の温度変化[℃]

時間[s]

理論式

ROM

0

20

19.99

80

36.6

36.92

160

53.3

53.84

240

70

70.79

320

86.6

87.71

24

(6.9)

2

表 6.2 を見ると、ROM による水の温度変化の再現は概ねできている。ここで、誤差要因としては、①低ランク近似による誤差と②展開係数の近似(内挿)の精度による誤差と考えることができる。①の誤差を小さくするためには、採用する特異値の数を増やせばよいが、扱う必要のあるデータが増加することに留意が必要である。また、②の誤差を小さくする場合、サンプル
データの数を増やすことにより近似の精度を上げることができること期待できる。一方、サン
プルデータの数を増やすためには、その分、シミュレーション計算を実施する必要があるため、計算時間の増加につながることに留意する必要がある。
ROM の大きな特徴は、入力パラメータとそれに対応した出力データのみで作成可能な点である。水の温度変化を理論式から算出するのに時間を要しないため、ROM の恩恵が少ないが、様々な物理現象を考慮したシミュレーションでは 1 サンプルの計算に長時間を要する場合がある。しかし、ROM を使えば少しのサンプル数で多数のサンプルの結果を再現することが可能となり、ROM の優位性が明らかになる。なお、予測する計算条件が ROM の作成に用いた条件から外れる場合、すなわち外挿になる場合、精度が悪化する可能性が高いことに留意が必要である。すなわち、ROM を作成する時点で、想定される入力データ範囲をカバーしておく必要がある。

6.2 例 2：一般的な時系列データ
本節では、より一般的なシミュレーション結果(時系列データ)を特異値分解し、その結果を用いて surrogate model を作成する方法について説明する。
シミュレーションは、様々な入力パラメータをユーザーが指定することで行う。ここで、乱数により入力データを変化させるランダムサンプリングを m 回実施し、解析コードで計算する。m 個の異なる結果が出力される。1 サンプル当たりの出力データ数(ここでは、時系列デー
タと仮定)は、t 個に統一する。m 回の計算から得られた出力データを行列に整形すると m 行
t 列の行列 A となる。

(6.10)

ここで、特異値分解を実施し行列 A を以下のように分解する。(𝑚𝑚 < 𝑡𝑡)

25

(6.11)

左特異ベクトル𝐔𝐔はm行m列の列直交行列：𝐔𝐔＝(𝑢𝑢�⃗1 , 𝑢𝑢�⃗2 , … 𝑢𝑢�⃗𝑚𝑚 )

特異値はm行t列： 𝚺𝚺

𝑇𝑇

𝑣𝑣1
����⃗
𝑇𝑇 ⎞
⎛
𝑣𝑣2
右特異ベクトルはt行t列の転置直交行列： 𝐓𝐓𝐕𝐕＝ ⎜ ����⃗
⋮ ⎟
⎝

𝑇𝑇

���⃗𝑡𝑡𝑣𝑣⎠

同一の解析結果は存在しないと仮定するのでrank𝐀𝐀 = m 、𝑢𝑢�⃗𝑖𝑖 はm行1列、𝑇𝑇𝑣𝑣⃗𝑖𝑖 は 1 行t列である。

ここで、低ランク近似を実施し、𝛼𝛼𝑘𝑘+1 ⋯ 𝛼𝛼𝑚𝑚 (1 ≤ 𝑘𝑘 ≤ 𝑚𝑚)が𝛼𝛼1 ⋯ 𝛼𝛼𝑘𝑘 と比べて小さければm − 𝑘𝑘個の項を無視する近似ができる。

(6.12)

左特異ベクトル𝐔𝐔′はm行k列の列直交行列：𝐔𝐔′＝(𝑢𝑢�⃗1 , 𝑢𝑢�⃗2 , … 𝑢𝑢�⃗𝑘𝑘 )
特異値はk行k列： 𝚺𝚺′

𝑇𝑇

𝑣𝑣1
����⃗
𝑇𝑇 ⎞
⎛
𝑣𝑣
����⃗
右特異ベクトルはk行t列の転置行直交行列： 𝐓𝐓𝐕𝐕′＝ ⎜ 2 ⎟
⋮
⎝

𝑇𝑇

����⃗
𝑣𝑣𝑘𝑘

⎠

ただし、rank𝐀𝐀 = k (1 ≤ 𝑘𝑘 ≪ 𝑚𝑚) 、𝑢𝑢�⃗𝑖𝑖 はm行1列、 𝑇𝑇𝑣𝑣⃗𝑖𝑖 は 1 行t列である。

この時、特異値に対応する右特異ベクトルの行成分は時系列データの基底と考えることができる。すると、あるサンプルの時系列データ 𝑥𝑥⃗ は右特異ベクトルの基底で展開できる。
𝑘𝑘

𝑥𝑥⃗ = � 𝛽𝛽𝑖𝑖 𝑇𝑇���⃗𝑣𝑣𝚤𝚤(6.13)

𝑖𝑖=0

基底ベクトル 𝑇𝑇���⃗𝑣𝑣𝚤𝚤 に対応する展開係数：𝛽𝛽𝑖𝑖(6.13)式における展開係数𝛽𝛽𝑖𝑖 は、左特異ベクトルと特異値の積から計算できる。展開係数𝛽𝛽𝑖𝑖 と入力パラメータの相関を用いて、入力パラメータから展開係数𝛽𝛽𝑖𝑖 の値を推定することによって未知の時系列データを再現することができる。

以上の要点をまとめると、以下のようになる。
26

①入力パラメータをランダムサンプリングし、シミュレーションコードを用いて計算する。
②得られた出力データを行列にまとめる。
③その行列に対して特異値分解を用いた低ランク近似を実施する。
④左特異ベクトルと特異値の積である展開係数と入力パラメータの相関を用いて ROM を作成する。
⑤ROM を用いて未知の入力パラメータの時の出力データを再現することができる。

27

7.入力データの次元削減
本章では、入力データ間の相関及び入力データの感度を用いた次元削減について説明する。

7.1 入力データ間の相関を用いた次元削減
本節では、入力データをランダムにサンプリングし、シミュレーションを実施、出力データを統計処理することで、入力データの不確かさが出力データに与える影響を定量評価するラン
ダムサンプリング法を例に取り、入力データ間の相関を用いた次元削減を説明する 3。

7.1.1 理論
(多変量正規乱数)
ここでは多変量正規乱数について説明する。まず、一変数の正規乱数について考える。平均0、分散 1(標準偏差 1)の正規分布を標準正規分布と呼び、これを N(0,1)で表すとすると、任意の正規分布に従う確率変数は N(0,1)に従う変数を用いて次式で表すことができる。
(7.1)

𝑥𝑥 = 𝜎𝜎 ∙ 𝑧𝑧 + 𝜇𝜇

ここで、σ は標準偏差、μ は平均値である。したがって、何らかの方法により標準正規分布に従う乱数(標準正規乱数)を作成すれば、任意の正規分布 N(μ, σ2)に従う乱数は(7.1)式の z に標準正規乱数を代入することによって作成することができる。
次にこれを多変量に拡張する。まず以下の記号を定義する。
𝑧𝑧1
𝜇𝜇1
𝑥𝑥1
𝑥𝑥2
𝑧𝑧2
𝜇𝜇2
𝑥𝑥⃗ = � … �, 𝑧𝑧⃗ = � … �, 𝜇𝜇⃗ = � … �
𝑥𝑥𝑛𝑛
𝑧𝑧𝑛𝑛
𝜇𝜇𝑛𝑛

(7.2)

ここで、𝑧𝑧⃗を n 次元の各要素(z1, z2,..., zn)が独立に標準正規分布 N(0,1)に従う確率変数であるベクトルとする。このとき、ベクトル𝑧𝑧⃗の従う確率分布のことを多変量正規分布の中でも特に標準多変量正規分布という。また、𝜇𝜇⃗の各要素(μ1, μ2,..., μn)は定数とする。このとき、標準多変量正
規分布に従うベクトル𝑧𝑧⃗を任意の行列 A により一次変換し、定数からなる n 次ベクトル𝜇𝜇⃗を足すことにより得られるベクトルを𝑥𝑥⃗とする。すなわち行列 A を任意の n × n 行列とすると、(7.3)

𝑥𝑥⃗ = 𝐀𝐀𝑧𝑧⃗ + 𝜇𝜇⃗と表すことができる。このとき、(7.3)式により得られる確率変数ベクトル𝑥𝑥⃗は多変量正規分布に従う。すなわち、ベクトル𝑥𝑥⃗の各要素はそれぞれ平均及び分散の異なる正規分布に従う確率変数となる。ここで、例として次のように 2 次元で考える。
𝑧𝑧1
𝜇𝜇1
𝑥𝑥1
𝑎𝑎
𝑥𝑥⃗ = �𝑥𝑥 � , 𝑧𝑧⃗ = �𝑧𝑧 � , 𝜇𝜇⃗ = �𝜇𝜇 � , 𝐀𝐀 = �
𝑐𝑐
2
2
2

このとき、(7.4)式の演算を行うと、

3

𝑏𝑏
�
𝑑𝑑

本節は、名古屋大学山本研究室の渡邊友章氏の修論をベースに作成した。
28

(7.4)

𝑥𝑥1
𝑎𝑎
�𝑥𝑥 � = �𝑐𝑐2

𝜇𝜇1
𝑎𝑎𝑧𝑧 + 𝑏𝑏𝑧𝑧2 + 𝜇𝜇1
𝑏𝑏 𝑧𝑧1� �𝑧𝑧 � + �𝜇𝜇 � = � 1
�𝑐𝑐𝑧𝑧1 + 𝑑𝑑𝑧𝑧2 + 𝜇𝜇2
𝑑𝑑 2
2

(7.5)

となる。(7.1)式を見ると分かるように、(7.5)式の𝑥𝑥1 と𝑥𝑥2 は正規分布に従う変数の和になっていることが分かる。そして、互いに独立した正規分布に従う変数の和は正規分布に従うという基本的性質から、𝑥𝑥1 と𝑥𝑥2 はともに正規分布に従う変数である。これは一般に n 次元の場合でも成���⃗の各要素は全てが正規分布に従う変数となる。
り立ち、ベクトル𝑥𝑥
1���⃗の期待値(平均)と分散(共分ここで、(7.3)式により得られる多変量正規分布に従うベクトル𝑥𝑥
1���⃗の期待値を(7.6)式の記号で表すとする。
散)を調べる。まずベクトル𝑥𝑥1
E[𝑥𝑥1 ]
E[𝑥𝑥2 ]
�𝐄𝐄[𝑥𝑥⃗] = �…
E[𝑥𝑥𝑛𝑛 ](7.6)

このとき、(7.3)式を用いて、行列𝐀𝐀およびベクトル𝜇𝜇⃗の各要素が定数であることと、ベクトル𝑧𝑧⃗の各要素の期待値は 0 であることに注意すると、𝐄𝐄[𝑥𝑥⃗]は次のように変形できる。
𝐄𝐄[𝑥𝑥⃗] = 𝐄𝐄[𝐀𝐀𝑧𝑧⃗ + 𝜇𝜇⃗]= 𝐄𝐄[𝐀𝐀𝑧𝑧⃗] + 𝐄𝐄[𝜇𝜇⃗]= 𝐀𝐀 ∙ 𝐄𝐄[𝑧𝑧⃗] + 𝜇𝜇⃗

(7.7)

= 𝜇𝜇⃗

このように、ベクトル𝑥𝑥⃗の各要素の期待値は、定数ベクトル𝜇𝜇⃗の各要素と等しいことが分かる。
次に、ベクトル𝑥𝑥⃗の分散共分散行列を var[𝑥𝑥⃗]とする。このとき、var[𝑥𝑥⃗]はその定義から次のよ

うに変形することができる。

var[𝑥𝑥⃗] = 𝐄𝐄[(𝑥𝑥⃗ − 𝐸𝐸[𝑥𝑥⃗])(𝑥𝑥⃗ − 𝐸𝐸[𝑥𝑥⃗])𝑇𝑇 ]= 𝐄𝐄[(𝐀𝐀𝑧𝑧⃗ + 𝜇𝜇⃗ − 𝐄𝐄[𝐀𝐀𝑧𝑧⃗ + 𝜇𝜇⃗])(𝐀𝐀𝑧𝑧⃗ + 𝜇𝜇⃗ − 𝐄𝐄[𝐀𝐀𝑧𝑧⃗ + 𝜇𝜇⃗])𝑇𝑇 ]= 𝐄𝐄[(𝐀𝐀𝑧𝑧⃗ + 𝜇𝜇⃗ − 𝜇𝜇⃗)(𝐀𝐀𝑧𝑧⃗ + 𝜇𝜇⃗ − 𝜇𝜇⃗)𝑇𝑇 ]= 𝐄𝐄[(𝐀𝐀𝑧𝑧⃗)(𝐀𝐀𝑧𝑧⃗)𝑇𝑇 ](7.8)= 𝐄𝐄[𝐀𝐀𝑧𝑧⃗𝑧𝑧⃗T 𝐀𝐀T ]= 𝐀𝐀 ∙ 𝐄𝐄[𝑧𝑧⃗𝑧𝑧⃗T ]𝐀𝐀T = 𝐀𝐀 ∙ 𝐯𝐯𝐯𝐯𝐯𝐯[𝑧𝑧⃗]𝐀𝐀T = 𝐀𝐀𝐀𝐀T

なお、(7.8)式の変形では、ベクトル𝑧𝑧⃗の各要素が独立な標準正規分布に従う変数であり、そのた
めその共分散行列が単位行列になることを用いている。(7.8)式より、ベクトル𝑥𝑥⃗の共分散行列
は行列𝐀𝐀とその転置行列 𝐓𝐓𝐀𝐀の積の行列となる。

さて、ここで任意の平均値ベクトルを𝜇𝜇⃗、共分散行列を𝚺𝚺として、それらを満足する多変量正

規分布に従う乱数、すなわち多変量正規乱数を作成することを考える。このとき、(7.1)式に示した一変数の正規乱数の場合と同様に考えると、(7.3)式の𝑧𝑧⃗の要素に互いに独立な標準正規乱数を代入することで、任意の多変量正規乱数が作成可能である。ただ(7.8)式より、(7.3)式で得られる多変量正規乱数の共分散行列は𝐀𝐀 T𝐀𝐀であるため、任意の共分散行列𝚺𝚺を満足する多変量
29

正規乱数を作成するためには、以下の式を満足する必要がある。
𝚺𝚺 = 𝐀𝐀 T𝐀𝐀(7.9)

したがって、(7.3)式の演算を行うためには、共分散行列𝚺𝚺から(7.9)式を満たす行列𝐀𝐀を求める必要がある。
以上より、多変量正規乱数に基づいたランダムサンプリングを実施するには、基本的には(7.3)式を用いればよいが、そのためには(7.9)式を満たす行列𝐀𝐀を求める必要がある。その方法はいくつか存在するが、特異値分解を用いることができる。
(特異値分解の適用)
ここでは、多変量正規乱数に基づくランダムサンプリングを実施するために用いた特異値分解法について改めて説明する。
行列𝚺𝚺を任意の n × p 行列としたとき、次式で表される行列の分解を特異値分解という。
𝚺𝚺 = 𝐔𝐔𝐔𝐔 T𝐕𝐕(7.10)

ここで𝐔𝐔は行列𝚺𝚺の左特異ベクトルを列に持つ n × n 行列、𝐒𝐒は対角成分に行列𝚺𝚺の特異値を持ち、非対角成分が 0 の n × p 行列、𝐕𝐕 T は行列𝚺𝚺の右特異ベクトルを行に持つ p × p 行列である。ここで、特異値および特異ベクトルは次のように言い換えることができる。
行列𝚺𝚺の非零特異値⇔ 行列𝚺𝚺 T𝚺𝚺(または T𝚺𝚺𝚺𝚺)の固有値の平方根

行列𝚺𝚺の左特異ベクトル ⇔ 行列𝚺𝚺 T𝚺𝚺の固有ベクトル行列𝚺𝚺の右特異ベクトル ⇔ 行列 T𝚺𝚺𝚺𝚺の固有ベクトル

したがって、特異値分解とは、ある行列をその特異値および特異ベクトルから成る行列に分解する手法であり、言い換えると、ある行列𝚺𝚺を行列𝚺𝚺 T𝚺𝚺及び T𝚺𝚺𝚺𝚺の固有値と固有ベクトルから成る行列に分解する手法である。一般に(7.10)式の行列𝐒𝐒の対角成分の特異値は大きい順に並べら

れ、それに伴い、𝐔𝐔と𝐕𝐕の特異ベクトルも特異値に対応した順に並べられる。また、特異ベクト
ル(固有ベクトル)の直交性から、𝐔𝐔と𝐕𝐕はどちらもユニタリ行列であり、以下の式が成り立つ。
T

ここで、𝐈𝐈は単位行列を表す。

(7.11)𝐔𝐔𝐔𝐔 = 𝐈𝐈T(7.12)𝐕𝐕𝐕𝐕 = 𝐈𝐈

以上は特異値分解に関する一般的な説明であり、ここからは本研究での利用方法に焦点を当てる。まず上記の説明において、任意の行列𝚺𝚺が共分散行列の場合を考える。このとき、行列𝚺𝚺は共分散行列の性質である対称行列(𝚺𝚺 = T𝚺𝚺)となり、以下が成り立つ。
𝚺𝚺 T𝚺𝚺 = T𝚺𝚺𝚺𝚺(7.13)

この性質により行列𝚺𝚺の左特異ベクトルと右特異ベクトルは一致し、それにより行列𝐔𝐔と行列𝐕𝐕は等しくなる。したがって、行列𝚺𝚺が共分散行列のとき、(7.10)式は次式のように書き直すこと

30

ができる。
𝚺𝚺 = 𝐔𝐔𝐔𝐔 T𝐔𝐔(7.14)

また、行列𝐒𝐒は行列𝚺𝚺の特異値を対角に持つ対角行列であり、特異値は行列𝚺𝚺𝚺𝚺 T(または𝚺𝚺 T 𝚺𝚺)の固

有値の平方根であるという性質から、行列𝐒𝐒の対角成分はすべて非負である。ここで、行列𝚺𝚺の1

特異値の平方根を対角にもつ対角行列を𝐒𝐒2 とすると、(7.14)式は次のように変形できる。
1 1

なお、

1

1

T

1

1

𝚺𝚺 = 𝐔𝐔𝐔𝐔 T𝐔𝐔 = 𝐔𝐔𝐒𝐒2 𝐒𝐒2 T𝐔𝐔 = �𝐔𝐔𝐒𝐒2� �𝐒𝐒2 T𝐔𝐔� = �𝐔𝐔𝐒𝐒2� �𝐒𝐒2 𝐔𝐔� = 𝐀𝐀 T𝐀𝐀
1(7.15)(7.16)

𝐀𝐀 = 𝐔𝐔𝐒𝐒2である。(7.15)式は(7.9)式と同じ形となっており、そして行列𝐀𝐀は(7.16)式で表される。つまり、共分散行列に対して特異値分解を行うことによって行列𝐀𝐀を求めることができる。
行列𝐒𝐒の対角成分は共分散行列𝚺𝚺を構成するデータに対して特殊な意味を持つ。共分散行列𝚺𝚺
に従うデータ点をランダムサンプリングし、プロットすると、ある広がりを持つ「雲」の様なデータ点の集合を作ることが出来る。この「雲」(データ点の集合)は様々な方向に広がりを持つことになる。ここで、行列𝐔𝐔の i 列目のデータからなる特異ベクトルを考える。このとき、この特異ベクトルの方向に対する「雲」の広がりが行列𝐒𝐒の i 番目の対角成分に対応している。具
体的には、特異ベクトル方向におけるデータ点の集合の分散が行列𝐒𝐒の i 番目の対角成分に、標1

準偏差が行列𝐒𝐒2 の i 番目の対角成分に対応する。この点については、7.1.2 節で具体例を用いて説明する。

最終的に、(7.16)式を(7.3)式に代入することで、次式を用いて多変量正規乱数に基づくランダ
ムサンプリングが実施可能である。
1

1(7.17)𝑥𝑥⃗ = 𝐔𝐔𝐒𝐒2 𝑧𝑧⃗ + 𝜇𝜇⃗

先に説明したように、行列𝐒𝐒2 の i 番目の対角要素は、行列𝐔𝐔の i 列目のデータからなる特異ベク
1

トル方向の分散(ばらつき)を示している。従って、行列𝐒𝐒2 の i 番目の対角要素が十分に小さけ

れば、対応する方向への摂動(ランダムサンプリング)は行う必要がない。特異値分解を行った場合の特異値は、降順に並んでいることから、例えば i 番目の対角要素が十分に小さければ、それ以降の特異値は考慮する必要が無い。言い換えると、最初から i-1 番目の方向についてのみ、ランダムサンプリングを行えば良いことになる。つまり、i 番目以降のデータ(次元)を削減することが出来る。
なお、本研究では特異値分解を用いたが、(7.15)式の分解の方法は一意ではない。すなわち、(7.15)式を満たす行列𝐀𝐀は他の分解方法によっても計算可能である。例えば、コレスキー分解という行列分解法を用いることもできる。

31

7.1.2 具体例
ここでは、特異値分解を用いたランダムサンプリング法について、具体例を用いて考察する。
簡単のため 2 次元のパラメータを考える。平均値ベクトル𝜇𝜇⃗は無視して(7.17)式を展開すると、次式のように書ける。

𝑥𝑥1𝜎𝜎
�𝑥𝑥 � = [����⃗𝑢𝑢2 ] � 1
𝑢𝑢1 ����⃗
0
2

𝑧𝑧1
0 𝑧𝑧1
𝑢𝑢1 + 𝜎𝜎2 𝑧𝑧2 ����⃗𝑢𝑢2
� �𝑧𝑧 � = [𝜎𝜎1����⃗𝑢𝑢1 𝜎𝜎2 ����⃗𝑢𝑢2 ] �𝑧𝑧 � = 𝜎𝜎1 𝑧𝑧1 ����⃗𝜎𝜎2 2
2(7.18)

𝑢𝑢1
����⃗が共分散行列𝚺𝚺の特異ベクトルで、
𝜎𝜎12 及び𝜎𝜎22 が共分散行列𝚺𝚺の特異値である。
ここで、����⃗及び𝑢𝑢2(7.18)式を見ると、ランダムサンプリングにより得られる変数ベクトル𝑥𝑥⃗は、特異ベクトルの線形和で表されることが分かる。つまり、ベクトル𝑥𝑥⃗の向きは特異ベクトルの方向の足し合わせ����⃗と𝑢𝑢によって決定する。また、特異値の平方根である σ1 と σ2 も同様に𝑢𝑢
1 ����⃗それぞれに掛かって2

いるが、これらは各方向の固有のばらつきの度合いを表しており、特異値の平方根はそれに対応する特異ベクトルの標準偏差に相当していると考えられる(特異値が分散に相当する)。以上をまとめると、以下のようになる。


特異ベクトル: ランダムサンプリングで与えられるばらつきの向きを決定する。



特異値: 各ばらつきの向き(特異ベクトル)のばらつきの固有の大きさ(変動量の平均)を決定する。

����⃗に掛かり、z����⃗に掛かっていることから、標
そして、標準正規乱数 z1 がベクトル𝑢𝑢
2 はベクトル𝑢𝑢
12

準正規乱数ベクトル𝑧𝑧⃗は各特異ベクトルの向きの変動量を決定しているといえる。

次に、以上の性質を簡単な例によって確認する。以下の二つの共分散行列をそれぞれ与えてランダムサンプリングを実施した。
1 0.8
1 −1
� , 𝚺𝚺𝟐𝟐 = ��𝚺𝚺𝟏𝟏 = �
0.8 1
−1 4

(7.19)

𝚺𝚺𝟏𝟏 は共分散(非対角成分)が正であるため正の相関があるケースで、𝚺𝚺𝟐𝟐 は逆に負の相関があるケースである。このときの𝚺𝚺𝟏𝟏 及び𝚺𝚺𝟐𝟐 の特異値分解の結果として得られる特異ベクトル行列𝐔𝐔𝟏𝟏及び
𝐔𝐔𝟐𝟐 と特異値行列𝐒𝐒𝟏𝟏 及び𝐒𝐒𝟐𝟐 を示す。

0.71 −0.71
1.80 0.00
� , 𝐒𝐒𝟏𝟏 = ��𝐔𝐔𝟏𝟏 = �0.71 0.71
0.00 0.20
2.60 0.00−0.53 0.85
� , 𝐒𝐒 = �
𝐔𝐔𝟐𝟐 = ��
0.00 0.38
0.85 0.53 𝟐𝟐(7.20)

それぞれの共分散行列を用いて、(7.17)式によりサンプル数 3000(標準正規乱数ベクトル z を3000 個作成)としてランダムサンプリングをおこなったときの𝑥𝑥⃗の散布図を図 7.1 に示す。なお、

標準正規乱数はプログラミング言語やアプリケーションによっては初めから実装されていることもあるが、本研究では Box-Muller 法という方法により一様乱数を正規乱数へ変換している。

図 7.1 において、プロットがサンプリングで得られた各サンプル点の位置を表しており、図����⃗及び𝑢𝑢����⃗の向きを表している。図7.1 よ内に引かれた直線 u1 及び u2 はそれぞれ特異ベクトル𝑢𝑢
1232

り、サンプル点は楕円形に集中しており、これは x1 と x2 に相関が与えられていることを意味����⃗となっており、𝑢𝑢����⃗はそしている。そして、サンプル点の最もばらつき(分散)が大きい方向が𝑢𝑢
12

れに直交する向きとなっている。これは、大きな特異値に対応する特異ベクトルの向きほどその向きのばらつきが大きいという、前述の性質を示している。より具体的に言うと、𝚺𝚺𝟏𝟏 の場合、
𝑢𝑢1 = (0.71, 0.71)方向のデータの分散が 1.80 であり、𝑢𝑢����⃗����⃗2 = (0.71, −0.71)方向のデータの分散が0.20 である。それぞれに対応する標準偏差は、√1.8 ≈ 1.34, √0.2 ≈ 0.45となり、これがデータのばらつきを示している。

また、特異ベクトルが直交しているのも重要な性質である。これにより、特異値分解によるランダムサンプリングでは、まず相関を無視できる方向への座標の変換が行われ、その方向に対して独立にサンプリングが行われる。ここでは 2 次元の場合を考えたが、これは多次元の場
合でも以上の性質は同様であり、特異値の大きい順に対応して、特異ベクトルが直交かつ分散の大きい方向を表すこととなる。
6

u2

4
2

2

0

0

-2

-2

-4

-4

-6

-6

-4

-2

u1
u2

4

x2

x2

6

u1

0
x1

2

4

-6

6

(a) 𝚺𝚺𝟏𝟏 の場合

-6

-4

-2

0
x1

2

4

6

(b) 𝚺𝚺𝟐𝟐 の場合

図 7.1 ランダムサンプリングの結果と特異ベクトル

7.2 入力データの感度係数を用いた次元削減
7.2.1 概要
本節では、感度解析を例にとり、入力データの次元削減について、説明する 4。具体的には、ROM の考え方を感度係数評価に適用することで、感度係数評価におけるモデルを単純化(低次元化)し、少ない解析数で近似的な感度係数を評価する手法について解説する 5,6。
4

本節の説明は、名古屋大学山本研究室の横井公洋氏の研究会資料をベースとして作成したものである。
5 Y. Bang, H. Abdel-Khalik, M. Jessee, et al. “Hybrid reduced order modeling for assembly
calculations,” Nucl. Eng. Des., 295, pp. 661–666,(2015).
6 R. Katano, T. Endo, A. Yamamoto, et al., " Estimation of sensitivity coefficients of core characteristics
based on reduced-order modeling using sensitivity matrix of assembly characteristics," J. Nucl. Sci.
33

𝑇𝑇

m 種類の核反応(m 箇所の放射化量)𝑅𝑅�⃗ = 𝑇𝑇(𝑅𝑅1 , … 𝑅𝑅𝑚𝑚 )と n 種類の巨視的核反応断面積Σ�⃗ =(Σ1 , … Σ𝑛𝑛 )について、 j 番目の巨視的核反応断面積に対する i 番目の放射化量の感度係数𝐺𝐺𝑖𝑖𝑖𝑖

を i 行 j 列の要素に持つ m 行 n 列の感度係数行列𝐆𝐆(𝑅𝑅 ← Σ)を用いることにより、放射化量の変化量は次の(7.21)式のように表すことができる。

�⃗∆𝑅𝑅�⃗ = 𝐆𝐆(𝑅𝑅 ← Σ)∆Σ

(7.21)

ここで、式中の変数は以下のものを示す。

∆𝑅𝑅�⃗: 放射化量の変動ベクトル、∆Σ�⃗: 巨視的核反応断面積の変動ベクトル
𝐺𝐺11 ⋯ 𝐺𝐺1𝑛𝑛⋱⋮ �
𝐆𝐆(𝑅𝑅 ← Σ) = � ⋮𝐺𝐺𝑚𝑚1 ⋯ 𝐺𝐺𝑚𝑚𝑚𝑚(7.22)

具体的な例として、巨視的核反応断面積が 3 種、放射化量が 5 種の場合における 5 つの感

1
1
1
2
2
2
3
3
3
度係数ベクトル𝑔𝑔������⃗������⃗������⃗������⃗𝑅𝑅1 = (𝜕𝜕Σ , 𝜕𝜕Σ , 𝜕𝜕Σ )、𝑔𝑔
𝑅𝑅2 = (𝜕𝜕Σ , 𝜕𝜕Σ , 𝜕𝜕Σ )、𝑔𝑔
𝑅𝑅3 = (𝜕𝜕Σ , 𝜕𝜕Σ , 𝜕𝜕Σ )、𝑔𝑔
𝑅𝑅4 =

𝜕𝜕𝑅𝑅

1

(

𝜕𝜕𝑅𝑅4 𝜕𝜕𝑅𝑅4 𝜕𝜕𝑅𝑅4,,

𝜕𝜕Σ1 𝜕𝜕Σ2 𝜕𝜕Σ3)、𝑔𝑔������⃗𝑅𝑅5 = (

𝜕𝜕𝑅𝑅

2

𝜕𝜕𝑅𝑅

3

𝜕𝜕𝑅𝑅5 𝜕𝜕𝑅𝑅5 𝜕𝜕𝑅𝑅5,,

𝜕𝜕Σ1 𝜕𝜕Σ2 𝜕𝜕Σ3

𝜕𝜕𝑅𝑅

1

𝜕𝜕𝑅𝑅

2

𝜕𝜕𝑅𝑅

3

𝜕𝜕𝑅𝑅

1

𝜕𝜕𝑅𝑅

2

𝜕𝜕𝑅𝑅

3

)について考え、ROM のイメージについて説明する。図 7.2が 5 つの感度係数ベクトルを図示したものである。

図 7.2

5 つの感度係数ベクトルのイメージ

仮に、5 つの感度係数ベクトルが青色の領域で示した 2 次元平面上におおよそ存在しているとすると、その 2 次元平面に直交する方向へ入力パラメータを摂動させても、感度係数との内積がほぼ 0 となるため、核特性量の変化は微小となる。感度解析は、出力パラメータに対する入力パラメータの感度が大きいものを特定することを目的としているため、核特性量の変化が微小である方向へ入力データを摂動させても感度係数の評価には意味が無い。そのため、この場合、入力データ(巨視的核反応断面積)は図 7.3 で示されるような 2 次元平面上の直Technol, 54(6), pp.637-647(2017).
34交座標𝑥𝑥1 , 𝑥𝑥2 の方向に摂動させるべきである。

’

’

’
’
’

図 7.3 2 次元平面上の直交座標𝒙𝒙𝟏𝟏 , 𝒙𝒙𝟐𝟐������⃗������⃗������⃗������⃗������⃗図 7.3 の𝑔𝑔
𝑅𝑅1 , 𝑔𝑔
𝑅𝑅2 , 𝑔𝑔
𝑅𝑅3 , 𝑔𝑔
𝑅𝑅4 , 𝑔𝑔
𝑅𝑅5 は、新たな直交座標𝑥𝑥1 , 𝑥𝑥2 を用いて、3 次元座標系におけ
′

′

′

′

′

������⃗,������⃗,������⃗,������⃗,������⃗を表したものである。2次元平面上の任意の感度係数ベクトルはる𝑔𝑔
𝑅𝑅1 𝑔𝑔
𝑅𝑅2 𝑔𝑔
𝑅𝑅3 𝑔𝑔
𝑅𝑅4 𝑔𝑔
𝑅𝑅5������⃗������⃗������⃗������⃗𝑥𝑥1 , 𝑥𝑥2 方向への 2 回の摂動計算により評価することができる。従って、𝑔𝑔
𝑅𝑅1 , 𝑔𝑔
𝑅𝑅2 , 𝑔𝑔
𝑅𝑅3 , 𝑔𝑔
𝑅𝑅4 ,
′

′

′

′

������⃗𝑔𝑔
𝑅𝑅5 はいずれも𝑥𝑥1 , 𝑥𝑥2 方向への 2 回の摂動計算により評価することができる。その後、図 6.3 の′������⃗������⃗������⃗������⃗������⃗ように 2 次元平面上における感度係数ベクトル𝑔𝑔
𝑅𝑅1 , 𝑔𝑔
𝑅𝑅2 , 𝑔𝑔
𝑅𝑅3 , 𝑔𝑔
𝑅𝑅4 , 𝑔𝑔
𝑅𝑅5 を元の 3 次元座標系
′

′

′

′

′

������⃗,������⃗,������⃗,������⃗,������⃗を評価することができる。
に座標変換することで、感度係数ベクトル𝑔𝑔
𝑅𝑅1 𝑔𝑔
𝑅𝑅2 𝑔𝑔
𝑅𝑅3 𝑔𝑔
𝑅𝑅4 𝑔𝑔
𝑅𝑅5������⃗,������⃗,������⃗,������⃗,������⃗なお、感度係数ベクトル𝑔𝑔
𝑅𝑅1 𝑔𝑔
𝑅𝑅2 𝑔𝑔
𝑅𝑅3 𝑔𝑔
𝑅𝑅4 𝑔𝑔
𝑅𝑅5 のうち、2 次元平面に直交する成分は、この手法では考慮できず(無視され)、最終的には、感度係数ベクトルの誤差となる。

摂動方向𝑥𝑥1 , 𝑥𝑥2 がわかっている場合、2 回の摂動計算で感度係数を評価することができる。

ただし、上述の 2 次元平面(以後、Active Subspace、AS と呼ぶ)というのは、もともと感度係������⃗,������⃗,������⃗,������⃗,������⃗が分かっていないと算出できない。すなわち、ROMを用い数ベクトル𝑔𝑔
𝑅𝑅1 𝑔𝑔
𝑅𝑅2 𝑔𝑔
𝑅𝑅3 𝑔𝑔
𝑅𝑅4 𝑔𝑔𝑅𝑅5て低コストで感度係数を評価するためには、目的とする感度係数ベクトルを評価しなければならないという自己撞着に陥る。そこで、一般的には、比較的低計算コストで実施できる簡易モデル等から、近似的な AS を構築する。

’

’

’

座標変換

’
’

図 7.4 座標変換のイメージ図
一般的な核特性解析においては、入出力パラメータ共に例題よりもはるかに次元が多いが、ROM の考え方は次元が増えても同じである。本例題における計算コストとしては、全種類の35

核特性量について感度係数を評価しようとすると、直接法では 3 回、摂動論では 4 回、ROM では 2 回の摂動計算+AS の構築のための計算を実施することとなる。AS の構築を低計算コストで実施することができれば、放射化解析へ ROM を適用することにより、感度解析の計算コストを低減させることが期待できる。

7.2.2 理論
本節では ROM を作成するための具体的な理論について述べる。核特性量の感度係数ベクトルの厳密な AS は、(7.22)式で示される感度係数𝐆𝐆(𝑅𝑅 ← Σ)を特異値分解することにより評価で

きる。

ここで、
𝐔𝐔: m 行 m 列の正規直交行列

𝐆𝐆(𝑅𝑅 ← Σ) = 𝐔𝐔𝐔𝐔 𝐓𝐓𝐕𝐕

(7.23)

𝐃𝐃: m 行 n 列の𝐆𝐆(𝑅𝑅 ← Σ)の特異値を対角成分に持つ対角行列

𝐕𝐕: n 行 n 列の正規直交行列

特異値分解で得られた𝐕𝐕の縦ベクトルが AS の正規直交基底に該当する。この基底方向への入力データの摂動∆Σ�⃗は、その方向への摂動量を成分としたベクトル∆𝛼𝛼⃗を用いて、(7.24)式のようにあらわすことができる。

ΔΣ1∆Σ�⃗ = �ΔΣ
𝑣𝑣1 + ∆𝛼𝛼2 ����⃗𝑣𝑣2 + ⋯ + ∆𝛼𝛼𝑘𝑘 ����⃗𝑣𝑣𝑘𝑘…2 � = ∆𝛼𝛼1 ����⃗
ΔΣ𝑛𝑛
∆𝛼𝛼1
∆𝛼𝛼
= (𝑣𝑣
����⃗,
𝑣𝑣2 … , ����⃗)𝑣𝑣𝑘𝑘 � …2 � = 𝐕𝐕∆𝛼𝛼⃗1 ����⃗,∆𝛼𝛼2(7.24)(7.23)式、(7.24)式を用いることで、放射化量の変動ベクトル∆𝑅𝑅�⃗は次のように変形することができる。ただし、 𝐓𝐓𝐕𝐕𝐕𝐕 = 𝐄𝐄であることを使用した(𝐄𝐄は単位行列)。
∆𝑅𝑅�⃗ = 𝐆𝐆(𝑅𝑅 ← Σ)∆Σ�⃗�⃗= 𝐔𝐔𝐔𝐔 𝐓𝐓𝐕𝐕∆Σ(7.25)= 𝐔𝐔𝐔𝐔 𝐓𝐓𝐕𝐕𝐕𝐕∆𝛼𝛼⃗= 𝐔𝐔𝐔𝐔∆𝛼𝛼⃗
ここで、対角行列𝐃𝐃について特異値が相対的に大きい上から r 個までの成分を残し、それ以外の相対的に十分小さな特異値の成分を 0 とした対角行列𝐃𝐃𝐫𝐫 を考えると、(7.25)式は(7.26)式のように変形できる。

∆𝑅𝑅�⃗ = 𝐔𝐔𝐔𝐔∆𝛼𝛼⃗(7.26)≈ 𝐔𝐔𝐃𝐃𝐫𝐫 ∆𝛼𝛼⃗(7.26)式より、∆𝛼𝛼⃗の最初の r 個の要素以外は放射化量に影響を与えない。つまり、正規直交行列𝐕𝐕の最初から r 列目までの基底の方向にのみ感度を考慮する形になる。なお、以降、正規直交行列𝐕𝐕の最初から r 列目までの基底を取り出した正規直交行列を𝐕𝐕𝐫𝐫 と表す。

入力データに対する感度係数行列𝐆𝐆(𝑅𝑅 ← Σ)を評価するには、基底方向への入力データの摂動36に対する感度係数行列𝐆𝐆(𝑅𝑅 ← 𝛼𝛼)を摂動計算により評価し、(7.27)式で示される座標変換を行う

ことにより、実施可能である。

𝐆𝐆(𝑅𝑅 ← Σ) = 𝐆𝐆(𝑅𝑅 ← 𝛼𝛼) 𝐓𝐓𝐕𝐕𝐫𝐫(7.27)

なお、(7.27)式の座標変換式は以下のような手順により導出される。まず、核特性量の変動ベ
クトル∆𝑅𝑅�⃗は(7.28)式のように変形できる∆𝑅𝑅�⃗ = 𝐆𝐆(𝑅𝑅 ← Σ)∆Σ�⃗≈ 𝐔𝐔𝐃𝐃𝐫𝐫 𝐓𝐓𝐕𝐕𝐫𝐫 ∆Σ�⃗
= 𝐔𝐔𝐃𝐃𝐫𝐫 𝐓𝐓𝐕𝐕𝐫𝐫 𝐕𝐕r ∆𝛼𝛼⃗
= 𝐔𝐔𝐃𝐃𝐫𝐫 ∆𝛼𝛼⃗(7.28)

また、放射化量の変動ベクトル∆𝑅𝑅�⃗は AS の基底方向への摂動量を成分としたベクトル∆𝛼𝛼⃗と、𝛼𝛼に対する感度係数𝐆𝐆(𝑅𝑅 ← 𝛼𝛼)を用いて(7.29)式のように示される∆𝑅𝑅�⃗ = 𝐆𝐆(𝑅𝑅 ← 𝛼𝛼)∆𝛼𝛼⃗(7.29)(7.28)式と(7.29)式より、𝐆𝐆(𝑅𝑅 ← 𝛼𝛼)は(7.30)式のように示される。
𝐆𝐆(𝑅𝑅 ← 𝛼𝛼) = 𝐔𝐔𝐫𝐫 𝐃𝐃𝐫𝐫

(7.30)

𝑮𝑮(𝑅𝑅 ← Σ) = 𝐔𝐔𝐫𝐫 𝐃𝐃𝐫𝐫 𝐓𝐓𝐕𝐕𝐫𝐫 = 𝐆𝐆(𝑅𝑅 ← 𝛼𝛼) 𝐓𝐓𝐕𝐕𝐫𝐫(7.31)

以上により、入力データに対する核特性量の感度係数ベクトル𝐆𝐆(𝑅𝑅 ← Σ)は、基底方向への入力データの摂動に対する核特性量の感度係数ベクトル𝐆𝐆(𝑅𝑅 ← 𝛼𝛼)を用いて(7.31)式のように変形することができる。

実際には ROM を実施する前に厳密な AS はわかっていないため、近似的な AS を𝐕𝐕として代用する。
7.2.3 具体例
本節では、ROM を簡単なモデルに適用する例について示す。対象とするモデルとしては、(7.32)式で示されるようなモデルを考える。
𝐸𝐸 = 1𝑎𝑎 + 5𝑏𝑏 + 2𝑐𝑐
� 𝐹𝐹 = 2𝑎𝑎 + 10𝑏𝑏 + 5𝑐𝑐
𝐺𝐺 = 5𝑎𝑎 + 20𝑏𝑏 + 8𝑐𝑐

ここで、(7.32)

𝑎𝑎, 𝑏𝑏, 𝑐𝑐:入力パラメータ、𝐸𝐸, 𝐹𝐹, 𝐺𝐺:出力パラメータ(7.1)式で示されるモデルの感度係数行列は(7.33)式となる。ここでは、この感度係数行列を

ROM により求めることを考える。

1 5 2
𝐆𝐆 = �2 10 5�
5 20 8(7.33)(7.33)式を特異値分解すると、Active Subspace の基底は式(7.34)のようになる。ここで、特異値は、順に 25.4365306, 0.9705102, 0.2025406 である。

−0.21432
0.518048
0.828066
𝑣𝑣1 = �−0.90074� , ����⃗����⃗
𝑣𝑣2 = � 0.223083 � , ����⃗𝑣𝑣3 = �−0.37269�
−0.3778−0.82575
0.418816(7.34)

����⃗のみを用いて、感度解析を行う。まず、(7.1)式において、入力データ(a,b, c)を
ここで、基底𝑣𝑣
137

(0,0,0)から(-0.21432, -0.90074, -0.3778)に変化させると、出力データ(E, F, G)は(0, 0, 0)から(5.47362, -11.32504, -22.1088)に変動する。従って、−5.47362
𝐆𝐆(𝑅𝑅 ← 𝛼𝛼) = �−11.32504�
(7.35)−22.1088(7.27)式より、基底𝑣𝑣
����⃗のみを用いた場合、
以下のように近似的な感度係数行列を評価できる。
1
𝐆𝐆(𝑅𝑅 ← 𝛼𝛼) 𝑇𝑇𝑣𝑣
����⃗1 ≈ 𝐆𝐆(𝑅𝑅 ← 𝛼𝛼) 𝐓𝐓𝐕𝐕𝐫𝐫 = 𝐆𝐆(𝑅𝑅 ← Σ)

(7.36)−5.47362
�−11.32504� (−0.21432 −0.90074 −0.3778)−22.1088
1.173100 4.930319 2.067935(7.37)= � 2.427169 10.200938 4.278604�
4.738331 19.914322 8.352711
1 5 2≈ �2 10 5�
5 20 8(7.27)式より、感度がおおよそ再現できていることがわかる。直接法では 3 回の感度解析が必要となるが、基底𝒗𝒗𝟏𝟏 が事前に何らかの方法でわかっているならば、ROM を用いることで 1回の感度解析でおおよその感度係数が評価可能であることがわかる。

38

8.エピローグ
Reduced Order Model の何たるかを理解したあなたは、さっそく上司から仰せつかった仕事に適用してみた。その結果、入力データの感度解析に ROM を適用することで、計算時間を大幅に短縮することに成功し、期限内に仕事を完遂することができたのであった。
「【いい経験をさせてくれた上司】さん、報告書がまとまりました。概要はこちらです。」「どれどれ・・・・。おお、これは素晴らしい。うまくはまると、ROM はこれだけの威力を発揮するのか！【よく ROM の勉強をしたね】君、大変ご苦労だった。良い経験になっただろう。そういえば、関連する国際学会がカリブ海のリゾート地で来年あるそうだ。顧客と連名で今回の成果を発表したらどうかね？」「【おお、太っ腹な上司】さん、いいですね。早速、担当の方にコンタクトしてみます。」
そして 10 ヶ月後・・・。

39