Habréの自動微分（HELL）については、
ここと
ここですでに説明さ
れてい
ます 。 この記事では、DelphiのAD（Embarcadero XE2、XE6でテスト済み）の実装と、非線形方程式f（x）= 0およびシステムF（X）= 0を解くための便利なニュートン法のクラスを提案します。しかし、スパースマトリックスを備えた優れたSLAUソルバーを除き、これは見つかりませんでした（以下を参照）。
Delphiの選択はレガシーコードによって決定されますが、C ++では次のように問題を解決できることを最初に注意しましょう。 まず、このタイプのニュートン法（基本的なニュートン法、ブロイデン法）には、
C ++の数値レシピがあり
ます 。 以前は、レシピは純粋なCのみであり、独自のクラスラッパーを作成する必要がありました。 次に、
Autodiff.orgリストからADライブラリの1つを取得できます。 私の経験では、CPPADの使用はFADBADおよびTrilinos :: Sacadoよりも約30％高速ですが、説明から判断すると、最も速いのは新しい
ADEPTライブラリです。 第三に、行列ベクトル演算では、
SLAEを
解決するために個別のライブラリを使用する場合（たとえば、直接法には
SuiteSparceまたは
Pardiso 、反復法には
ITL ）、
uBlas 、時間
テスト済み 、または新しい高速の
Armadilloおよび
blaze-libを使用できます。 統合線形代数ライブラリとSLAE Eigen、
MTL 、
PETScのソルバーを使用することも魅力的
です （
Cには
ニュートンのソルバーもあります）。 私の知る限り、ヘッダーからのトライアド全体は、
Trilinosでのみ完全に実装されてい
ます 。
数値微分の適用について
デリバティブは、分析的および数値的に計算できます。 分析方法には、選択したプログラミング言語の手段で表現された、手動による差別化、シンボリック（Maple、Wolframなど）および直接的な自動差別化が含まれます。
血圧の使用に関する現在の傾向は、1つの単純な理由によって正当化されます-この手法を使用すると、コードの冗長性と重複が排除されます。 もう1つの議論は、たとえば、グリッド法によって非線形微分方程式（システム）を解くとき、F（X）を計算する方法自体は
重要なタスクであるということです。 実際の問題では、不一致F（X）は異なるプログラムレイヤーからの関数呼び出しの重ね合わせによって表され、手動による区別はその可視性を失います。 また、非定常プロセスをモデル化する場合、F（X）が各タイムステップで変化するため、未知のXのベクトル自体も変化する可能性があることに注意してください。ADを使用すると、F（X）ヤコビアンdF（X）/ dX。 事実、残差を計算するとき、中間変数の型を標準のdoubleからHELLの型に変更することを忘れることがあります（そして多くのライブラリーはHELL型の暗黙的なキャストをdoubleに持っています）。 この場合の問題は、微分の式に誤差があったとしても、反復回数が増えてもニュートン法が収束する可能性があることです。 これは、一部の初期データでは知覚できない場合があり、他のプロセスとのプロセスの相違につながる可能性があります。
したがって、どのような分析の微分df / dxが選択されたとしても、数値微分（f（x + h）-f（x））/ hとの比較で補完することが非常に望ましいです。そうでなければ、コードの正確性について常に疑問が生じます。 当然、数値微分は正しい分析の精度と一致することはありませんが、それでも次の単体テスト手法を推奨できます。 X、F（X）ベクトルおよびdF（X）/ dXマトリックスをテキストファイルにアップロードし、SVNに配置できます。 次に、dF（X）/ dXを数値で取得し、既存のファイルの上にファイルを保存してから、ファイルを相互に視覚的に比較します。 ここでは、値がどのように変化したかを常に確認でき、大きな偏差（分数ではない）でマトリックス要素の座標をローカライズできます-これが分析微分の誤差です。
BPの実装
XE5より前のEmbarcadero Delphiでは、クラスの算術演算をオーバーロードする方法はありませんが、レコード構造
（リンク）にはこのような可能性があります。
TAutoDiff = packed record public class operator Equal(a, b: TAutoDiff): Boolean; class operator Negative(v: TAutoDiff): TAutoDiff; class operator Add(a, b: TAutoDiff): TAutoDiff;
通常、C ++ ADでは、導関数ベクトルの次元は変数であり、コンストラクターで初期化されます。 Delphiでは、構造体の代入演算子をオーバーロード
することはできません（
回避する試みがあります ）。ビット単位のコピーに関連して、導関数のベクトルは固定長で作成する必要があります。 対応する定数TAutoDiff.n_allは、最初に手動で設定する必要があります。
例1
procedure TestAutoDiff_1; var i: integer; x, dy: double; x_ad, y_ad: TAutoDiff; begin x := 4;
現時点では、三角関数とブールシフト演算を除き、ほとんどすべての二項演算子と単項演算子がライブラリにオーバーロードされています。 欠落している関数は、個別に簡単に変更できます。
スパース行列のBP実装
一方では、n_allの固定値は重要な制限です。これは、ベクトルの次元が外部から来る可能性があるためです。 一方、一部のタスクでは、弱めることができます。 事実は、連続体力学の方程式を解くための数値的方法について話す場合、それらで生じる行列はまばらな構造を持っています-古典的な例は、座標（i、j）を持つノードの方程式でのラプラス演算子の「クロススキーム」です2D）5つのノードのみが含まれます：（i、j）、（i-1、j）、（i + 1、j）、（i、j-1）、（i、j + 1）。 アイデアを要約すると、この特定のタスクについて次のことを行う必要があります。
const n_juncs = 5;
解決する問題のノードの総数をNとします。次に、ヤコビ行列にN_all = N * n_junc_vars列があり、そのうちn_allのみが非ゼロです。 TAutoDiff構造内に整数ベクトルn_juncsを挿入し、その各要素がノードのグローバルインデックスを0からN-1に定義する場合、範囲[0、n_all-1]のローカルインデックスを持つdx関数は、範囲[ 0、N_all-1]。 例2は、このようなインターフェイスの使用の詳細を示しています。その利点は、以下のニュートン法を実装するときに最も完全に明らかになります。
例2
procedure TestAutoDiff_2; const N = 100;
次のパートでは、ニュートン型のメソッドのクラスと、スパースSLAEソルバーを選択する問題について検討する予定です。 それまでの間、私は読者に以下を提供します。
- 動作のセマンティクスを使用して、C ++ 11でADを記述してみてください。1）非常に高速に動作するはずです。 2）式テンプレートなしで演算子のオーバーロードを行うことができます。 3）これは初めてです。
- Roslyn CTPでのADの実装に関するコースのアイデア：F（X）の算術式に関するすべての情報を含む構文ツリーをすぐに使用できます。
AutoDiff.pasライブラリファイル unit AutoDiff; interface const SMALL = 1e-12; const n_juncs = 5; type TAutoDiffJuncVector = array[0..n_juncs - 1] of integer; TAutoDiff = packed record const n_junc_vars = 10; const n_all = n_juncs * n_junc_vars; private juncs_offset: TAutoDiffJuncVector;