Реализация алгоритма Левенберга-Марквардта в MATLAB, код проблемы не сходится.

Я пытаюсь реализовать алгоритм Левенберга-Маркуарда в MATLAB, но мой код не сходим и ошибки становятся всё больше и больше после первой итерации. Я проверил матрицу Якоби, обратную функцию и так далее, но не смог понять, в чем проблема.

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

Код находит низкую остаточную ошибку в начале, но после вычисления матрицы Якоби и нахождения значения delta_value и изменения весов на второй или третьей итерации норма остаточных значений увеличивается.

Есть ли что-то, что я не вижу? Спасибо за вашу помощь, кстати.

clear; clc;
% Определить гиперпараметры
epoch = 10;            % Максимальное количество эпох
tol = 0.02;             % Допуск для сходимости
lambda = 0.01;          % Начальный параметр регуляризации
layer_number = 3;       % количество слоев (входной и выходной включены)
input_size= 1;          % Количество входных признаков
hidden_size = 15;       % Количество нейронов в скрытом слое
output_size = 1;        % Количество выходов
num_runs = 5;           % Количество запусков для расчета среднего
number_of_sample=46;    % Количество выборок

% Создание обучающих данных для задачи один Синусоидальная волна y = 1/2 + 1/4sin(37pix)
input_vals = linspace(-1, 1, number_of_sample)';        % 40 входов, распределенных на интервале [-1, 1]
target_vals = 0.5 + 0.25 * sin(3 * pi * input_vals);    % Целевое значение на основе данной формулы

% Запуск алгоритма 5 раз для расчета средней ошибки для построения графика
errors_of_epochs_for_average_5_run = zeros(epoch,1); 
for run = 1:num_runs

    W1 = rand(hidden_size,input_size);
    b1 = rand(hidden_size,1);
    W2 = rand(output_size,hidden_size);
    b2 = rand(output_size,1);

    parameterVector= [W1(:);b1(:);W2(:);b2];                       % вектор параметров для матрицы Якоби;

    errors_of_epochs = zeros(epoch,1);                          % для хранения значений ошибок каждой эпохи
    for epoch_number = 1:epoch
        % Инициализация матрицы для хранения векторов ошибок для каждого запуска
        output_vals=zeros(number_of_sample,output_size);        % выходы для всех выборок
        errors_for_each=zeros(number_of_sample,output_size);    % значения ошибок для каждой входной-выходной пары

        % цикл прямого хода
        for i=1:number_of_sample
            output_vals(i,output_size) = forward_pass(input_vals(i,input_size), W1, b1, W2, b2);
        end
       errors_for_each = target_vals - output_vals;
       error_norm = norm(errors_for_each);
       errors_of_epochs(epoch_number,1) = error_norm;        % расчет средней суммы квадратов для каждой эпохи

       Jaco=jacobian(parameterVector, target_vals, input_vals,hidden_size); % расчет матрицы Якоби
       while true               % для повторного расчета delta, когда ошибка больше предыдущей

           JTJ=Jaco'*Jaco;

           I=eye(size(JTJ));
           JTJ_Lambda_Eye=JTJ+lambda*I;
           JTe=Jaco'*errors_for_each;
           try
                [R,Flag]=chol(JTJ_Lambda_Eye);
                if Flag
                    disp('Матрица должна быть симметричной и положительно определенной.');

                end

                % JTJLEyeInverse=inv(JTJ_Lambda_Eye);

                L = chol(JTJ_Lambda_Eye, 'lower');
                n=size(JTJ_Lambda_Eye,1);
                I=eye(n);
                Y=zeros(n);
                % Решение для Y с использованием прямой подстановки
                for i = 1:n
                    Y(:, i) = L \ I(:, i);
                end
                % Теперь решите L^T * X = Y для X
                JTJLEyeInverse = L' \ Y;
           catch
                disp('Разложение Холецкого не удалось.');
                break;
           end
           deltaValue=JTJLEyeInverse*JTe;
           parameterVectorTest=parameterVector;
           parameterVectorTest=parameterVectorTest-deltaValue;
           % проверить, создаст ли новая ошибка лучшее значение.
           W1_test = reshape(parameterVectorTest(1:hidden_size * input_size), hidden_size, input_size);
           b1_test = reshape(parameterVectorTest(hidden_size * input_size + 1:hidden_size * input_size + hidden_size),hidden_size, input_size);
           W2_test = reshape(parameterVectorTest(hidden_size * input_size + hidden_size + 1:end - 1), output_size, hidden_size);
           b2_test = parameterVectorTest(end);

           % тестировать изменение delta, слишком ли они малы для изменения
           if norm(deltaValue) errors_of_epochs(epoch_number,1)
               lambda=lambda*10;
               % если lambda>1000
               %     lambda=1000;
               % end
           else
               break;
           end
       end % заканчивается точка больше ошибки и увеличение lambda

       if test_error_norm <= errors_of_epochs(epoch_number,1)
           W1=W1_test;
           b1=b1_test;
           W2=W2_test;
           b2=b2_test;
           lambda=max(lambda/10, 1e-7);
       end

       if test_error_norm <= tol
           disp(['Сошлось на эпохе ', num2str(epoch_number), ' в запуске ', num2str(run)]);
           break;
       end

    end

    errors_of_epochs_for_average_5_run = errors_of_epochs_for_average_5_run + errors_of_epochs; % для расчета среднего из 5 запусков

end

% Построение результата для исследования
% Расчет средней ошибки по всем запускам
avg_epoch_errors = mean(errors_of_epochs_for_average_5_run, 1);

% Построение графика средней ошибки против эпох
plot_target=zeros(number_of_sample ,1);
 for i_test = 1:input_size
    X_test = input_vals(i_test);
    % Прямой проход для получения текущего выходного значения
    plot_target(i_test) = forward_pass(X_test, W1, b1, W2, b2);
    % Расчет текущей ошибки  
end
figure;
hold on;
%plot(1:biggest_last_epoch_number, avg_epoch_errors(1:biggest_last_epoch_number), '-o');
plot(input_vals, target_vals, 'r-o');
plot(input_vals, plot_target, 'b-o');
xlabel('Эпоха');
ylabel('Средняя ошибка');
title('Средняя ошибка за эпоху за 5 запусков');
grid on;
hold off;

% Определите функцию прямого прохода
function [outputf,A1] = forward_pass(X, W1, b1, W2, b2)
    Z1 = (W1 * X) + b1;                 % Линейная комбинация для скрытого слоя
    A1 = sigmoid(Z1);                  % Сигмоидальная активация для скрытого слоя
    Z2 = (W2 * A1) + b2;                % Линейная комбинация для выходного слоя
    outputf = Z2;                      % Выходной слой
end
% Сигмоидальная активационная функция
function cikti = sigmoid(x)
    cikti = 1 ./ (1 + exp(-x));
end

% Производная сигмоидальной функции
function cikti = sigmoidDerivative(x)
    cikti = sigmoid(x) .* (1 - sigmoid(x));
end

% Определите функцию для вычисления Якобиана
function jaco_return=jacobian(jparameterVector, jtarget_vals, jinput_vals,jhidden_size)
    jaco_return=zeros(length(jtarget_vals), length(jparameterVector));
    W1 = reshape(jparameterVector(1:jhidden_size * 1), jhidden_size, 1);
    b1 = reshape(jparameterVector(jhidden_size * 1 + 1:jhidden_size * 1 + jhidden_size),jhidden_size, 1);
    W2 = reshape(jparameterVector(jhidden_size * 1 + jhidden_size + 1:end - 1), 1, jhidden_size);
    b2 = jparameterVector(end);
    for i=1:length(jtarget_vals)            % создание каждой строки Якобиана для каждой входной-выходной пары
        [~ , jA1] = forward_pass(jinput_vals(i),W1,b1,W2,b2);
        derivW1 = -1 * W2'.* sigmoidDerivative(jA1) *jinput_vals(i);    % частные производные W1
        derivb1 = -1 * W2'.* sigmoidDerivative(jA1);                    % частные производные b1
        derivW2 = -1 * jA1;                                             % частные производные W2
        derivb2 = -1;                                                   % частные производные b2 линейная активационная функция
        jaco_return(i,:)= [derivW1(:);derivb1(:);derivW2(:);derivb2]; % каждая строка матрицы Якоби      
    end
    epsilon=1e-10;
    jaco_test=zeros(length(jtarget_vals), length(jparameterVector));
    jparameterVector_pus=jparameterVector;

    for jsatir=1:length(jtarget_vals)
        for jsutun=1:length(jparameterVector)
            outpus_test_normal=forward_pass(jinput_vals(jsatir),W1,b1,W2,b2);

            jparameterVector_pus(jsutun)=jparameterVector(jsutun)+epsilon;

            W1_test = reshape(jparameterVector_pus(1:jhidden_size * 1), jhidden_size, 1);
            b1_test = reshape(jparameterVector_pus(jhidden_size * 1 + 1:jhidden_size * 1 + jhidden_size),jhidden_size, 1);
            W2_test = reshape(jparameterVector_pus(jhidden_size * 1 + jhidden_size + 1:end - 1), 1, jhidden_size);
            b2_test = jparameterVector_pus(end);
            jparameterVector_pus(jsutun)=jparameterVector(jsutun);
            output_vals_test_plus = forward_pass(jinput_vals(jsatir),W1_test,b1_test,W2_test,b2_test);
            jaco_test(jsatir, jsutun)=((output_vals_test_plus- outpus_test_normal)/epsilon)*-1;
        end
    end
    jaco_return=jaco_test;
end