Matlab实现Levenberg-Marquardt算法不收敛问题求助

问题描述 #

我在Matlab中实现Levenberg-Marquardt算法训练神经网络时遇到不收敛问题：首次迭代后误差持续增大。已检查Jacobian矩阵（尝试过数值微分和解析公式）、逆函数模块（测试过内置逆函数与Cholesky分解），还通过Excel验证了矩阵和偏置计算值。算法初始残差较低，但计算Jacobian、求解delta并更新权重后，残差范数从第二次或第三次迭代开始持续增大，求排查可能的问题。

实现代码 #

clear; clc;
% Define hyperparameters
epoch = 10;            % Maximum number of epochs
tol = 0.02;             % Tolerance for convergence
lambda = 0.01;          % Initial regularization parameter
layer_number = 3;       % number of layer(input and output are included)
input_size= 1;          % Number of input features
hidden_size = 15;       % Number of neurons in the hidden layer
output_size = 1;        % Number of outputs
num_runs = 5;           % Number of runs to calculate average
number_of_sample=46;    % Number of sample

% Create training data for problem one Sinus wave y = 1/2 + 1/4sin(37pix)
input_vals = linspace(-1, 1, number_of_sample)';        % 40 inputs scattered in the interval [-1, 1]
target_vals = 0.5 + 0.25 * sin(3 * pi * input_vals);    % Target output according to given formula

% Run the algorithm 5 times for calculating average error to create plot
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];                       %parameter vector for jacobian matrix;
    

    errors_of_epochs = zeros(epoch,1);                          % to store error values of each epoch
    for epoch_number = 1:epoch
        % Initialize matrix to store error vectors for each run
        output_vals=zeros(number_of_sample,output_size);        %outputs for all samples
        errors_for_each=zeros(number_of_sample,output_size);    %error values for each input output value
        
        %forwardpass loop
        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;        % calculate average sum of square for each epoch
       parameterVector= [W1(:);b1(:);W2(:);b2];              % recreate paramVector with new weights
       Jaco=jacobian(parameterVector, target_vals, input_vals,hidden_size); % calculate jacobian

       while true               % to calculate delta again when error is bigger than previous one
           
           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('Matrix must be symmetric and positive definite.');
                    
                end
                   
                % JTJLEyeInverse=inv(JTJ_Lambda_Eye);
    
                L = chol(JTJ_Lambda_Eye, 'lower');
                n=size(JTJ_Lambda_Eye,1);
                I=eye(n);
                Y=zeros(n);
                % Solve for Y using forward substitution
                for i = 1:n
                    Y(:, i) = L \ I(:, i);
                end
                % Now solve L^T * X = Y for X
                JTJLEyeInverse = L' \ Y;
           catch
                disp('Cholesky factorization failed.');
                break;
           end
           deltaValue=JTJLEyeInverse*JTe;
           parameterVectorTest=parameterVector;
           parameterVectorTest=parameterVectorTest-deltaValue;
           % check wheather the new error create better value.
           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);

           %test delta value whether they are too small for change or not
           if norm(deltaValue)<tol;
               disp("Parameter changes are too small to continue");
               break;
           end
           %test feedforward
           test_output_vals=zeros(number_of_sample,output_size);
           for i=1:number_of_sample
               test_output_vals(i,output_size) = forward_pass(input_vals(i,input_size), W1_test, b1_test, W2_test, b2_test);
           end

           test_errors_for_each = target_vals - test_output_vals;
           test_error_norm = norm(test_errors_for_each);
           
           if test_error_norm > errors_of_epochs(epoch_number,1)
               lambda=lambda*10;
               % if lambda>1000
               %     lambda=1000;
               % end
           else
               break;
           end
       end %finis point of bigger error and increasing 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(['Converged at epoch ', num2str(epoch_number), ' in run ', num2str(run)]);
           break;
       end

    end

    errors_of_epochs_for_average_5_run = errors_of_epochs_for_average_5_run + errors_of_epochs; %to calculate average of 5 run

end

% Plotting the result for investigating
% Calculate the average error across all runs
avg_epoch_errors = mean(errors_of_epochs_for_average_5_run, 1);

% Plot average error vs. epochs
plot_target=zeros(number_of_sample ,1);
 for i_test = 1:input_size
    X_test = input_vals(i_test);
    % Forward pass to get the current output
    plot_target(i_test) = forward_pass(X_test, W1, b1, W2, b2);
    % Calculate the current error  
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('Epoch');
ylabel('Average Error');
title('Average Error per Epoch across 5 Runs');
grid on;
hold off;


% Define forward pass function
function [outputf,A1] = forward_pass(X, W1, b1, W2, b2)
    Z1 = (W1 * X) + b1;                 % Linear combination for hidden layer
    A1 = sigmoid(Z1);                  % Sigmoid activation for hidden layer
    Z2 = (W2 * A1) + b2;                % Linear combination for output layer
    outputf = Z2;                      % Output layer
end
% Sigmoid activation function
function cikti = sigmoid(x)
    cikti = 1 ./ (1 + exp(-x));
end

% Derivative of the sigmoid function
function cikti = sigmoidDerivative(x)
    cikti = sigmoid(x) .* (1 - sigmoid(x));
end

% Define function for calculating Jacobian
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)            % to creat each line of jacobian for each input output pair
        [~ , jA1] = forward_pass(jinput_vals(i),W1,b1,W2,b2);
        derivW1 = -1 * W2'.* sigmoidDerivative(jA1) *jinput_vals(i);    %W1 partial derivatives
        derivb1 = -1 * W2'.* sigmoidDerivative(jA1);                    %b1 partial derivatives
        derivW2 = -1 * jA1;                                             %W2 partial derivatives
        derivb2 = -1;                                                   %b2 partial derivatives linear activation function
        jaco_return(i,:)= [derivW1(:);derivb1(:);derivW2(:);derivb2]; %each line of jacobian matrix       
    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

核心问题与修复方案 #

1. 参数更新符号错误（最关键） #

Levenberg-Marquardt的参数更新逻辑应为：参数 = 参数 + Δ，其中Δ = (J^T J + λI)^{-1} J^T e（e = target - output为残差）。你的代码错误地使用减法更新参数，直接导致误差反向增大：

parameterVectorTest=parameterVectorTest-deltaValue;

修复：改为加法更新：

parameterVectorTest=parameterVectorTest+deltaValue;

2. b1参数reshape维度错误 #

b1是隐藏层偏置，维度应为(hidden_size, 1)，但代码中错误地指定为hidden_size, input_size：

b1_test = reshape(..., hidden_size, input_size);

虽然当前input_size=1时表面无异常，但逻辑错误，后续修改input_size会直接报错。
修复：修正reshape目标维度：

b1_test = reshape(parameterVectorTest(hidden_size * input_size + 1:hidden_size * input_size + hidden_size), hidden_size, 1);

3. Delta过小判断的位置错误 #

代码在计算delta后立刻判断norm(deltaValue)<tol并终止迭代，此时尚未验证新参数的误差是否更优，会导致提前停止有效更新。
修复：将该判断移到参数更新逻辑之后，仅在确认更新有效时检查：

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);
    % 新增：参数变化过小则终止
    if norm(deltaValue)<tol
        disp("Parameter changes are too small to continue");
        break;
    end
end

4. Lambda无上限导致步长过小 #

当测试误差增大时，lambda持续乘以10但无上限，会导致步长趋近于0，无法有效更新参数。
修复：添加lambda上限：

if test_error_norm > errors_of_epochs(epoch_number,1)
    lambda=lambda*10;
    lambda=min(lambda, 1e6); % 限制lambda最大值
else
    break;
end

5. 绘图逻辑错误 #

绘图循环仅执行input_size次（input_size=1），导致只绘制第一个样本的预测值，无法展示整体拟合效果。
修复：修改循环范围为样本总数：

plot_target=zeros(number_of_sample ,1);
for i_test = 1:number_of_sample % 改为遍历所有样本
    X_test = input_vals(i_test);
    plot_target(i_test) = forward_pass(X_test, W1, b1, W2, b2);
end

内容的提问来源于stack exchange，提问作者SeAk