Метод levenberg-marquardt для решения нелинейных уравнений

Question

Я попытался реализовать метод levenberg-marquardt для решения нелинейных уравнений на Джулии на основе презентации Numerical Optimization using the Levenberg-Marquardt Algorithm. Этот мой код:

function get_J(ArrOfFunc,X,delta) 
    N = length(ArrOfFunc) 
    J = zeros(Float64,N,N) 
    for i = 1:N 
    for j=1:N 
     Temp = copy(X); 
     Temp[j]=Temp[j]+delta; 
     J[i,j] = (ArrOfFunc[i](Temp)-ArrOfFunc[i](X))/delta; 
    end 
    end 
    return J 
end 

function get_resudial(ArrOfFunc,Arg) 
    return map((x)->x(Arg),ArrOfFunc) 
end 

function lm_solve(Funcs,Init) 
    X = copy(Init) 
    delta = 0.01; 
    Lambda = 0.01; 
    Factor = 2; 
    J = get_J(Funcs,X,delta) 
    R = get_resudial(Funcs,X) 
    N = 5 
    for t = 1:N 

    G = J'*J+Lambda.*eye(length(X)) 
    dC = J'*R 
    C = sum(R.*R)/2; 
    Xnew = X-(inv(G)\dC); 
    Rnew = get_resudial(Funcs,Xnew) 
    Cnew = sum(Rnew.*Rnew)/2; 
    if (Cnew < C) 
     X = Xnew; 
     R = Rnew; 
     Lambda = Lambda/Factor; 
     J = get_J(Funcs,X,delta) 
    else 
     Lambda = Lambda*Factor; 
    end 
    if(maximum(abs(Rnew)) < 0.001) 
     return X 
    end 
    end 
    return X 
end 

function test() 
    ArrOfFunc = [ 
    (X)->X[1]+X[2]-2; 
    (X)->X[1]-X[2] 
    ]; 

    X = lm_solve(ArrOfFunc,Float64[3;3]) 
    println(X) 
    return X 
end 

Но с любой отправной точкой шаг не принимается. Что я делаю неправильно? Любая помощь будет оценена по достоинству.

Answer 1

У меня на данный момент нет способа проверить это, но одна строка не имеет смысла математически:

При вычислении Xnew она должна быть либо inv(G)*dC или G\dC, но не смесь обоих. Предпочтительно второе, поскольку для решения линейной системы не требуется вычисление обратной матрицы.

С этим неправильным вычислением в центре итерации траектория вычисления почти наверняка сбивается с пути.