logoChat01.ai

обнаружил, что более простая форрмула градиеного с...

생성일: 2026년 3월 4일

답변에 사용된 모델 GPT-5.2 Thinking by Chat01

질문

обнаружил, что более простая форрмула градиеного спуска в локальной доводке даёт лучший результат, при коде, котороый сейчас лучшая точность и скорость, код локальной доводки сейчас: " auto evalAt = [&](const float t) -> float {
map.map01ToPoint(t, q_local);
float f = cost(q_local, x, y);
if (f < fmaf(bestF, adaptive_coeff, 0.0f)) {
// Параметры Armijo
const float c1 = 1e-4f; // достаточное убывание
const float tau = 0.5f; // коэффициент уменьшения шага
const int max_outer_iters = 50.0f * (1.0f + 0.65f * p); // максимум итераций градиентного спуска
const int max_backtrack = 20.0f * (1.0f + 0.65f * p); // максимум дроблений шага

text
		float eta = 1.0f;                  // начальный шаг

		int outer = 0;
		while (outer < max_outer_iters) {
			// ---- Вычисляем градиент в текущей точке q_local ----
			float acc = 0.0f;
			int ii = 0;

#pragma loop ivdep
while (ii < n) {
acc = fmaf(q_local[ii], 1.0f, acc);
phi[ii] = acc;
++ii;
}
FABE13_SINCOS(phi, s_arr, c_arr, n);

text
			float as = 0.0f;
			float ac = 0.0f;
			int k = n - 1;
			while (k >= 0) {
				const float Lk = cost.variableLen ? q_local[n + k] : 1.0f;
				as = fmaf(Lk, s_arr[k], as);
				ac = fmaf(Lk, c_arr[k], ac);
				sum_s[k] = as;
				sum_c[k] = ac;
				--k;
			}

			const float dx = fmaf(x, 1.0f, -cost.targetX);
			const float dy = fmaf(y, 1.0f, -cost.targetY);
			const float dist = sqrtf(fmaf(dx, dx, dy * dy));
			const float inv_dist = 1.0f / dist;

			// Массив градиентов и норма
			float grad[64];
			float grad_norm2 = 0.0f;

			// Градиент по углам
			int i = 0;

#pragma loop ivdep
while (i < n) {
float gpen = 0.0f;
// Производная штрафа за превышение minTheta
{
const float ai = fabsf(q_local[i]);
const float v = fmaf(ai, 1.0f, -cost.minTheta);
if (v > 0.0f) {
const float scale_arg = fmaf(2.0f / cost.minTheta, v * 0.69314718f, 0.0f);
const float exp_val = fmaf(scale_arg,
fmaf(scale_arg,
fmaf(scale_arg,
fmaf(scale_arg, 0.00833333377f, 0.0416666679f),
0.16666667f),
0.5f),
1.0f);
const float dpen = fmaf(cost.sharpW, exp_val * (1.38629436f / cost.minTheta), 0.0f);
gpen = fmaf(dpen, copysignf(1.0f, q_local[i]), gpen);
}
}
// Производная арк-штрафа
{
const float tsg = fmaf(-q_local[i], cost.archBiasK, 0.0f);
const float exp_arg = -tsg;
const float exp_val = fmaf(exp_arg,
fmaf(exp_arg,
fmaf(exp_arg,
fmaf(exp_arg, 0.00833333377f, 0.0416666679f),
0.16666667f),
0.5f),
1.0f);
const float sig = 1.0f / (exp_val + 1.0f);
gpen = fmaf(-cost.archBiasW * cost.archBiasK, sig, gpen);
}

text
				const float g_main = fmaf(dx, -sum_s[i], dy * sum_c[i]) * inv_dist;
				float gi = g_main + gpen;
				grad[i] = gi;
				grad_norm2 = fmaf(gi, gi, grad_norm2);
				++i;
			}

			// Градиент по длинам
			if (cost.variableLen) {
				int j = 0;

#pragma loop ivdep
while (j < n) {
const float gi = fmaf(dx, c_arr[j], dy * s_arr[j]) * inv_dist;
grad[n + j] = gi;
grad_norm2 = fmaf(gi, gi, grad_norm2);
++j;
}
}

text
			if (grad_norm2 < 1e-12f) break;

			// ---- Линейный поиск с условием Армихо ----
			float eta_trial = eta;
			float f_current = f;
			bool found = false;
			int backtrack = 0;
			while (backtrack < max_backtrack) {
				// Пробная точка
				i = 0;

#pragma loop ivdep
while (i < n) {
q_try[i] = fmaf(-eta_trial, grad[i], q_local[i]);
const float lo = (i == 0) ? -1.0471975511965977f : -2.6179938779914944f;
const float hi = 2.6179938779914944f;
if (q_try[i] < lo) q_try[i] = lo;
else if (q_try[i] > hi) q_try[i] = hi;
++i;
}
if (cost.variableLen) {
i = 0;
#pragma loop ivdep
while (i < n) {
q_try[n + i] = fmaf(-eta_trial, grad[n + i], q_local[n + i]);
if (q_try[n + i] < 0.5f) q_try[n + i] = 0.5f;
else if (q_try[n + i] > 2.0f) q_try[n + i] = 2.0f;
++i;
}
}

text
				float x2, y2;
				float f_try = cost(q_try, x2, y2);

				// Условие Армихо
				if (f_try <= f_current - c1 * eta_trial * grad_norm2) {
					// Шаг принят
					memcpy(q_local, q_try, static_cast<size_t>(dim) * sizeof(float));
					f = f_try;
					x = x2;
					y = y2;
					eta = eta_trial;        // запоминаем удачный шаг
					found = true;
					break;
				}

				eta_trial *= tau;
				++backtrack;
			}

			if (!found) break;   // не удалось сделать шаг — выход

			++outer;
		}

		// Обновляем лучшую точку
		if (f < bestF) {
			bestF = f;
			bestQ.assign(q_local, q_local + dim);
			bestX = x;
			bestY = y;
			no_improve = 0;
		}
		else {
			++no_improve;
		}
	}
	return f;
	};", но в этой доводке нет L-BFGS, может ли он быть полезен? тот ли это механизм, который точно не помешает при локальной доводке? ранее локальная доводка выглядела так: "auto evalAt = [&](const float t) -> float {
	map.map01ToPoint(t, q_local);
	float f = cost(q_local, x, y);
	if (f < fmaf(bestF, adaptive_coeff, 0.0f)) {
		// ----- Вычисляем градиент для оценки нормы и выбора метода -----
		float acc_g = 0.0f;
		int ii_g = 0;

#pragma loop ivdep
while (ii_g < n) {
acc_g = fmaf(q_local[ii_g], 1.0f, acc_g);
phi[ii_g] = acc_g;
++ii_g;
}
FABE13_SINCOS(phi, s_arr, c_arr, n);
float as_g = 0.0f;
float ac_g = 0.0f;
int kk_g = n - 1;
while (kk_g >= 0) {
const float Lk_g = cost.variableLen ? q_local[n + kk_g] : 1.0f;
as_g = fmaf(Lk_g, s_arr[kk_g], as_g);
ac_g = fmaf(Lk_g, c_arr[kk_g], ac_g);
sum_s[kk_g] = as_g;
sum_c[kk_g] = ac_g;
--kk_g;
}
const float dx_g = fmaf(x, 1.0f, -cost.targetX);
const float dy_g = fmaf(y, 1.0f, -cost.targetY);
const float dist_g = sqrtf(fmaf(dx_g, dx_g, dy_g * dy_g));
const float inv_dist_g = 1.0f / dist_g;
float grad_norm2_init = 0.0f;
int i_g = 0;
#pragma loop ivdep
while (i_g < n) {
float gpen_g = 0.0f;
// Производная штрафа за превышение minTheta
{
const float ai = fabsf(q_local[i_g]);
const float v = fmaf(ai, 1.0f, -cost.minTheta);
if (v > 0.0f) {
const float scale_arg = fmaf(2.0f / cost.minTheta, v * 0.69314718f, 0.0f);
const float exp_val = fmaf(scale_arg,
fmaf(scale_arg,
fmaf(scale_arg,
fmaf(scale_arg, 0.00833333377f, 0.0416666679f),
0.16666667f),
0.5f),
1.0f);
const float dpen = fmaf(cost.sharpW, exp_val * (1.38629436f / cost.minTheta), 0.0f);
gpen_g = fmaf(dpen, copysignf(1.0f, q_local[i_g]), gpen_g);
}
}
// Производная арк-штрафа
{
const float tsg = fmaf(-q_local[i_g], cost.archBiasK, 0.0f);
const float exp_arg = -tsg;
const float exp_val = fmaf(exp_arg,
fmaf(exp_arg,
fmaf(exp_arg,
fmaf(exp_arg, 0.00833333377f, 0.0416666679f),
0.16666667f),
0.5f),
1.0f);
const float sig = 1.0f / (exp_val + 1.0f);
gpen_g = fmaf(-cost.archBiasW * cost.archBiasK, sig, gpen_g);
}
const float g_main_g = fmaf(dx_g, -sum_s[i_g], dy_g * sum_c[i_g]) * inv_dist_g;
float gi_g = g_main_g + gpen_g;
grad_norm2_init = fmaf(gi_g, gi_g, grad_norm2_init);
++i_g;
}
if (cost.variableLen) {
int j_g = 0;
#pragma loop ivdep
while (j_g < n) {
const float gi_g = fmaf(dx_g, c_arr[j_g], dy_g * s_arr[j_g]) * inv_dist_g;
grad_norm2_init = fmaf(gi_g, gi_g, grad_norm2_init);
++j_g;
}
}
// ----- Ваши эмпирические параметры (общие) -----
const float c1 = 1e-4f; // Armijo
const float c2 = 0.9f; // curvature (strong Wolfe)
const float tau = 0.5f; // уменьшение шага
const float expansion = 2.0f; // увеличение шага
const float adaptive_window = fmaf(1.0f / fmaf(bestF, adaptive_coeff, 0.0f),
fmaf(bestF, adaptive_coeff, -f), -1.0f);
const float k = fmaf(adaptive_window,
fmaf(adaptive_window,
fmaf(adaptive_window,
fmaf(adaptive_window, 0.00833333377f, 0.0416666679f),
0.16666667f),
0.5f),
1.0f);
const int max_outer_iters = static_cast<int>(fmaf(k, fmaf(sqrtf(dim_f), 30.0f, 0.0f), 0.0f));
const int max_zoom_iters = max_outer_iters >> 1;

text
		// ----- Выбор метода локальной доводки на основе глобального прогресса p -----
		const float grad_threshold = 1e-3f;
		bool use_lbfgs = (p > 0.9f && grad_norm2_init < grad_threshold);

		// ----- Состояние для Adam (используется только если use_lbfgs == false) -----
		float m[64] = { 0.0f };
		float v[64] = { 0.0f };
		int t_adam = 0;

		// ----- Состояние для L-BFGS (используется только если use_lbfgs == true) -----
		const int m_lbfgs = 5;
		float s_storage[32][5];
		float y_storage[32][5];
		float rho_storage[5];
		int lbfgs_iter = 0;

		int outer = 0;
		while (outer < max_outer_iters) {
			// ---------- Вычисляем градиент в текущей точке q_local ----------
			float acc = 0.0f;
			int ii = 0;

#pragma loop ivdep
while (ii < n) {
acc = fmaf(q_local[ii], 1.0f, acc);
phi[ii] = acc;
++ii;
}
FABE13_SINCOS(phi, s_arr, c_arr, n);

text
			float as = 0.0f;
			float ac = 0.0f;
			int kk = n - 1;
			while (kk >= 0) {
				const float Lk = cost.variableLen ? q_local[n + kk] : 1.0f;
				as = fmaf(Lk, s_arr[kk], as);
				ac = fmaf(Lk, c_arr[kk], ac);
				sum_s[kk] = as;
				sum_c[kk] = ac;
				--kk;
			}

			const float dx = fmaf(x, 1.0f, -cost.targetX);
			const float dy = fmaf(y, 1.0f, -cost.targetY);
			const float dist = sqrtf(fmaf(dx, dx, dy * dy));
			const float inv_dist = 1.0f / dist;

			float grad[64];
			float grad_norm2 = 0.0f;

			// Градиент по углам
			int i = 0;

#pragma loop ivdep
while (i < n) {
float gpen = 0.0f;
// Производная штрафа за превышение minTheta
{
const float ai = fabsf(q_local[i]);
const float v = fmaf(ai, 1.0f, -cost.minTheta);
if (v > 0.0f) {
const float scale_arg = fmaf(2.0f / cost.minTheta, v * 0.69314718f, 0.0f);
const float exp_val = fmaf(scale_arg,
fmaf(scale_arg,
fmaf(scale_arg,
fmaf(scale_arg, 0.00833333377f, 0.0416666679f),
0.16666667f),
0.5f),
1.0f);
const float dpen = fmaf(cost.sharpW, exp_val * (1.38629436f / cost.minTheta), 0.0f);
gpen = fmaf(dpen, copysignf(1.0f, q_local[i]), gpen);
}
}
// Производная арк-штрафа
{
const float tsg = fmaf(-q_local[i], cost.archBiasK, 0.0f);
const float exp_arg = -tsg;
const float exp_val = fmaf(exp_arg,
fmaf(exp_arg,
fmaf(exp_arg,
fmaf(exp_arg, 0.00833333377f, 0.0416666679f),
0.16666667f),
0.5f),
1.0f);
const float sig = 1.0f / (exp_val + 1.0f);
gpen = fmaf(-cost.archBiasW * cost.archBiasK, sig, gpen);
}

text
				const float g_main = fmaf(dx, -sum_s[i], dy * sum_c[i]) * inv_dist;
				float gi = g_main + gpen;
				grad[i] = gi;
				grad_norm2 = fmaf(gi, gi, grad_norm2);
				++i;
			}

			// Градиент по длинам
			if (cost.variableLen) {
				int j = 0;

#pragma loop ivdep
while (j < n) {
const float gi = fmaf(dx, c_arr[j], dy * s_arr[j]) * inv_dist;
grad[n + j] = gi;
grad_norm2 = fmaf(gi, gi, grad_norm2);
++j;
}
}

text
			if (grad_norm2 < 1e-12f) break;

			// ----- Выбор направления спуска в зависимости от выбранного метода -----
			float direction[64];
			if (!use_lbfgs) {
				// ---------- Adam ----------
				t_adam++;
				float pow_beta1_t = 1.0f;
				float pow_beta2_t = 1.0f;
				for (int p = 0; p < t_adam; ++p) {
					pow_beta1_t *= 0.9f;
					pow_beta2_t *= 0.999f;
				}
				float one_minus_beta1_t = 1.0f - pow_beta1_t;
				float one_minus_beta2_t = 1.0f - pow_beta2_t;

				i = 0;

#pragma loop ivdep
while (i < dim) {
float g = grad[i];
m[i] = fmaf(0.9f, m[i], 0.1f * g);
v[i] = fmaf(0.999f, v[i], 0.001f * g * g);

text
					float m_hat = m[i] / one_minus_beta1_t;
					float v_hat = v[i] / one_minus_beta2_t;

					direction[i] = -m_hat / (sqrtf(v_hat) + 1e-8f);
					++i;
				}
			}
			else {
				// ---------- L-BFGS ----------
				if (lbfgs_iter == 0) {
					i = 0;

#pragma loop ivdep
while (i < dim) {
direction[i] = -grad[i];
++i;
}
}
else {
float q[64];
i = 0;
#pragma loop ivdep
while (i < dim) {
q[i] = grad[i];
++i;
}

text
					int jj = lbfgs_iter - 1;
					while (jj >= 0) {
						float alpha_j = 0.0f;
						i = 0;

#pragma loop ivdep
while (i < dim) {
alpha_j = fmaf(s_storage[i][jj], q[i], alpha_j);
++i;
}
alpha_j *= rho_storage[jj];
i = 0;
#pragma loop ivdep
while (i < dim) {
q[i] = fmaf(-alpha_j, y_storage[i][jj], q[i]);
++i;
}
--jj;
}

text
					float scale = 1.0f;
					if (lbfgs_iter > 0) {
						float ys = 0.0f, yy = 0.0f;
						i = 0;

#pragma loop ivdep
while (i < dim) {
ys = fmaf(y_storage[i][lbfgs_iter - 1], s_storage[i][lbfgs_iter - 1], ys);
yy = fmaf(y_storage[i][lbfgs_iter - 1], y_storage[i][lbfgs_iter - 1], yy);
++i;
}
if (yy > 0.0f) scale = ys / yy;
}

					i = 0;

#pragma loop ivdep
while (i < dim) {
direction[i] = scale * q[i];
++i;
}

text
					jj = 0;
					while (jj < lbfgs_iter) {
						float beta = 0.0f;
						i = 0;

#pragma loop ivdep
while (i < dim) {
beta = fmaf(y_storage[i][jj], direction[i], beta);
++i;
}
beta *= rho_storage[jj];
i = 0;
#pragma loop ivdep
while (i < dim) {
direction[i] = fmaf(s_storage[i][jj], -beta, direction[i]);
++i;
}
++jj;
}

					i = 0;

#pragma loop ivdep
while (i < dim) {
direction[i] = -direction[i];
++i;
}
}
}

text
			// Проверка, что направление спуска
			float dir_dot_g = 0.0f;
			i = 0;

#pragma loop ivdep
while (i < dim) {
dir_dot_g = fmaf(direction[i], grad[i], dir_dot_g);
++i;
}
if (dir_dot_g >= 0.0f) {
i = 0;
#pragma loop ivdep
while (i < dim) {
direction[i] = -grad[i];
++i;
}
dir_dot_g = 0.0f;
i = 0;
#pragma loop ivdep
while (i < dim) {
dir_dot_g = fmaf(direction[i], grad[i], dir_dot_g);
++i;
}
}
float directional_deriv = dir_dot_g;

text
			// ---- Сохраняем текущую точку для линейного поиска ----
			float q_local_save[32];
			memcpy(q_local_save, q_local, dim * sizeof(float));
			float x_save = x, y_save = y;
			float f_save = f;

			// ---- Начальный шаг alpha из вашей эмпирической формулы ----
			float alpha = fmaf(sqrtf(fmaf(-dim_f, 2.0f, 23.0f)),
				fmaf(fmaf(adaptive_window + 1.0f,
					fmaf(adaptive_window + 1.0f,
						fmaf(adaptive_window + 1.0f,
							fmaf(adaptive_window + 1.0f, 0.00833333377f, 0.0416666679f),
							0.16666667f),
						0.5f),
					1.0f),
					0.215f, 0.7f),
				0.0f);
			float alpha_lo = 0.0f;
			float alpha_hi = 0.0f;
			float f_lo = f_save;
			bool found = false;
			int bracket_iter = 0;
			const int max_bracket = 20;

			// Фаза 1: Bracketing
			while (bracket_iter < max_bracket) {
				i = 0;

#pragma loop ivdep
while (i < n) {
q_try[i] = fmaf(alpha, direction[i], q_local_save[i]);
const float lo = (i == 0) ? -1.0471975511965977f : -2.6179938779914944f;
const float hi = 2.6179938779914944f;
if (q_try[i] < lo) q_try[i] = lo;
else if (q_try[i] > hi) q_try[i] = hi;
++i;
}
if (cost.variableLen) {
i = 0;
#pragma loop ivdep
while (i < n) {
q_try[n + i] = fmaf(alpha, direction[n + i], q_local_save[n + i]);
if (q_try[n + i] < 0.5f) q_try[n + i] = 0.5f;
else if (q_try[n + i] > 2.0f) q_try[n + i] = 2.0f;
++i;
}
}
float x_try, y_try;
float f_try = cost(q_try, x_try, y_try);

text
				if (f_try <= f_save + c1 * alpha * directional_deriv) {
					// Armijo выполнено, вычисляем градиент в пробной точке
					float acc_try = 0.0f;
					float phi_try[32];
					int ii_try = 0;

#pragma loop ivdep
while (ii_try < n) {
acc_try = fmaf(q_try[ii_try], 1.0f, acc_try);
phi_try[ii_try] = acc_try;
++ii_try;
}
float s_try[32], c_try[32];
FABE13_SINCOS(phi_try, s_try, c_try, n);

text
					float as_try = 0.0f, ac_try = 0.0f;
					float sum_s_try[32], sum_c_try[32];
					int k_try = n - 1;
					while (k_try >= 0) {
						const float Lk = cost.variableLen ? q_try[n + k_try] : 1.0f;
						as_try = fmaf(Lk, s_try[k_try], as_try);
						ac_try = fmaf(Lk, c_try[k_try], ac_try);
						sum_s_try[k_try] = as_try;
						sum_c_try[k_try] = ac_try;
						--k_try;
					}

					const float dx_try = fmaf(x_try, 1.0f, -cost.targetX);
					const float dy_try = fmaf(y_try, 1.0f, -cost.targetY);
					const float dist_try = sqrtf(fmaf(dx_try, dx_try, dy_try * dy_try));
					const float inv_dist_try = 1.0f / dist_try;

					float dot = 0.0f;
					i = 0;

#pragma loop ivdep
while (i < n) {
float gpen_try = 0.0f;
{
const float ai = fabsf(q_try[i]);
const float v = fmaf(ai, 1.0f, -cost.minTheta);
if (v > 0.0f) {
const float scale_arg = fmaf(2.0f / cost.minTheta, v * 0.69314718f, 0.0f);
const float exp_val = fmaf(scale_arg,
fmaf(scale_arg,
fmaf(scale_arg,
fmaf(scale_arg, 0.00833333377f, 0.0416666679f),
0.16666667f),
0.5f),
1.0f);
const float dpen = fmaf(cost.sharpW, exp_val * (1.38629436f / cost.minTheta), 0.0f);
gpen_try = fmaf(dpen, copysignf(1.0f, q_try[i]), gpen_try);
}
}
{
const float tsg = fmaf(-q_try[i], cost.archBiasK, 0.0f);
const float exp_arg = -tsg;
const float exp_val = fmaf(exp_arg,
fmaf(exp_arg,
fmaf(exp_arg,
fmaf(exp_arg, 0.00833333377f, 0.0416666679f),
0.16666667f),
0.5f),
1.0f);
const float sig = 1.0f / (exp_val + 1.0f);
gpen_try = fmaf(-cost.archBiasW * cost.archBiasK, sig, gpen_try);
}

text
						const float g_main_try = fmaf(dx_try, -sum_s_try[i], dy_try * sum_c_try[i]) * inv_dist_try;
						float gi_try = g_main_try + gpen_try;
						dot = fmaf(gi_try, direction[i], dot);
						++i;
					}
					if (cost.variableLen) {
						int j = 0;

#pragma loop ivdep
while (j < n) {
const float gi_try = fmaf(dx_try, c_try[j], dy_try * s_try[j]) * inv_dist_try;
dot = fmaf(gi_try, direction[n + j], dot);
++j;
}
}

text
					if (fabsf(dot) <= c2 * fabsf(directional_deriv)) {
						memcpy(q_local, q_try, dim * sizeof(float));
						f = f_try;
						x = x_try;
						y = y_try;
						found = true;
						break;
					}
					else {
						alpha_lo = alpha;
						f_lo = f_try;
						alpha *= expansion;
						++bracket_iter;
						continue;
					}
				}
				else {
					alpha_hi = alpha;
					alpha *= tau;
					++bracket_iter;
					continue;
				}
			}

			if (found) {
				// Обновление истории для L-BFGS (если используется)
				if (use_lbfgs) {
					// Вычисляем новый градиент в точке q_local
					float new_grad[64];
					float acc2 = 0.0f;
					int i2 = 0;

#pragma loop ivdep
while (i2 < n) {
acc2 = fmaf(q_local[i2], 1.0f, acc2);
phi[i2] = acc2;
++i2;
}
FABE13_SINCOS(phi, s_arr, c_arr, n);
float as2 = 0.0f, ac2 = 0.0f;
int k2 = n - 1;
while (k2 >= 0) {
const float Lk = cost.variableLen ? q_local[n + k2] : 1.0f;
as2 = fmaf(Lk, s_arr[k2], as2);
ac2 = fmaf(Lk, c_arr[k2], ac2);
sum_s[k2] = as2;
sum_c[k2] = ac2;
--k2;
}
const float dx2 = fmaf(x, 1.0f, -cost.targetX);
const float dy2 = fmaf(y, 1.0f, -cost.targetY);
const float dist2 = sqrtf(fmaf(dx2, dx2, dy2 * dy2));
const float inv_dist2 = 1.0f / dist2;
i2 = 0;
#pragma loop ivdep
while (i2 < n) {
float gpen2 = 0.0f;
{
const float ai = fabsf(q_local[i2]);
const float v = fmaf(ai, 1.0f, -cost.minTheta);
if (v > 0.0f) {
const float scale_arg = fmaf(2.0f / cost.minTheta, v * 0.69314718f, 0.0f);
const float exp_val = fmaf(scale_arg,
fmaf(scale_arg,
fmaf(scale_arg,
fmaf(scale_arg, 0.00833333377f, 0.0416666679f),
0.16666667f),
0.5f),
1.0f);
const float dpen = fmaf(cost.sharpW, exp_val * (1.38629436f / cost.minTheta), 0.0f);
gpen2 = fmaf(dpen, copysignf(1.0f, q_local[i2]), gpen2);
}
}
{
const float tsg = fmaf(-q_local[i2], cost.archBiasK, 0.0f);
const float exp_arg = -tsg;
const float exp_val = fmaf(exp_arg,
fmaf(exp_arg,
fmaf(exp_arg,
fmaf(exp_arg, 0.00833333377f, 0.0416666679f),
0.16666667f),
0.5f),
1.0f);
const float sig = 1.0f / (exp_val + 1.0f);
gpen2 = fmaf(-cost.archBiasW * cost.archBiasK, sig, gpen2);
}
const float g_main2 = fmaf(dx2, -sum_s[i2], dy2 * sum_c[i2]) * inv_dist2;
new_grad[i2] = g_main2 + gpen2;
++i2;
}
if (cost.variableLen) {
int j2 = 0;
#pragma loop ivdep
while (j2 < n) {
new_grad[n + j2] = fmaf(dx2, c_arr[j2], dy2 * s_arr[j2]) * inv_dist2;
++j2;
}
}

text
					float s_new[32], y_new[32];
					i = 0;

#pragma loop ivdep
while (i < dim) {
s_new[i] = q_local[i] - q_local_save[i];
y_new[i] = new_grad[i] - grad[i];
++i;
}

text
					float ys = 0.0f;
					i = 0;

#pragma loop ivdep
while (i < dim) {
ys = fmaf(y_new[i], s_new[i], ys);
++i;
}

text
					const float eps_lbfgs = 1e-8f;
					if (ys > eps_lbfgs) {
						int new_idx;
						if (lbfgs_iter < m_lbfgs) {
							new_idx = lbfgs_iter;
							lbfgs_iter++;
						}
						else {
							for (int idx = 0; idx < m_lbfgs - 1; ++idx) {
								for (int d = 0; d < dim; ++d) {
									s_storage[d][idx] = s_storage[d][idx + 1];
									y_storage[d][idx] = y_storage[d][idx + 1];
								}
								rho_storage[idx] = rho_storage[idx + 1];
							}
							new_idx = m_lbfgs - 1;
						}
						for (int d = 0; d < dim; ++d) {
							s_storage[d][new_idx] = s_new[d];
							y_storage[d][new_idx] = y_new[d];
						}
						rho_storage[new_idx] = 1.0f / ys;
					}
					else {
						lbfgs_iter = 0;
					}
				}
				++outer;
				continue;
			}

			// Если не нашли шаг, пробуем zoom
			if (alpha_hi == 0.0f) break;

			int zoom_iter = 0;
			while (zoom_iter < max_zoom_iters) {
				float alpha_new = (alpha_lo + alpha_hi) * 0.5f;

				i = 0;

#pragma loop ivdep
while (i < n) {
q_try[i] = fmaf(alpha_new, direction[i], q_local_save[i]);
const float lo = (i == 0) ? -1.0471975511965977f : -2.6179938779914944f;
const float hi = 2.6179938779914944f;
if (q_try[i] < lo) q_try[i] = lo;
else if (q_try[i] > hi) q_try[i] = hi;
++i;
}
if (cost.variableLen) {
i = 0;
#pragma loop ivdep
while (i < n) {
q_try[n + i] = fmaf(alpha_new, direction[n + i], q_local_save[n + i]);
if (q_try[n + i] < 0.5f) q_try[n + i] = 0.5f;
else if (q_try[n + i] > 2.0f) q_try[n + i] = 2.0f;
++i;
}
}
float x_try, y_try;
float f_try = cost(q_try, x_try, y_try);

text
				if (f_try <= f_save + c1 * alpha_new * directional_deriv) {
					// Armijo выполнено, вычисляем градиент в пробной точке (аналогично)
					float acc_try = 0.0f;
					float phi_try[32];
					int ii_try = 0;

#pragma loop ivdep
while (ii_try < n) {
acc_try = fmaf(q_try[ii_try], 1.0f, acc_try);
phi_try[ii_try] = acc_try;
++ii_try;
}
float s_try[32], c_try[32];
FABE13_SINCOS(phi_try, s_try, c_try, n);

text
					float as_try = 0.0f, ac_try = 0.0f;
					float sum_s_try[32], sum_c_try[32];
					int k_try = n - 1;
					while (k_try >= 0) {
						const float Lk = cost.variableLen ? q_try[n + k_try] : 1.0f;
						as_try = fmaf(Lk, s_try[k_try], as_try);
						ac_try = fmaf(Lk, c_try[k_try], ac_try);
						sum_s_try[k_try] = as_try;
						sum_c_try[k_try] = ac_try;
						--k_try;
					}

					const float dx_try = fmaf(x_try, 1.0f, -cost.targetX);
					const float dy_try = fmaf(y_try, 1.0f, -cost.targetY);
					const float dist_try = sqrtf(fmaf(dx_try, dx_try, dy_try * dy_try));
					const float inv_dist_try = 1.0f / dist_try;

					float dot = 0.0f;
					i = 0;

#pragma loop ivdep
while (i < n) {
float gpen_try = 0.0f;
{
const float ai = fabsf(q_try[i]);
const float v = fmaf(ai, 1.0f, -cost.minTheta);
if (v > 0.0f) {
const float scale_arg = fmaf(2.0f / cost.minTheta, v * 0.69314718f, 0.0f);
const float exp_val = fmaf(scale_arg,
fmaf(scale_arg,
fmaf(scale_arg,
fmaf(scale_arg, 0.00833333377f, 0.0416666679f),
0.16666667f),
0.5f),
1.0f);
const float dpen = fmaf(cost.sharpW, exp_val * (1.38629436f / cost.minTheta), 0.0f);
gpen_try = fmaf(dpen, copysignf(1.0f, q_try[i]), gpen_try);
}
}
{
const float tsg = fmaf(-q_try[i], cost.archBiasK, 0.0f);
const float exp_arg = -tsg;
const float exp_val = fmaf(exp_arg,
fmaf(exp_arg,
fmaf(exp_arg,
fmaf(exp_arg, 0.00833333377f, 0.0416666679f),
0.16666667f),
0.5f),
1.0f);
const float sig = 1.0f / (exp_val + 1.0f);
gpen_try = fmaf(-cost.archBiasW * cost.archBiasK, sig, gpen_try);
}

text
						const float g_main_try = fmaf(dx_try, -sum_s_try[i], dy_try * sum_c_try[i]) * inv_dist_try;
						float gi_try = g_main_try + gpen_try;
						dot = fmaf(gi_try, direction[i], dot);
						++i;
					}
					if (cost.variableLen) {
						int j = 0;

#pragma loop ivdep
while (j < n) {
const float gi_try = fmaf(dx_try, c_try[j], dy_try * s_try[j]) * inv_dist_try;
dot = fmaf(gi_try, direction[n + j], dot);
++j;
}
}

text
					if (fabsf(dot) <= c2 * fabsf(directional_deriv)) {
						memcpy(q_local, q_try, dim * sizeof(float));
						f = f_try;
						x = x_try;
						y = y_try;
						found = true;
						break;
					}
					else {
						alpha_lo = alpha_new;
						f_lo = f_try;
					}
				}
				else {
					alpha_hi = alpha_new;
				}
				++zoom_iter;
			}

			if (!found) break;

			// Обновление истории L-BFGS после zoom (если используется)
			if (use_lbfgs) {
				// Вычисляем новый градиент в точке q_local
				float new_grad[64];
				float acc2 = 0.0f;
				int i2 = 0;

#pragma loop ivdep
while (i2 < n) {
acc2 = fmaf(q_local[i2], 1.0f, acc2);
phi[i2] = acc2;
++i2;
}
FABE13_SINCOS(phi, s_arr, c_arr, n);
float as2 = 0.0f, ac2 = 0.0f;
int k2 = n - 1;
while (k2 >= 0) {
const float Lk = cost.variableLen ? q_local[n + k2] : 1.0f;
as2 = fmaf(Lk, s_arr[k2], as2);
ac2 = fmaf(Lk, c_arr[k2], ac2);
sum_s[k2] = as2;
sum_c[k2] = ac2;
--k2;
}
const float dx2 = fmaf(x, 1.0f, -cost.targetX);
const float dy2 = fmaf(y, 1.0f, -cost.targetY);
const float dist2 = sqrtf(fmaf(dx2, dx2, dy2 * dy2));
const float inv_dist2 = 1.0f / dist2;
i2 = 0;
#pragma loop ivdep
while (i2 < n) {
float gpen2 = 0.0f;
{
const float ai = fabsf(q_local[i2]);
const float v = fmaf(ai, 1.0f, -cost.minTheta);
if (v > 0.0f) {
const float scale_arg = fmaf(2.0f / cost.minTheta, v * 0.69314718f, 0.0f);
const float exp_val = fmaf(scale_arg,
fmaf(scale_arg,
fmaf(scale_arg,
fmaf(scale_arg, 0.00833333377f, 0.0416666679f),
0.16666667f),
0.5f),
1.0f);
const float dpen = fmaf(cost.sharpW, exp_val * (1.38629436f / cost.minTheta), 0.0f);
gpen2 = fmaf(dpen, copysignf(1.0f, q_local[i2]), gpen2);
}
}
{
const float tsg = fmaf(-q_local[i2], cost.archBiasK, 0.0f);
const float exp_arg = -tsg;
const float exp_val = fmaf(exp_arg,
fmaf(exp_arg,
fmaf(exp_arg,
fmaf(exp_arg, 0.00833333377f, 0.0416666679f),
0.16666667f),
0.5f),
1.0f);
const float sig = 1.0f / (exp_val + 1.0f);
gpen2 = fmaf(-cost.archBiasW * cost.archBiasK, sig, gpen2);
}
const float g_main2 = fmaf(dx2, -sum_s[i2], dy2 * sum_c[i2]) * inv_dist2;
new_grad[i2] = g_main2 + gpen2;
++i2;
}
if (cost.variableLen) {
int j2 = 0;
#pragma loop ivdep
while (j2 < n) {
new_grad[n + j2] = fmaf(dx2, c_arr[j2], dy2 * s_arr[j2]) * inv_dist2;
++j2;
}
}

text
				float s_new[32], y_new[32];
				i = 0;

#pragma loop ivdep
while (i < dim) {
s_new[i] = q_local[i] - q_local_save[i];
y_new[i] = new_grad[i] - grad[i];
++i;
}

text
				float ys = 0.0f;
				i = 0;

#pragma loop ivdep
while (i < dim) {
ys = fmaf(y_new[i], s_new[i], ys);
++i;
}

text
				const float eps_lbfgs = 1e-8f;
				if (ys > eps_lbfgs) {
					int new_idx;
					if (lbfgs_iter < m_lbfgs) {
						new_idx = lbfgs_iter;
						lbfgs_iter++;
					}
					else {
						for (int idx = 0; idx < m_lbfgs - 1; ++idx) {
							for (int d = 0; d < dim; ++d) {
								s_storage[d][idx] = s_storage[d][idx + 1];
								y_storage[d][idx] = y_storage[d][idx + 1];
							}
							rho_storage[idx] = rho_storage[idx + 1];
						}
						new_idx = m_lbfgs - 1;
					}
					for (int d = 0; d < dim; ++d) {
						s_storage[d][new_idx] = s_new[d];
						y_storage[d][new_idx] = y_new[d];
					}
					rho_storage[new_idx] = 1.0f / ys;
				}
				else {
					lbfgs_iter = 0;
				}
			}

			++outer;
		}

		if (f < bestF) {
			bestF = f;
			bestQ.assign(q_local, q_local + dim);
			bestX = x;
			bestY = y;
			no_improve = 0;
		}
		else {
			++no_improve;
		}
	}
	return f;
	};", возьми L-BFGS механизм из этой доводки и примени его  к текущей, добавь в простую доводку L-BFGS для очень перспективных шагов градиентного спуска, то есть логика должна быть следующая, после каждого шага градиентного спуска проверяется относительное улучшение относительно точки изнаально поданной на градиентный спуск, и если процентное улучшение велико например больше 70 процентов - то запускать сразу же из этой точки L-BFGS, если же во время работы L-BFGS наблюдаться будет расходимость или увеличение функции то также предусмотреть механизм возврата из L-BFGS к градиентному спуску, приведи полный код реализации без пропусков, обрати внимание что код градиентного спуска точно правильный но код L-BFGS может быть неверным так как все эксперименты с L-BFGS дают худшие результаты - так что возможно нужно корректно переписать L-BFGS, но идея должна остаться той же - градентный спуск в начале, затем из него L-BFGS для очень удачных точек, с механизмом возврата к градиентному спуску в случае расходимости или увеличения функции, условие попадание в L-BFGS тоже нужно оставить таким же - по 70ти процентному относительному улучшению

사고

Integrating L-BFGS into GD refinement

The user wants to integrate L-BFGS into a simple GD local refinement process controlled by a relative improvement measure. I need to apply L-BFGS when the improvement surpasses a certain threshold, ensuring fallback in case of divergence. The task requires maintaining precise gradient computations, respecting array sizes and other specifics. I'll implement the L-BFGS with a two-loop recursion, fitting it into the existing code. I also need to correct any potential issues in the old L-BFGS code.

Implementing L-BFGS with GD refinement

The user wants to introduce L-BFGS after a promising step in the GD process. After each accepted GD step, I’ll compute the relative improvement, and if it exceeds 70%, I’ll switch to L-BFGS. I'll set up history arrays and establish conditions to handle failure, divergence, and fallback to the last GD point if needed.

Handling divergence and reversion to GD from L-BFGS

The user wants to revert to GD from L-BFGS if divergence occurs. I'll implement a mechanism to save a snapshot before L-BFGS runs and restore it when needed. If L-BFGS shows divergence or failure, we'll return to GD using the snapshot, possibly lowering the step size. I'll also create a gradient function to compute gradients without duplication and ensure they’re consistent with the current state.

Handling gradient dependencies and potential issues with arrays

The gradient computation relies on x and y values returned by the cost function. I'll make sure the gradient function uses the correct x, y, and q values while calculating related arrays. There's a risk in the old code where array indexing for s_storage was incorrect. If dim exceeds 32 due to variableLen, this could break the code. I'll ensure the arrays are sized properly to handle this.

Implementing circular buffer and L-BFGS updates

We’ll allocate s_hist and y_hist as circular buffers, and implement the two-loop recursion with backtracking Armijo. Using a circular buffer will help avoid unnecessary array shifting, though we may shift when m = 5. For line search, we’ll use the Armijo condition and check directional derivatives with clamping. Even with potential inaccuracies from clamping, I'll apply practices similar to the original GD since their code works. We'll keep the approach minimal and readable.

Computing updates and managing L-BFGS triggers

After accepting a step, we need to update the variables and compute the gradient at the new point. If the relative improvement is significant (≥70%) and conditions are met, we can trigger L-BFGS once per evaluation. We'll check if L-BFGS diverged and decide to either revert to GD or continue with L-BFGS, depending on its success. We'll remove redundant checks and manage this process efficiently, using a cooldown approach.

Proposing strategy for integrating GD and L-BFGS

I propose starting with gradient descent (GD), and if there’s substantial improvement, we'll run L-BFGS for up to the max iterations. If L-BFGS succeeds, we'll end GD, but if it fails, we'll revert to GD. The condition for L-BFGS to trigger is based on relative improvement. We'll also track the best solution and adjust the step size during the process, using backtracking to refine performance.

Handling LBFGS divergence and full code structure

For LBFGS, I’ll detect divergence if the cost becomes NaN/Inf, if line search fails, or if step size reduces too low. If the step leads to an increase in function value, I'll restore a snapshot and revert to GD. I'll respect constraints and preserve the original gradient computation, even if it means duplicating some code. I'll ensure the full code is compiled and ready, including needed arrays and helper functions.

Tracking previous q values for s calculation

To calculate s, I need a separate q_prev array, like q_prev[64]. Their code uses q_local and q_try arrays preallocated, possibly to size 64. I’ll allocate the necessary local arrays inside the evalAt function, making sure to handle the dimension correctly with memcpy.

Crafting evaluation code and gradient descent

I’m constructing code for the evaluation function, defining parameters for line search, LBFGS, and gradient computation. The history arrays will store past iterations, and I’ll implement line search with a backtracking mechanism. The goal is to trigger LBFGS after a specific improvement, and control the iteration loop, adjusting parameters like eta and alpha.

Implementing LBFGS with safeguards and gradient updates

I’m implementing the LBFGS process, focusing on direction computation, gradient checks, and step acceptance. I’ll ensure that the computeGrad function respects the existing arrays and includes necessary safeguards against potential inf or NaN values. I'll also be careful to replicate certain behavior such as updates to alpha_j and indexing adjustments.

Thought for 2m 32s

답변

Ниже — полная версия твоей текущей простой локальной доводки (GD + Armijo backtracking), в которую добавлен корректный L-BFGS (двухпетлевая рекурсия) и логика:

  • После каждого принятого шага GD считаем относительное улучшение относительно f_start (значение функции в точке, с которой запустили локальную доводку).
  • Если улучшение >= 70% → сразу переключаемся в L-BFGS из текущей точки.
  • В L-BFGS используется Armijo backtracking по общему правилу f(x+αd) <= f(x) + c1*α*(gᵀd).
  • Если в L-BFGS:
    • не найден шаг (backtracking провалился), или
    • f стала хуже, или
    • появились NaN/Inf,
    то срабатывает fallback: возвращаемся к градиентному спуску (причём сохраняем лучшую точку, найденную в L-BFGS, если она лучше, чем точка входа в L-BFGS).

Важное: я не трогал твою формулу/код вычисления градиента и текущий GD-механизм (он перенесён 1-в-1), а L-BFGS переписан в более “каноничном” и устойчивом виде (история s,y, двухпетлевая рекурсия, Armijo по gᵀd, проверка убывания/валидности).

cpp
auto evalAt = [&](const float t) -> float {
    map.map01ToPoint(t, q_local);
    float f = cost(q_local, x, y);

    if (f < fmaf(bestF, adaptive_coeff, 0.0f)) {

        // --- значение функции в точке старта локальной доводки (для триггера 70%) ---
        const float f_start = f;

        // Параметры Armijo
        const float c1 = 1e-4f;                // достаточное убывание
        const float tau = 0.5f;                // коэффициент уменьшения шага
        const int max_outer_iters = (int)(50.0f * (1.0f + 0.65f * p));
        const int max_backtrack  = (int)(20.0f * (1.0f + 0.65f * p));

        // --- L-BFGS параметры ---
        const float lbfgs_trigger = 0.70f;     // 70% относительное улучшение
        const int   m_lbfgs = 5;               // память
        const int   max_lbfgs_iters = (int)(25.0f * (1.0f + 0.65f * p));
        const float eps_lbfgs_curv = 1e-8f;    // порог кривизны ys
        const float eps_descent = 1e-12f;

        float eta = 1.0f;                      // шаг для GD (и как стартовый для line-search)

        // ------------------------------------------------------------
        // Вспомогательное: clamp в допустимые границы (как у тебя)
        // ------------------------------------------------------------
        auto clampPoint = [&](float* q) {
            int i = 0;
#pragma loop ivdep
            while (i < n) {
                const float lo = (i == 0) ? -1.0471975511965977f : -2.6179938779914944f;
                const float hi =  2.6179938779914944f;
                if (q[i] < lo) q[i] = lo;
                else if (q[i] > hi) q[i] = hi;
                ++i;
            }
            if (cost.variableLen) {
                i = 0;
#pragma loop ivdep
                while (i < n) {
                    if (q[n + i] < 0.5f) q[n + i] = 0.5f;
                    else if (q[n + i] > 2.0f) q[n + i] = 2.0f;
                    ++i;
                }
            }
        };

        // ------------------------------------------------------------
        // Вспомогательное: вычислить градиент (код 1-в-1 как у тебя)
        // Вход: q_in, x_in, y_in
        // Выход: grad_out[dim], grad_norm2_out
        // ------------------------------------------------------------
        auto computeGrad = [&](const float* q_in, const float x_in, const float y_in,
                               float* grad_out, float& grad_norm2_out) {
            // ---- Вычисляем градиент в текущей точке q_in ----
            float acc = 0.0f;
            int ii = 0;
#pragma loop ivdep
            while (ii < n) {
                acc = fmaf(q_in[ii], 1.0f, acc);
                phi[ii] = acc;
                ++ii;
            }
            FABE13_SINCOS(phi, s_arr, c_arr, n);

            float as = 0.0f;
            float ac = 0.0f;
            int k = n - 1;
            while (k >= 0) {
                const float Lk = cost.variableLen ? q_in[n + k] : 1.0f;
                as = fmaf(Lk, s_arr[k], as);
                ac = fmaf(Lk, c_arr[k], ac);
                sum_s[k] = as;
                sum_c[k] = ac;
                --k;
            }

            const float dx = fmaf(x_in, 1.0f, -cost.targetX);
            const float dy = fmaf(y_in, 1.0f, -cost.targetY);
            const float dist = sqrtf(fmaf(dx, dx, dy * dy));
            const float inv_dist = 1.0f / dist;

            grad_norm2_out = 0.0f;

            // Градиент по углам
            int i = 0;
#pragma loop ivdep
            while (i < n) {
                float gpen = 0.0f;
                // Производная штрафа за превышение minTheta
                {
                    const float ai = fabsf(q_in[i]);
                    const float v = fmaf(ai, 1.0f, -cost.minTheta);
                    if (v > 0.0f) {
                        const float scale_arg = fmaf(2.0f / cost.minTheta, v * 0.69314718f, 0.0f);
                        const float exp_val = fmaf(scale_arg,
                            fmaf(scale_arg,
                                fmaf(scale_arg,
                                    fmaf(scale_arg, 0.00833333377f, 0.0416666679f),
                                        0.16666667f),
                                    0.5f),
                                1.0f);
                        const float dpen = fmaf(cost.sharpW, exp_val * (1.38629436f / cost.minTheta), 0.0f);
                        gpen = fmaf(dpen, copysignf(1.0f, q_in[i]), gpen);
                    }
                }
                // Производная арк-штрафа
                {
                    const float tsg = fmaf(-q_in[i], cost.archBiasK, 0.0f);
                    const float exp_arg = -tsg;
                    const float exp_val = fmaf(exp_arg,
                        fmaf(exp_arg,
                            fmaf(exp_arg,
                                fmaf(exp_arg, 0.00833333377f, 0.0416666679f),
                                    0.16666667f),
                                0.5f),
                            1.0f);
                    const float sig = 1.0f / (exp_val + 1.0f);
                    gpen = fmaf(-cost.archBiasW * cost.archBiasK, sig, gpen);
                }

                const float g_main = fmaf(dx, -sum_s[i], dy * sum_c[i]) * inv_dist;
                float gi = g_main + gpen;
                grad_out[i] = gi;
                grad_norm2_out = fmaf(gi, gi, grad_norm2_out);
                ++i;
            }

            // Градиент по длинам
            if (cost.variableLen) {
                int j = 0;
#pragma loop ivdep
                while (j < n) {
                    const float gi = fmaf(dx, c_arr[j], dy * s_arr[j]) * inv_dist;
                    grad_out[n + j] = gi;
                    grad_norm2_out = fmaf(gi, gi, grad_norm2_out);
                    ++j;
                }
            }
        };

        // ------------------------------------------------------------
        // Вспомогательное: общий Armijo backtracking для направления dir
        // Условие: f_new <= f_cur + c1 * alpha * (g^T d)
        // ------------------------------------------------------------
        auto armijoLineSearch = [&](const float* q_base,
                                    const float f_base,
                                    const float x_base, const float y_base,
                                    const float* grad_base,
                                    const float* dir,
                                    const float gtd,
                                    float& alpha_io,
                                    float* q_out,
                                    float& f_out,
                                    float& x_out,
                                    float& y_out) -> bool {
            float alpha = alpha_io;
            int backtrack = 0;

            while (backtrack < max_backtrack) {
                int i = 0;
#pragma loop ivdep
                while (i < dim) {
                    q_out[i] = fmaf(alpha, dir[i], q_base[i]);
                    ++i;
                }

                // clamp как у тебя
                clampPoint(q_out);

                float x2, y2;
                float f_try = cost(q_out, x2, y2);

                // проверки на NaN/Inf
                if (!(f_try == f_try)) { // NaN
                    alpha *= tau;
                    ++backtrack;
                    continue;
                }

                // Armijo
                if (f_try <= f_base + c1 * alpha * gtd) {
                    alpha_io = alpha;
                    f_out = f_try;
                    x_out = x2;
                    y_out = y2;
                    return true;
                }

                alpha *= tau;
                ++backtrack;
            }

            return false;
        };

        // ------------------------------------------------------------
        // 1) Градиентный спуск (твой текущий, без изменений)
        // 2) При улучшении >= 70% -> L-BFGS, с fallback обратно в GD
        // ------------------------------------------------------------
        bool lbfgs_already_tried = false;

        int outer = 0;
        while (outer < max_outer_iters) {

            // ---- Градиент в текущей точке ----
            float grad[64];
            float grad_norm2 = 0.0f;
            computeGrad(q_local, x, y, grad, grad_norm2);

            if (!(grad_norm2 == grad_norm2) || grad_norm2 < 1e-12f) break;

            // ---- Направление GD: d = -g ----
            float dir_gd[64];
            int i = 0;
#pragma loop ivdep
            while (i < dim) {
                dir_gd[i] = -grad[i];
                ++i;
            }

            // Для d = -g: g^T d = -||g||^2
            const float gtd_gd = -grad_norm2;

            // ---- Line-search (Armijo) ----
            float eta_trial = eta;
            float f_new = f, x_new = x, y_new = y;

            // пробная точка в q_try
            bool found = armijoLineSearch(q_local, f, x, y, grad, dir_gd, gtd_gd,
                                          eta_trial, q_try, f_new, x_new, y_new);

            if (!found) break; // не удалось сделать шаг — выход

            // Шаг принят
            memcpy(q_local, q_try, static_cast<size_t>(dim) * sizeof(float));
            f = f_new;
            x = x_new;
            y = y_new;
            eta = eta_trial; // запоминаем удачный шаг

            // --------------------------------------------------------
            // Триггер: относительное улучшение от f_start
            // improvement = (f_start - f) / f_start
            // --------------------------------------------------------
            const float denom = fmaxf(f_start, 1e-12f);
            const float rel_impr = (f_start - f) / denom;

            if (!lbfgs_already_tried && rel_impr >= lbfgs_trigger) {
                lbfgs_already_tried = true;

                // --- Сохраняем точку, чтобы уметь вернуться в GD ---
                float q_resume[64];
                memcpy(q_resume, q_local, static_cast<size_t>(dim) * sizeof(float));
                float f_resume = f;
                float x_resume = x;
                float y_resume = y;
                float eta_resume = eta;

                // --- Лучшая точка, найденная внутри L-BFGS ---
                float q_best_lbfgs[64];
                memcpy(q_best_lbfgs, q_local, static_cast<size_t>(dim) * sizeof(float));
                float f_best_lbfgs = f;
                float x_best_lbfgs = x;
                float y_best_lbfgs = y;

                // --- История L-BFGS ---
                float s_hist[m_lbfgs][64];
                float y_hist[m_lbfgs][64];
                float rho_hist[m_lbfgs];
                float alpha_hist[m_lbfgs];
                int hist_size = 0;

                // --- Текущий градиент для L-BFGS ---
                float gk[64];
                float gk_norm2 = 0.0f;
                computeGrad(q_local, x, y, gk, gk_norm2);

                bool lbfgs_ok = true;

                // стартовый alpha для L-BFGS (можно взять eta, но чаще лучше 1.0)
                float alpha_k = 1.0f;

                int it = 0;
                while (it < max_lbfgs_iters) {
                    if (!(gk_norm2 == gk_norm2) || gk_norm2 < 1e-12f) break;

                    // ---- Считаем направление L-BFGS ----
                    float dir[64];

                    if (hist_size == 0) {
                        int d = 0;
#pragma loop ivdep
                        while (d < dim) {
                            dir[d] = -gk[d];
                            ++d;
                        }
                    } else {
                        // two-loop recursion
                        float q_vec[64];
                        int d = 0;
#pragma loop ivdep
                        while (d < dim) {
                            q_vec[d] = gk[d];
                            ++d;
                        }

                        for (int jj = hist_size - 1; jj >= 0; --jj) {
                            float dot_sq = 0.0f;
                            d = 0;
#pragma loop ivdep
                            while (d < dim) {
                                dot_sq = fmaf(s_hist[jj][d], q_vec[d], dot_sq);
                                ++d;
                            }
                            const float a = dot_sq * rho_hist[jj];
                            alpha_hist[jj] = a;

                            d = 0;
#pragma loop ivdep
                            while (d < dim) {
                                q_vec[d] = fmaf(-a, y_hist[jj][d], q_vec[d]);
                                ++d;
                            }
                        }

                        // gamma = (s_{k-1}^T y_{k-1}) / (y_{k-1}^T y_{k-1})
                        float gamma = 1.0f;
                        {
                            const int last = hist_size - 1;
                            float yy = 0.0f;
                            d = 0;
#pragma loop ivdep
                            while (d < dim) {
                                yy = fmaf(y_hist[last][d], y_hist[last][d], yy);
                                ++d;
                            }
                            // ys = 1/rho
                            const float ys = 1.0f / rho_hist[last];
                            if (yy > 0.0f) gamma = ys / yy;
                        }

                        float r_vec[64];
                        d = 0;
#pragma loop ivdep
                        while (d < dim) {
                            r_vec[d] = gamma * q_vec[d];
                            ++d;
                        }

                        for (int jj = 0; jj < hist_size; ++jj) {
                            float dot_yr = 0.0f;
                            d = 0;
#pragma loop ivdep
                            while (d < dim) {
                                dot_yr = fmaf(y_hist[jj][d], r_vec[d], dot_yr);
                                ++d;
                            }
                            const float b = dot_yr * rho_hist[jj];
                            const float coeff = alpha_hist[jj] - b;

                            d = 0;
#pragma loop ivdep
                            while (d < dim) {
                                r_vec[d] = fmaf(coeff, s_hist[jj][d], r_vec[d]);
                                ++d;
                            }
                        }

                        d = 0;
#pragma loop ivdep
                        while (d < dim) {
                            dir[d] = -r_vec[d];
                            ++d;
                        }
                    }

                    // ---- Проверка, что направление действительно спусковое ----
                    float gtd = 0.0f;
                    int d = 0;
#pragma loop ivdep
                    while (d < dim) {
                        gtd = fmaf(gk[d], dir[d], gtd);
                        ++d;
                    }
                    if (!(gtd == gtd) || gtd >= -eps_descent) {
                        // fallback на -grad
                        d = 0;
#pragma loop ivdep
                        while (d < dim) {
                            dir[d] = -gk[d];
                            ++d;
                        }
                        gtd = -gk_norm2;
                        if (gtd >= -eps_descent) { lbfgs_ok = false; break; }
                    }

                    // ---- Сохраняем старую точку/градиент для s,y ----
                    float q_old[64];
                    memcpy(q_old, q_local, static_cast<size_t>(dim) * sizeof(float));
                    float g_old[64];
                    memcpy(g_old, gk, static_cast<size_t>(dim) * sizeof(float));
                    const float f_old = f;
                    const float x_old = x;
                    const float y_old = y;

                    // ---- Line-search Armijo по L-BFGS направлению ----
                    float alpha_try = alpha_k;
                    float f_try, x_try, y_try;

                    bool step_ok = armijoLineSearch(q_local, f, x, y, gk, dir, gtd,
                                                   alpha_try, q_try, f_try, x_try, y_try);

                    if (!step_ok) { lbfgs_ok = false; break; }

                    // ---- Принимаем шаг ----
                    memcpy(q_local, q_try, static_cast<size_t>(dim) * sizeof(float));
                    f = f_try;
                    x = x_try;
                    y = y_try;
                    alpha_k = alpha_try;

                    // защитные проверки (на всякий случай)
                    if (!(f == f) || f > f_old + 1e-6f) { lbfgs_ok = false; break; }

                    // ---- Новый градиент ----
                    computeGrad(q_local, x, y, gk, gk_norm2);
                    if (!(gk_norm2 == gk_norm2)) { lbfgs_ok = false; break; }

                    // ---- Обновляем лучшую точку внутри L-BFGS ----
                    if (f < f_best_lbfgs) {
                        f_best_lbfgs = f;
                        memcpy(q_best_lbfgs, q_local, static_cast<size_t>(dim) * sizeof(float));
                        x_best_lbfgs = x;
                        y_best_lbfgs = y;
                    }

                    // ---- s = x_{k+1}-x_k, y = g_{k+1}-g_k ----
                    float ys = 0.0f;
                    float s_new[64];
                    float y_new[64];

                    d = 0;
#pragma loop ivdep
                    while (d < dim) {
                        const float sd = q_local[d] - q_old[d];
                        const float yd = gk[d] - g_old[d];
                        s_new[d] = sd;
                        y_new[d] = yd;
                        ys = fmaf(yd, sd, ys);
                        ++d;
                    }

                    if (!(ys == ys) || ys <= eps_lbfgs_curv) {
                        // плохая кривизна -> сброс истории
                        hist_size = 0;
                    } else {
                        // пушим в историю (сдвигом; m маленькое)
                        if (hist_size < m_lbfgs) {
                            const int idx = hist_size;
                            rho_hist[idx] = 1.0f / ys;
                            d = 0;
#pragma loop ivdep
                            while (d < dim) {
                                s_hist[idx][d] = s_new[d];
                                y_hist[idx][d] = y_new[d];
                                ++d;
                            }
                            ++hist_size;
                        } else {
                            // shift 0..m-2 <- 1..m-1
                            for (int jj = 0; jj < m_lbfgs - 1; ++jj) {
                                rho_hist[jj] = rho_hist[jj + 1];
                                d = 0;
#pragma loop ivdep
                                while (d < dim) {
                                    s_hist[jj][d] = s_hist[jj + 1][d];
                                    y_hist[jj][d] = y_hist[jj + 1][d];
                                    ++d;
                                }
                            }
                            const int idx = m_lbfgs - 1;
                            rho_hist[idx] = 1.0f / ys;
                            d = 0;
#pragma loop ivdep
                            while (d < dim) {
                                s_hist[idx][d] = s_new[d];
                                y_hist[idx][d] = y_new[d];
                                ++d;
                            }
                            hist_size = m_lbfgs;
                        }
                    }

                    // ---- критерий выхода ----
                    if (gk_norm2 < 1e-12f) break;

                    ++it;
                }

                if (lbfgs_ok) {
                    // L-BFGS успешно: фиксируем лучшую точку из L-BFGS и завершаем доводку
                    memcpy(q_local, q_best_lbfgs, static_cast<size_t>(dim) * sizeof(float));
                    f = f_best_lbfgs;
                    x = x_best_lbfgs;
                    y = y_best_lbfgs;
                    break; // выходим из GD-цикла: L-BFGS — финальная стадия
                } else {
                    // L-BFGS провалился: откатываемся (но если L-BFGS успел найти лучше — берем лучше)
                    if (f_best_lbfgs < f_resume) {
                        memcpy(q_local, q_best_lbfgs, static_cast<size_t>(dim) * sizeof(float));
                        f = f_best_lbfgs;
                        x = x_best_lbfgs;
                        y = y_best_lbfgs;
                    } else {
                        memcpy(q_local, q_resume, static_cast<size_t>(dim) * sizeof(float));
                        f = f_resume;
                        x = x_resume;
                        y = y_resume;
                    }

                    // возвращаем шаг GD и продолжаем GD
                    eta = eta_resume;
                }
            }

            ++outer;
        }

        // Обновляем лучшую точку
        if (f < bestF) {
            bestF = f;
            bestQ.assign(q_local, q_local + dim);
            bestX = x;
            bestY = y;
            no_improve = 0;
        } else {
            ++no_improve;
        }
    }

    return f;
};

Если захочешь усилить “не мешает”/безопасность ещё сильнее (без ухудшения твоего GD), самый практичный тюнинг обычно такой (я не внедрял, чтобы не менять задумку):

  • ограничить L-BFGS по времени/итерациям ещё жёстче (например 10–15),
  • запускать L-BFGS только если не в клипе (много клипа ломает квазиньютоновские аппроксимации),
  • и/или сбрасывать историю L-BFGS при частом клипе или при слишком малом фактическом шаге.

Но текущая версия уже делает главное: L-BFGS включается только на “суперудачных” точках и не может сломать результат — есть откат обратно в GD.

이 Q&A 공유하기