在Numerical Optimization using the Levenberg-Marquardt Algorithm表示的基础上,实现了求解朱莉娅非线性方程组的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但从任何一个起点来看,这一步都不被接受。我做错什么了?任何帮助都将不胜感激。
发布于 2016-09-30 19:40:16
我目前无法检验这一点,但有一行在数学上没有意义:
在计算Xnew时,应该是inv(G)*dC或G\dC,而不是两者兼而有之。最好是第二种,因为线性系统的解不需要求逆矩阵。
由于这一错误的计算处于迭代的中心,计算的轨迹几乎肯定会误入歧途。
https://stackoverflow.com/questions/39780570
复制相似问题