技 術 論 文

微分値の高精度計算法

High-Precision calculation methods for differential value
林 健太郎＊ 小堀 壮彦＊
Kentaro Hayashi, Takehiko Kobori

数値計算を行うプログラムにおいて微分値を計算する場合、数値微分では桁落ちによる精度劣化を招く。導関数をコード化するのは、独立変数の数や微分の階数が増えるほど誤りを生じやすい。これ
らの欠点を解消する方法として、記号微分と自動微分を紹介する。これらの方法を簡単に利用できるようにするため、適切なインタフェースを持つライブラリを作成することは意義がある。本稿では、
記号微分による方法を試作した C++クラス･ライブラリについても述べる。
In the case of calculating differential values in a numerical calculation program, degradation of
precision caused by cancellation often occurs. The more independent variables or higher differentialorder we have, the more mistakes we tend to make through coding derived functions. In order tosolve these problems, symbolic differentiation and automatic differentiation are introduced. It issignificant to make up some libraries which have appropriate interface to use these methods easily.
This paper also refers to a prototype C++ class library of symbolic differentiation.

1．まえがき桁落ちによる精度劣化は免れない。盪では、人手もしくは数式処理ツールによる計算により導関数を求め、それコンピュータは、もともと数値計算を高速に実行するをプログラムに書き込む手間が必要であり、誤りも生じために開発された機械であった。最初の高級言語であるやすい。独立変数の数や微分の階数が多くなるほど、導FORTRANも、数値計算を用途として、計算したい数関数を表す式は多く複雑になるため、その傾向は顕著と式を記述できるように仕様が決定された。

なる。

現在では、世に出回るプログラムのうち、数値計算用

本稿では、これらの欠点を解消する微分値の計算方法のものは極めて少数となってしまった。コンピュータをとして、記号微分による方法と自動微分による方法を紹用いて数値計算をさせる場合でも、Mathematicaや介する。つまり、この２つの方法は、MATLABなどの市販ツールを利用した方が便利なこと・桁落ちによる精度劣化を生じない。

が多い。しかし、大規模なシミュレーションや高速性が・導関数をプログラムに書き込む必要がない。

要求されるシミュレーション，あるいは高度な制御理論という特長を持つ。もちろん方法自体はプログラムとしを応用したリアルタイム処理システムにおいては、数値て実装する必要があるが、微分値計算のユーザ自身が実計算プログラムの価値は失われていない。

装しなくてすむよう、適切なインタフェースを持つライ

さて、数値計算の入門書では、関数の積分や常微分方ブラリ化を行うことは可能である。

程式の解法を必ず扱っているものであるが、微分値の計

2．記号微分算について詳しく述べられているものは少ない。これは、次のような直感的な方法で間に合わせることが多いためである。

2.1

記号微分の使用例

記号微分（symbolic differentiation）（1）（2）とは、言わ盧 数値微分：例えば、ある刻み幅での関数値の差分をその刻み幅で割る。
盪 独立変数の式として表された関数であれば、導関数をプログラムに書き込む。

ば数式処理による計算方法である。実際のコード例を示したほうがわかりやすいと思うので、まずは筆者らの試作したC++ によるクラス･ライブラリの使用例（図１）を示して説明する。

盧の数値微分は、いろいろ工夫はあるけれど、同程度のこのクラス･ライブラリの使用には、C++言語につい数値の差をとって小さい値で割ることに変わりはなく、ての詳細を知る必要はない。C++は「Cよりも優れてい24＊つくば事業部第一技術部

るC」（3）でもあり、オブジェクト指向技術に深く立ち入

x*sin（y）と表せるとき、その構造は図２に示すように表ることなく使用できる。

現することができる。ここで、木構造の葉を三角形で、（a）,（b）は、x, yが記号微分の独立変数であることノードを円で表している。図２からわかるように、木構の宣言である。（）内の文字列は、デバッグ用のもので造の葉には独立変数または定数が相当し、ノードには加あり、本質的ではない。（c）でuは記号微分の変数であ

法，乗法などの演算や関数が相当する。ノードは、演算り、x*sin（y）と表せることを宣言する。SinはSDV-や関数の引数の数だけの枝をもつ。独立変数に値が与えariableを引数とし、SDVariableを戻り値とする関数でられたとき、付随した木構造に従って計算していけば、ある。また、多重定義した演算子*は、SDVariable同士最終的にその変数の値を得ることができる。の乗算を行う。独立変数も変数として扱えるよう、記号微分における微分操作とは、この木構造を別の木

SDIndependentVariableからSDVariableへの暗黙の型変構造に変形することに他ならない。この変形の方法は、換を使用できるように設計しているので、（c）の右辺の微分公式に応じて次のとおり再帰的に定めればよい。

ように書くことができる。（d）は変数dudxが、uをxで盧 ノードにおいては、その演算や関数の微分に応偏微分したものであることの宣言である。Differentiate

じて変形する。図３にノードの変形の例を示す。

はクラスSDVariableのメンバ関数であり、SDIndepen-

四角形で囲まれた部分は、それより下（葉に近dentVariableを引数とする。（e）はさらにdudxをyで偏い部分）の木構造全体を表している。箱の中の微分した変数をdu2dxdyと宣言していることを示す。

微分記号は、その箱の中の木構造に対して微分（f）,（g）で独立変数x, yに対し、それぞれdouble型操作による変形を施すことを意味している。の値x0, y0を設定する。SetValue はクラスSDIndepen-盪 葉における微分操作では、定数は全て定数０にdentVariable のメンバ関数であり、doubleを引数とする。（h）〜（j）で、クラスSDVariable のメンバ関数

uCalculateValue の戻り値として、x0, y0の値に応じた変+数u, dudx, du2dxdy の値を得ることができる。つまり、,，た値が入る。u＝x*sin(y)の偏導関数

x2を計算しsin

xy,をコードとして記述しなくても、偏微分係

図２数を計算することができるのである。
2.2**u̲x0y0, dudx̲x0y0, du2dxdy̲x0y0 にはそれぞれ+

記号微分のデータ構造とアルゴリズム

,u,xv

u

値を精度よく計算することが目的なので、より簡単な数+,,x

本節では、試作した記号微分の設計結果として、そのデータ構造とアルゴリズムについて述べる。微分係数の式２*x＋x*sin(y)の木構造,u,v+,x,xu + v式に変形する簡約化などの本格的な数式処理を行うことは、設計では考慮していない。

,u,xx.SetValue（ x0）;
y.SetValue（ y0）;
double u̲x0y0= u.CalculateValue（ ）;double dudx̲x0y0= dudx.CalculateValue（ ）;double du2dxdy̲x0y0 = du2dxdy.CalculateValue（ ）;//（f）//（g）//（h）//（i）//（j）図１

試作したC＋＋クラス・ライブラリにおける記号微分の使用例

,v,xu

v,u,v
*v + u*,x,xu*v//（a）//（b）//（c）//（d）//（e）**vuる。例えば、変数uが独立変数x, y によって u = 2*x +SDIndependentVariable x（ "x"）;SDIndependentVariable y（ "y"）;SDVariable u= x*Sin（ y）;
SDVariable dudx
= u.Differentiate（ x）;
SDVariable du2dxdy = dudx.Differentiate（ y）;

+,,x

*

独立変数の式として表された変数に対し、その式を情報として付加しておく。式の情報は木構造として表現す,v,x,,x

sin*,u,xcosuusin (u)図３,u* cos (u),x微分操作におけるノードの変形の例

MSS技報・Vol.1825

技 術 論 文て簡単に述べる。

置き換える。独立変数は、もしそれが微分する

例えば、変数u, v に対して u*v を計算するときは、独立変数と一致していれば定数１に、そうでな次の（1）〜（3）の左辺の値を全て計算し、その計算式い場合は定数０に置き換える。

は右辺による。

これにより、例えば図２に示した木構造の微分操作によ（1）る変形は、図４のように表される。このように偏微分して得られる変数においても、独立変数の式として表す情（2）報を持っているので、その値を計算したり、再度微分したりすることができる。
この設計では、変数ごとに独立変数の式を情報として（3）持っているので、メモリの使用量が多い。図４に示すように、０を掛けられる式の情報など、無駄な情報も持っている。また、微分の計算では木構造を変形したり、値また、sin（u）を計算するときは、次の（4）〜（6）の左の計算ではいちいち独立変数の値まで戻って計算するな

辺の値を全て計算し、その計算式は右辺による。
（4）

ど、実行時間もかかる。改良の余地はあるが、数値微分や３章で示す自動微分など、他の微分値の計算方法に比（5）べて、計算機資源を浪費することに変わりはない。
しかしながら、プログラム中に書いた式を何度でも気軽に微分できるのが、記号微分の魅力である。数値計算（6）プログラムにおいて、局所的に微分値の計算を組み込む必要がある場合に有用であろう。例えば、ちょっと込み入った式で表されるような、非線型方程式をNewton-

自動微分においては、独立変数x, yやその値x0, y0をRaphson法4）で解きたいとか、曲線や曲面の曲率を求め意識する必要はない。もちろん、独立変数を計算元としたいといった場合がその例である。て選択することは可能で、その場合はとすればよい。

3．自動微分
3.1

ここで説明したのは、２階以上の微分係数は必要でない場合の方法である。より高階の微分係数まで要する場自動微分の計算方法

合は、演算や関数計算の各段階において、その階数まで自動微分（automatic differentiation）（1）（6）（7）では、各変数について、その値と必要な階数までの微分値を組の微分値の計算を要する。高階の微分値の計算式はにしたデータとして扱う。つまり、加算や乗算などの演（2）,（3）,（5）,（6）のように簡単ではない。階数が算やsin, cosなどの関数計算を行う場合には、その値だ高くなるほど、また独立変数の数が増えるほど、必要なけでなく微分値まで計算するのである。したがって計算微分係数の数も増え、微分値の計算は複雑になる。

元の変数（例えば入力変数）に対しては、その値だけで3.2なく必要とする階数までの微分値も設定する必要がある。

自動微分の使用例

自動微分の使用例として、変数の数及び微分の階数を以下、２変数の１階微分係数まで必要とする場合につい,u,xu++**2xxsin+,,x*0y+*x2*11*

sinyx*

0cosy

図４26

図２に示した木構造の微分操作による変形，

ADVariable u;u.SetValue（ u̲x0y0）;
u.SetDerivativeX（dudx̲x0y0）;
u.SetDerivativeY（dudy̲x0y0）;
ADVariable v;
v.SetValue（ v̲x0y0）;
v.SetDerivativeX（dvdx̲x0y0）;
v.SetDerivativeY（dvdy̲x0y0）;//（a）//（b）//（c）//（d）//（e）//（f）//（g）//（h）

ADVariable w= u*Sin（ v）;
double w̲x0y0= w.GetValue（ ）;double dwdx̲x0y0 = w.GetDerivativeX（）;
double dwdy̲x0y0 = w.GetDerivativeY（）;//（i）//（j）//（k）//（l）

合に適用できる、最適化問題（非線型計画法）が挙げられる。最適化問題の数値解法の多くは、目的関数および制約関数の独立変数に関する偏微分係数が必要であり、数値微分で間に合わせることもあるが、自動微分による偏微分係数を利用すれば、桁落ちによる精度劣化の悪影響を避けることができる。
自動微分の計算を高速化するための発展的な方法として、関数値を計算するプログラムから偏導関数の値を計算するプログラムを生成することも行われている （8）（9）。

図５

自動微分の使用例「自動微分」、"automatic differentiation" でインターネット上を検索すれば、本稿より詳しい情報が得られるの２変数、１階まで限定する場合について述べる。このよ

で、興味を持った方はお試しになられたい。

うに簡単な場合は、変数の値および１階微分係数の値を4．むすびデータ･メンバとして持つクラス（図５のADVariable）を１つだけ作ればよい。演算や関数計算、例えば乗算や桁落ちによる精度劣化を生じない微分値の計算方法とSin関数においては、式（1）〜（3）や式（4）〜（6）して、記号微分による方法と自動微分による方法を紹介に基づいてデータ・メンバの値を計算する。使用例は図した。

５のようになる。

記号微分や自動微分のように、数式に付随して有益な（a）はu が自動微分を行う変数であることの宣言で情報を付け加えていく方法は、数値計算プログラムの新ある。（b）によりuの値としてdouble型のu̲x0y0を設定しい可能性を示すものである。このような方法として、する。（c）,（d）でuをそれぞれx, yで偏微分した微分係

他には「精度保証つき数値計算」（1）（あるいは「数値計数としてdouble型のdudx̲x0y0, dudy̲x0y0を設定する。

算の品質保証」（9））がある。

SetValue , SetDerivativeX , SetDerivativeYは、クラス

これらの方法では、図１のSDVariableや図５のADVariable のメンバ関数である。（e）〜（h）で変数vADVariableのように、組み込み型でないクラスを使用について同じことを行う。

する。このようなユーザ定義型クラスに対しても、図１（i）で変数wが、u*sin（v）と表せることを宣言する。の（c）や図５の（i）のように計算したい数式になるべ

本節の例において実質的な微分値の計算を行うのは（i）

く沿って書けることが、プログラムの可読性を高めるたの右辺だけである。多重定義した演算子 * は、めに重要である。したがって、これらの方法を実装する

ADVariable同士の乗算として式（1）〜（3）の計算をプログラム言語としては、演算子の多重定義を利用でき実行する。また、SinはADVariableを引数としADVariableる方が有利である。

を戻り値とする関数であり、式（4）〜（6）の計算を実参考文献・参考サイト

行する。
（j）〜（l）で、クラスADVariable のメンバ関数
GetValue , GetDerivativeX , GetDerivativeY の戻り値として、wの値，wのxによる偏微分係数の値，wのyに計算機で微分方程式の厳密解を求める 大石進一
http://www.sci.waseda.ac.jp/journal/vol1/no1/

よる偏微分係数の値がそれぞれ得られる。

oishi/oish00.htm

ここで示した使用例では、独立変数の数と微分の階数盧盪

平野拓一

http://www-antenna.ee.titech.ac.jp/を、それぞれ２と１に固定することを前提にしてクラス

˜hira/hobby/symbolic/index.html

ADVariableを設計した。独立変数の数と微分の階数を固定しない自動微分の設計も考えられるが、インタフェ

記号数式処理蘯Bjarne Stroustrup著，長尾高弘訳 「プログラミンースに工夫が必要となると思われる。高階の微分値の計

グ言語 C++ 第３版」アジソン･ウェスレイ･パブリ算も必要なため、実装も難しい。

ッシャーズ･ジャパン，ISBN4-7561-1895-X, 1998,
pp.51

シミュレーションなどの計算量の多い数値計算では、計算機資源を浪費する記号微分よりも自動微分の方が有盻

William H. Press, Saul A. Teukolsky, William T.

用である。利用方法としては、例えば、衛星の遷移軌道

Vetterling , Brian P. Flannery 著，丹慶勝市・奥村を消費推進薬がなるべく少なくなるように作成したい場

晴 彦 ・ 佐 藤 俊 郎 ・ 小 林 誠 訳 「 NUMERICAL

MSS技報・Vol.18

27

技 術 論 文

RECIPES in C ［日本語版］」，技術評論社，ISBN487408-560-1, 1993，pp.264
眈同，pp.307

眇M. Berz, "Differential Algebraic Description ofBeam Dynamics to Very High Order", Particle
Accelerators, Vol.24, pp.109-124, 1989

眄Wikipedia
http://en.wikipedia.org/wiki/Automatic̲differentiation

眩

久保田光一，伊理正夫共著「アルゴリズムの自動微分と応用」現代非線形科学シリーズ３，コロナ社
ISBN4-339-02602-6, 1998

眤

福井義成，野寺隆志，久保田光一，戸川隼人著「新数値計算」インターネット時代の数学シリーズ２，
共立出版，ISBN4-320-01641-6，1999，Part 328