关于优化例程/使用限制的建议
我正在尝试做一些数字计算，但很难确定解决问题的适当方法并寻求一些反馈。
到目前为止，我已经在 Mathematica 中完成了所有工作，但是，我相信现在是我需要对算法进行更多控制的时候了。
我还无法发布图片，所以这里有一个 链接 到它们 其中 H 是亥维赛阶跃函数。 C(k) 只是 C(r) 和 m=4 的 FT。在我的例子中，N 是 2000，因此您可以看到 omega 几乎是大量指数的总和。 rho 只是密度。 C(r) 如您所见，因为 m=4 具有不同的 a 系数。 IRISM 最终是 a 系数的函数。
我认为这三个方程在 Mathematica 中工作正常，但是我试图最小化 IRISM 并找到 4 个 a 值。我遇到的问题是，由于显而易见的原因，当积分中的对数等于零时会出现不连续性。我似乎找不到修改 Mathematica 算法的方法（它们是黑盒，这是正确的术语吗？）以便检查试验 a 值。我使用 Nelder-Meade 和差分进化并尝试不同的约束。然而，我似乎只能得到想象中的结果，显然是从负对数中得到的，或者如果我限制得足够好以避免明显仅局部最小值，因为我的结果与“正确”结果不匹配。我尝试了几次使用梯度的最小化算法，但运气不佳。
我认为前进的最佳方法是从头开始编写一个最小化例程，或者修改其他代码，这样我就可以在集成之前检查 IRISM 是否存在不连续性。我已经阅读了一些关于惩罚函数、对数屏障等的内容，但对于编程有些陌生，希望有人能够让我知道什么是好的开始方法。我认为最重要的是，关于优化的信息几乎太多了，我发现很难知道从哪里开始。
编辑：这是原始输入。如果我需要以不同的方式发布，请告诉我。
OverHat[c][a1_, a2_, a3_, a4_, k_] := (a1*(4*Pi*(Sin[k] - k*Cos[k])))/k^3 +
(a2*(4*Pi*(k*Sin[k] + 2*Cos[k] - 2)))/k^4 +
(a3*(8*Pi*(2*k - 3*Sin[k] + k*Cos[k])))/k^5 +
(a4*(-(24*Pi*(k^2 + k*Sin[k] + 4*Cos[k] - 4))))/k^6
Subscript[OverHat[\[Omega]], \[Alpha]\[Gamma]][k_, \[Alpha]\[Gamma]_, n_] :=
Exp[(-k^2)*\[Alpha]\[Gamma]*((n - \[Alpha]\[Gamma])/(6*n))]
OverHat[\[Omega]][k_] := Sum[Subscript[OverHat[\[Omega]], \[Alpha]\[Gamma]][k, \[Alpha]\[Gamma], n],
{\[Alpha]\[Gamma], 1, n}] /. n -> 2000
IRISM[a1_, a2_, a3_, a4_, \[Rho]_, kmax_] :=
\[Rho]^2*(1/15)*(20*a1 - 5*a2 + 2*a3 - a4)*Pi -
(1/(8*Pi^3))*NIntegrate[(\[Rho]*OverHat[\[Omega]][k]*OverHat[c][a1, a2, a3, a4, k] +
Log[1 - \[Rho]*OverHat[\[Omega]][k]*OverHat[c][a1, a2, a3, a4, k]])*4*Pi*k^2,
{k, 0, kmax}, WorkingPrecision -> 80]
NMinimize[IRISM[a1, a2, a3, a4, 0.9, 30], {a1, a2, a3, a4},
Method -> "DifferentialEvolution"]
I am trying to do some numerics and having a difficult time determining the appropriate way to solve a problem and looking for some feedback.
So far I have done all my work in Mathematica, however, I believe that the time has come where I need more control over my algorithms.
I can't post Images yet so here is a link to them
Where H is the heaviside stepfunciton. C(k) is just the FT of C(r) and m=4. N in my case is 2000 so you can see the omega is mearly the sum of a large number of exponentials. rho is just the densitiy. C(r) as you can see because m=4 has for different a coefficients. IRISM is ultimately a function of those for a coefficients.
I have these three equations working correctly I think within Mathematica however I am trying to minimize IRISM and find the 4 a values. The problem I am having is that, for obvious reasons, there is a discontinuity when the log with in the integral is equal to zero. I cannot seem to find a way to modify the Mathematica Algorithm (they are blackbox is that the right term?) so as to check the trial a values. I was using Nelder-Meade and Differential Evolution and attempting different constraints. However, I seemed to only get either imaginary results, obviously from a negative Log, or if I constrained well enough to avoid obviously only local minimum as my results did not match the "correct" results. I tried a few times with minimization algorithms that used gradients however I did not have much luck.
I think my best way to move forward is to just write a minimization routine from scratch, or modify other code, so as I can check IRISM ahead of integration for discontinuity. I have read up some on penalty function, log-barrier etc. but being somewhat new to programming was hoping someone might be able to let me know what a good approach would be to start off with. I think more than anything there is almost too much information out there on optimization and I am finding it difficult to know where to begin.
Edit: Here is the raw input. If I need to post it in a different way please let me know.
OverHat[c][a1_, a2_, a3_, a4_, k_] := (a1*(4*Pi*(Sin[k] - k*Cos[k])))/k^3 +
(a2*(4*Pi*(k*Sin[k] + 2*Cos[k] - 2)))/k^4 +
(a3*(8*Pi*(2*k - 3*Sin[k] + k*Cos[k])))/k^5 +
(a4*(-(24*Pi*(k^2 + k*Sin[k] + 4*Cos[k] - 4))))/k^6
Subscript[OverHat[\[Omega]], \[Alpha]\[Gamma]][k_, \[Alpha]\[Gamma]_, n_] :=
Exp[(-k^2)*\[Alpha]\[Gamma]*((n - \[Alpha]\[Gamma])/(6*n))]
OverHat[\[Omega]][k_] := Sum[Subscript[OverHat[\[Omega]], \[Alpha]\[Gamma]][k, \[Alpha]\[Gamma], n],
{\[Alpha]\[Gamma], 1, n}] /. n -> 2000
IRISM[a1_, a2_, a3_, a4_, \[Rho]_, kmax_] :=
\[Rho]^2*(1/15)*(20*a1 - 5*a2 + 2*a3 - a4)*Pi -
(1/(8*Pi^3))*NIntegrate[(\[Rho]*OverHat[\[Omega]][k]*OverHat[c][a1, a2, a3, a4, k] +
Log[1 - \[Rho]*OverHat[\[Omega]][k]*OverHat[c][a1, a2, a3, a4, k]])*4*Pi*k^2,
{k, 0, kmax}, WorkingPrecision -> 80]
NMinimize[IRISM[a1, a2, a3, a4, 0.9, 30], {a1, a2, a3, a4},
Method -> "DifferentialEvolution"]
如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱，如果其他用户或者作者回复了您的评论，将不能在第一时间通知您！
发布评论
评论（1）
如果 Mathematica 的
FindMinimum发现虚数，则会中止。即使您的目标在约束内是实值，这种情况也可能发生，因为对于默认的 Barrier 方法，它的精度控制很差，并且偶尔会超出范围。最简单的方法是将您的目标包装在Re中。如果您发布完整的代码，您可能会得到更好的答案。一些一般建议：
尝试简化 Mathematica 的目标比重新实现优化算法更容易。原因是一种算法失败通常意味着它是一个难题，其他算法也会失败。
下方绘制目标表面时，它是有意义的
我曾经遇到过一个问题，
FindMinimum发出警告，但未能收敛到正确的最小值，我可以通过分析确定该问题，它是用不同的方法发生的，当我在（来源：yaroslavvb.com）
在此在这种情况下，您可以看到问题在最小值处的条件非常糟糕（几乎是一个平台），并且最小值很难定位。
当存在不等式约束时，默认方法是 Barrier 方法，该方法成本较高且精度控制较差。非常低效的做法是将等式约束指定为不等式对，即用
a>=b和a<=b< 代替a=b/代码>。这可能会慢 3-10 倍，而且在数值上也更差——a 和 b 的结果可能只是近似相等。理想情况下，目标是得到一个凸问题，没有任何不等式约束并且条件良好。
Mathematica's
FindMinimumaborts if it sees an imaginary number. This can happen even if your objective is real-valued inside the constraints because for default Barrier method because it poor accuracy control and can occasionally step out of bounds. Simplest way around it is to wrap your objective insideRe. You may get better answers if you post complete code.Some general advice:
It's easier to try to simplify your objective for Mathematica than re-implement optimization algorithms. The reason is that one algorithm failing often means it's a difficult problem, and other algorithms will fail as well.
I once had a problem where
FindMinimumgave warnings and failed to converge to correct minimum, which I could determine analytically, it happened with different methods, and it made sense when I plotted the objective surface, below(source: yaroslavvb.com)
In this case, you can see the problem is very badly conditioned at the minimum (it's almost a plateau) and minimum is hard to localize.
When you have inequality constraints, default method is Barrier method, which is expensive and offers poor precision control. Very inefficient thing to do is to specify equality constraints as pairs of inequalities, ie instead of
a=b, havea>=banda<=b. This can be 3-10 times slower, and also numerically worse -- a and b might be only approximately equal in the result.Ideally the goal is to get a problem which is convex, doesn't have any inequality constraints and is well conditioned.