diff --git a/doc/img/fminbnd_forward.png b/doc/img/fminbnd_forward.png new file mode 100644 index 00000000..0a5f38d1 Binary files /dev/null and b/doc/img/fminbnd_forward.png differ diff --git a/doc/img/fminbnd_forward_2.png b/doc/img/fminbnd_forward_2.png new file mode 100644 index 00000000..d240fba3 Binary files /dev/null and b/doc/img/fminbnd_forward_2.png differ diff --git a/doc/img/fminbnd_principle.png b/doc/img/fminbnd_principle.png new file mode 100644 index 00000000..fb9f51ba Binary files /dev/null and b/doc/img/fminbnd_principle.png differ diff --git a/doc/img/ls-eg1.png b/doc/img/ls-eg1.png index ba714002..c42cc489 100644 Binary files a/doc/img/ls-eg1.png and b/doc/img/ls-eg1.png differ diff --git a/doc/img/ls-eg3.png b/doc/img/ls-eg3.png index 581a5368..d401d00d 100644 Binary files a/doc/img/ls-eg3.png and b/doc/img/ls-eg3.png differ diff --git a/doc/img/ls-equations.png b/doc/img/ls-equations.png index 2a8e92ce..c9de8179 100644 Binary files a/doc/img/ls-equations.png and b/doc/img/ls-equations.png differ diff --git a/doc/img/ls-line-vertical.png b/doc/img/ls-line-vertical.png index 7f807bdd..8f38b240 100644 Binary files a/doc/img/ls-line-vertical.png and b/doc/img/ls-line-vertical.png differ diff --git a/doc/img/ls-vector-vertical.png b/doc/img/ls-vector-vertical.png index 6440dd18..83e5f32a 100644 Binary files a/doc/img/ls-vector-vertical.png and b/doc/img/ls-vector-vertical.png differ diff --git a/doc/img/region_reduce.png b/doc/img/region_reduce.png new file mode 100644 index 00000000..2080e5d0 Binary files /dev/null and b/doc/img/region_reduce.png differ diff --git a/doc/img/upper_base.png b/doc/img/upper_base.png index f99f674c..3631af25 100644 Binary files a/doc/img/upper_base.png and b/doc/img/upper_base.png differ diff --git a/doc/rmvl.bib b/doc/rmvl.bib index 59e4407c..d4491085 100644 --- a/doc/rmvl.bib +++ b/doc/rmvl.bib @@ -1,8 +1,9 @@ -@misc{opencv, +@software{opencv, author = {opencv}, title = {opencv}, - year = {2001}, - url = {https://github.com/matheecs/ICRA-2019-DJI-RoboMaster} + license = {Apache-2.0}, + url = {https://github.com/matheecs/ICRA-2019-DJI-RoboMaster}, + year = {2001} } @misc{icra2019, author = {matheecs}, @@ -19,13 +20,13 @@ @article{kalman pages = {35--45}, year = {1960} } -@misc{onnx-rt, +@software{onnx-rt, author = {microsoft}, title = {onnx runtime inference}, year = {2021}, url = {https://github.com/microsoft/onnxruntime} } -@misc{apriltag, +@software{apriltag, author = {AprilRobotics}, title = {AprilTag 3}, year = {2023}, @@ -62,4 +63,50 @@ @book{numerical_analysis year = {2008}, publisher = {华南理工大学出版社}, address = {广州} -} \ No newline at end of file +} +@software{Agarwal_Ceres_Solver_2022, + author = {Agarwal, Sameer and Mierle, Keir and The Ceres Solver Team}, + title = {{Ceres Solver}}, + license = {Apache-2.0}, + url = {https://github.com/ceres-solver/ceres-solver}, + version = {2.2}, + year = {2023}, + month = {10} +} +@article{ConjGrad, + author = {Hestenes, Magnus R. and Stiefel, Eduard}, + title = {Methods of Conjugate Gradients for Solving Linear Systems}, + journal = {Journal of Research of the National Bureau of Standards}, + volume = {49}, + number = {6}, + pages = {409}, + year = {1952}, + month = {12}, + doi = {10.6028/jres.049.044}, + url = {https://nvlpubs.nist.gov/nistpubs/jres/049/jresv49n6p409_A1b.pdf} +} +@article{Ridders, + author = {C. J. F. Ridders}, + title = {Accurate computation of F'(x) and F'(x) F"(x)}, + journal = {Advances in Engineering Software}, + volume = {4}, + number = {2}, + pages = {75--76}, + year = {1978} +} +@article{NelderMead, + author = {Nelder, J. A. and Mead, R.}, + title = "{A Simplex Method for Function Minimization}", + journal = {The Computer Journal}, + volume = {7}, + number = {4}, + pages = {308-313}, + year = {1965}, + month = {01}, + abstract = "{A method is described for the minimization of a function of n variables, which depends on the comparison of function values at the (n + 1) vertices of a general simplex, followed by the replacement of the vertex with the highest value by another point. The simplex adapts itself to the local landscape, and contracts on to the final minimum. The method is shown to be effective and computationally compact. A procedure is given for the estimation of the Hessian matrix in the neighbourhood of the minimum, needed in statistical estimation problems.}", + issn = {0010-4620}, + doi = {10.1093/comjnl/7.4.308}, + url = {https://doi.org/10.1093/comjnl/7.4.308}, + eprint = {https://academic.oup.com/comjnl/article-pdf/7/4/308/1013182/7-4-308.pdf}, +} + diff --git a/doc/tutorials/modules/algorithm/auto_differential.md b/doc/tutorials/modules/algorithm/auto_differential.md new file mode 100644 index 00000000..b50c61ac --- /dev/null +++ b/doc/tutorials/modules/algorithm/auto_differential.md @@ -0,0 +1,168 @@ +自动求导、数值微分 {#tutorial_modules_auto_differential} +============ + +@author 赵曦 +@date 2024/05/06 +@version 1.0 +@brief 包含中心差商以及 Richardson 外推原理的介绍 + +@prev_tutorial{tutorial_modules_runge_kutta} + +@next_tutorial{tutorial_modules_fminbnd} + +@tableofcontents + +------ + +\f[ +\def\red#1{\color{red}{#1}} +\def\teal#1{\color{teal}{#1}} +\def\green#1{\color{green}{#1}} +\def\transparent#1{\color{transparent}{#1}} +\def\orange#1{\color{orange}{#1}} +\def\fml#1{\text{(#1)}} +\f] + +### 1. 基于 Taylor 公式的数值微分 + +利用以下两个 2 阶 Taylor 公式： + +\f[\begin{align}f(x+h)&=f(x)+f'(x)h+\frac{f''(x)}2h^2+\frac{f'''(x)}{3!}h^3\tag{1-1a}\\ +f(x-h)&=f(x)-f'(x)h+\frac{f''(x)}2h^2-\frac{f'''(x)}{3!}h^3\tag{1-1b}\end{align}\f] + +可以得到两个数值微分公式及其对应的截断误差（余项）\f$R\f$： + +\f[\begin{align}f'(x)&\approx\frac{f(x+h)-f(x)}h,\qquad R=-\frac{f''(\xi)}2h+o(h^2)\tag{1-2a}\\ +f'(x)&\approx\frac{f(x)-f(x-h)}h,\qquad R=\frac{f''(\eta)}2h+o(h^2)\tag{1-2b}\end{align}\f] + +其中\f$\fml{1-2a}\f$称为前向差商，\f$\fml{1-2b}\f$称为后向差商。由于截断误差的存在，差商法的精度不高，因此我们可以考虑使用更高阶的差商法。我们将\f$\fml{1-1a}\f$和\f$\fml{1-1b}\f$相减，得到： + +\f[f(x+h)-f(x-h)=2f'(x)h+\frac2{3!}f'''(\xi)h^3+o(h^3)\tag{1-3a}\f] + +整理得到 + +\f[f'(x)=\frac{f(x+h)-f(x-h)}{2h},\qquad R=-\frac{f'''(\xi)}{3!}h^2+o(h^2)\tag{1-3b}\f] + +如果导数\f$f'(x)\f$用\f$\fml{1-3b}\f$计算，那么截断误差为\f$o(h^2)\f$，比\f$\fml{1-2a}\f$和\f$\fml{1-2b}\f$的截断误差\f$O(h)\f$更小。这种方法称为中心差商。 + +### 2. Richardson 外推原理 + +#### 2.1 公式推导 + +将公式\f$\fml{1-3b}\f$扩展，可以得到更完全的写法： + +\f[f'(x)=\green{\frac{f(x+h)-f(x-h)}{2h}}-\left[\frac1{3!}f^{(3)}(x)h^2+\frac1{5!}f^{(5)}(x)h^4+\frac1{7!}f^{(7)}(x)h^6+\cdots\right]\tag{2-1}\f] + +一般的，可以记作 + +\f[L=\green{\varphi(h)}+\left[a_2h^2+a_4h^4+a_6h^6+\cdots\right]\tag{2-2}\f] + +正如上一节所说，如果\f$L\f$用\f$\varphi(h)\f$来近似表示，则截断误差是按\f$h^2\f$展开的幂级数，误差阶为\f$o\left(h^2\right)\f$。 + +事实上，对于公式\f$\fml{2-2}\f$，求解\f$L\f$的过程还可以继续向前推进（外推），不妨以\f$\frac h2\f$替换\f$\fml{2-2}\f$中的\f$h\f$，可得 + +\f[L=\varphi\left(\frac h2\right)+\left[a_2\frac{h^2}4+a_4\frac{h^4}{16}+a_6\frac{h^6}{64}+\cdots\right]\tag{2-3}\f] + +由\f$4\times\fml{2-3}-\fml{2-2}\f$可以得到 + +\f[L=\left[\frac43\varphi\left(\frac h2\right)-\frac13\varphi(h)\right]-\left[a_4\frac{h^4}4+5a_6\frac{h^6}{16}+\cdots\right]\tag{2-4}\f] + +上式表达了外推过程的第 1 步，说明\f$\varphi(h)\f$和\f$\varphi\left(\frac h2\right)\f$的一个线性组合提供了导数即\f$L\f$的新的计算公式，其精度提高至了\f$o\left(h^4\right)\f$ + +继续进行第 2 步，令 + +\f[\psi(h)=\frac43\varphi\left(\frac h2\right)-\frac13\varphi(h)=\frac1{4^1-1}\left[4^1\varphi\left(\frac h2\right)-\varphi(h)\right]\tag{2-5}\f] + +那么\f$\fml{2-4}\f$可以改写成 + +\f[L=\psi(h)+\left[b_4h^4+b_6h^6+\cdots\right]\tag{2-6}\f] + +与第 1 步的外推过程类似，以\f$\frac h2\f$替换\f$\fml{2-6}\f$中的\f$h\f$，可得 + +\f[L=\psi\left(\frac h2\right)+\left[b_4\frac{h^4}{16}+b_6\frac{h^6}{64}+\cdots\right]\tag{2-7}\f] + +由\f$4^2\times\fml{2-7}-\fml{2-6}\f$可以得到 + +\f[L=\left[\frac{16}{15}\psi\left(\frac h2\right)-\frac1{15}\psi(h)\right]-\left[b_6\frac{h^6}{20}+\cdots\right]\tag{2-8}\f] + +此式表达了外推过程的第 2 步，说明\f$\psi(h)\f$和\f$\psi\left(\frac h2\right)\f$的一个线性组合提供了导数即\f$L\f$的新的计算公式，其精度提高至了\f$o\left(h^6\right)\f$ + +如果我们令 + +\f[\theta(h)=\frac{16}{15}\psi\left(\frac h2\right)-\frac1{15}\psi(h)\tag{2-9}\f] + +不难猜出，新的导数计算公式为 + +\f[L\approx\frac{64}{63}\theta\left(\frac h2\right)-\frac1{63}\theta(h)\tag{2-10}\f] + +读者可以自行验证，按照这个思路可以进一步外推，可执行任意多步得到精度不断提高的新公式，这就是 Richardson 外推原理。 + +#### 2.2 外推算法总结 + +令\f$\varphi(h)=\frac{f(x+h)-f(x-h)}{2h}\f$，外推\f$M\f$步，则外推算法的步骤如下 + +① 选取一个适当的\f$h\f$值，计算\f$M+1\f$个数，记为\f$T(*,*)\f$ + +\f[T(n,0)=\varphi\left(\frac h{2^n}\right),\qquad n=0,1,2,\cdots,M\tag{2-11}\f] + +② 按公式计算 + +\f[T(n,k)=\frac1{4^k-1}\left[4^kT(n,k-1)-T(n-1,k-1)\right],\qquad\left\{\begin{align} +k&=1,2,\cdots,M\\n&=k,k+1,\cdots,M +\end{align}\right.\tag{2-12}\f] + +按照以上步骤计算得到的\f$T(n,k)\f$，满足等式 + +\f[L=T(n,k)+o(h^{2k+2})\qquad(h\to0)\tag{2-13}\f] + +即\f$T(n,k)\f$具备\f$2k+2\f$阶的精度。 + +### 3. 如何使用 + +RMVL 中提供了一元函数以及多元函数的微分工具，求解一元函数导数时可使用 rm::derivative ，求解多元函数梯度是需要使用 rm::grad 。 + +以下展示了使用自动求导、数值微分的例子。 + +#### 3.1 创建项目 + +1. 添加源文件 `main.cpp` + ```cpp + #include + #include + + // 自定义函数 f(x)=x²+4x-3 + inline double quadratic(double x) { return x * x + 4 * x - 3; } + + int main() + { + double dydx = rm::derivative(quadratic, 1); + printf("f'(1) = %f\n", dydx); + } + ``` + +2. 添加 `CMakeLists.txt` + ```cmake + cmake_minimum_required(VERSION 3.10) + project(DerivativeDemo) + find_package(RMVL COMPONENTS core REQUIRED) + add_executable(demo main.cpp) + target_link_libraries(demo PRIVATE rmvl_core) + ``` + +#### 3.2 构建、运行 + +在项目根目录打开终端，输入 + +```bash +mkdir build +cd build +cmake .. +make -j2 +./demo +``` + +可以看到运行结果 + +``` +f'(1) = 6 +``` diff --git a/doc/tutorials/modules/algorithm/ew_topsis.md b/doc/tutorials/modules/algorithm/ew_topsis.md index 26f78573..e999da75 100644 --- a/doc/tutorials/modules/algorithm/ew_topsis.md +++ b/doc/tutorials/modules/algorithm/ew_topsis.md @@ -4,7 +4,7 @@ @author RoboMaster Vision Community @date 2023/01/11 -@prev_tutorial{tutorial_modules_runge_kutta} +@prev_tutorial{tutorial_modules_fminunc} @next_tutorial{tutorial_modules_kalman} diff --git a/doc/tutorials/modules/algorithm/fminbnd.md b/doc/tutorials/modules/algorithm/fminbnd.md new file mode 100644 index 00000000..4b1be81f --- /dev/null +++ b/doc/tutorials/modules/algorithm/fminbnd.md @@ -0,0 +1,158 @@ +一维最优化方法 {#tutorial_modules_fminbnd} +============ + +@author 赵曦 +@date 2024/05/06 +@version 1.0 +@brief 包含区间搜索以及黄金分割法的介绍 + +@prev_tutorial{tutorial_modules_auto_differential} + +@next_tutorial{tutorial_modules_fminunc} + +@tableofcontents + +------ + +\f[ +\def\red#1{\color{red}{#1}} +\def\teal#1{\color{teal}{#1}} +\def\green#1{\color{green}{#1}} +\def\transparent#1{\color{transparent}{#1}} +\def\orange#1{\color{orange}{#1}} +\def\fml#1{\text{(#1)}} +\f] + +### 1. 区间搜索 + +在使用一维优化方法搜索目标函数的极小值点时，首先要确定搜索区间\f$[a,b]\f$，这个搜索区间要包含目标函数的极小值点，而且是单峰区间，即在这个区间内，目标函数只有一个极小值点。 + +在单峰区间内，函数值具有“高—低—高”的特点。根据这一特点，可以采用进退法来找到搜索区间。 + +进退法分为 + +- 初始探察确定进退 +- 前进或者后退寻查 + +两步，其步骤如下 + +
+ +![图 1-1 向前一步](fminbnd_forward.png) + +![图 1-2 向前两步](fminbnd_forward.png) + +
+ +1. 选择一个初始点\f$\alpha_1\f$和一个初始步长\f$h\f$。 +2. 如图 1-1 所示，计算点\f$\alpha_1\f$和点\f$\alpha_1+h\f$对应的函数值\f$f(\alpha_1)\f$和\f$f(\alpha_1+h)\f$，令\f[f_1=f(\alpha_1),\quad f_2=f(\alpha_1+h)\f] +3. 比较\f$f_1\f$和\f$f_2\f$，若\f$f_1>f_2\f$，则执行前进运算，将步长加大\f$k\f$倍（如加大 2 倍），取新点\f$\alpha_1+3h\f$，如图 1-2 所示，计算其函数值，并令\f[f_1=f(\alpha_1+h),\quad f_2=f(\alpha_1+3h)\f]若\f$f_1< f_2\f$，则初始搜索区间端点为\f[a=\alpha_1,\quad b=\alpha_1+3h\f]若\f$f_1=f_2\f$，则初始搜索区间端点为\f[a=\alpha_1+h,\quad b=\alpha_1+3h\f]若\f$f_1>f_2\f$，则要继续做前进运算，且步长再加大两倍，取第 4 个点\f$\alpha_1+7h\f$，再比较第 3 和第 4 个点处的函数值……如此反复循环，直到在连续的 3 个点的函数值出现“两头大、中间小”的情况为止。 +4. 如果在 **步骤 3** 中出现\f$f_1< f_2\f$的情况，则执行后退运算，将步长变为负值，取新点\f$\alpha_1-h\f$，计算函数值，令\f[f_1=f(\alpha_1-h),\quad f_2=f(\alpha_1)\f]若\f$f_1>f_2\f$，则初始搜索区间端点为\f[a=\alpha_1-h,\quad b=\alpha_1+h\f]若\f$f_1=f_2\f$，则初始搜索区间端点为\f[a=\alpha_1-h,\quad b=\alpha_1\f]若\f$f_1< f_2\f$，则要继续做后退运算，且步长再加大两倍，取第 4 个点\f$\alpha_1-3h\f$，再比较第 3 和第 4 个点处的函数值……如此反复循环，直到在连续的 3 个点的函数值出现“两头大、中间小”的情况为止。 + +**示例** + +试用进退法确定目标函数\f$f(\alpha)=\alpha^2-5\alpha+8\f$的一维优化初始搜索区间\f$[a,b]\f$。设初始点\f$\alpha_1=0\f$，初始步长\f$h=1\f$ + +**解：** 由初始点\f$\alpha_1=0\f$，初始步长\f$h=1\f$，得\f[f_1=f(\alpha_1)=8,\quad f_2=f(\alpha_1+h)=4\f]因为\f$f_1>f_2\f$，所以执行前进运算，将步长加大 2 倍，取新点\f$\alpha_1+3h=3\f$，令\f[f_1=f(\alpha_1+h)=4,\quad f_2=f(\alpha_1+3h)=2\f]因为\f$f_1>f_2\f$，应继续做前进运算，且步长再加大 2 倍，取第 4 个点\f$\alpha_1+7h=7\f$。令\f[f_1=f(\alpha_1+3h)=2,\quad f_2=f(\alpha_1+7h)=22\f]此时\f$f_1< f_2\f$，在连续的 3 个点\f$\alpha_1+h\f$，\f$\alpha_1+3h\f$，\f$\alpha_1+7h\f$的函数值出现了“两头大、中间小”的情况，则初始搜索区间端点为\f[a=\alpha_1+h=1,\quad b=\alpha_1+7h=7\f] + +### 2. 黄金分割法 + +#### 2.1 区间消去 + +黄金分割法是利用区间消去法的原理，通过不断缩小单峰区间长度，即每次迭代都消去一部分不含极小值点的区间，使搜索区间不断缩小，从而逐渐逼近目标函数极小值点的一种优化方法。黄金分割法是直接寻优法，通过直接比较区间上点的函数值的大小来判断区间的取舍。这种方法具有计算简单、收敛速度快等优点。 + +如图1-2所示，在已确定的单峰区间[a,b]内任取\f$\alpha_1\f$，\f$\alpha_2\f$两点，计算并比较两点处的函数值\f$f(\alpha_1)\f$和\f$f(\alpha_2)\f$,可能出现3种情况： + +1. \f$f(\alpha_1)< f(\alpha_2)\f$，因为函数是单峰的，所以极小值点必定位于点\f$\alpha_2\f$的左侧，即\f$\alpha^*\in[a,\alpha_2]\f$，搜索区间可以缩小为\f$[a,\alpha_2]\f$，如图 2-1(a) 所示； +2. \f$f(\alpha_1)>f(\alpha_2)\f$，极小值点必定位于点\f$\alpha_1\f$的右侧，即\f$\alpha^*\in[\alpha_1,b]\f$，搜索区间可以缩小为\f$[\alpha_1,b]\f$，如图 2-1(b) 所示； +3. \f$f(\alpha_1)=f(\alpha_2)\f$，极小值点必定位于点\f$\alpha_1\f$和\f$\alpha_2\f$之间，即\f$\alpha^*\in[\alpha_1,\alpha_2]\f$，搜索区间可缩小为\f$[\alpha_1,\alpha_2]\f$，如图 2-1(c) 所示。 + +
+ +![图 2-1 搜索区间缩小示意图](region_reduce.png) + +
+ +根据上述方法，可在新搜索区间里再取两个新点比较函数值来继续缩小区间，但这样做效率较低，应该考虑利用已经计算过的区间内剩下的那个点。对于以上的 (1)、(2) 两种情况，可以在新搜索区间内取一点\f$x_3\f$和该区间内剩下的那个点（第 (1) 种情况的\f$\alpha_1\f$点或第 (2) 种情况的\f$\alpha_2\f$点）进行函数值的比较来继续缩短搜索区间。而第 (3) 种情况则要取两个新点进行比较，为统一起见，将前面 3 种情况综合为以下两种情况: + +1. 若\f$f(\alpha_1)< f(\alpha_2)\f$，则将搜索区间缩小为\f$[a,\alpha_2]\f$ +2. 若\f$f(\alpha_2)\ge f(\alpha_2)\f$，则将搜索区间缩小为\f$[\alpha_1,b]\f$。 + +#### 2.2 黄金分割法原理 + +黄金分割法就是基于上述原理来选择区间内计算点的位置的，它有以下要求： +1. 点\f$\alpha_1\f$和\f$\alpha_2\f$相对区间\f$[a,b]\f$的边界要对称布置，即区间\f$[a,\alpha_1)\f$的大小与区间\f$(\alpha_2,b]\f$的大小相等； +2. 每次计算一个新点，要求保留的区间长度\f$l\f$与原区间长度\f$L\f$之比等于被消去的区间长度\f$L-l\f$与保留区间长度\f$l\f$之比，即要求下式成立：\f[\frac lL=\frac{L-l}l\tag{2-1}\f] + +令\f[\lambda=\frac lL\f]将上式代入式\f$\fml{2-1}\f$，有\f[\lambda^2+\lambda-1=0\tag{2-2}\f]求解式\f$\fml{2-2}\f$，得\f[\lambda=\frac{\sqrt5-1}2\approx0.618\f] + +该方法保证每次迭代都以同一比率缩小区间，缩短率为\f$0.618\f$，故黄金分割法又称为\f$0.618\f$法。保留的区间长度为整个区间长度的\f$0.618\f$倍，消去的区间长度为整个区间长度的\f$0.382\f$倍。 + +黄金分割法的计算步骤如下： + +
+ +![图 2-2 黄金分割法原理](fminbnd_principle.png) + +
+ +1. 在\f$[a,b]\f$内取两点\f$\alpha_1\f$，\f$\alpha_2\f$，如图 2-2(a) 所示，使\f[\begin{align}\alpha_1&=a+0.382(b-a)\\\alpha_2&=a+0.618(b-a)\end{align}\f]令\f$f_1=f(\alpha_1)\f$，\f$f_2=f(\alpha_2)\f$。 +2. 比较\f$f_1\f$和\f$f_2\f$。当\f$f_1< f_2\f$时，消去区间\f$(\alpha_2,b]\f$。做置换\f$b=\alpha_2\f$，\f$\alpha_2=\alpha_1\f$并另取新点\f$\alpha_1=a+0.382(b-a)\f$，如图 2-2(b) 所示，令\f$f_1=f(\alpha_1)\f$。当\f$f_1\ge f_2\f$时，消去区间\f$[a,\alpha_1)\f$。做置换\f$a=\alpha_1\f$，\f$\alpha_1=\alpha_2\f$，\f$f_1=f_2\f$，并另取新点\f$\alpha_2=a+0.618(b-a)\f$，如图 2-2(c) 所示，令 \f$f_2=f(\alpha_2)\f$。 +3. 检查终止条件。若\f$b-a\le\epsilon\f$，则输出最优解\f$\alpha^*=\frac12(a+b)\f$和最优值\f$f(\alpha^*)\f$；否则转第 (2) 步。 + +### 3. 如何使用 + +RMVL 中提供了一维寻优的函数 rm::fminbnd ，以下展示了一维寻优的例子。 + +#### 3.1 创建项目 + +1. 添加源文件 `main.cpp` + ```cpp + #include + #include + + // 自定义函数 f(x)=x²+4x-3 + inline double quadratic(double x) { return x * x + 4 * x - 3; } + + int main() + { + // 确定搜索区间 + auto [x1, x2] = rm::region(quadratic, 0); + // 添加优化选项（容许误差和最大迭代次数） + rm::OptimalOptions options; + options.max_iter = 100; // 最大迭代次数是 100，不加以设置默认是 1000 + options.tol = 1e-4; // 容许误差是 1e-4，不加以设置默认是 1e-6 + // 一维寻优 + auto [x, fval] = rm::fminbnd(quadratic, x1, x2, options); + printf("x = %f\n", x); + printf("fval = %f\n", fval); + } + ``` + +2. 添加 `CMakeLists.txt` + ```cmake + cmake_minimum_required(VERSION 3.10) + project(FminbndDemo) + find_package(RMVL COMPONENTS core REQUIRED) + add_executable(demo main.cpp) + target_link_libraries(demo PRIVATE rmvl_core) + ``` + +#### 3.2 构建、运行 + +在项目根目录打开终端，输入 + +```bash +mkdir build +cd build +cmake .. +make -j2 +./demo +``` + +可以看到运行结果 + +``` +x = -2.000000 +fval = -7.000000 +``` \ No newline at end of file diff --git a/doc/tutorials/modules/algorithm/fminunc.md b/doc/tutorials/modules/algorithm/fminunc.md new file mode 100644 index 00000000..31ec157c --- /dev/null +++ b/doc/tutorials/modules/algorithm/fminunc.md @@ -0,0 +1,17 @@ +多维无约束最优化方法 {#tutorial_modules_fminunc} +============ + +@author 赵曦 +@date 2024/05/06 +@version 1.0 +@brief 包含最速下降、共轭梯度、Nelder-Mead 单纯形法的介绍 + +@prev_tutorial{tutorial_modules_fminbnd} + +@next_tutorial{tutorial_modules_ew_topsis} + +@tableofcontents + +------ + +@warning 未完成 \ No newline at end of file diff --git a/doc/tutorials/modules/algorithm/least_square.md b/doc/tutorials/modules/algorithm/least_square.md index ab5c7f12..97769e52 100644 --- a/doc/tutorials/modules/algorithm/least_square.md +++ b/doc/tutorials/modules/algorithm/least_square.md @@ -19,9 +19,9 @@ 最小二乘法诞生于统计学，其最初的目标是使用一条直线拟合一系列离散的点，那么我们该如何寻找这一直线？最小二乘法的原理是使得拟合直线的误差的平方和最小。这里的误差定义成每一个离散的点\f$x_i\f$与拟合直线的 **距离** 。在这一最初的用法中，点到直线的距离能够很直观的描述成点到直线的 **垂线** 的长度。
-![ls-line-vertical](ls-line-vertical.png) -图 1-1 +![图 1-1](ls-line-vertical.png) +
如图 1-1 所示，令直线外一点为\f$P\f$，直线为\f$l\f$，垂线为\f$l_0\f$，垂足为点\f$O\f$。最小二乘法就是构建一条直线，使得这些点\f$P_i\f$到这条直线的距离平方和最小，此时\f$O_i\f$点被称为 **最小二乘解** 。 @@ -57,9 +57,9 @@ 我们先用一个最简单的例子。
-![ls-vector-vertical](ls-vector-vertical.png) -图 1-2 +![图 1-2](ls-vector-vertical.png) +
对于图 1-2 中使用向量表示的方式，最小二乘解可以表示为向量\f$OP\f$的坐标，在这种情况下，\f$BP\f$与\f$OA\f$垂直，即： @@ -86,16 +86,14 @@ 这里 \f$A=\begin{bmatrix}2&1\\1&-1\\1&1\end{bmatrix}\f$，\f$\pmb b=(1,0,2)^T\f$，在图 1-3 中表示如下。
-![ls-equations](ls-equations.png) -图 1-3 +![图 1-3](ls-equations.png) +
一般的，对于一个系数矩阵 \f$A=(a_{ij})_{s\times{n}}\f$，\f$\pmb b=(b_1,b_2,\cdots,b_s)^T\f$，\f$\pmb x=(x_1,x_2,\cdots,x_s)^T\f$，若满足 \f$\text{rank}(A)\leq s\f$ 时，线性方程组没有数值解，但我们希望找到一组 **最优解** ，衡量此最优解的方法仍然可以采用最小二乘法。设法找出一组解 \f$\hat{\pmb x}=(x_1^0,\ x_2^0,\ x_3^0,\ \cdots,\ x_n^0)\f$ 使得每一项的误差 \f$\delta_i\f$ 平方和最小，如何定量这个误差？能否继续采用最小二乘法的点到直线的最短距离作为出发点？答案是肯定的，这里先给出误差平方和的表达式。 -\f[ -\delta^2=\sum_{i=0}^{s-1}(a_{i1}x_1+a_{i2}x_2+\cdots+a_{in}x_n-b_i)^2\tag{2-2} -\f] +\f[\delta^2=\sum_{i=0}^{s-1}(a_{i1}x_1+a_{i2}x_2+\cdots+a_{in}x_n-b_i)^2\tag{2-2}\f] 上式也可以写成 @@ -174,9 +172,9 @@ y&=x_1 这是一个参数方程，消去 \f$x_1\f$ 可以得到 \f$y=x\f$，这就是系数矩阵 \f$A=(1,\ 1)^T\f$ 所生成的列空间。
-![ls-eg1](ls-eg1.png) -图 2-1 +![图 2-1](ls-eg1.png) +
相当于现在需要在 \f$y=x\f$ 上找到一个点，使得其到 \f$\pmb b\f$ 的距离最短，这是初等几何的内容，做垂线即可，最终得到交点的坐标 @@ -196,9 +194,9 @@ y&=x_1 列空间外向量（右端项）：\f$\pmb b=(1,0,2)^T\f$
-![ls-eg2](ls-eg2.png) -图 2-2 +![图 2-2](ls-eg2.png) +
需要满足 \f$(\pmb y-\pmb b)\perp\alpha\f$，则需要分别满足 \f$(\pmb y-\pmb b)\perp\pmb\alpha_1\f$ 和 \f$(\pmb y-\pmb b)\perp\pmb\alpha_2\f$，下面使用几何法对最小二乘法的原理进行验证。 @@ -408,9 +406,9 @@ d_1=(\phi_1,f)&=\sum_{i=0}^3\phi_1(t_i)f(t_i)=\sum_{i=0}^3t_i\times f(t_i)\\&=0+ 解得：\f$\left\{\begin{align}a_0&=-0.5\\a_1&=0.8\end{align}\right.\f$，即拟合曲线为 \f$y=-0.5+0.8t\f$，使用数学绘图软件验证，结果正确。
-![ls-eg3](ls-eg3.png) -图 4-1 +![图 4-1](ls-eg3.png) +
### 5. 部署使用 diff --git a/doc/tutorials/modules/tools/opcua.md b/doc/tutorials/modules/tools/opcua.md index c2642bb3..00c05e81 100644 --- a/doc/tutorials/modules/tools/opcua.md +++ b/doc/tutorials/modules/tools/opcua.md @@ -35,9 +35,9 @@ OPC UA 的设计目标是建立一种通用的、独立于厂商和平台的通 #### 地址空间
-![opcua](opcua.svg) -图 1-1 OPC UA 地址空间模型 +![图 1-1 OPC UA 地址空间模型](opcua.svg) +
1. 对象类型节点 rm::ObjectType ：提供对象的定义，即对象的抽象，与类相当，且子类可以继承父类的特征，方便模型的扩充。该节点包括对象的各种数据类型，数据的语义，以及控制方式。OPC UA 命名空间 `0` 中规定了多个基础的对象类型节点。如使用最广的 BaseObjectType（在 RMVL 中表示为 `rm::nodeBaseObjectType`），所有对象类型节点都需要继承该节点再进行扩充。在对具体设备建模的过程中，应该将设备组成的各部分分解为不同的对象分别建模，再用引用节点将各部分按照实际设备中的关系相关联，从而得到完整设备的对象类型节点。 diff --git a/doc/tutorials/modules/tutorial_modules.md b/doc/tutorials/modules/tutorial_modules.md index 15073309..e0117a32 100644 --- a/doc/tutorials/modules/tutorial_modules.md +++ b/doc/tutorials/modules/tutorial_modules.md @@ -33,6 +33,9 @@ - @subpage tutorial_modules_least_square - @subpage tutorial_modules_func_iteration - @subpage tutorial_modules_runge_kutta +- @subpage tutorial_modules_auto_differential +- @subpage tutorial_modules_fminbnd +- @subpage tutorial_modules_fminunc #### 数据与信号处理 diff --git a/modules/core/CMakeLists.txt b/modules/core/CMakeLists.txt index 2f2772f8..30215d24 100644 --- a/modules/core/CMakeLists.txt +++ b/modules/core/CMakeLists.txt @@ -39,6 +39,6 @@ if(BUILD_PERF_TESTS) rmvl_add_test( core Performance DEPENDS core - DEPEND_TESTS benchmark::benchmark_main + DEPEND_TESTS benchmark::benchmark_main opencv_video ) endif(BUILD_PERF_TESTS) diff --git a/modules/core/include/rmvl/core/math.hpp b/modules/core/include/rmvl/core/math.hpp index 696aa5db..3858f930 100644 --- a/modules/core/include/rmvl/core/math.hpp +++ b/modules/core/include/rmvl/core/math.hpp @@ -121,31 +121,31 @@ template constexpr Tp rad2deg(Tp rad) { return rad * static_cast(180) / static_cast(PI); } /** - * @brief Point类型转换为Matx类型 + * @brief Point 类型转换为 Matx 类型 * * @tparam Tp 数据类型 - * @param[in] point Point类型变量 - * @return Matx类型变量 + * @param[in] point Point 类型变量 + * @return Matx 类型变量 */ template constexpr cv::Matx point2matx(cv::Point3_ point) { return cv::Matx(point.x, point.y, point.z); } /** - * @brief Matx类型转换为Point类型 + * @brief Matx 类型转换为 Point 类型 * * @tparam Tp 数据类型 - * @param[in] matx Matx类型变量 - * @return Point类型变量 + * @param[in] matx Matx 类型变量 + * @return Point 类型变量 */ template constexpr cv::Point3_ matx2point(cv::Matx matx) { return cv::Point3_(matx(0), matx(1), matx(2)); } /** - * @brief Matx类型转换为Vec类型 + * @brief Matx 类型转换为 Vec 类型 * * @tparam Tp 数据类型 - * @param[in] matx Matx类型变量 - * @return Vec类型变量 + * @param[in] matx Matx 类型变量 + * @return Vec 类型变量 */ template constexpr cv::Vec matx2vec(cv::Matx matx) { return cv::Vec(matx(0), matx(1), matx(2)); } @@ -421,6 +421,156 @@ typename ForwardIterator::value_type calculateModeNum(ForwardIterator first, For ->first; } +/** + * @brief 使用 `std::vector` 表示的向量加法 + * + * @tparam T 数据类型 + * @param[in] vec1 向量 1 + * @param[in] vec2 向量 2 + * @return 和向量 + */ +template +inline std::vector operator+(const std::vector &vec1, const std::vector &vec2) +{ + std::vector retval(vec1.size()); + std::transform(vec1.cbegin(), vec1.cend(), vec2.cbegin(), retval.begin(), std::plus()); + return retval; +} + +/** + * @brief 使用 `std::vector` 表示的向量减法 + * + * @tparam T 数据类型 + * @param[in] vec1 向量 1 + * @param[in] vec2 向量 2 + * @return 差向量 + */ +template +inline std::vector operator-(const std::vector &vec1, const std::vector &vec2) +{ + std::vector retval(vec1.size()); + std::transform(vec1.cbegin(), vec1.cend(), vec2.cbegin(), retval.begin(), std::minus()); + return retval; +} + +/** + * @brief 使用 `std::vector` 表示的向量自加 + * + * @tparam T 数据类型 + * @param[in] vec1 向量 1 + * @param[in] vec2 向量 2 + * @return 和向量 + */ +template +inline std::vector &operator+=(std::vector &vec1, const std::vector &vec2) +{ + std::transform(vec1.cbegin(), vec1.cend(), vec2.cbegin(), vec1.begin(), std::plus()); + return vec1; +} + +/** + * @brief 使用 `std::vector` 表示的向量自减 + * + * @tparam T 数据类型 + * @param[in] vec1 向量 1 + * @param[in] vec2 向量 2 + * @return 差向量 + */ +template +inline std::vector &operator-=(std::vector &vec1, const std::vector &vec2) +{ + std::transform(vec1.cbegin(), vec1.cend(), vec2.cbegin(), vec1.begin(), std::minus()); + return vec1; +} + +/** + * @brief 使用 `std::vector` 表示的向量取反 + * + * @tparam T 数据类型 + * @param[in] vec 向量 + * @return 向量取反后的结果 + */ +template +inline std::vector operator-(const std::vector &vec) +{ + std::vector retval(vec.size()); + std::transform(vec.cbegin(), vec.cend(), retval.begin(), std::negate()); + return retval; +} + +/** + * @brief 使用 `std::vector` 表示的向量乘法（数乘） + * + * @tparam T 数据类型 + * @param[in] vec 向量 + * @param[in] val 数乘因子 + * @return 向量乘法（数乘） + */ +template +inline std::vector operator*(const std::vector &vec, T val) +{ + std::vector retval(vec.size()); + std::transform(vec.cbegin(), vec.cend(), retval.begin(), [val](const T &x) { return x * val; }); + return retval; +} + +/** + * @brief 使用 `std::vector` 表示的向量乘法（数乘） + * + * @tparam T 数据类型 + * @param[in] val 数乘因子 + * @param[in] vec 向量 + * @return 向量乘法（数乘） + */ +template +inline std::vector operator*(T val, const std::vector &vec) { return vec * val; } + +/** + * @brief 使用 `std::vector` 表示的向量乘法（数乘） + * + * @tparam T 数据类型 + * @param[in] vec 向量 + * @param[in] val 数乘因子 + * @return 向量乘法（数乘） + */ +template +inline std::vector &operator*=(std::vector &vec, T val) +{ + std::transform(vec.cbegin(), vec.cend(), vec.begin(), [val](const T &x) { return x * val; }); + return vec; +} + +/** + * @brief 使用 `std::vector` 表示的向量除法（数乘） + * + * @tparam T 数据类型 + * @param[in] vec 向量 + * @param[in] val 除数 + * @return 向量除法（数乘） + */ +template +inline std::vector operator/(const std::vector &vec, T val) +{ + std::vector retval(vec.size()); + std::transform(vec.cbegin(), vec.cend(), retval.begin(), [val](const T &x) { return x / val; }); + return retval; +} + +/** + * @brief 使用 `std::vector` 表示的向量除法（数乘） + * + * @tparam T 数据类型 + * @param[in] val 除数 + * @param[in] vec 向量 + * @return 向量除法（数乘） + */ +template +inline std::vector &operator/=(std::vector &vec, T val) +{ + std::transform(vec.cbegin(), vec.cend(), vec.begin(), [val](const T &x) { return x / val; }); + return vec; +} + //! @} core } // namespace rm diff --git a/modules/core/include/rmvl/core/numcal.hpp b/modules/core/include/rmvl/core/numcal.hpp index dd2ffe8d..83188bcf 100644 --- a/modules/core/include/rmvl/core/numcal.hpp +++ b/modules/core/include/rmvl/core/numcal.hpp @@ -11,11 +11,11 @@ #pragma once -#include -#include #include #include #include +#include +#include //! @addtogroup core //! @{ @@ -23,6 +23,10 @@ //! @{ //! @brief 包含函数插值、曲线拟合、非线性方程（组）数值解、常微分方程数值解等数值计算算法 //! @} core_numcal +//! @defgroup core_optimal 最优化算法库 +//! @{ +//! @brief 包含一维函数最小值搜索、无约束多维函数最小值搜索等最优化算法 +//! @} core_optimal //! @} core namespace rm @@ -59,12 +63,7 @@ class Polynomial * @param[in] x 指定点的 x 坐标 * @return 多项式在指定点的函数值 */ - inline double operator()(double x) const noexcept - { - double y = _coeffs.back(); - std::for_each(_coeffs.rbegin() + 1, _coeffs.rend(), [&](double val) { y = y * x + val; }); - return y; - } + double operator()(double x) const noexcept; }; ///////////////////// 函数插值 ///////////////////// @@ -305,4 +304,92 @@ class RungeKutta4 : public RungeKutta //! @} core_numcal +//! @addtogroup core_optimal +//! @{ + +//! 一维函数 +using Func1d = std::function; + +//! 多维函数 +using FuncNd = std::function &)>; + +//! 梯度/导数计算模式 +enum DiffMode : uint8_t +{ + Diff_Central, //!< 中心差商 + Diff_Ridders, //!< Richardson 外推 +}; + +//! 多维函数最优化模式 +enum OptimizeMode : uint8_t +{ + Optm_ConjGrad, //!< 共轭梯度法 + Optm_Simplex, //!< 单纯形法 +}; + +//! 优化选项 +struct OptimalOptions +{ + DiffMode diff_mode{}; //!< 梯度计算模式，默认为中心差商 `Central` + OptimizeMode optm_mode{}; //!< 优化模式，默认为共轭梯度法 `ConjGrad` + int max_iter{1000}; //!< 最大迭代次数 + double dx{1e-2}; //!< 求解步长 + double tol{1e-6}; //!< 误差容限 +}; + +/** + * @brief 计算一元函数的导数 + * + * @param[in] func 多元函数 + * @param[in] x 指定位置的自变量 + * @param[in] mode 导数计算模式，默认为中心差商 `Diff_Central` + * @param[in] dx 坐标的微小增量，默认为 `1e-3` + * @return 函数在指定点的导数 + */ +double derivative(Func1d func, double x, DiffMode mode = Diff_Central, double dx = 1e-3); + +/** + * @brief 计算多元函数的梯度 + * + * @param[in] func 多元函数 + * @param[in] x 指定位置的自变量 + * @param[in] mode 梯度计算模式，默认为中心差商 `Diff_Central` + * @param[in] dx 计算偏导数时，坐标的微小增量，默认为 `1e-3` + * @return 函数在指定点的梯度向量 + */ +std::vector grad(FuncNd func, const std::vector &x, DiffMode mode = Diff_Central, double dx = 1e-3); + +/** + * @brief 采用进退法确定搜索区间 + * + * @param[in] func 一维函数 + * @param[in] x0 初始点 + * @param[in] delta 搜索步长 + * @return 搜索区间 + */ +std::pair region(Func1d func, double x0, double delta = 1); + +/** + * @brief 一维函数最小值搜索 + * + * @param[in] func 一维函数 + * @param[in] x1 搜索区间左端点 + * @param[in] x2 搜索区间右端点 + * @param[in] options 优化选项，默认为 `1e-4` 的误差容限和 `100` 的最大迭代次数 + * @return `[x, fval]` 最小值点和最小值 + */ +std::pair fminbnd(Func1d func, double x1, double x2, const OptimalOptions &options = {}); + +/** + * @brief 无约束多维函数最小值搜索 \cite ConjGrad \cite NelderMead + * + * @param[in] func 多维函数 + * @param[in] x0 初始点 + * @param[in] options 优化选项，默认为 `1e-4` 的误差容限和 `100` 的最大迭代次数 + * @return `[x, fval]` 最小值点和最小值 + */ +std::pair, double> fminunc(FuncNd func, const std::vector &x0, const OptimalOptions &options = {}); + +//! @} core_optimal + } // namespace rm diff --git a/modules/core/perf/perf_kalman.cpp b/modules/core/perf/perf_kalman.cpp index 72ade041..8e98bdbb 100644 --- a/modules/core/perf/perf_kalman.cpp +++ b/modules/core/perf/perf_kalman.cpp @@ -10,11 +10,10 @@ */ #include - +#include #include #include "rmvl/core/kalman.hpp" -#include namespace rm_test { @@ -76,7 +75,7 @@ void kalman42_opencv(benchmark::State &state) } } -BENCHMARK(kalman42_rmvl)->Name("KalmanFilter42 - RMVL")->Iterations(10); -BENCHMARK(kalman42_opencv)->Name("KalmanFilter42 - OpenCV")->Iterations(10); +BENCHMARK(kalman42_rmvl)->Name("kf (x_dim: 4, z_dim: 2) - by rmvl ")->Iterations(10); +BENCHMARK(kalman42_opencv)->Name("kf (x_dim: 4, z_dim: 2) - by opencv")->Iterations(10); } // namespace rm_test diff --git a/modules/core/perf/perf_optimal.cpp b/modules/core/perf/perf_optimal.cpp new file mode 100644 index 00000000..59732f98 --- /dev/null +++ b/modules/core/perf/perf_optimal.cpp @@ -0,0 +1,105 @@ +/** + * @file perf_optimal.cpp + * @author zhaoxi (535394140@qq.com) + * @brief 最优化算法库基准测试 + * @version 1.0 + * @date 2024-05-04 + * + * @copyright Copyright 2024 (c), zhaoxi + * + */ + +#include +#include + +#include "rmvl/core/numcal.hpp" + +namespace rm_test +{ + +/////////////////////// Quadratic Function /////////////////////// + +static double quadraticFunc(const std::vector &x) +{ + const auto &x1 = x[0], &x2 = x[1]; + return 60 - 10 * x1 - 4 * x2 + x1 * x1 + x2 * x2 - x1 * x2; +}; + +static void cg_quadratic_rmvl(benchmark::State &state) +{ + for (auto _ : state) + rm::fminunc(quadraticFunc, {0, 0}); +} + +class Quadratic : public cv::MinProblemSolver::Function +{ +public: + int getDims() const override { return 2; } + + double calc(const double *x) const override + { + const auto &x1 = x[0], &x2 = x[1]; + return 60 - 10 * x1 - 4 * x2 + x1 * x1 + x2 * x2 - x1 * x2; + } +}; + +static void cg_quadratic_cv(benchmark::State &state) +{ + for (auto _ : state) + { + auto quadratic = cv::makePtr(); + auto solver = cv::ConjGradSolver::create(quadratic); + cv::Vec2d x; + solver->minimize(x); + } +} + +BENCHMARK(cg_quadratic_rmvl)->Name("fminunc (conj_grad, quadratic) - by rmvl ")->Iterations(50); +BENCHMARK(cg_quadratic_cv)->Name("fminunc (conj_grad, quadratic) - by opencv")->Iterations(50); + +/////////////////////// Rosenbrock Function /////////////////////// + +static double rosenbrockFunc(const std::vector &x) +{ + const auto &x1 = x[0], &x2 = x[1]; + return 100 * (x2 - x1 * x1) * (x2 - x1 * x1) + (1 - x1) * (1 - x1); +}; + +void splx_rosenbrock_rmvl(benchmark::State &state) +{ + for (auto _ : state) + { + rm::OptimalOptions options; + options.optm_mode = rm::Optm_Simplex; + rm::fminunc(rosenbrockFunc, {5, 3}, options); + } +} + +class Rosenbrock : public cv::DownhillSolver::Function +{ +public: + int getDims() const override { return 2; } + + double calc(const double *x) const override + { + const auto &x1 = x[0], &x2 = x[1]; + return 100 * (x2 - x1 * x1) * (x2 - x1 * x1) + (1 - x1) * (1 - x1); + } +}; + +void splx_rosenbrock_cv(benchmark::State &state) +{ + for (auto _ : state) + { + auto rosenbrock = cv::makePtr(); + auto solver = cv::DownhillSolver::create(rosenbrock); + cv::Vec2d x(5, 3); + solver->setInitStep(x); + solver->minimize(x); + } +} + +BENCHMARK(splx_rosenbrock_rmvl)->Name("fminunc (simplex, rosenbrock) - by rmvl ")->Iterations(50); +BENCHMARK(splx_rosenbrock_cv)->Name("fminunc (simplex, rosenbrock) - by opencv")->Iterations(50); + +} // namespace rm_test diff --git a/modules/core/src/numcal.cpp b/modules/core/src/numcal.cpp index e90dd90b..11945799 100644 --- a/modules/core/src/numcal.cpp +++ b/modules/core/src/numcal.cpp @@ -9,10 +9,9 @@ * */ -#include - #include +#include "rmvl/core/math.hpp" #include "rmvl/core/numcal.hpp" #include "rmvl/core/util.hpp" @@ -21,6 +20,13 @@ namespace rm { +double Polynomial::operator()(double x) const noexcept +{ + double y = _coeffs.back(); + std::for_each(_coeffs.rbegin() + 1, _coeffs.rend(), [&](double val) { y = y * x + val; }); + return y; +} + ///////////////////// 函数插值 ///////////////////// Interpolator::Interpolator(const std::vector &xs, const std::vector &ys) @@ -139,36 +145,6 @@ RungeKutta::RungeKutta( RMVL_Error(RMVL_StsBadArg, "\"r[i].size()\" must be greater than or equal to the \"i\"."); } -// 加法 -template -inline static std::vector operator+(const std::vector &vec1, const std::vector &vec2) -{ - std::vector retval(vec1.size()); - std::transform(vec1.cbegin(), vec1.cend(), vec2.cbegin(), retval.begin(), std::plus()); - return retval; -} - -// 自加 -template -inline static std::vector &operator+=(std::vector &vec1, const std::vector &vec2) -{ - std::transform(vec1.cbegin(), vec1.cend(), vec2.cbegin(), vec1.begin(), std::plus()); - return vec1; -} - -// 数乘 -template -inline static std::vector operator*(const std::vector &vec, T val) -{ - std::vector retval(vec.size()); - std::transform(vec.cbegin(), vec.cend(), retval.begin(), [val](T x) { return x * val; }); - return retval; -} - -// 数乘 -template -inline static std::vector operator*(T val, const std::vector &vec) { return vec * val; } - /** * @brief 计算数值解 * diff --git a/modules/core/src/optimal.cpp b/modules/core/src/optimal.cpp new file mode 100644 index 00000000..cb08a3bf --- /dev/null +++ b/modules/core/src/optimal.cpp @@ -0,0 +1,286 @@ +/** + * @file optimal.cpp + * @author zhaoxi (535394140@qq.com) + * @brief 最优化算法库 + * @version 1.0 + * @date 2024-05-02 + * + * @copyright Copyright 2024 (c), zhaoxi + * + */ + +#include + +#include + +#include "rmvl/core/math.hpp" +#include "rmvl/core/numcal.hpp" + +namespace rm +{ + +// 中心差商计算一元函数导数 +static inline double partial(Func1d func, double x_dx, double dx) +{ + x_dx += dx; + double f1 = func(x_dx); + x_dx -= 2 * dx; + double f2 = func(x_dx); + x_dx += dx; + return (f1 - f2) / (2 * dx); +} + +// 中心差商计算多元函数偏导数 +static inline double partial(FuncNd func, std::vector &x_dx, std::size_t idx, double dx) +{ + x_dx[idx] += dx; + double f1 = func(x_dx); + x_dx[idx] -= 2 * dx; + double f2 = func(x_dx); + x_dx[idx] += dx; + return (f1 - f2) / (2 * dx); +} + +double derivative(Func1d func, double x, DiffMode mode, double dx) +{ + double x_dx{x}; + if (mode == Diff_Ridders) + { + double T00{partial(func, x_dx, 2 * dx)}; + double T10{partial(func, x_dx, dx)}; + double T20{partial(func, x_dx, dx / 2)}; + double T11{(4. * T10 - T00) / 3.}; + double T21{(4. * T20 - T10) / 3.}; + return (16. * T21 - T11) / 15.; + } + else + return partial(func, x_dx, dx); +} + +/** + * @brief 计算多元函数的梯度（无返回值） + * + * @param[in] func 多元函数 + * @param[in] x 指定位置的自变量 + * @param[out] xgrad 函数在指定点的梯度向量 + * @param[in] mode 梯度计算模式 + * @param[in] dx 计算偏导数时的步长 + * @return 函数在指定点的梯度向量 + */ +static void calcGrad(FuncNd func, const std::vector &x, std::vector &xgrad, DiffMode mode, double dx) +{ + std::vector x_dx{x}; + if (mode == Diff_Ridders) + for (std::size_t i = 0; i < x_dx.size(); ++i) + { + double T00{partial(func, x_dx, i, 2 * dx)}; + double T10{partial(func, x_dx, i, dx)}; + double T20{partial(func, x_dx, i, dx / 2)}; + double T11{(4. * T10 - T00) / 3.}; + double T21{(4. * T20 - T10) / 3.}; + xgrad[i] = (16. * T21 - T11) / 15.; + } + else + for (std::size_t i = 0; i < x_dx.size(); ++i) + xgrad[i] = partial(func, x_dx, i, dx); +} + +std::vector grad(FuncNd func, const std::vector &x, DiffMode mode, double dx) +{ + std::vector ret(x.size()); + calcGrad(func, x, ret, mode, dx); + return ret; +} + +/** + * @brief 计算向量的二范数 + * + * @param[in] x 向量 + * @return 二范数 + */ +static inline double norm(const std::vector &x) +{ + double retval{}; + for (auto &&v : x) + retval += v * v; + return std::sqrt(retval); +} + +std::pair region(Func1d func, double x0, double delta) +{ + double f1{func(x0)}, f2{func(x0 + delta)}; + if (f1 > f2) + { + f1 = f2; + f2 = func(x0 + 3 * delta); + return f1 < f2 ? std::make_pair(x0, x0 + 3 * delta) : (f1 == f2 ? std::make_pair(x0 + delta, x0 + 3 * delta) : region(func, x0 + delta, 2 * delta)); + } + else if (f1 < f2) + { + f2 = f1; + f1 = func(x0 - delta); + return f1 > f2 ? std::make_pair(x0 - delta, x0 + delta) : (f1 == f2 ? std::make_pair(x0 - delta, x0) : region(func, x0 - 3 * delta, 2 * delta)); + } + else + return {x0, x0 + delta}; +} + +std::pair fminbnd(Func1d func, double x1, double x2, const OptimalOptions &options) +{ + constexpr double phi = 0.618033988749895; + double a1{x1 + (1.0 - phi) * (x2 - x1)}, a2{x1 + phi * (x2 - x1)}; + double f1{func(a1)}, f2{func(a2)}; + for (int i = 0; i < options.max_iter; ++i) + { + if (std::abs(x1 - x2) < options.tol) + break; + if (f1 < f2) // 替换右端点 + { + x2 = a2, a2 = a1, f2 = f1; + a1 = x1 + (1.0 - phi) * (x2 - x1); + f1 = func(a1); + } + else // 替换左端点 + { + x1 = a1, a1 = a2, f1 = f2; + a2 = x1 + phi * (x2 - x1); + f2 = func(a2); + } + } + return {0.5 * (x1 + x2), func(0.5 * (x1 + x2))}; +} + +// 共轭梯度法 +static double fminunc_cg(FuncNd func, std::vector &xk, const OptimalOptions &options) +{ + std::vector s = -xk; + std::vector xk_grad = grad(func, xk, options.diff_mode, options.dx), xk2_grad(xk_grad.size()); + // 判断是否收敛 + double nbl_xk = norm(xk_grad); + if (nbl_xk < options.tol) + return func(xk); + // 一维搜索函数 + auto func_alpha = [&](double alpha) { return func(xk + alpha * s); }; + double retfval{}; + for (int i = 0; i < options.max_iter; ++i) + { + // 一维搜索 alpha + auto [a, b] = region(func_alpha, 1); + RMVL_Assert(b > 0); + auto [alpha, fval] = fminbnd(func_alpha, a < 0. ? 0. : a, b, options); + // 更新 xk，并计算对应的梯度 + for (std::size_t j = 0; j < xk.size(); ++j) + xk[j] += alpha * s[j]; + retfval = fval; + calcGrad(func, xk, xk2_grad, options.diff_mode, options.dx); + auto nbl_xk2 = norm(xk2_grad); + if (nbl_xk2 < options.tol) + break; + // 计算 beta + double nbl_xk_xk2 = std::inner_product(xk_grad.begin(), xk_grad.end(), xk2_grad.begin(), 0.0); + double beta = std::max((nbl_xk2 * nbl_xk2 - nbl_xk_xk2) / (nbl_xk * nbl_xk), 0.0); + // 更新搜索方向 + for (std::size_t j = 0; j < s.size(); ++j) + s[j] = -xk2_grad[j] + beta * s[j]; + xk_grad = xk2_grad; + } + return retfval; +} + +// 单纯形法 +static double fminunc_splx(FuncNd func, std::vector &xk, const OptimalOptions &options) +{ + const std::size_t dim = xk.size(); // 维度 + const std::size_t N = dim + 1; // 单纯形点的个数 + std::vector, double>> splx(N); + for (auto &p : splx) + p = {xk, 0}; + for (std::size_t i = 0; i < dim; ++i) + splx[i + 1].first[i] += 100 * options.dx; + for (std::size_t i = 0; i < N; ++i) + splx[i].second = func(splx[i].first); + // 单纯形迭代 + for (int i = 0; i < options.max_iter; ++i) + { + std::sort(splx.begin(), splx.end(), [](const auto &a, const auto &b) { return a.second < b.second; }); + // 求出除最大值点外的所有点的中心 + std::vector xc(dim); + std::for_each(splx.begin(), splx.end() - 1, [&](const auto &vp) { xc += vp.first; }); + xc /= static_cast(N - 1); + // 反射 xr = xc + alpha * (xc - xh)，其中 alpha = 1 + std::vector xr(dim); + for (std::size_t j = 0; j < dim; ++j) + xr[j] = 2 * xc[j] - splx.back().first[j]; + + double fxr = func(xr); + // f(xn-1) <= f(xr)，反射点函数值大于最差点，则重新计算反射点（反压缩） + if (splx.back().second <= fxr) + { + // 压缩 xs = xc - beta * (xc - xh)，其中 beta = 0.5 + xr = xc + splx.back().first; + xr /= 2.0; + fxr = func(xr); + } + + // f(xr) < f(x0)，反射点函数值小于最优点 + if (fxr < splx[0].second) + { + // 扩展 xe = xc + gamma * (xc - xh)，其中 gamma = 2 + std::vector xe(dim); + for (std::size_t j = 0; j < dim; ++j) + xe[j] = 3 * xc[j] - 2 * splx.back().first[j]; + double fxe = func(xe); + splx.back() = fxe < fxr ? std::make_pair(xe, fxe) : std::make_pair(xr, fxr); + } + // f(x0) <= f(xr) < f(xn-2)，反射点函数值大于最优点，小于次差点 + else if (splx[0].second <= fxr && fxr < splx[N - 2].second) + splx.back() = {xr, fxr}; + // f(xn-2) <= f(xr) < f(xn)，反射点函数值大于次差点，小于最差点 + else + { + // 压缩 xs = xc + beta * (xc - xh)，其中 beta = 0.5 + std::vector xs(dim); + for (std::size_t j = 0; j < dim; ++j) + xs[j] = 1.5 * xc[j] - 0.5 * splx.back().first[j]; + double fxs = func(xs); + if (fxs < splx.back().second) + splx.back() = {xs, fxs}; + else + { + for (std::size_t i = 1; i < N; ++i) + { + for (std::size_t j = 0; j < dim; ++j) + splx[i].first[j] = 0.5 * (splx[i].first[j] + splx[0].first[j]); + splx[i].second = func(splx[i].first); + } + } + } + // 判断是否收敛 + if (std::abs(splx.back().second - splx.front().second) < options.tol) + break; + } + xk = std::move(splx[0].first); + return splx[0].second; +} + +std::pair, double> fminunc(FuncNd func, const std::vector &x0, const OptimalOptions &options) +{ + if (x0.empty()) + RMVL_Error(RMVL_StsBadArg, "x0 is empty"); + std::vector xk{x0}; + double fval{}; + switch (options.optm_mode) + { + case Optm_ConjGrad: + fval = fminunc_cg(func, xk, options); + break; + case Optm_Simplex: + fval = fminunc_splx(func, xk, options); + break; + default: + break; + } + return {xk, fval}; +} + +} // namespace rm diff --git a/modules/core/test/test_optimal.cpp b/modules/core/test/test_optimal.cpp new file mode 100644 index 00000000..0ee7af9a --- /dev/null +++ b/modules/core/test/test_optimal.cpp @@ -0,0 +1,70 @@ +/** + * @file test_optimal.cpp + * @author zhaoxi (535394140@qq.com) + * @brief 最优化算法库测试 + * @version 1.0 + * @date 2024-05-04 + * + * @copyright Copyright 2024 (c), zhaoxi + * + */ + +#include +#include + +#include + +namespace rm_test +{ + +static inline double f1d(double x) { return x * x - 4 * x + 7; } +static inline double rosenbrock(const std::vector &x) +{ + const auto &x1 = x[0], &x2 = x[1]; + return 100 * (x2 - x1 * x1) * (x2 - x1 * x1) + (1 - x1) * (1 - x1); +} +static inline double quadratic(const std::vector &x) +{ + const auto &x1 = x[0], &x2 = x[1]; + return 60 - 10 * x1 - 4 * x2 + x1 * x1 + x2 * x2 - x1 * x2; +} + +TEST(Optimal, derivative) +{ + auto func = [](double x) { return std::pow(x, 7) - std::pow(x, 5); }; + auto val_central = rm::derivative(func, 1); + EXPECT_NEAR(val_central, 2, 1e-3); + auto val_ridders = rm::derivative(func, 1, rm::Diff_Ridders); + EXPECT_NEAR(val_ridders, 2, 1e-12); +} + +TEST(Optimal, fminbnd) +{ + auto [x1, x2] = rm::region(f1d, -5); + EXPECT_LE(x1, 2); + EXPECT_GE(x2, 2); + auto [x, fval] = rm::fminbnd(f1d, x1, x2); + EXPECT_NEAR(x, 2, 1e-4); + EXPECT_NEAR(fval, 3, 1e-4); +} + +TEST(Optimal, fminunc_simplex) +{ + rm::OptimalOptions options; + options.optm_mode = rm::Optm_Simplex; + auto [x, fval] = rm::fminunc(rosenbrock, {1, -2}, options); + EXPECT_NEAR(x[0], 1, 1e-2); + EXPECT_NEAR(x[1], 1, 1e-2); + EXPECT_NEAR(fval, 0, 1e-2); +} + +TEST(Optimal, fminunc_conjgrad) +{ + rm::OptimalOptions options; + auto [x, fval] = rm::fminunc(quadratic, {0, 0}, options); + EXPECT_NEAR(x[0], 8, 1e-4); + EXPECT_NEAR(x[1], 6, 1e-4); + EXPECT_NEAR(fval, 8, 1e-4); +} + +} // namespace rm_test