Реализация алгоритма усечённого метода Ньютона с сопряжёнными градиентами в MQL5

Введение

При минимизации вещественнозначных функций разработчики MQL5 часто используют утилиты, предоставляемые библиотекой ALGLIB. Модуль minLBFGS библиотеки ALGLIB представляет собой надёжную реализацию метода Бройдена — Флетчера — Гольдфарба — Шанно с ограниченным объёмом памяти (LBFGS). К сожалению, в MQL5-версии библиотеки ALGLIB отсутствует версия этого метода для задач с ограничениями. Это упущение зачастую вынуждает пользователей прибегать к другим нелинейным оптимизаторам из библиотеки, что в некоторых случаях может оказаться излишним для простых задач с покомпонентными ограничениями на переменные. Нам нужен практичный решатель для задач с покомпонентными ограничениями на переменные, который можно было бы так же легко интегрировать в проект, как minLBFGS из библиотеки ALGLIB, и который поддерживал бы как аналитические градиенты, так и устойчивые численные приближения вблизи границ с учётом ограничений типа «ящик».

В данной статье мы представляем реализацию алгоритма усечённого метода Ньютона с сопряжёнными градиентами (TNC). Помимо рассмотрения теоретических основ, мы проверим способность данного метода находить глобальный минимум функции Розенброка, а также продемонстрируем его практическое применение, использовав его для реализации логистической регрессии в качестве альтернативы оптимизатору LBFGS.

Что такое TNC и чем он отличается от LBFGS?

Усечённый метод Ньютона с сопряжёнными градиентами (TNC) представляет собой алгоритм оптимизации второго порядка, аппроксимирующий работу метода Ньютона. В контексте оптимизации метод Ньютона представляет собой процедуру, используемую для нахождения минимума функции. В общем виде метод начинается с начального приближения, затем вычисляются градиент функции и матрица Гессе, после чего обновляется точка/вектор параметров. Этот процесс повторяется до достижения сходимости.

Градиент функции (первая производная) содержит информацию о наклоне, указывая направление наиболее крутого возрастания. Матрица Гессе (вторая производная) даёт информацию о форме, а именно о кривизне функции в заданной точке. Если матрица Гессе демонстрирует положительную кривизну во всех направлениях, это указывает на локальный минимум (нижнюю точку). Если кривизна отрицательна во всех направлениях, то эта точка является локальным максимумом (вершиной). Если матрица Гессе показывает изменение направления кривизны — сначала в одну сторону, затем в другую, — это указывает на седловую точку.

 

При использовании чистого метода Ньютона оптимизатору требуется вычисление матрицы Гессе, содержащей производные второго порядка, описывающие кривизну функции. В случае задач с большим количеством переменных вычисление и хранение этой матрицы требуют значительных вычислительных затрат. TNC решает эту проблему за счёт усечения ньютоновского процесса. Вместо точного решения уравнения Ньютона в нём используется внутренний цикл алгоритма сопряжённых градиентов (CG) для нахождения приближённого направления поиска. В этом внутреннем цикле требуются только произведения матрицы Гессе на вектор, которые можно эффективно оценить с помощью конечных разностей градиента. Это позволяет алгоритму учитывать информацию о кривизне, ни разу явно не формируя большую матрицу Гессе.

LBFGS также является квазиньютоновским методом, но основан на принципиально иной концепции. LBFGS хранит обновляемую историю градиентов и позиций за последние несколько итераций, чтобы построить неявную аппроксимацию обратной матрицы Гессе. Поскольку LBFGS повторно использует информацию из предыдущих шагов, он часто работает чрезвычайно быстро и требует небольших вычислительных затрат на каждой итерации, что делает его отраслевым стандартом для гладких функций без патологий в машинном обучении.

Обзор API TNC

Код, реализующий оптимизатор TNC, можно интегрировать в любую программу на MQL5 посредством включения заголовочного файла tnc.mqh. Использование API TNC состоит из четырёх этапов:

• Сначала мы задаём целевую функцию, её якобиан и все покомпонентные ограничения на переменные.
• Во-вторых, мы настраиваем сам оптимизатор.
• На третьем этапе мы запускаем оптимизатор.
• И последний шаг — получение результатов.

Интерфейс API состоит из двух классов. Класс CFunctor отвечает за все аспекты работы целевой функции, а CTruncNewtonCG представляет собой оптимизатор TNC.

Описание целевой функции выполняется путём определения класса, унаследованного от базового класса CFunctor. Целевую функцию предполагается определять путём переопределения метода orig_fun класса CFunctor. Если градиент или якобиан известны, их можно явно задать, переопределив метод 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, который принимает матрицу, где каждая строка соответствует одному измерению, а столбцы представляют нижнюю и верхнюю границы именно в таком порядке. Если метод grad_fun был переопределён, его необходимо включить, вызвав setGradOption с опцией GRAD_POINT_CALLABLE. В противном случае градиент будет вычисляться методом конечных разностей. Параметры, доступные в перечислении ENUM_DIFF_POINTS, позволяют пользователям выбирать метод конечных разностей, применяемый в случаях, когда функция градиента либо не задана, либо не включена. Более подробная информация об этих вариантах приведена в последующих разделах данной статьи.

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 для управления выводом в терминал: от одной строки на итерацию до подробной информации о масштабировании или только итоговой причины завершения.
• Метод SetMaxCGit устанавливает ограничение на количество внутренних итераций метода сопряжённых градиентов, чтобы внутренний цикл оставался эффективным.
• Метод SetMaxFunCalls задаёт жёсткое ограничение на общее количество вычислений значения целевой функции.
• Методы SetFtol и SetXtol устанавливают пороги сходимости для относительных изменений значения функции и вектора переменных соответственно. Параметр ftol (допуск по функции) отслеживает относительное изменение значения целевой функции между последовательными итерациями. Если уменьшение значения функции становится настолько малым, что падает ниже этого порога, оптимизатор считает, что достиг плато, где дальнейшие итерации дадут лишь пренебрежимо малые улучшения. Если установить это значение слишком высоким, возникнет «преждевременная сходимость»: алгоритм остановится до достижения истинного минимума. Слишком низкое значение, напротив, может привести к лишним вычислениям, которые фактически сводятся к борьбе с численным шумом. Параметр xtol (допуск по шагу) ориентирован на перемещение в пространстве поиска, а не на значение функции. Он измеряет евклидово расстояние между текущим и предыдущим положением; если оптимизатор вносит лишь микроскопические корректировки в переменные, срабатывает остановка процесса. Это особенно полезно в «плоских» областях функции, где её значение почти не меняется, но иначе алгоритм мог бы бесцельно блуждать в течение многих итераций, фактически не улучшая решение.
• Метод SetPGtol определяет пороговое значение для величины проецированного градиента. Алгоритм останавливается, если градиент опускается ниже этого уровня. pgtol (допуск по проецированному градиенту) — это, пожалуй, наиболее строгое с математической точки зрения условие остановки. Он измеряет крутизну локального наклона с учётом ограничений; нулевой градиент указывает на стационарную точку, что является теоретическим условием локального минимума. Очень малое значение pgtol заставляет алгоритм точно находить нижнюю точку долины, тогда как большее значение позволяет решателю завершить работу, как только он окажется достаточно близко к области оптимума.
• Метод SetAccuracy задаёт машинную точность или оценочную погрешность функции для выбора шага при вычислении градиента.
• Метод SetEta принимает значение от 0 до 1 для поиска по направлению, задавая компромисс между уменьшением значения функции и эффективностью шага.

Запуск минимизатора осуществляется путём вызова метода Minimize класса CTruncNewtonCG.

int return_code = tnc_minim.Minimize(sf);

В качестве единственного аргумента ему требуется экземпляр объекта класса, производного от CFunctor. Результаты процесса оптимизации затем можно получить с помощью следующих методов  CTruncNewtonCG:

• Метод Solution предоставляет переменные, соответствующие глобальному минимуму.
• Метод ObjectiveResult возвращает значение целевой функции в точке глобального минимума.
• Метод ObjectiveGradient возвращает значения градиента в точке глобального минимума.
• NumFeval возвращает целочисленное значение — количество раз, когда вычислялось значение целевой функции.
• Метод NumIters возвращает целочисленное значение, соответствующее количеству итераций процесса оптимизации.

В следующем разделе мы проверим нашу реализацию на практике.

Тест оптимизации с ограничениями

В классической литературе по оптимизации функция Розенброка является распространённой тестовой функцией для оценки производительности алгоритмов оптимизации. Она считается сложной задачей минимизации, поскольку имеет плоскую узкую долину, содержащую глобальный минимум, что может затруднять эффективную сходимость многих алгоритмов. Функции с подобной геометрией широко распространены в таких областях, как экономика и машинное обучение.

Мы можем визуализировать эту функцию в тестере стратегий MetaTrader 5 с помощью возможности вычисления математических функций. Советник (Expert Advisor, 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);
  }
//+------------------------------------------------------------------+

Запуск советника в тестере показывает, что глобальный минимум находится в точке [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 из-за заданных допусков. Чтобы продемонстрировать более практический пример применения, теперь мы заменим минимизатор LBFGS в реализации логистической регрессии на наш решатель TNC.

Реализация логистической регрессии на основе TNC

Класс CLogit для логистической регрессии ранее полагался на реализацию LBFGS из ALGLIB как на единственный механизм оптимизации. В этом разделе мы добавляем в класс оптимизатор 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
  };

Кроме того, в заголовочный файл добавлен класс CLogitFunctor, производный от CFunctor. Две закрытые функции-члены — fit_lbfgs и 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 проверяет работу класса CLogit на классическом наборе данных Iris. Программа подбирает модель для части данных с использованием обоих оптимизаторов, а полученные параметры модели выводятся в журнал терминала.

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, и универсального вспомогательного компонента, который определяет обёртки для целевой функции, её производной и матрицы Гессе. Сначала мы сосредоточимся на файле tnc.mqh, поскольку в нём реализован основной алгоритм. Код в файле 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 выступают в качестве связующего звена между низкоуровневым решателем и состоянием 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 &gtest1, double &gtest2, 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 &gtest1, double &gtest2, 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;
           }
        }
     }

Обёртка целевой функции

Второй компонент реализации обрабатывает целевую функцию и её производные первого и второго порядка, которые определены в отдельном заголовочном файле 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);
  } 

Заключение

В данной статье представлена реализация алгоритма оптимизации — усечённого метода Ньютона с сопряжёнными градиентами (TNC) — на языке MQL5. Представленная реализация поддерживает решение задач минимизации как с ограничениями типа «ящик», так и без них, предлагая разработчикам MQL5 универсальный инструмент. Компонент, отвечающий за вычисление целевой функции, отличается гибкостью и поддерживает как явно заданные функции градиента, так и автоматическое численное дифференцирование. Мы проверили работоспособность реализации, применив её к сложной функции Розенброка, что продемонстрировало её способность работать со сложными математическими поверхностями. В заключение мы продемонстрировали практическое применение, интегрировав оптимизатор TNC в модель логистической регрессии в качестве альтернативы методу LBFGS. Весь код, упомянутый в статье, прилагается и приведён ниже. Читатели получают готовый к использованию оптимизатор: достаточно подключить основные заголовочные файлы tnc.mqh и num_diff.mqh.

Файл	Описание
MQL5/experts/RosenBrock.mq5	Советник (EA) для вычисления значения функции Розенброка в тестере стратегий.
MQL5/files/iris.csv	Набор данных Iris, используемый в скрипте LogisticRegression.
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 на функции Розенброка.