我尝试基于 Numerical Optimization using the Levenberg-Marquardt Algorithm 在 Julia 上实现求解非线性方程的 levenberg-marquardt 方法推介会。这是我的代码:
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
但从任何起点开始,该步骤都不被接受。我做错了什么? 任何帮助将不胜感激。
最佳答案
我目前无法测试这一点,但是一行在数学上没有意义:
在计算Xnew
时,它应该是inv(G)*dC
或G\dC
,但不能是两者的混合。最好是第二个,因为线性系统的解不需要计算逆矩阵。
由于迭代中心的这一错误计算,计算的轨迹几乎肯定会误入歧途。
关于algorithm - 求解非线性方程的 levenberg-marquardt 方法,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/39780570/