DownhillSimplex[]関数の能力が、デグレードしているように見えます
目標関数として正定値な二次形式(下に凸な放物面)を使った場合の例は以下の通り
N = 1と 2のケースで明らかな劣化を確認しました
(* from Packages/optimize.n 1.42 *)
SimplexNEW=Class[{},{},{CurrentSimplex,Status="Stop",
limit,n,lo,tol,per,minq,exp,cont,quick,exl,cl,rr={1,-1}},
Amoebar[p_,f_,r_,kr_]:=Module[
{pr,fac1=(1-r)/(Length[p]-1),fr,tp=Take[p[[,2]],rr]},
! {pr,fac1=(1-r)/(L),fr,tp=Take[p[[,2]],rr]},
tp=(Plus@@tp)/Length[tp];
pr=(1-r)*tp+p[[kr,2]]*r;
{pr,fr=f[pr]+
(Plus@@MapThread[(#<#2||#>#3)&,{pr,limit[[1]],limit[[2]]}])*
Abs[p[[-1,1]]-p[[1,1]]]*10,
If[fr<p[[kr,1]],
Sort[ReplacePart[p,{fr,pr},kr]],
p]}];
Amoeba[p_,f_]:=Module[
{pr,pr1,pl,fr,fr1,prr,n,kr,fr0,pr0,pt},
Do[
Status=StandardForm["Reflect "//kr];
{pt,fr,pr}=Amoebar[p,f,-1,-kr];
(* Print["dh1 ",{fr,p[[exl,1]]}]; *)
If[kr==1,fr0=fr;pr0=pr];
If[fr <= p[[exl,1]],
Status="Extend";
{pt,fr1,pr1}=Amoebar[p,f,-exp,-kr];
Return[If[fr1 <= fr,
pr1,
pr]]],
(* If[fr1<=pr[[-1,1]],
Sort[ReplacePart[pr,{fr1,pt},-1]],
pr]]]], *)
{kr,cont}];
If[cont<>1,fr=fr0;pr=pr0];
If[Length[p]>2 && fr < p[[-1-cl,1]],
pr,
If[fr < p[[-1,1]],
Status="Contract+",
Status="Contract-"];
{pt,fr1,pr1}=Amoebar[pr,f,0.5,-1];
If[fr1 < pr[[-1-cl,1]],
pr1,
pl=pr[[1,2]];
prr=Thread[(pl+Thread[Drop[pr,1][[,2]]])/2];
prr=(MapThread[(If[#2<#<#3,(#>#4)*(#3-#2)+#2,#])&,
{#,pl-minq,pl+minq,pl}]&/@prr);
Do[
Status=StandardForm["Shrink "//n];
fv=f[prr[[n-1]]];
pr[[n]]={fv,prr[[n-1]]};
If[quick && (fv<pr[[n-1,1]] || fv>pr[[-1,1]]),
Break[]],
{n,Length[pr],2,-1}];
If[per<>1,
pr[[1]]={pr[[1,1]]*per+f[pl]*(1-per),pl}];
Sort[pr]]]];
DownhillSimplex[p0_,f_,opt___]:=Module[
{def={VariableRange:>{Table[Max[],{Length[p0[[1,2]]]}],
Table[Min[],{Length[p0[[1,2]]]}]},
MaxIteration:>Max[100,10*Length[p0]],
ExtendLevel->1,ContractLevel->1,ReflectRange->{1,-1},
Output->6,MinimumSize->0,ExpansionRatio->2,TryReflect->1,
Tolerance->1e-6,Persistence->1,QuickShrink->0},
lo,p=Sort[p0],nc=Floor[(PageWidth-3)/22]},
{limit,n,lo,tol,per,minq,exp,cont,quick,exl,cl,rr}=
{VariableRange,MaxIteration,Output,
Tolerance,Persistence,MinimumSize,ExpansionRatio,TryReflect,QuickShrink,
ExtendLevel,ContractLevel,ReflectRange}/.{opt}/.def;
CurrentSimplex=p;
exl=Restrict[exl,1,Length[p]];
cont=Restrict[cont,1,Length[p]];
Check[
Catch[
Do[CurrentSimplex=p=Amoeba[p,f];
If[lo > 0, StandardForm[
$FORM='10.6';
Write[lo,Drop[p,{nc+1,Length[p0]-nc}][[,1]]]]];
If[(p[[-1,1]]-p[[1,1]])/(1e-100+Abs[p[[-1,1]]]+Abs[p[[1,1]]]) < tol,
Break[]],{n}]],];
Status="Stop";
p];
];
(* from Packages/optimize.n 1.41 *)
SimplexOLD=Class[{},{},{CurrentSimplex,Status="Stop"},
Amoebar[p_,f_,r_,limit_,kr_]:=Module[
{pr,fac1=(1-r)/(Length[p]-1),fr},
pr=Plus@@p[[,2]]*fac1+p[[kr,2]]*(r-fac1);
{pr,fr=f[pr]+
(Plus@@MapThread[(#<#2||#>#3)&,{pr,limit[[1]],limit[[2]]}])*
Abs[p[[-1,1]]-p[[1,1]]]*10,
If[fr<p[[kr,1]],
Sort[ReplacePart[p,{fr,pr},kr]],
p]}];
Amoeba[p_,f_,limit_,per_,minq_,exp_,cont_,quick_,exl_]:=Module[
{pr,pr1,pl,fr,fr1,prr,n,kr,fr0,pr0,pt},
Do[
Status=StandardForm["Reflect "//kr];
{pt,fr,pr}=Amoebar[p,f,-1,limit,-kr];
(* Print["dh1 ",{fr,p[[exl,1]]}]; *)
If[kr==1,fr0=fr;pr0=pr];
If[fr <= p[[exl,1]],
Status="Extend";
{pt,fr1,pr1}=Amoebar[p,f,-exp,limit,-kr];
Return[If[fr1 <= fr,
pr1,
pr]]],
(* If[fr1<=pr[[-1,1]],
Sort[ReplacePart[pr,{fr1,pt},-1]],
pr]]]], *)
{kr,cont}];
If[cont<>1,fr=fr0;pr=pr0];
If[Length[p]>2 && fr < p[[-2,1]],
pr,
If[fr < p[[-1,1]],
Status="Contract+",
Status="Contract-"];
{pt,fr1,pr1}=Amoebar[pr,f,0.5,limit,-1];
If[fr1 < pr[[-1,1]],
pr1,
pl=pr[[1,2]];
prr=Thread[(pl+Thread[Drop[pr,1][[,2]]])/2];
prr=(MapThread[(If[#2<#<#3,(#>#4)*(#3-#2)+#2,#])&,
{#,pl-minq,pl+minq,pl}]&/@prr);
Do[
Status=StandardForm["Shrink "//n];
fv=f[prr[[n-1]]];
pr[[n]]={fv,prr[[n-1]]};
If[quick && (fv<pr[[n-1,1]] || fv>pr[[-1,1]]),
Break[]],
{n,Length[pr],2,-1}];
If[per<>1,
pr[[1]]={pr[[1,1]]*per+f[pl]*(1-per),pl}];
Sort[pr]]]];
DownhillSimplex[p0_,f_,opt___]:=Module[
{def={VariableRange:>{Table[Max[],{Length[p0[[1,2]]]}],
Table[Min[],{Length[p0[[1,2]]]}]},
MaxIteration:>Max[100,10*Length[p0]],
ExtendLevel->1,
Output->6,MinimumSize->0,ExpansionRatio->2,TryReflect->1,
Tolerance->1e-6,Persistence->1,QuickShrink->0},
limit,n,tol,per,minq,exp,cont,quick,exl,
lo,p=Sort[p0],nc=Floor[(PageWidth-3)/22]},
{limit,n,lo,tol,per,minq,exp,cont,quick,exl}=
{VariableRange,MaxIteration,Output,
Tolerance,Persistence,MinimumSize,ExpansionRatio,TryReflect,QuickShrink,
ExtendLevel}/.{opt}/.def;
CurrentSimplex=p;
exl=Restrict[exl,1,Length[p]];
cont=Restrict[cont,1,Length[p]];
Check[
Catch[
Do[CurrentSimplex=p=Amoeba[p,f,limit,per,minq,exp,cont,quick,exl];
If[lo > 0, StandardForm[
$FORM='10.6';
Write[lo,Drop[p,{nc+1,Length[p0]-nc}][[,1]]]]];
If[(p[[-1,1]]-p[[1,1]])/(1e-100+Abs[p[[-1,1]]]+Abs[p[[1,1]]]) < tol,
Break[]],{n}]],];
Status="Stop";
p];
];
The$SimplexOLD = SimplexOLD[];
The$SimplexNEW = SimplexNEW[];
func = With[{l = Range[Length[#]]},
Plus@@(Reverse[l] * (# + l)^2)]&;
f1 = func[{#1}]&;
f2 = func[{#1, #2}]&;
N = 1;
size = 1e-1;
x0 = Table[0, {N}];
p = {func[x0 + #], x0 + #}&/@Append[size * IdentityMatrix[N], Table[0, {N}]];
l = Transpose[Table[2 * N * {-1, 1}, {N}]];
opt = {VariableRange->l, Tolerance->1e-12, MaxIteration->1000};
rO = The$SimplexOLD@DownhillSimplex[p, func, Null@@opt];
rN = The$SimplexNEW@DownhillSimplex[p, func, Null@@opt];
Print["\nOLD ->\n", rO]; Print["\nNEW ->\n", rN];
Switch[N,
1, Plot[f1[z], {z, -2, 2}]; TkWait[],
2, ListContourPlot[Table[func[{i, j}], {i, -4, 4, 0.1}, {j, -4, 4, 0.1}],
MeshRange->{{-4, 4}, {-4, 4}}, Contours->30]; TkWait[],
_, (* Nothing *)];