MQL5における切断ニュートン共役勾配(TNC)アルゴリズムの実装
はじめに
実数値関数を最適化する際、MQL5開発者はALGLIBライブラリが提供する最適化機能を利用することが一般的です。ALGLIBのminLBFGSモジュールは、L-BFGS (Limited-memory Broyden-Fletcher-Goldfarb-Shanno)法を堅牢に実装しています。しかしながら、この手法の制約付きバージョンは、MQL5へ移植されたALGLIBには含まれていません。この欠落により、単純なボックス制約付き問題であっても、ユーザーはライブラリ内の他の非線形最適化アルゴリズムを使用せざるを得ない場合があり、それらはこのような用途には過剰な場合があります。 そのため、ALGLIBのminLBFGSと同じように容易にプロジェクトへ組み込むことができ、解析的勾配と境界付近でも安定した数値微分の両方をサポートし、さらにボックス制約にも対応した、実用的なソルバーが求められています。
本記事では、この要求を満たす切断ニュートン共役勾配法(TNC, Truncated Newton Conjugate-Gradient)アルゴリズムの実装を紹介します。アルゴリズムの基礎となる理論を解説するとともに、ローゼンブロック関数の大域最小値を探索する性能を評価します。また、実用例として、TNCをL-BFGSオプティマイザの代替として用いたロジスティック回帰の実装方法についても紹介します。
TNCとは何か?L-BFGSとは何が異なるのか?
切断ニュートン共役勾配法(TNC, Truncated Newton Conjugate-Gradient)法は、ニュートン法の挙動を近似する二次最適化アルゴリズムです。最適化の文脈において、ニュートン法は関数の最小値を求めるために用いられる手法です。一般的な手順では、まず初期推定値を設定し、その後、目的関数の勾配とヘッセ行列を計算して変数を更新します。この処理を収束するまで繰り返します。
関数の勾配(一次微分)は、その関数の傾きに関する情報を提供し、最も急激に増加する方向を示します。一方、ヘッセ行列(二次微分)は、関数の形状、すなわちその点における曲率に関する情報を提供します。ヘッセ行列がすべての方向に対して正の曲率を示す場合、その点は局所最小値(谷底)であることを意味します。すべての方向に対して負の曲率を示す場合、その点は局所最大値(山頂)です。また、ある方向では正、別の方向では負の曲率を示す場合、その点は鞍点であることを意味します。
通常のニュートン法では、関数の曲率を表す二階微分から構成されるヘッセ行列を計算する必要があります。しかし、変数の数が多い問題では、この行列を計算して保持することは計算コストの面で非常に高価になります。TNCは、この問題に対してニュートン法を「切断(truncate)」することで対応します。ニュートン方程式を厳密に解く代わりに、内部ループで共役勾配法(CG, Conjugate Gradient)を用いて探索方向の近似解を求めます。この内部CGループでは、必要となるのはヘッセ行列そのものではなく、ヘッセ行列と任意のベクトルとの積だけです。この積は、勾配の有限差分を用いて効率的に近似できます。その結果、大規模なヘッセ行列を明示的に構築することなく、関数の曲率情報を利用できるようになります。

L-BFGSも準ニュートン法の一種ですが、その基本的な考え方はTNCとは異なります。L-BFGSは、直近数回の反復で得られた勾配とパラメータ位置の履歴を保持し、それらを利用して逆ヘッセ行列を暗黙的に近似します。過去の反復で得られた情報を再利用するため、L-BFGSは1回の反復当たりの計算コストが非常に低く、滑らかで性質の良い目的関数に対しては非常に高速に動作します。このため、L-BFGSは機械学習において事実上の標準的な最適化アルゴリズムとなっています。

TNC APIの概要
TNCオプティマイザを実装したコードは、tnc.mqhヘッダーファイルをインクルードすることで、任意のMQL5プログラムに統合できます。TNC APIの利用は、以下の4つの手順で進めます。
- 目的関数、そのヤコビ行列、および必要に応じてボックス制約を設定する
- オプティマイザを構成する
- オプティマイザを実行する
- 最適化結果を取得する
APIは2つのクラスで構成されています。 CFunctorクラスは目的関数に関するすべての処理を担当し、CTruncNewtonCGクラスはTNCオプティマイザを表します。
目的関数は、CFunctorを基底クラスとする派生クラスを定義して指定します。目的関数そのものは、CFunctorのorig_funメソッドをオーバーライドして実装します。勾配またはヤコビ行列が既知である場合は、grad_funメソッドをオーバーライドして明示的に定義できます。
class CScalarFunc:public CFunctor { public: CScalarFunc(void) { } ~CScalarFunc(void) { } virtual double orig_fun(vector& x) { return np::rosen(x); } virtual vector grad_fun(vector& x) { return np::rosen_gradient(x); } };
CFunctor派生クラスのインスタンスは、目的関数およびその導関数に関する各種設定をおこなうためのインターフェイスを提供します。これには、setBoundsメソッドによるボックス制約の指定が含まれます。このメソッドには行列を渡します。各行は1つの次元を表し、第1列に下限、第2列に上限を、この順序で指定します。grad_funメソッドをオーバーライドした場合は、setGradOptionをGRAD_POINT_CALLABLEオプションで呼び出して有効にしてください。有効にされなかった場合、勾配は有限差分法によって数値的に計算されます。ENUM_DIFF_POINTS列挙体で提供されるオプションを使用すると、grad_funが実装されていない場合、または有効化されていない場合に適用される有限差分法を選択できます。これらのオプションの詳細については、本稿の後半で説明します。
CScalarFunc sf;
matrix bd = {{-2, 2}, {-1, 3}};
sf.setBounds(bd);
sf.setGradOption(GRAD_POINT_CALLABLE);
if(!sf.initialize(init_params))
return;目的関数の設定は、initializeメソッドを呼び出すことで完了します。このメソッドには解の初期推定値を格納したベクトルを渡します。実行時にはtrueが返されることを確認してください。このメソッドがfalseを返した場合は、最適化処理を続行してはいけません。これは、初期設定が正しく完了していないことを示しています。
tnc::CTruncNewtonCG tnc_minim; tnc_minim.SetLoglevel(tnc::TNC_MSG_ALL); tnc_minim.SetMaxCGit(in_maxCGit); tnc_minim.SetMaxFunCalls(in_maxnfeval); tnc_minim.SetEta(in_eta); tnc_minim.SetStepMax(in_stepmx); tnc_minim.SetAccuracy(in_accuracy); tnc_minim.SetFmin(in_fmin); tnc_minim.SetFtol(in_ftol); tnc_minim.SetXtol(in_xtol); tnc_minim.SetPGtol(in_pgtol); tnc_minim.SetRescaleFactor(in_rescale);
続いて、TNCオプティマイザを構成します。ここでは、TNCミニマイザ(最小化器)のさまざまなパラメータを調整できます。TNCオプティマイザの設定に使用する主なメソッドは以下のとおりです。
- SetLoglevelメソッドは、ENUM_TNC_MESSAGE列挙体の値を受け取り、ターミナルへの出力内容を制御します。反復ごとに1行だけ出力する設定から、詳細なスケーリング情報を表示する設定、あるいは最終的な終了理由のみを表示する設定まで選択できます。
- SetMaxCGitメソッドは、内部で実行される共役勾配法の反復回数の上限を設定し、内部ループの効率を維持します。
- SetMaxFunCallsメソッドは、目的関数の評価回数に対する上限を設定します。
- SetFtolメソッドおよびSetXtolメソッドは、それぞれ目的関数値および位置ベクトルの相対変化量に対する収束判定の閾値を設定します。ftol(関数許容誤差)は、連続する反復間における目的関数値の相対変化を監視します。目的関数値の減少量がこの閾値を下回ると、オプティマイザは、これ以上計算を続けても改善はごくわずかであると判断して処理を終了します。この値を大きく設定しすぎると、真の最小値に到達する前に終了してしまう「早期収束」が発生する可能性があります。一方、小さく設定しすぎると、数値ノイズの影響しか残らない状況でも不要な計算を続けることになります。xtol(ステップ許容誤差)は、目的関数値ではなく探索空間における移動量に基づいて収束を判定します。現在の位置と前回の位置とのユークリッド距離を評価し、その移動量が十分に小さくなると最適化を終了します。これは、目的関数値がほとんど変化しない平坦な領域で、アルゴリズムが実質的な改善なく微小な移動を繰り返すことを防ぐために有効です。
- SetPGtolメソッドは、射影勾配の大きさに対する閾値を設定します。勾配がこの値を下回ると、アルゴリズムは収束したと判断します。pgtol(射影勾配許容誤差)は、数学的に最も厳密な収束判定条件です。制約条件を考慮した後の勾配の大きさを評価し、勾配がゼロであれば停留点、すなわち局所最小値の理論的条件を満たします。pgtolを非常に小さく設定すると、アルゴリズムは谷底に近い最小点をより高い精度で探索します。一方、大きめに設定すると、最適解近傍に到達した時点で処理を終了できます。
- SetAccuracyメソッドは、勾配計算時のステップ幅を決定するために使用する機械精度、または目的関数の推定誤差を設定します。
- SetEtaメソッドは、0~1の範囲の値を受け取り、ラインサーチにおいて目的関数の減少量とステップ効率のバランスを調整します。
最小化処理は、CTruncNewtonCGのMinimizeメソッドを呼び出して実行します。
int return_code = tnc_minim.Minimize(sf);このメソッドは、唯一の引数としてCFunctor派生オブジェクトを受け取ります。最適化が完了した後は、以下のCTruncNewtonCGメソッドを使用して結果を取得できます。
- Solutionメソッドは、得られた最小解に対応する変数を返す
- ObjectiveResultメソッドは、得られた最小解における目的関数の値を返す
- ObjectiveGradientメソッドは、得られた最小解における目的関数の勾配を返す
- NumFevalメソッドは、目的関数が評価された回数を返す
- NumItersメソッドは、最適化処理で実行された反復回数を返す
次のセクションでは、本実装を実際に適用し、その動作を検証します。
制約付き最適化テスト
最適化分野では、ローゼンブロック関数は最適化アルゴリズムの性能を評価するための代表的なベンチマーク関数として広く用いられています。この関数は、得られた最小解を含む細く平坦な谷を持つため、多くのアルゴリズムでは効率的な収束が難しく、最小化問題として高い難易度を持つことで知られています。同様の幾何学的特徴を持つ関数は、経済学や機械学習などの分野でも数多く見られます。
MetaTrader 5ストラテジテスターの数式評価機能を使用すると、この関数を可視化できます。この評価処理は、エキスパートアドバイザー(EA) RosenBrock.mq5で実装されています。
//+------------------------------------------------------------------+ //| RosenBrock.mq5 | //| Copyright 2025, MetaQuotes Ltd. | //| https://www.mql5.com | //+------------------------------------------------------------------+ #property copyright "Copyright 2025, MetaQuotes Ltd." #property link "https://www.mql5.com" #property version "1.00" //--- input parameters input double x=-1.2;// start=-2, step=0.01, stop=2 input double y=1.0;// start=-1, step=0.01, stop=3 //+------------------------------------------------------------------+ //| Tester function | //+------------------------------------------------------------------+ double OnTester() { //--- double ret= pow(1.-x,2.0)+100.*pow(y-pow(x,2.),2.0); //--- return(ret); } //+------------------------------------------------------------------+
EAをテスターで実行すると、[1,1]に得られた最小解があることがわかります。

三次元可視化により、最小値を隠している特徴的な「狭い谷」が強調されます。この性質のため、最適化アルゴリズムは振動を繰り返したり、収束が停滞したりすることがよくあります。

TNCミニマイザが真の大域最小値にどれほど近づけるかを確認するために、TestTNC.mq5スクリプトを使用します。
//+------------------------------------------------------------------+ //| TestTNC.mq5 | //| Copyright 2025, MetaQuotes Ltd. | //| https://www.mql5.com | //+------------------------------------------------------------------+ #property copyright "Copyright 2025, MetaQuotes Ltd." #property link "https://www.mql5.com" #property version "1.00" #property script_show_inputs #include<tnc/tnc.mqh> #include<np.mqh> //--- input double in_eta = -1.0; input double in_stepmx = 0.0; input double in_accuracy = 0.0; input double in_fmin = 0.0; input double in_ftol = -1.0; input double in_xtol = -1.0; input double in_pgtol = -1.0; input double in_rescale = -1.0; input int in_maxCGit = -1; input int in_maxnfeval = 100; //+------------------------------------------------------------------+ //| Script program start function | //+------------------------------------------------------------------+ void OnStart() { //--- vector init_params = {-1.2, 1.0};//vector::Zeros(2); CScalarFunc sf; matrix bd = {{-2, 2}, {-1, 3}}; sf.setBounds(bd); sf.setGradOption(GRAD_POINT_CALLABLE); if(!sf.initialize(init_params)) return; tnc::CTruncNewtonCG tnc_minim; tnc_minim.SetLoglevel(tnc::TNC_MSG_ALL); tnc_minim.SetMaxCGit(in_maxCGit); tnc_minim.SetMaxFunCalls(in_maxnfeval); tnc_minim.SetEta(in_eta); tnc_minim.SetStepMax(in_stepmx); tnc_minim.SetAccuracy(in_accuracy); tnc_minim.SetFmin(in_fmin); tnc_minim.SetFtol(in_ftol); tnc_minim.SetXtol(in_xtol); tnc_minim.SetPGtol(in_pgtol); tnc_minim.SetRescaleFactor(in_rescale); int return_code = tnc_minim.Minimize(sf); Print(" Optimization return code ", tnc::tnc_rc_string[return_code]); Print(" optimization solution ", tnc_minim.Solution()); } //+------------------------------------------------------------------+ class CScalarFunc:public CFunctor { public: CScalarFunc(void) { } ~CScalarFunc(void) { } virtual double orig_fun(vector& x) { return np::rosen(x); } virtual vector grad_fun(vector& x) { return np::rosen_gradient(x); } }; //+------------------------------------------------------------------+
スクリプトを実行すると、以下の出力が得られます。
QO 0 19:06:52.209 TestTNC (GBPUSD,D1) NIT NF F GTG GJ 0 19:06:52.209 TestTNC (GBPUSD,D1) 0 1 2.420000000000000E+01 5.42273600E+04 MQ 0 19:06:52.209 TestTNC (GBPUSD,D1) tnc: fscale = 0.00107357 DE 0 19:06:52.209 TestTNC (GBPUSD,D1) 1 3 4.567781791360708E+00 9.57591179E+02 RH 0 19:06:52.209 TestTNC (GBPUSD,D1) 2 5 4.127793942693432E+00 3.78238779E+00 LO 0 19:06:52.209 TestTNC (GBPUSD,D1) tnc: fscale = 0.128546 GK 0 19:06:52.209 TestTNC (GBPUSD,D1) 3 7 4.116988513105550E+00 1.81606916E+01 MO 0 19:06:52.209 TestTNC (GBPUSD,D1) 4 16 3.317053801700717E+00 2.94434047E+02 MS 0 19:06:52.209 TestTNC (GBPUSD,D1) 5 20 3.172067249923803E+00 5.34629683E+02 HD 0 19:06:52.209 TestTNC (GBPUSD,D1) 6 26 1.770939721007689E+00 6.64565256E+01 LI 0 19:06:52.209 TestTNC (GBPUSD,D1) 7 28 1.651272391609239E+00 4.99232721E+00 GR 0 19:06:52.210 TestTNC (GBPUSD,D1) 8 38 1.269069338899967E+00 4.60590202E+01 JF 0 19:06:52.210 TestTNC (GBPUSD,D1) 9 42 1.110240253054296E+00 7.40799184E+01 MJ 0 19:06:52.210 TestTNC (GBPUSD,D1) 10 44 7.993902318805092E-01 3.13860496E+00 LM 0 19:06:52.210 TestTNC (GBPUSD,D1) 11 48 6.601238685636209E-01 3.15085702E+01 EP 0 19:06:52.210 TestTNC (GBPUSD,D1) 12 53 3.391820074593478E-01 5.23994236E+01 DD 0 19:06:52.210 TestTNC (GBPUSD,D1) 13 55 2.644309665366499E-01 7.39869479E-01 NO 0 19:06:52.210 TestTNC (GBPUSD,D1) 14 61 1.675328075835125E-01 1.16689434E+01 ER 0 19:06:52.210 TestTNC (GBPUSD,D1) 15 65 1.224744553458436E-01 3.21610219E+01 LE 0 19:06:52.210 TestTNC (GBPUSD,D1) 16 67 6.925326451315299E-02 6.87199782E-01 JI 0 19:06:52.210 TestTNC (GBPUSD,D1) 17 73 2.930633981879897E-02 1.06326643E+01 GL 0 19:06:52.210 TestTNC (GBPUSD,D1) 18 77 1.341804114201537E-03 2.64590552E+00 DG 0 19:06:52.210 TestTNC (GBPUSD,D1) 19 79 2.447262820598820E-04 6.40221664E-04 DE 0 19:06:52.210 TestTNC (GBPUSD,D1) tnc: fscale = 9.88041 FR 0 19:06:52.210 TestTNC (GBPUSD,D1) 20 81 2.442825206410282E-04 2.70915012E-04 EE 0 19:06:52.210 TestTNC (GBPUSD,D1) 21 83 5.876599419779381E-06 1.17530118E-02 FI 0 19:06:52.210 TestTNC (GBPUSD,D1) 22 85 3.385708330669465E-11 1.17585414E-10 QS 0 19:06:52.210 TestTNC (GBPUSD,D1) tnc: fscale = 23054.9 NG 0 19:06:52.210 TestTNC (GBPUSD,D1) tnc: |fn-fn-1] = 7.61917e-14 -> convergence FN 0 19:06:52.210 TestTNC (GBPUSD,D1) 23 87 3.378089160219990E-11 3.50223471E-11 DN 0 19:06:52.210 TestTNC (GBPUSD,D1) tnc: Converged (|f_n-f_(n-1)| ~= 0) MK 0 19:06:52.210 TestTNC (GBPUSD,D1) Optimization return code 1 LS 0 19:06:52.210 TestTNC (GBPUSD,D1) optimization solution [0.9999941953936243,0.9999883612512912]
結果から、オプティマイザは大域最小値にかなり近い解を得ることができることが分かります。ただし、指定された許容誤差の影響により、最終的な値は必ずしも正確に1.0にはなりません。より実用的な利用例を示すために、ここではロジスティック回帰の実装において、L-BFGSミニマイザをTNCソルバーに置き換えます。
TNCに基づくロジスティック回帰の実装
これまでロジスティック回帰用のCLogitクラスでは、最適化エンジンとしてALGLIBのL-BFGS実装のみを利用していました。このセクションでは、ロジスティック回帰モデルを構築するための、堅牢な代替手段としてTNCオプティマイザを追加します。これは、既存のFitメソッドを変更し、使用するオプティマイザを選択できるようにすることで実現します。
//+------------------------------------------------------------------+ //| fit a model | //+------------------------------------------------------------------+ bool Clogit::Fit(matrix &predictors, vector &targets, ENUM_MINIM_METHOD minimizer=MINIM_LBFGS, tnc::ENUM_TNC_MESSAGE in_display = tnc::TNC_MSG_NONE, double in_eta = -1.0, double in_stepmx = 0.0, double in_accuracy = 0.0, double in_fmin = 0.0, double in_ftol = -1.0, double in_xtol = -1.0, double in_pgtol = -1.0, double in_rescale = -1.0, int in_maxCGit = -1, int in_maxnfeval = 100) { switch(minimizer) { case MINIM_LBFGS: return fit_lbfgs(predictors,targets); case MINIM_TNC: return fit_tnc(predictors,targets,in_display,in_eta,in_stepmx,in_accuracy,in_fmin,in_ftol,in_xtol,in_pgtol,in_rescale,in_maxCGit,in_maxnfeval); } return false; }
これを可能にするために、logistic.mqh内のコードを更新し、ENUM_MINIM_METHOD列挙型を定義します。この列挙型は、利用可能な最適化アルゴリズムをまとめて管理するものです。
//+------------------------------------------------------------------+ //|minimizer used | //+------------------------------------------------------------------+ enum ENUM_MINIM_METHOD { MINIM_LBFGS=0,//LBGS MINIM_TNC//TNC };
さらに、CFunctorクラスを継承したCLogitFunctorクラスがヘッダーファイルに追加されました。また、2つのprivateメンバー関数であるfit_L-BFGSとfit_tncが、それぞれ対応する最適化アルゴリズムの実行処理を担当します。
bool fit_lbfgs(matrix &predictors, vector &targets); bool fit_tnc(matrix &predictors, vector &targets, tnc::ENUM_TNC_MESSAGE in_display = tnc::TNC_MSG_NONE, double in_eta = -1.0, double in_stepmx = 0.0, double in_accuracy = 0.0, double in_fmin = 0.0, double in_ftol = -1.0, double in_xtol = -1.0, double in_pgtol = -1.0, double in_rescale = -1.0, int in_maxCGit = -1, int in_maxnfeval = 100);
スクリプトLogisticRegression.mq5では、古典的なIrisデータセットを用いてCLogitクラスの動作を検証します。このプログラムでは、データの一部を使用して2種類のオプティマイザでモデルを学習させ、その結果得られたモデルパラメータをターミナルの操作ログに出力します。
void OnStart() { //--- CHighQualityRandStateShell rngstate; CHighQualityRand::HQRndSeed(Random_Seed,Random_Seed+Random_Seed,rngstate.GetInnerObj()); //--- Print(iris_data); string lines[],cells[]; int ncells,nlines; //--- nlines = StringSplit(iris_data,StringGetCharacter("\n",0),lines); matrix data = matrix::Zeros(0,0); for(int i = 1; i<nlines-2; ++i) { ncells = StringSplit(lines[i],StringGetCharacter(",",0),cells); if(!data.Rows()) data.Resize(nlines-2,ncells-1); for(int k = 1; k<(ncells); ++k) data[i-1,k-1] = StringToDouble(cells[k]); } //--- COneHotEncoder enc; ulong colum[1] = {4}; if(!enc.fit(data,colum)) { Print(" failed to encode data "); return; } //--- data = enc.transform(data); //--- long rindices[],trainset[],testset[]; np::arange(rindices,int(data.Rows())); //--- np::shuffleArray(rindices,GetPointer(rngstate)); ArrayCopy(trainset,rindices,0,0,int(ceil(Tra_Test_Split*rindices.Size()))); ArraySort(trainset); //--- CSortedSet<long> test_set(rindices); //--- test_set.ExceptWith(trainset); //--- test_set.CopyTo(testset); //--- matrix testdata = np::selectMatrixRows(data,testset); matrix test_predictors = np::sliceMatrixCols(testdata,0,4); vector test_targets = testdata.Col(4); matrix traindata = np::selectMatrixRows(data,trainset); matrix tra_preditors = np::sliceMatrixCols(traindata,0,4); vector tra_targets = traindata.Col(4); //--- logistic::Clogit logit; if(!logit.Fit(tra_preditors,tra_targets)) { Print(" failed to fit data with lbfgs"); return; } //--- Print(" LBFGS logit results "); Print(" coefs ", logit.Get_Coefs()); Print(" bias ", logit.Get_Bias()); //--- if(!logit.Fit(tra_preditors,tra_targets,logistic::MINIM_TNC,0,_eta_,_stepmx_,_accuracy_,_fmin_,_ftol_,_xtol_,_pgtol_,_rescale_,_maxCGit_,_maxnfeval_)) { Print(" failed to fit data with tnc "); return; } //--- Print(" TNC logit results "); Print(" coefs ", logit.Get_Coefs()); Print(" bias ", logit.Get_Bias()); }
結果から、両方のオプティマイザがほぼ同じパラメータに収束することが確認できます。これは、TNCがロジスティック回帰における代替のオプティマイザとして利用可能であることを示しています。
CG 0 19:11:15.136 LogisticRegression (GBPUSD,D1) index,sl,sw,pl,pw,target GS 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 0,5.1,3.5,1.4,0.2,0 HK 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 1,4.9,3.0,1.4,0.2,0 LR 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 2,4.7,3.2,1.3,0.2,0 CM 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 3,4.6,3.1,1.5,0.2,0 ED 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 4,5.0,3.6,1.4,0.2,0 DL 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 5,5.4,3.9,1.7,0.4,0 OG 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 6,4.6,3.4,1.4,0.3,0 IN 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 7,5.0,3.4,1.5,0.2,0 HI 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 8,4.4,2.9,1.4,0.2,0 QQ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 9,4.9,3.1,1.5,0.1,0 NK 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 10,5.4,3.7,1.5,0.2,0 HP 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 11,4.8,3.4,1.6,0.2,0 PI 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 12,4.8,3.0,1.4,0.1,0 CF 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 13,4.3,3.0,1.1,0.1,0 OO 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 14,5.8,4.0,1.2,0.2,0 HD 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 15,5.7,4.4,1.5,0.4,0 JL 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 16,5.4,3.9,1.3,0.4,0 PE 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 17,5.1,3.5,1.4,0.3,0 OR 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 18,5.7,3.8,1.7,0.3,0 NK 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 19,5.1,3.8,1.5,0.3,0 HP 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 20,5.4,3.4,1.7,0.2,0 GI 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 21,5.1,3.7,1.5,0.4,0 HF 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 22,4.6,3.6,1.0,0.2,0 HO 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 23,5.1,3.3,1.7,0.5,0 MG 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 24,4.8,3.4,1.9,0.2,0 NL 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 25,5.0,3.0,1.6,0.2,0 OE 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 26,5.0,3.4,1.6,0.4,0 LR 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 27,5.2,3.5,1.5,0.2,0 KK 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 28,5.2,3.4,1.4,0.2,0 NP 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 29,4.7,3.2,1.6,0.2,0 LI 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 30,4.8,3.1,1.6,0.2,0 JF 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 31,5.4,3.4,1.5,0.4,0 NN 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 32,5.2,4.1,1.5,0.1,0 QG 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 33,5.5,4.2,1.4,0.2,0 JL 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 34,4.9,3.1,1.5,0.2,0 KE 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 35,5.0,3.2,1.2,0.2,0 IR 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 36,5.5,3.5,1.3,0.2,0 RK 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 37,4.9,3.6,1.4,0.1,0 LP 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 38,4.4,3.0,1.3,0.2,0 OI 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 39,5.1,3.4,1.5,0.2,0 NQ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 40,5.0,3.5,1.3,0.3,0 NN 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 41,4.5,2.3,1.3,0.3,0 QG 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 42,4.4,3.2,1.3,0.2,0 QL 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 43,5.0,3.5,1.6,0.6,0 IE 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 44,5.1,3.8,1.9,0.4,0 RR 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 45,4.8,3.0,1.4,0.3,0 LK 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 46,5.1,3.8,1.6,0.2,0 GP 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 47,4.6,3.2,1.4,0.2,0 HH 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 48,5.3,3.7,1.5,0.2,0 CQ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 49,5.0,3.3,1.4,0.2,0 PN 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 50,7.0,3.2,4.7,1.4,1 KG 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 51,6.4,3.2,4.5,1.5,1 JL 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 52,6.9,3.1,4.9,1.5,1 JE 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 53,5.5,2.3,4.0,1.3,1 QR 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 54,6.5,2.8,4.6,1.5,1 LK 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 55,5.7,2.8,4.5,1.3,1 QS 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 56,6.3,3.3,4.7,1.6,1 IH 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 57,4.9,2.4,3.3,1.0,1 QQ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 58,6.6,2.9,4.6,1.3,1 NN 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 59,5.2,2.7,3.9,1.4,1 EG 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 60,5.0,2.0,3.5,1.0,1 OL 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 61,5.9,3.0,4.2,1.5,1 JE 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 62,6.0,2.2,4.0,1.0,1 PR 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 63,6.1,2.9,4.7,1.4,1 DJ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 64,5.6,2.9,3.6,1.3,1 NS 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 65,6.7,3.1,4.4,1.4,1 PH 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 66,5.6,3.0,4.5,1.5,1 JQ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 67,5.8,2.7,4.1,1.0,1 NN 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 68,6.2,2.2,4.5,1.5,1 JG 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 69,5.6,2.5,3.9,1.1,1 RL 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 70,5.9,3.2,4.8,1.8,1 RE 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 71,6.1,2.8,4.0,1.3,1 EM 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 72,6.3,2.5,4.9,1.5,1 JJ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 73,6.1,2.8,4.7,1.2,1 JS 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 74,6.4,2.9,4.3,1.3,1 QH 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 75,6.6,3.0,4.4,1.4,1 KQ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 76,6.8,2.8,4.8,1.4,1 JN 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 77,6.7,3.0,5.0,1.7,1 JG 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 78,6.0,2.9,4.5,1.5,1 PL 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 79,5.7,2.6,3.5,1.0,1 LD 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 80,5.5,2.4,3.8,1.1,1 GM 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 81,5.5,2.4,3.7,1.0,1 FJ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 82,5.8,2.7,3.9,1.2,1 DS 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 83,6.0,2.7,5.1,1.6,1 RH 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 84,5.4,3.0,4.5,1.5,1 OQ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 85,6.0,3.4,4.5,1.6,1 MN 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 86,6.7,3.1,4.7,1.5,1 NG 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 87,6.3,2.3,4.4,1.3,1 NO 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 88,5.6,3.0,4.1,1.3,1 QD 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 89,5.5,2.5,4.0,1.3,1 IM 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 90,5.5,2.6,4.4,1.2,1 DJ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 91,6.1,3.0,4.6,1.4,1 DS 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 92,5.8,2.6,4.0,1.2,1 LH 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 93,5.0,2.3,3.3,1.0,1 RQ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 94,5.6,2.7,4.2,1.3,1 MN 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 95,5.7,3.0,4.2,1.2,1 QF 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 96,5.7,2.9,4.2,1.3,1 MO 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 97,6.2,2.9,4.3,1.3,1 LD 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 98,5.1,2.5,3.0,1.1,1 HM 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 99,5.7,2.8,4.1,1.3,1 EI 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 100,6.3,3.3,6.0,2.5,2 FP 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 101,5.8,2.7,5.1,1.9,2 QK 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 102,7.1,3.0,5.9,2.1,2 DR 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 103,6.3,2.9,5.6,1.8,2 RM 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 104,6.5,3.0,5.8,2.2,2 ED 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 105,7.6,3.0,6.6,2.1,2 PO 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 106,4.9,2.5,4.5,1.7,2 EF 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 107,7.3,2.9,6.3,1.8,2 OQ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 108,6.7,2.5,5.8,1.8,2 RH 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 109,7.2,3.6,6.1,2.5,2 FS 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 110,6.5,3.2,5.1,2.0,2 RJ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 111,6.4,2.7,5.3,1.9,2 NE 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 112,6.8,3.0,5.5,2.1,2 IL 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 113,5.7,2.5,5.0,2.0,2 IG 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 114,5.8,2.8,5.1,2.4,2 IN 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 115,6.4,3.2,5.3,2.3,2 KI 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 116,6.5,3.0,5.5,1.8,2 MP 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 117,7.7,3.8,6.7,2.2,2 HK 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 118,7.7,2.6,6.9,2.3,2 LR 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 119,6.0,2.2,5.0,1.5,2 LM 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 120,6.9,3.2,5.7,2.3,2 FD 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 121,5.6,2.8,4.9,2.0,2 RO 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 122,7.7,2.8,6.7,2.0,2 FF 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 123,6.3,2.7,4.9,1.8,2 CQ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 124,6.7,3.3,5.7,2.1,2 MH 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 125,7.2,3.2,6.0,1.8,2 LS 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 126,6.2,2.8,4.8,1.8,2 FJ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 127,6.1,3.0,4.9,1.8,2 GE 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 128,6.4,2.8,5.6,2.1,2 DL 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 129,7.2,3.0,5.8,1.6,2 FG 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 130,7.4,2.8,6.1,1.9,2 DN 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 131,7.9,3.8,6.4,2.0,2 OI 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 132,6.4,2.8,5.6,2.2,2 LP 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 133,6.3,2.8,5.1,1.5,2 GK 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 134,6.1,2.6,5.6,1.4,2 HR 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 135,7.7,3.0,6.1,2.3,2 GM 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 136,6.3,3.4,5.6,2.4,2 ND 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 137,6.4,3.1,5.5,1.8,2 PO 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 138,6.0,3.0,4.8,1.8,2 FF 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 139,6.9,3.1,5.4,2.1,2 GQ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 140,6.7,3.1,5.6,2.4,2 NH 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 141,6.9,3.1,5.1,2.3,2 GS 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 142,5.8,2.7,5.1,1.9,2 FJ 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 143,6.8,3.2,5.9,2.3,2 EE 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 144,6.7,3.3,5.7,2.5,2 FL 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 145,6.7,3.0,5.2,2.3,2 NG 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 146,6.3,2.5,5.0,1.9,2 CN 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 147,6.5,3.0,5.2,2.0,2 LI 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 148,6.2,3.4,5.4,2.3,2 LP 0 19:11:15.136 LogisticRegression (GBPUSD,D1) 149,5.9,3.0,5.1,1.8,2 LI 0 19:11:15.136 LogisticRegression (GBPUSD,D1) EN 0 19:11:15.138 LogisticRegression (GBPUSD,D1) LBFGS logit results QF 0 19:11:15.138 LogisticRegression (GBPUSD,D1) coefs [[-0.3544914154418004,0.7640895553785291,-2.029702076841132,-0.8297514969607693]] EK 0 19:11:15.138 LogisticRegression (GBPUSD,D1) bias [5.587065638186691] RJ 0 19:11:15.141 LogisticRegression (GBPUSD,D1) TNC logit results FR 0 19:11:15.141 LogisticRegression (GBPUSD,D1) coefs [[-0.3543682771949164,0.7643878803868314,-2.029620459724486,-0.8298337657974659]] IN 0 19:11:15.141 LogisticRegression (GBPUSD,D1) bias [5.585350403781237]
この記事の以降のセクションでは、TNCオプティマイザの実装について、より詳細に説明します。具体的には、コードがどのように動作するのか、また特定の変数がどのように定義されているのかについて解説します。
TNCのコア実装
TNC最適化の実装は、tnc.mqhに含まれるコアロジック(主要なアルゴリズム処理)と目的関数、その導関数、およびヘッセ行列のラッパーを定義する汎用的な補助コンポーネントの2つの部分から構成されています。まず、tnc.mqhに焦点を当てます。これはTNCアルゴリズムの中核となる処理を担当しているためです。tnc.mqhのコードは、内部的な情報交換をおこなうための構造的な基盤となる、いくつかの列挙型を定義することから始まります。
enum ENUM_TNC_MESSAGE { TNC_MSG_NONE = 0, /* No messages */ TNC_MSG_ITER = 1, /* One line per iteration */ TNC_MSG_INFO = 2, /* Informational messages */ TNC_MSG_EXIT = 8, /* Exit reasons */ TNC_MSG_ALL = TNC_MSG_ITER | TNC_MSG_INFO | TNC_MSG_EXIT /* All messages */ }; enum ENUM_TNC_RC { TNC_MINRC = -3, /* Constant to add to get the rc_string */ TNC_ENOMEM = -3, /* Memory allocation failed */ TNC_EINVAL = -2, /* Invalid parameters (n<0) */ TNC_INFEASIBLE = -1, /* Infeasible (low bound > up bound) */ TNC_LOCALMINIMUM = 0, /* Local minima reach (|pg| ~= 0) */ TNC_FCONVERGED = 1, /* Converged (|f_n-f_(n-1)| ~= 0) */ TNC_XCONVERGED = 2, /* Converged (|x_n-x_(n-1)| ~= 0) */ TNC_MAXFUN = 3, /* Max. number of function evaluations reach */ TNC_LSFAIL = 4, /* Linear search failed */ TNC_CONSTANT = 5, /* All lower bounds are equal to the upper bounds */ TNC_NOPROGRESS = 6, /* Unable to progress */ TNC_USERABORT = 7 /* User requested end of minization */ };
- ENUM_TNC_MESSAGE列挙型は、出力の詳細度を制御するためのビットマスクフラグを定義する
- ENUM_TNC_RC列挙型は、終了理由を示すリターンコードのセットを提供する
- tnc_rc_string変数は、これらのリターンコードを人が理解しやすいメッセージに対応付ける文字列配列
const string tnc_rc_string[11] = { "Memory allocation failed", "Invalid parameters (n<0)", "Infeasible (low bound > up bound)", "Local minima reach (|pg| ~= 0)", "Converged (|f_n-f_(n-1)| ~= 0)", "Converged (|x_n-x_(n-1)| ~= 0)", "Maximum number of function evaluations reached", "Linear search failed", "All lower bounds are equal to the upper bounds", "Unable to progress", "User requested end of minimization" };
さらに、線形探索サブルーチンおよび点探索サブルーチンの戻り状態を定義する専用の列挙型も用意されています。これにより、アルゴリズムは各入れ子ループ内での処理の進行状況を正確に追跡できます。
/* * getptc return codes */ enum ENUM_GETPTC_RC { GETPTC_OK = 0, /* Suitable point found */ GETPTC_EVAL = 1, /* Function evaluation required */ GETPTC_EINVAL = 2, /* Bad input values */ GETPTC_FAIL = 3 /* No suitable point found */ }; /* * linearSearch return codes */ enum ENUM_LS_RC { LS_OK = 0, /* Suitable point found */ LS_MAXFUN = 1, /* Max. number of function evaluations reach */ LS_FAIL = 2, /* No suitable point found */ LS_USERABORT = 3, /* User requested end of minimization */ LS_ENOMEM = 4 /* Memory allocation failed */ };
CObjectiveクラスは、最適化対象となる関数の状態を管理する役割を担います。このクラスは、問題の次元数、失敗の発生有無、実際の目的関数の処理を指すポインタなどのメタデータを保持します。
//+------------------------------------------------------------------+ //|objective function state | //+------------------------------------------------------------------+ class CObjective { protected: bool m_failed,m_fupdated,m_gupdated,m_hupdated; double m_lowestx,m_lowest_f; ulong m_size; IObjective *m_objective; public: CObjective(void):m_failed(false), m_size(0) { } ~CObjective(void) { } void set_objective(IObjective* &fun_obj) { m_objective = fun_obj; } void has_failed(bool yes) { m_failed = yes; } void set_dim(ulong size) { m_size = size; } bool failed(void) { return m_failed; } bool abort(void) { return IsStopped(); } ulong size(void) { return m_size; } ObjReturn objective(vector& x) { return m_objective.fun_and_grad(x); } virtual int callback(vector& x) { return 0; } }
このクラスはラッパーとして機能することで、TNCアルゴリズムが関数評価および勾配計算を要求できるようにすると同時に、最適化の進行状況に関する内部フラグを維持します。この抽象化レイヤーにより、コアソルバーは、最小化対象となる関数固有の数学的詳細から切り離された状態を保つことができます。クラス定義に続いて、コードでは目的関数およびコールバック関数のシグネチャを標準化するために、typedefによる関数ポインタが定義されています。
//+------------------------------------------------------------------+ //| function pointers | //+------------------------------------------------------------------+ typedef int(*tnc_function)(vector&, double&, vector&, CObjective&); typedef void(*tnc_callback)(vector&, CObjective&); //+------------------------------------------------------------------+ //| fpointer | //+------------------------------------------------------------------+ int func(vector& x, double& f, vector& g, CObjective& state) { ulong n = state.size(); vector x_data,g_data; if(state.abort()) return 2; if(state.failed()) return 1; x_data = x; ObjReturn fg = state.objective(x_data); f = fg.f; if(fg.g.Size()!=n) { printf(" tnc: gradient must have length ", n); return 2; } g = fg.g; return 0; } //+------------------------------------------------------------------+ //| callback | //+------------------------------------------------------------------+ void call_back(vector& x, CObjective& state) { ulong n = state.size(); vector x_data = x; if(state.callback(x_data)) state.has_failed(true); else state.has_failed(false); return;
func関数とcall_back関数およびcall_back関数は、低レベルソルバーとCObjectiveの状態管理機構との間を仲介するブリッジロジックとして機能します。これらの関数は、ベクトルからのデータ抽出、勾配の次元数の検証、ユーザーによって開始された中断処理の確認などをおこなった後、その結果をソルバーへ返します。これにより、コアアルゴリズムは呼び出し側で使用される特定のデータ構造から分離された状態を維持できます。
//+------------------------------------------------------------------+ //|Optimization results | //+------------------------------------------------------------------+ struct OptimizeResult { int return_code; int nfeval; int niter; vector solution; vector objective_result; vector objective_gradient; OptimizeResult(void) { return_code = WRONG_VALUE; nfeval = niter = 0; solution = objective_result = objective_gradient = vector::Zeros(0); } OptimizeResult(int rc,int feval,int iter,vector &x, vector& f, vector& g) { return_code = rc; nfeval = feval; niter = iter; solution = x; objective_result = f; objective_gradient = g; } OptimizeResult(OptimizeResult& other) { return_code = other.return_code; nfeval = other.nfeval; niter = other.niter; solution = other.solution; objective_result = other.objective_result; objective_gradient = other.objective_gradient; } void operator=(OptimizeResult& other) { return_code = other.return_code; nfeval = other.nfeval; niter = other.niter; solution = other.solution; objective_result = other.objective_result; objective_gradient = other.objective_gradient; } }
OptimizeResult構造体は、最適化処理の最終結果を格納するために設計された軽量なデータコンテナです。この構造体には、リターンコード、関数評価回数、実行された総反復回数など、重要なメタデータが保持されます。さらに、解を表す最終ベクトル、目的関数の評価結果、および勾配ベクトルも格納します。ソルバーの処理完了後に、アプリケーション内の異なる部分間でデータを円滑に受け渡せるように、この構造体には複数のコンストラクタと代入演算子が実装されています。
//+------------------------------------------------------------------+ //|class encapsulating TNC minimizer | //+------------------------------------------------------------------+ class CTruncNewtonCG: public CObject { private: int tnc(int n, vector& x, double &f, vector& g, tnc_function& function, CObjective &state, vector& low, vector& up, vector& scale, vector& offset, int messages, int maxCGit, int maxnfeval, double eta, double stepmx, double accuracy, double _fmin, double ftol, double xtol, double pgtol, double rescale, int &nfeval, int &niter, tnc_callback &callback); void coercex(int n, vector& x, const vector& low, const vector& up); void unscalex(int n, vector& x, const vector& xscale, const vector& xoffset); void scalex(int n, vector& x, const vector& xscale, const vector& xoffset); void scaleg(int n, vector& g, const vector& xscale, double& fscale); void setConstraints(int n, vector& x, int& pivot[], vector& xscale, vector& xoffset, vector& low, vector& up); ENUM_TNC_RC minize_tnc(int n, vector& x, double &f, vector& gfull, tnc_function function, CObjective &state, vector& xscale, vector& xoffset, double &fscale, vector& low, vector& up, ENUM_TNC_MESSAGE messages, int maxCGit, int maxnfeval, int &nfeval, int &niter, double eta, double stepmx, double accuracy, double _fmin, double ftol, double xtol, double pgtol, double rescale, tnc_callback& callback); void printCurrentIteration(int n, double f, vector& g, int niter, int nfeval, int &pivot[]); void project(int n, vector& x, const int &pivot[]); void projectConstants(int n, vector& x, const vector& xscale); double stepMax(double step, int n, vector& x, vector& dir, int &pivot[], vector& low, vector& up, vector& xscale, vector& xoffset); bool addConstraint(int n, vector& x, vector& p, int &pivot[], vector& low, vector& up, vector& xscale, vector& xoffset); bool removeConstraint(double gtpnew, double gnorm, double pgtolfs, double f, double fLastConstraint, vector& g, int &pivot[], int n); */ int tnc_direction(vector &zsol, vector &diagb, vector &x, vector& g, int n, int maxCGit, int maxnfeval, int &nfeval, bool upd1, double yksk, double yrsr, vector &sk, vector &yk, vector &sr, vector &yr, bool lreset, tnc_function & function, CObjective &state, vector& xscale, vector& xoffset, double fscale, int &pivot[], double accuracy, double gnorm, double xnorm, vector& low, vector& up); void diagonalScaling(int n, vector& e, vector& v, vector& gv, vector& r); double initialStep(double fnew, double _fmin, double gtp, double smax); int hessianTimesVector(vector& v, vector& gv, int n, vector& x, vector& g, tnc_function & function, CObjective &state, vector& xscale, vector& xoffset, double fscale, double accuracy, double xnorm, vector& low, vector& up); int msolve(vector& g, vector& y, int n, vector& sk, vector& yk, vector& diagb, vector& sr, vector& yr, bool upd1, double yksk, double yrsr, bool lreset); void ssbfgs(int n, double gamma, vector& sj, vector& hjv, vector& hjyj, double yjsj, double yjhyj, double vsj, double vhyj, vector& hjp1v); int initPreconditioner(vector& diagb, vector& emat, int n, bool lreset, double yksk, double yrsr, vector& sk, vector& yk, vector& sr, vector& yr, bool upd1); ENUM_LS_RC linearSearch(int n, tnc_function & function, CObjective &state, vector& low, vector& up, vector& xscale, vector& xoffset, double fscale, int &pivot[], double eta, double ftol, double xbnd, vector& p, vector& x, double &f, double &alpha, vector& gfull, int maxnfeval, int &nfeval); ENUM_GETPTC_RC getptcInit(double &reltol, double &abstol, double tnytol, double eta, double rmu, double xbnd, double &u, double &fu, double &gu, double &xmin, double &_fmin, double &gmin, double &xw, double &fw, double &gw, double &a, double &b, double &oldf, double &b1, double &scxbnd, double &e, double &step, double &factor, bool & braktd, double >est1, double >est2, double &tol); ENUM_GETPTC_RC getptcIter(double big, double rtsmll, double &reltol, double &abstol, double tnytol, double fpresn, double xbnd, double &u, double &fu, double &gu, double &xmin, double &_fmin, double &gmin, double &xw, double &fw, double &gw, double &a, double &b, double &oldf, double &b1, double &scxbnd, double &e, double &step, double &factor, bool & braktd, double >est1, double >est2, double &tol); void dxpy1(int n, const vector& dx, vector& dy); void daxpy1(int n, double da, const vector& dx, vector& dy); void dcopy1(int n, const vector& dx, vector& dy); void dneg1(int n, vector& v); double ddot1(int n, const vector& dx, const vector& dy); double dnrm21(int n, const vector& dx); OptimizeResult tnc_minimize(CFunctor &fungrad,vector& scale, vector& offset, int messages, int maxCGit, int maxfun, double eta, double stepmax, double accuracy, double fmin_, double ftol, double xtol,double pgtol, double rescale); vector m_scale,m_offset; int m_messages,m_maxCGit, m_maxfun; double m_eta, m_stepmax, m_accuracy, m_fmin,m_ftol,m_xtol,m_pgtol, m_rescale; OptimizeResult m_result; public: CTruncNewtonCG(void); ~CTruncNewtonCG(void); void SetScale(vector& scale); void SetOffset(vector& offset); void SetLoglevel(ENUM_TNC_MESSAGE messages); void SetMaxCGit(int maxCGit); void SetMaxFunCalls(int maxfun); void SetEta(double eta); void SetStepMax(double stepmax); void SetAccuracy(double accuracy); void SetFmin(double f_min); void SetFtol(double ftol); void SetXtol(double xtol); void SetPGtol(double pgtol); void SetRescaleFactor(double rescale); int Minimize(CFunctor &fungrad); vector Solution(void); double ObjectiveResult(void); vector ObjectiveGradient(void); int NumFevals(void); int NumIters(void); };
CTruncNewtonCGクラスには、最適化アルゴリズムの中核となるロジックが含まれています。このクラスは、初期パラメータの検証、変数スケーリングの処理、許容誤差およびステップサイズのデフォルト値を設定するtncエントリメソッドを備えています。
int tnc(int n, vector& x, double &f, vector& g, tnc_function& function, CObjective &state, vector& low, vector& up, vector& scale, vector& offset, int messages, int maxCGit, int maxnfeval, double eta, double stepmx, double accuracy, double _fmin, double ftol, double xtol, double pgtol, double rescale, int &nfeval, int &niter, tnc_callback &callback) { int rc, frc, i, nc, nfeval_local, free_low = TNC_FALSE,free_up = TNC_FALSE, free_g = TNC_FALSE; double fscale, rteps; vector xscale,xoffset; nfeval = nfeval_local = 0; /* Check for errors in the input parameters */ if(n == 0) { rc = TNC_CONSTANT; nfeval = (nfeval==0)?nfeval_local:nfeval; if(bool(messages & TNC_MSG_EXIT)) printf(" tnc: %s", tnc_rc_string[rc - TNC_MINRC]); return rc; } if(n < 0) { rc = TNC_EINVAL; nfeval = (nfeval==0)?nfeval_local:nfeval; if(bool(messages & TNC_MSG_EXIT)) printf(" tnc: %s", tnc_rc_string[rc - TNC_MINRC]); return rc; } /* Check bounds arrays */ if(!low.Size()) { low = vector::Zeros(n); free_low = TNC_TRUE; for(i = 0; i < n; i++) { low[i] = -HUGE_VAL; } } if(!up.Size()) { up = vector::Zeros(n); if(up.Size()==0) { rc = TNC_ENOMEM; nfeval = (nfeval==0)?nfeval_local:nfeval; if(bool(messages & TNC_MSG_EXIT)) printf(" tnc: %s", tnc_rc_string[rc - TNC_MINRC]); return rc; } free_up = TNC_TRUE; for(i = 0; i < n; i++) { up[i] = HUGE_VAL; } } /* Coherency check */ for(i = 0; i < n; i++) { if(low[i] > up[i]) { rc = TNC_INFEASIBLE; nfeval = (nfeval==0)?nfeval_local:nfeval; if(bool(messages & TNC_MSG_EXIT)) printf(" tnc: %s", tnc_rc_string[rc - TNC_MINRC]); return rc; } } /* Coerce x into bounds */ coercex(n, x, low, up); if(maxnfeval < 1) { rc = TNC_MAXFUN; nfeval = (nfeval==0)?nfeval_local:nfeval; if(bool(messages & TNC_MSG_EXIT)) printf(" tnc: %s", tnc_rc_string[rc - TNC_MINRC]); return rc; } /* Allocate g if necessary */ if(g.Size()==0) { g = vector::Zeros(n); if(g.Size()==0) { rc = TNC_ENOMEM; nfeval = (nfeval==0)?nfeval_local:nfeval; if(bool(messages & TNC_MSG_EXIT)) printf(" tnc: %s", tnc_rc_string[rc - TNC_MINRC]); return rc; } free_g = TNC_TRUE; } /* Initial function evaluation */ //Print(__FUNCTION__," - ", x); frc = function(x, f, g, state); (nfeval)++; if(frc) { rc = TNC_USERABORT; nfeval = (nfeval==0)?nfeval_local:nfeval; if(bool(messages & TNC_MSG_EXIT)) printf(" tnc: %s", tnc_rc_string[rc - TNC_MINRC]); return rc; } /* Constant problem ? */ for(nc = 0, i = 0; i < n; i++) { if((low[i] == up[i]) || (scale.Size() != 0 && scale[i] == 0.0)) { nc++; } } if(nc == n) { rc = TNC_CONSTANT; nfeval = (nfeval==0)?nfeval_local:nfeval; if(bool(messages & TNC_MSG_EXIT)) printf(" tnc: %s", tnc_rc_string[rc - TNC_MINRC]); return rc; } /* Scaling parameters */ xscale = vector::Zeros(n); if(xscale.Size()==0) { rc = TNC_ENOMEM; nfeval = (nfeval==0)?nfeval_local:nfeval; if(bool(messages & TNC_MSG_EXIT)) printf(" tnc: %s", tnc_rc_string[rc - TNC_MINRC]); return rc;; } xoffset = vector::Zeros(n); if(xoffset.Size()==0) { rc = TNC_ENOMEM; nfeval = (nfeval==0)?nfeval_local:nfeval; if(bool(messages & TNC_MSG_EXIT)) printf(" tnc: %s", tnc_rc_string[rc - TNC_MINRC]); return rc; } fscale = 1.0; for(i = 0; i < n; i++) { if(scale.Size() != NULL) { xscale[i] = fabs(scale[i]); if(xscale[i] == 0.0) { xoffset[i] = low[i] = up[i] = x[i]; } } else if(low[i] != -HUGE_VAL && up[i] != HUGE_VAL) { xscale[i] = up[i] - low[i]; xoffset[i] = (up[i] + low[i]) * 0.5; } else { xscale[i] = 1.0 + fabs(x[i]); xoffset[i] = x[i]; } if(offset.Size() != NULL) { xoffset[i] = offset[i]; } } /* Default values for parameters */ rteps = sqrt(DBL_EPSILON); if(stepmx < rteps * 10.0) { stepmx = 1.0e1; } if(eta < 0.0 || eta >= 1.0) { eta = 0.25; } if(rescale < 0) { rescale = 1.3; } if(maxCGit < 0) /* maxCGit == 0 is valid */ { maxCGit = n / 2; if(maxCGit < 1) { maxCGit = 1; } else if(maxCGit > 50) { maxCGit = 50; } } if(maxCGit > n) { maxCGit = n; } if(accuracy <= DBL_EPSILON) { accuracy = rteps; } if(ftol < 0.0) { ftol = accuracy; } if(pgtol < 0.0) { pgtol = 1e-2 * sqrt(accuracy); } if(xtol < 0.0) { xtol = rteps; } /* Optimisation */ rc = minize_tnc(n, x, f, g, function, state, xscale, xoffset, fscale, low, up,(ENUM_TNC_MESSAGE)messages, maxCGit, maxnfeval, nfeval, niter, eta, stepmx, accuracy, _fmin, ftol, xtol, pgtol, rescale, callback); if(bool(messages & TNC_MSG_EXIT)) printf(" tnc: %s", tnc_rc_string[rc - TNC_MINRC]); return rc; }
このセクションでは、変数を定義された範囲内に収める処理と、スケーリングパラメータの初期化を担当します。これらのパラメータにより、アルゴリズムは、桁違いに大きさの異なる変数を扱う場合でも、数値的な安定性を維持できます。
ENUM_TNC_RC minize_tnc(int n, vector& x, double &f, vector& gfull, tnc_function function, CObjective &state, vector& xscale, vector& xoffset, double &fscale, vector& low, vector& up, ENUM_TNC_MESSAGE messages, int maxCGit, int maxnfeval, int &nfeval, int &niter, double eta, double stepmx, double accuracy, double _fmin, double ftol, double xtol, double pgtol, double rescale, tnc_callback& callback) { double fLastReset, difnew, epsred, oldgtp, difold, oldf, xnorm, newscale, gnorm, ustpmax, fLastConstraint, spe, yrsr, yksk; vector temp,sk,yk,diagb,sr,yr,oldg,pk,g; double alpha = 0.0; /* Default unused value */ int i, icycle, oldnfeval, frc; bool lreset, newcon, upd1, remcon; ENUM_TNC_RC rc = TNC_ENOMEM; /* Default error */ niter = 0; int pivot[]; /* Allocate temporary vectors */ oldg = vector::Zeros(n); if(oldg.Size()==0) { return rc; } g = vector::Zeros(n); if(g.Size()==0) { return rc; } temp = vector::Zeros(n); if(temp.Size()==0) { return rc; } diagb = vector::Zeros(n); if(diagb.Size()==0) { return rc; } pk = vector::Zeros(n); if(pk.Size()==0) { return rc; } sk = vector::Zeros(n); if(sk.Size()==0) { return rc; } yk = vector::Zeros(n); if(yk.Size()==0) { return rc; } sr = vector::Zeros(n); if(sr.Size()==0) { return rc; } yr = vector::Zeros(n); if(yr.Size()==0) { return rc; } ArrayResize(pivot,n); if(pivot.Size()==0) { return rc; } /* Initialize variables */ difnew = 0.0; epsred = 0.05; upd1 = TNC_TRUE; icycle = n - 1; newcon = TNC_TRUE; /* Unneeded initialisations */ lreset = TNC_FALSE; yrsr = 0.0; yksk = 0.0; /* Initial scaling */ scalex(n, x, xscale, xoffset); f *=fscale; /* initial pivot calculation */ setConstraints(n, x, pivot, xscale, xoffset, low, up); dcopy1(n, gfull, g); scaleg(n, g, xscale, fscale); /* Test the lagrange multipliers to see if they are non-negative. */ for(i = 0; i < n; i++) { if(-pivot[i] * g[i] < 0.0) { pivot[i] = 0; } } project(n, g, pivot); /* Set initial values to other parameters */ gnorm = dnrm21(n, g); fLastConstraint = f; /* Value at last constraint */ fLastReset = f; /* Value at last reset */ if(bool(messages & TNC_MSG_ITER)) { printf(" NIT NF F GTG"); } if(bool(messages & TNC_MSG_ITER)) { printCurrentIteration(n, f / fscale, gfull, niter, nfeval, pivot); } /* Set the diagonal of the approximate hessian to unity. */ diagb = vector::Ones(diagb.Size()); /*for(i = 0; i < n; i++) { diagb[i] = 1.0; } */ /* Start of main iterative loop */ while(TNC_TRUE && !IsStopped()) { /* Local minimum test */ if(dnrm21(n, g) <= pgtol * (fscale)) { /* |PG| == 0.0 => local minimum */ dcopy1(n, gfull, g); project(n, g, pivot); if(bool(messages & TNC_MSG_INFO)) { printf("tnc: |pg| = %g -> local minimum", dnrm21(n,g)/(fscale)); } rc = TNC_LOCALMINIMUM; break; } /* Terminate if more than maxnfeval evaluations have been made */ if(nfeval >= maxnfeval) { rc = TNC_MAXFUN; break; } /* Rescale function if necessary */ newscale = dnrm21(n, g); if((newscale > DBL_EPSILON) && (fabs(log10(newscale)) > rescale)) { newscale = 1.0 / newscale; f *= newscale; fscale *= newscale; gnorm *= newscale; fLastConstraint *= newscale; fLastReset *= newscale; difnew *= newscale; g*=newscale; diagb = vector::Ones(diagb.Size()); upd1 = TNC_TRUE; icycle = n - 1; newcon = TNC_TRUE; if(bool(messages & TNC_MSG_INFO)) { printf("tnc: fscale = %g", fscale); } } dcopy1(n, x, temp); project(n, temp, pivot); xnorm = dnrm21(n, temp); oldnfeval = nfeval; /* Compute the new search direction */ frc = tnc_direction(pk, diagb, x, g, n, maxCGit, maxnfeval, nfeval, upd1, yksk, yrsr, sk, yk, sr, yr, lreset, function, state, xscale, xoffset, fscale, pivot, accuracy, gnorm, xnorm, low, up); if(frc == -1) { rc = TNC_ENOMEM; break; } if(frc) { rc = TNC_USERABORT; break; } if(!newcon) { if(!lreset) { /* Compute the accumulated step and its corresponding gradient difference. */ dxpy1(n, sk, sr); dxpy1(n, yk, yr); icycle++; } else { /* Initialize the sum of all the changes */ dcopy1(n, sk, sr); dcopy1(n, yk, yr); fLastReset = f; icycle = 1; } } dcopy1(n, g, oldg); oldf = f; oldgtp = ddot1(n, pk, g); /* Maximum unconstrained step length */ ustpmax = stepmx / (dnrm21(n, pk) + DBL_EPSILON); /* Maximum constrained step length */ spe = stepMax(ustpmax, n, x, pk, pivot, low, up, xscale, xoffset); if(spe > 0.0) { ENUM_LS_RC lsrc; /* Set the initial step length */ alpha = initialStep(f, _fmin / (fscale), oldgtp, spe); /* Perform the linear search */ lsrc = linearSearch(n, function, state, low, up, xscale, xoffset, fscale, pivot, eta, ftol, spe, pk, x, f, alpha, gfull, maxnfeval, nfeval); if(lsrc == LS_ENOMEM) { rc = TNC_ENOMEM; break; } if(lsrc == LS_USERABORT) { rc = TNC_USERABORT; break; } if(lsrc == LS_FAIL) { rc = TNC_LSFAIL; break; } /* If we went up to the maximum unconstrained step, increase it */ if(alpha >= 0.9 * ustpmax) { stepmx *= 1e2; if(bool(messages & TNC_MSG_INFO)) { printf("tnc: stepmx = %g", stepmx); } } /* If we went up to the maximum constrained step, a new constraint was encountered */ if(alpha - spe >= -DBL_EPSILON * 10.0) { newcon = TNC_TRUE; } else { /* Break if the linear search has failed to find a lower point */ if(lsrc != LS_OK) { if(lsrc == LS_MAXFUN) { rc = TNC_MAXFUN; } else { rc = TNC_LSFAIL; } break; } newcon = TNC_FALSE; } } else { /* Maximum constrained step == 0.0 => new constraint */ newcon = TNC_TRUE; } if(newcon) { if(!addConstraint(n, x, pk, pivot, low, up, xscale, xoffset)) { if(nfeval == oldnfeval) { rc = TNC_NOPROGRESS; break; } } fLastConstraint = f; } (niter)++; /* Invoke the callback function */ if(callback) { dcopy1(n, x, temp); unscalex(n, temp, xscale, xoffset); callback(temp, state); } /* Set up parameters used in convergence and resetting tests */ difold = difnew; difnew = oldf - f; /* If this is the first iteration of a new cycle, compute the percentage reduction factor for the resetting test */ if(icycle == 1) { if(difnew > difold * 2.0) { epsred += epsred; } if(difnew < difold * 0.5) { epsred *= 0.5; } } dcopy1(n, gfull, g); scaleg(n, g, xscale, fscale); dcopy1(n, g, temp); project(n, temp, pivot); gnorm = dnrm21(n, temp); /* Reset pivot */ remcon = removeConstraint(oldgtp, gnorm, pgtol * (fscale), f, fLastConstraint, g, pivot, n); /* If a constraint is removed */ if(remcon) { /* Recalculate gnorm and reset fLastConstraint */ dcopy1(n, g, temp); project(n, temp, pivot); gnorm = dnrm21(n, temp); fLastConstraint = f; } if(!remcon && !newcon) { /* No constraint removed & no new constraint : tests for convergence */ if(fabs(difnew) <= ftol * (fscale)) { if(bool(messages & TNC_MSG_INFO)) { printf( "tnc: |fn-fn-1] = %g -> convergence", fabs(difnew) / (fscale)); } rc = TNC_FCONVERGED; break; } if(alpha * dnrm21(n, pk) <= xtol) { if(bool(messages & TNC_MSG_INFO)) { printf( "tnc: |xn-xn-1] = %g -> convergence", alpha * dnrm21(n, pk)); } rc = TNC_XCONVERGED; break; } } project(n, g, pivot); if(bool(messages & TNC_MSG_ITER)) { printCurrentIteration(n, f / fscale, gfull, niter, nfeval, pivot); } /* Compute the change in the iterates and the corresponding change in the gradients */ if(!newcon) { yk = g - oldg; sk = alpha*pk; /* Set up parameters used in updating the preconditioning strategy */ yksk = ddot1(n, yk, sk); if(icycle == (n - 1) || difnew < epsred * (fLastReset - f)) { lreset = TNC_TRUE; } else { yrsr = ddot1(n, yr, sr); if(yrsr <= 0.0) { lreset = TNC_TRUE; } else { lreset = TNC_FALSE; } } upd1 = TNC_FALSE; } } if(bool(messages & TNC_MSG_ITER)) { printCurrentIteration(n, f / fscale, gfull, niter, nfeval, pivot); } /* Unscaling */ unscalex(n, x, xscale, xoffset); coercex(n, x, low, up); (f) /= fscale; return rc;
minimize_tncメソッドは、アルゴリズムの主要な反復ループを実装しています。このメソッドでは、探索方向の計算、最適なステップ長を決定するためのラインサーチの実行、さらに上限または下限の境界に到達した変数を追跡するための処理を管理します。また、このメソッドには、制約条件を考慮するための勾配の射影処理や、収束速度を向上させるための前処理戦略の更新機能も含まれています。このループは、勾配ノルムが特定のしきい値を下回るなどの収束条件が満たされるか、または関数評価回数の最大値に達するまで継続されます。
コードの最後のセクションには、ベクトル操作および制約管理に使用されるさまざまなユーティリティメソッドが含まれています。scalex、unscalex、scalegなどのメソッドは、ユーザーが扱う座標空間と内部で使用されるスケーリング済み空間との間でデータを変換する役割を担います。
/* Unscale x */ void unscalex(int n, vector& x, const vector& xscale, const vector& xoffset) { x = x*xscale+xoffset; } /* Scale x */ void scalex(int n, vector& x, const vector& xscale, const vector& xoffset) { x = (x-xoffset)/xscale; } /* Scale g */ void scaleg(int n, vector& g, const vector& xscale, double& fscale) { g*=xscale*fscale; }
projectメソッドおよびaddConstraintメソッドは、問題における境界条件を処理します。これらのメソッドは、勾配成分を手動でゼロに設定したり、許容範囲外へ移動しようとする変数を調整したりする役割を担っています。これにより、最適化処理全体を通じて解が常に実行可能領域内に保たれます。
void project(int n, vector& x, const int &pivot[]) { int i; for(i = 0; i < n; i++) { if(pivot[i] != 0) { x[i] = 0.0; } } } /* * Set x[i] = 0.0 if direction i is constant */ void projectConstants(int n, vector& x, const vector& xscale) { int i; for(i = 0; i < n; i++) { if(xscale[i] == 0.0) { x[i] = 0.0; } } }
目的関数ラッパー
実装の第2のコンポーネントでは、目的関数とその1次および2次導関数を処理します。これらはすべて、別途用意されたヘッダーファイルnum_diff.mqh内で定義されています。このヘッダーは、数値微分処理の設定項目として機能する複数の列挙型から始まります。これらの列挙型(例:ENUM_SCHEME_DIRECTION、ENUM_DIFF_POINTS、ENUM_HESS_DIFF_POINTS)を使用することで、ユーザーは片側または両側の探索方向スキームを選択したり、勾配およびヘッセ行列を推定するために使用する点の数を指定したりできます。これらの設定値によって、最終的には導関数近似における計算速度と数学的精度のバランスが決定されます。
//+------------------------------------------------------------------+ //| directional scheme options | //+------------------------------------------------------------------+ enum ENUM_SCHEME_DIRECTION { SCHEME_1=0,//1 sided SCHEME_2//2 sided }; //+------------------------------------------------------------------+ //| num points of evaluation | //+------------------------------------------------------------------+ enum ENUM_DIFF_POINTS { GRAD_POINT_2=0,//2-point GRAD_POINT_3,//3-point GRAD_POINT_CS,//complex GRAD_POINT_CALLABLE//callable }; //+------------------------------------------------------------------+ //| num points of evaluation | //+------------------------------------------------------------------+ enum ENUM_HESS_DIFF_POINTS { HESS_POINT_2=0,//2-point HESS_POINT_3,//3-point HESS_POINT_CS,//complex HESS_POINT_HESS_STRATEGY,//hessian update strategy HESS_POINT_CALLABLE//callable };
続いて、コードではObjReturn構造体とIObjectiveインターフェースが導入されており、最適化問題に対する標準化されたフレームワークを構築しています。ObjReturn構造体は、関数値と勾配ベクトルを同時に保持するために設計された単純なコンテナです。これにより、関数値と勾配を取得するために個別の冗長な関数呼び出しをおこなう必要がなくなります。IObjectiveインターフェースは、ソルバーに渡される任意の目的関数が一貫した構造に従うことを保証します。このインターフェースでは、目的関数値の直接計算と、関数値および勾配を組み合わせて返す処理の両方を実装することが要求されます。
//+------------------------------------------------------------------+ //|struct objective function return | //+------------------------------------------------------------------+ struct ObjReturn { double f; vector g; ObjReturn(void) { f = double(0); g = vector::Zeros(0); } ObjReturn(ObjReturn& other) { f = other.f; g = other.g; } void operator=(ObjReturn& other) { f = other.f; g = other.g; } }; //+------------------------------------------------------------------+ //|IObjective provides the base interface for an objective function | //|that will be provided to a minimizer routine | //+------------------------------------------------------------------+ interface IObjective { //---the objective function vector objective_function(vector& x); ObjReturn fun_and_grad(vector& x); };
GradDiffOptions構造体およびHessDiffOptions構造体は、有限差分計算に必要なパラメータを格納する設定用パッケージとして機能します。これらの構造体には、選択された推定手法、相対および絶対ステップサイズ、さらに境界制約が保持されます。これらの変数をまとめて管理することで、関数シグネチャを複雑化させることなく、微分計算の設定情報をソルバーの各レイヤー間で受け渡すことができます。
//+------------------------------------------------------------------+ //|differentiation options | //+------------------------------------------------------------------+ struct GradDiffOptions { ENUM_DIFF_POINTS method; vector rel_step; vector abs_step; matrix bounds; GradDiffOptions(void) { method = WRONG_VALUE; rel_step = abs_step = vector::Zeros(0); bounds = matrix::Zeros(0,0); } GradDiffOptions(ENUM_DIFF_POINTS m, vector& relstep, vector& absstep, matrix& bnds) { method = m; rel_step = relstep; abs_step = absstep; bounds = bnds; } GradDiffOptions(GradDiffOptions& other) { method = other.method; rel_step = other.rel_step; abs_step = other.abs_step; bounds = other.bounds; } void operator=(GradDiffOptions& other) { method = other.method; rel_step = other.rel_step; abs_step = other.abs_step; bounds = other.bounds; } }; //+------------------------------------------------------------------+ //|differentiation options | //+------------------------------------------------------------------+ struct HessDiffOptions { bool as_linear_operator; ENUM_HESS_DIFF_POINTS method; vector rel_step; vector abs_step; HessDiffOptions(void) { method = WRONG_VALUE; as_linear_operator = false; rel_step = abs_step = vector::Zeros(0); } HessDiffOptions(ENUM_HESS_DIFF_POINTS m, vector& relstep, vector& absstep, bool aslinearoperator) { method = m; rel_step = relstep; abs_step = absstep; as_linear_operator = aslinearoperator; } HessDiffOptions(HessDiffOptions& other) { method = other.method; rel_step = other.rel_step; abs_step = other.abs_step; as_linear_operator = other.as_linear_operator; } void operator=(HessDiffOptions& other) { method = other.method; rel_step = other.rel_step; abs_step = other.abs_step; as_linear_operator = other.as_linear_operator; } }
CFunctorクラスは、最適化状態を管理する中心的なクラスであり、関数値、勾配、およびヘッセ行列の計算とキャッシュ処理のロジックを実装しています。このクラスには、現在の位置に対して既に計算済みの値が存在するかを確認する内部メソッドが含まれており、コストの高い重複計算を回避します。また、このクラスは、これまでに見つかった最良の関数値とその位置を追跡します。さらに、ユーザーが提供した導関数と、自動的に実行される有限差分による推定処理を切り替えることができるラッパーとして機能します。
//+---------------------------------------------------------------------------+ //|function objective representing the objective function and its derivatives.| //+---------------------------------------------------------------------------+ class CFunctor:public IObjective { protected: vector m_xp; vector m_x; ulong m_n; matrix m_H; int m_nfev,m_ngev,m_nhev; bool m_fupdated,m_gupdated,m_hupdated; double m_lowest_f,m_f; vector m_lowest_x,m_g; GradDiffOptions m_grad_options; HessDiffOptions m_hess_options; void update_fun(void) { if(!m_fupdated) { double fx = wrapped_fun(m_x); if(fx<m_lowest_f) { m_lowest_f = fx; m_lowest_x = m_x; } m_f = fx; m_fupdated = true; } } void update_grad(void) { if(!m_gupdated) { if(m_grad_options.method!=GRAD_POINT_CALLABLE) update_fun(); vector ff(1); ff[0] = m_f; m_g = wrapped_grad(m_x,ff); m_gupdated = true; } } void update_hess(void) { if(!m_hupdated) { if(m_hess_options.method != HESS_POINT_CALLABLE) { update_grad(); m_H = wrapped_hess(m_x,m_g); } else { vector a = vector::Zeros(0); m_H = wrapped_hess(m_x,a); } m_hupdated = true; } } void update_x(vector& x) { m_x = x; m_fupdated = m_hupdated = m_gupdated = false; } public: CFunctor(void) { m_fupdated = m_hupdated = m_gupdated = false; m_lowest_f = DBL_MAX; m_lowest_x = m_g = vector::Zeros(0); m_H = matrix::Zeros(0,0); m_nfev = m_ngev = m_nhev = 0; m_grad_options.method = GRAD_POINT_2; m_hess_options.method = HESS_POINT_CALLABLE; m_hess_options.as_linear_operator = true; } ~CFunctor(void) { } void setGradOption(ENUM_DIFF_POINTS grad) { m_grad_options.method = grad; } void setAbsoluteStep(vector& epsilon) { m_grad_options.abs_step = epsilon; m_hess_options.abs_step = epsilon; } void setBounds(matrix& finite_bounds) { m_grad_options.bounds = finite_bounds; } void setRelativeStep(vector& finite_diff_rel_step) { m_grad_options.rel_step = finite_diff_rel_step; m_hess_options.rel_step = finite_diff_rel_step; } void setHessOption(ENUM_HESS_DIFF_POINTS hess) { m_hess_options.method = hess; } bool initialize(vector& x) { m_x = x; m_xp = m_x; m_n = x.Size(); if(m_grad_options.method != GRAD_POINT_CALLABLE && m_hess_options.method != HESS_POINT_CALLABLE) { Print(__FUNCTION__, "Whenever the gradient is estimated via " "finite-differences, it is required that" " the Hessian " "be estimated using one of the " "quasi-Newton strategies."); return false; } double check = orig_fun(m_x); if(MathClassify(check)!=FP_NORMAL) { Print(__FUNCTION__," check the implementation of the objective function, currently evaluates to an invalid number "); return false; } if(m_grad_options.method == GRAD_POINT_CALLABLE) { vector a = grad_fun(m_x); if(!a.Size()) { Print(__FUNCTION__, " check the implementation of the overriden gradient function, currently evaluates to an empty vector "); return false; } } update_fun(); update_grad(); if(m_hess_options.method == HESS_POINT_CALLABLE) { vector a = vector::Zeros(0); m_H = wrapped_hess(x,a); m_hupdated = true; } return true; } double wrapped_fun(vector& x) { m_nfev += 1; vector copy = x; return orig_fun(copy); } vector objective_function(vector& x) { vector r(1); r[0] = wrapped_fun(x); return r; } vector wrapped_grad(vector& x,vector& f0) { m_ngev += 1; vector copy = x; if(m_grad_options.method == GRAD_POINT_CALLABLE) return grad_fun(copy); IObjective* objective = GetPointer(this); matrix ad = approx_derivative(objective,copy,f0,m_grad_options.method,m_grad_options.rel_step,m_grad_options.abs_step,m_grad_options.bounds); //Print(__FUNCTION__, " - ", x, " - ", f0, " -> ", ad.Row(0), " | ", m_ngev); return ad.Row(0); } matrix wrapped_hess(vector& x, vector& f0) { m_nhev += 1; vector copy = x; if(m_hess_options.method == HESS_POINT_CALLABLE) return hess_fun(copy); IObjective* objective = GetPointer(this); return approx_derivative(objective,x,f0,m_grad_options.method,m_grad_options.rel_step,m_grad_options.abs_step,m_grad_options.bounds); } vector lower_bounds(void) { return m_grad_options.bounds.Col(0); } vector upper_bounds(void) { return m_grad_options.bounds.Col(1); } vector initial_params(void) { return m_xp; } virtual double orig_fun(vector& x) { return double("nan"); } virtual vector grad_fun(vector& x) { return vector::Zeros(0); } virtual matrix hess_fun(vector& x) { return matrix::Zeros(0,0); } ObjReturn fun_and_grad(vector& x) { vector dif = MathAbs(m_x - x); if(dif.Sum() >= DBL_EPSILON || dif.HasNan()) { update_x(x); } update_fun(); update_grad(); ObjReturn out; out.f = m_f; out.g = m_g; return out; } }
num_diff.mqhにおける次の主要なセクションでは、数値微分における最適なステップサイズを計算するユーティリティ関数eps_for_methodおよびcompute_absolute_stepが定義されています。コンピュータでは数値精度に限界があるため、ステップサイズが小さすぎる場合は丸め誤差が発生し、一方で大きすぎる場合は打ち切り誤差が発生します。これらの関数は、マシンイプシロンと使用する微分手法に基づいて、関数の傾きを評価する際に変数を移動させるための数学的に適切な距離を決定します。
//+----------------------------------------------------------------------------------+ //|Calculates relative EPS step to use for a given data type and numdiff step method.| //+----------------------------------------------------------------------------------+ double eps_for_method(ENUM_DIFF_POINTS method) { switch(method) { case GRAD_POINT_2: case GRAD_POINT_CS: return pow(2.220446049250313e-16,0.5); case GRAD_POINT_3: return pow(2.220446049250313e-16,(1./3.)); }; return DBL_EPSILON; } //+---------------------------------------------------------------------------------+ //|Computes an absolute step from a relative step for finite difference calculation.| //+---------------------------------------------------------------------------------+ vector compute_absolute_step(vector& rel_step,vector& x0, vector& f0, ENUM_DIFF_POINTS method) { vector signx0 = x0; vector abs_step = signx0; for(ulong i = 0; i<signx0.Size(); ++i) { if(x0[i] >= 0.) signx0[i]=1.*2-1; else signx0[i] = 0.0*2-1; } double rstep = eps_for_method(method); if(rel_step.Size()==0) for(ulong i = 0; i<abs_step.Size(); ++i) abs_step[i] = rstep*signx0[i]*MathMax(1.,fabs(x0[i])); else { abs_step = rstep*signx0*MathAbs(x0); vector dx = ((x0+abs_step) - x0); for(ulong i = 0; i<abs_step.Size(); ++i) if(dx[i] == 0.0) abs_step[i] = rstep*signx0[i]*MathMax(1.,fabs(x0[i])); } return abs_step; }
adjust_scheme_to_bounds関数は、数値的な評価点が許可された変数境界の外側へ移動することを防止します。アルゴリズムが境界付近にある場合、この関数は有限差分ステップの方向を自動的に反転させるか、中心差分スキームから片側差分スキームへ切り替えます。これにより、目的関数が無効な点で評価されることを防ぎ、制約付き最適化における安定性を確保します。
//+------------------------------------------------------------------+ //|Adjust final difference scheme to the presence of bounds. | //+------------------------------------------------------------------+ vector adjust_scheme_to_bounds(vector& x0, vector& h, int num_steps,ENUM_SCHEME_DIRECTION scheme, vector& lb, vector& ub, vector &one_sided) { switch(scheme) { case SCHEME_1: one_sided = vector::Ones(h.Size()); break; case SCHEME_2: one_sided = vector::Ones(h.Size()); h = MathAbs(h); break; } bool all_true = true; for(ulong i = 0; i<x0.Size(); ++i) if(lb[i] != -double("inf") || ub[i] != double("inf")) { all_true = false; break; } if(all_true) return h; vector h_total = h * double(num_steps); vector h_adjusted = h; vector lower_dist = x0 - lb; vector upper_dist = ub - x0; int forward,backward,fitting,violated,central, adjusted_central; forward = backward = violated = fitting = central = false; double x = 0.; double min_dist = 0.; switch(scheme) { case SCHEME_1: { for(ulong i = 0; i<h.Size(); ++i) { x = x0[i] + h_total[i]; violated = int(x<lb[i]|x>ub[i]); fitting = int(fabs(h_total[i])<=MathMax(lower_dist[i],upper_dist[i])); if(bool(violated & fitting)) h_adjusted[i]*=-1.; forward = int((upper_dist[i] >= lower_dist[i]) & ~fitting); if(forward) h_adjusted[i] = upper_dist[i]/double(num_steps); backward = int((upper_dist[i]<lower_dist[i]) & ~fitting); if(backward) h_adjusted[i] = -lower_dist[i]/double(num_steps); } } break; case SCHEME_2: { for(ulong i = 0; i<h.Size(); ++i) { central = int(((lower_dist[i]>=h_total[i]) & (upper_dist[i] >= h_total[i]))); forward = int(((upper_dist[i]>=lower_dist[i]) & ~central)); if(forward) { h_adjusted[i] = MathMin(h[i],0.5*upper_dist[i]/double(num_steps)); one_sided[i] = 1.; } backward = int(((upper_dist[i]<lower_dist[i]) & ~central)); if(backward) { h_adjusted[i] = -1.* MathMin(h[i],0.5*lower_dist[i]/double(num_steps)); one_sided[i] = 1.0; } min_dist = MathMin(upper_dist[i],lower_dist[i])/double(num_steps); adjusted_central = int((~central & (fabs(h_adjusted[i])<=min_dist))); if(adjusted_central) { h_adjusted[i] = min_dist; one_sided[i] = 0.; } } } break; } return h_adjusted; }
最後のセクションには、ヤコビ行列または勾配行列の構築に必要な数値計算を実行するdense_difference関数およびapprox_derivative関数が含まれています。dense_difference関数は、問題の各次元を順番に処理し、入力ベクトルに摂動を加え、その結果として生じる出力の変化を測定します。approx_derivativeメソッドは、高レベルの調整役として機能し、入力値の検証、ステップサイズの調整管理をおこなった上で、最終的な導関数行列を最適化ソルバーへ返します。
//+------------------------------------------------------------------+ //|dense difference | //+------------------------------------------------------------------+ matrix dense_difference(IObjective* fun, vector& x0, vector& f0, vector& h, vector& use_one_sided, ENUM_DIFF_POINTS method) { ulong m = f0.Size(); ulong n = x0.Size(); matrix j_transposed = matrix::Zeros(n,m); vector x1 = x0; vector x2 = x0; vector df = vector::Zeros(x0.Size()); for(ulong i = 0; i<h.Size(); ++i) { double dx = 1.e-12; if(method == GRAD_POINT_2) { x1[i] += h[i]; dx = x1[i] - x0[i]; df = fun.objective_function(x1) - f0; } else if(method == GRAD_POINT_3 && use_one_sided[i]!=0.0) { x1[i] += h[i]; x2[i] += 2. * h[i]; dx = x2[i] - x0[i]; df = -3.0 * f0 + 4 * fun.objective_function(x1) - fun.objective_function(x2); } else if(method == GRAD_POINT_3 && use_one_sided[i]==0.0) { x1[i] -= h[i]; x2[i] += h[i]; dx = x2[i] - x1[i]; df = fun.objective_function(x2) - fun.objective_function(x1); } j_transposed.Row(df/dx,i); x1[i] = x2[i] = x0[i]; } return j_transposed.Transpose(); } //+---------------------------------------------------------------------------------------+ //|Compute finite difference approximation of the derivatives of a vector-valued function.| //+---------------------------------------------------------------------------------------+ matrix approx_derivative(IObjective* fun,vector& x0,vector &f0,ENUM_DIFF_POINTS method, vector &rel_step,vector &abs_step, matrix& bounds/*sparsity,as linear_operator*/) { if(CheckPointer(fun)==POINTER_INVALID) { Print(__FUNCTION__, " fun variable is an invalid pointer "); return matrix::Zeros(0,0); } vector lb,ub; lb = bounds.Col(0); ub = bounds.Col(1); if(lb.Size()!=x0.Size() ||ub.Size()!=x0.Size()) { Print(__FUNCTION__, " inconsistent shaptes between bounds and x0 "); return matrix::Zeros(0,0); } if(!f0.Size()) f0 = fun.objective_function(x0); for(ulong i = 0; i<x0.Size(); ++i) if(x0[i]<lb[i] || x0[i]>ub[i]) { Print(__FUNCTION__, " x0 violates bound constraints "); return matrix::Zeros(0,0); } vector h; if(!abs_step.Size()) h = compute_absolute_step(rel_step,x0,f0,method); else { h = abs_step; vector signx0 = vector::Zeros(x0.Size()); for(ulong i = 0; i<x0.Size(); ++i) { if(x0[i]>=0.0) signx0[i] = 1.0*2.-1.; else signx0[i] = 0.0*2.-1.; if(((x0[i]+h[i]) - x0[i]) == 0.0) h[i] = eps_for_method(method)*signx0[i]*MathMax(1.,fabs(x0[i])); } } vector use_one_sided; switch(method) { case GRAD_POINT_2: h = adjust_scheme_to_bounds(x0,h,1,SCHEME_1,lb,ub,use_one_sided); break; case GRAD_POINT_3: h = adjust_scheme_to_bounds(x0,h,1,SCHEME_2,lb,ub,use_one_sided); break; case GRAD_POINT_CS: use_one_sided = vector::Zeros(x0.Size()); break; } return dense_difference(fun,x0,f0,h,use_one_sided,method); }
結論
本記事では、MQL5における切断ニュートン共役勾配法(TNC, Truncated Newton Conjugate-Gradient)最適化アルゴリズムの実装について解説しました。提供された実装は、ボックス制約を持つ最適化問題と、制約を持たない最適化問題の両方の最小化に対応しており、MQL5開発者にとって汎用的なツールとして利用できます。目的関数の計算を管理するコンポーネントは柔軟に設計されており、明示的に定義された勾配関数と、自動的な数値微分による勾配推定の両方をサポートしています。実装の検証では、複雑な数学的曲面を持つことで知られるローゼンブロック関数に適用し、困難な最適化問題においても探索能力を発揮できることを確認しました。さらに、L-BFGSの代替手法としてTNCオプティマイザをロジスティック回帰モデルへ統合する実用例も紹介しました。記事内で参照したすべてのコードは、以下に添付および一覧として示されています。読者は、主要なヘッダーファイルであるtnc.mqhおよびnum_diff.mqhをインクルードするだけで、すぐに利用可能な最適化ソルバーを使用できます。
| ファイル | 説明 |
|---|---|
| MQL5/experts/RosenBrock.mq5 | ストラテジーテスター上でローゼンブロック関数を評価するために使用されるEA |
| MQL5/files/iris.csv | LogisticRegressionスクリプトで使用されるirisデータセット |
| MQL5/include/tnc | tnc.mqhヘッダーファイルを含むフォルダ |
| MQL5/include/Regression | logistic.mqhヘッダーファイルを含むフォルダ |
| MQL5/include/np.mqh | 各種ベクトルおよび行列操作ユーティリティを含むヘッダーファイル |
| MQL5/include/num_diff.mqh | 数値微分の実装に使用されるユーティリティを含むヘッダーファイル |
| MQL5/scripts/LogisticRegression.mq5 | TNCソルバーを使用したロジスティック回帰の実装例を示すスクリプト |
| MQL5/scripts/TestTNC.mq5 | ローゼンブロック関数を用いTNCソルバーを評価するスクリプト |
MetaQuotes Ltdにより英語から翻訳されました。
元の記事: https://www.mql5.com/en/articles/21589
警告: これらの資料についてのすべての権利はMetaQuotes Ltd.が保有しています。これらの資料の全部または一部の複製や再プリントは禁じられています。
この記事はサイトのユーザーによって執筆されたものであり、著者の個人的な見解を反映しています。MetaQuotes Ltdは、提示された情報の正確性や、記載されているソリューション、戦略、または推奨事項の使用によって生じたいかなる結果についても責任を負いません。
グラフ理論:取引における深さ優先探索(DFS)の応用
MQL5取引ツール(第23回):カメラ制御対応DirectX 3Dグラフによる二項分布分析
MQL5経済指標カレンダーを用いたニュースフィルタリング(第2回):ニュースリリース中に管理ポジションを停止する
ラリー・ウィリアムズの『市場の秘密』(第14回):カスタムインジケータで隠れスマッシュデー反転を検出する
- 無料取引アプリ
- 8千を超えるシグナルをコピー
- 金融ニュースで金融マーケットを探索