\[
I_c=\int_0^1\frac{\cos x}{\sqrt{x}}\mathrm{d}x
\]

\[
I_s=\int_0^1\frac{\sin x}{\sqrt{x}}\mathrm{d}x
​\]

 

使用梯形公式或辛普森公式进行近似计算.

 

(a)  Use the composite trapezoidal rule with $n$ intervals of $[0,1]$.  "ignoring" the singularity at $x=0$.

(b)  Use the composite trapezoidal rule over the interval $[h,1]$ where $h=1/n$ in combination with a weighted Newton-Cotes rule .

证明 $x=\cos x$ 在 $(0,1)$ 内存在惟一的实根. 并求这个根的近似值.

设 $f(x)=x^5-4x^3-5$, 则 $f'(x)=5x^4-12x^2$. 令

\[
x_{n+1}=x_n-\frac{f(x_n)}{f'(x_n)}
\]

>> :mode polyn
Switch into polynomial mode.

>> diff(x^5-4x^3-5)
out> 5x^4-12x^2

------------------------

>> (x^5-4x^3-5)/(5x^4-12x^2)
in> (x^5-4x^3-5)/(5x^4-12x^2)

out>
 quotient> q(x) = 1|5x
remainder> r(x) = -8|5x^3-5

1|5x^1
------------------------

 

---

>> printRecursiveSeries(4/5*x_n+(8/5*x_n^3+5)/(5*x_n^4-12*x_n^2),x_n,1,100,\n,linenumber)
[1]     1
[2]     -0.1428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571428571429
[3]     -20.68684146042636608674344523401127174712080372457730948296986032835089438863023768684146042636608673092
[4]     -16.565023688617589109177537826433863432338260623780573240602345316225300420691945262445633989624803938736
[5]     -13.2714938151392127327086618324645366416606984746301979545605997238256662746354987884645019492889088259888
[6]     -10.64160729106856435253004280224457254617986844856116673479583215765055870684575566382988446238459617659104
[7]     -8.543927904421098082103266940566651336467393310841794835891510921669912488167846412120893765286900283972832
[8]     -6.8736750146567822836834747098409137278847103761224278023001309485647375068881413772811833775034090294782656
[9]     -5.54751385154459757117959853911308426694533117603318730701466364274880256207172022811686144943490408568261248
[10]    -4.499428383849649530577659730552688681850614839831792743978496600619449792627229845215124665110187212646089984
...  ...
[95]    -1.7523882613158145756642021439965730423292471867922368214477915180260352744599353088675504105795106405000000039951342135619708337669405749916685322147717924146834921482445434074138489452088000512
[96]    -1.75238826131581457566420214399657304232924718679223682144779151802603527445993530886755041057951064050000000319610737084957666701355245999333482577181743393174679371859563472593107915616704004096
[97]    -1.752388261315814575664202143996573042329247186792236821447791518026035274459935308867550410579510640500000002556885896679661333610841967994667860617453947145397434974876507780744863324933632032768
[98]    -1.7523882613158145756642021439965730423292471867922368214477915180260352744599353088675504105795106405000000020455087173437290668886735743957342884939631577163179479799012062245958906599469056262144
[99]    -1.75238826131581457566420214399657304232924718679223682144779151802603527445993530886755041057951064050000000163640697387498325351093885951658743079517052617305435838392096497967671252795752450097152
[100]   -1.752388261315814575664202143996573042329247186792236821447791518026035274459935308867550410579510640500000001309125579099986602808751087613269944636136420938443486707136771983741370022366019600777216

 

 

>> x=-1.752388261315814575664202143996573042329247186792236821447791518026035274459935308867550410579510640500000001309125579099986602808751087613269944636136420938443486707136771983741370022366019600777216
----------------------------
 type: double
 name: x
value: -0
--------------------
>> (x)^5-4*(x)^3-5
in> (-1.752388261315814575664202143996573042329247186792236821447791518026035274459935308867550410579510640500000001309125579099986602808751087613269944636136420938443486707136771983741370022366019600777216)^5-4*(-1.752388261315814575664202143996573042329247186792236821447791518026035274459935308867550410579510640500000001309125579099986602808751087613269944636136420938443486707136771983741370022366019600777216)^3-5

out> 0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000100246145459086898585802895870188576837095998950884918298641158723559282188775705998581874986989258205424303464763848840790191389632164176548802253105448816681594005165936343156722773516711893034676728450651453570649817816469002947814009596245777511620864799261738696376838799747784277730251509134065093393799300694265837562131421519346343974503745560305120305614651107308855715555811229082159622948908292166108381373370499489108074739974508281096479430654358092447625637625382610938780421530534478487833251837308685695313740517610641168059282723155215195454469839074609502737089508716321411361434807371548709509131679916412796742165441812738587766584562366496393623176677562279651152507508747918974460653536570760308781440856424720183327725525430635283070792546407203728655490798266409749217358117945684603185503201696956723823148920764841223782621217436808496957051984681847936038599655424
------------------------

 

From graph to manifold Laplacians: The convergence rate

Author:  A. Singer

Appl. Comput. Harmon. Anal. 21(2006)  128--134

 

以下是论文的读书笔记.

---

 

1. 介绍

图拉普拉斯在机器学习中被广泛用到, 主要用于维数约减、半监督学习(semi-supervised learning)和谱聚类(spectral clustering)(参考 [3--9] 及其所列文献). 

设 $\mathcal{M}\subset\mathbb{R}^m$ 是维数为 $d < m$ 的黎曼流形, 给定 $N$ 个数据点 $\{x_1, x_2, \ldots, x_N\}\subset\mathcal{M}$. 这些点被视为所在空间 $\mathbb{R}^m$ 中的向量. 任务是找到未知的底层流形 $\mathcal{M}$ 以及该流形的几何和低维表示.

谱方法的出发点是从如下所示的一个合适的半正定核(semi-positive kernel) $k$ 中抽取出一个 $N\times N$ 的权重矩阵 $W$,

\[
W_{ij}=k\Bigl(\|x_i-x_j\|^2/(2\varepsilon)\Bigr),
\]

这里 $\|\cdot\|$ 是外围空间 $\mathbb{R}^m$ 中的欧氏距离, $\varepsilon > 0$ 是核的带宽(bandwidth of the kernel). 关于核的一个流行选择是选用指数核 $k(x)=e^{-x}$, 当然选用其他核也是可以的.

加权矩阵 $W$ 被归一化为行随机矩阵，方法是将其除以对角矩阵 $D$，其对角线元素是 $W$ 相应行的行和.

\[
D=(D_{ij})_{N\times N},\quad D_{ij}=0, \text{若}\ i\neq j;\quad D_{ii}=\sum_{j=1}^{N}W_{ij}.
\]

并且（负）图拉普拉斯((negative defined) graph Laplacian) $L$ 定义如下:

\[
L=D^{-1}W-I,
\]

其中 $I$ 是 $N\times N$ 单位阵.

 

当数据点 $\{x_i\}_{i=1}^{N}$ 在流形 $\mathcal{M}$ 上是独立一致分布的(independently uniformly distributed), 图拉普拉斯收敛到流形上的连续 Laplace-Beltrami 算子 $\Delta_M$.  这种说法有两种表现形式. 首先, 若 $f:\mathcal{M}\rightarrow\mathbb{R}$ 是一光滑函数, 则关于偏差项(bias term)的以下结果已经由不同的作者([1, 5, 6] 等)建立.

\[
\frac{1}{\varepsilon}\lim_{N\rightarrow\infty}\sum_{j=1}^{N}L_{ij}f(x_j)=\frac{1}{2}\Delta_M f(x_i)+O(\varepsilon^{1/2}).\tag{1.4}
\]

其次, Hein 等[1] 建立了关于 $N$ 和 $\varepsilon$(方差项, the variance term)误差的一致估计, 即当 $N\rightarrow\infty$ 时按 $\frac{1}{\sqrt{N}}$ 递减; 但当 $\varepsilon\rightarrow 0$ 时以 $\frac{1}{\varepsilon^{1+d/4}}$ 递增,

\[
\frac{1}{\varepsilon}\sum_{j=1}^{N}L_{ij}f(x_j)=\frac{1}{2}\Delta_M f(x_i)+O\biggl(\frac{1}{N^{1/2}\varepsilon^{1+d/4}},\varepsilon^{1/2}\biggr).\tag{1.5}
\]

---

Footnote: 这里记号 $O(\cdot,\cdot)$ 意思是存在(不依赖于 $N$ 和 $\varepsilon$ 的)正常数 $C_1$, $C_2$, 使得

\[
\Biggl|O\biggl(\frac{1}{N^{1/2}\varepsilon^{1+d/4}},\varepsilon^{1/2}\biggr)\Biggr|\leqslant C_1\frac{1}{N^{1/2}\varepsilon^{1+d/4}}+C_2\varepsilon^{1/2},
\]

对充分大的 $N$ 和足够小的 $\varepsilon$ 成立.

---

在实际中, 我们不能取非常小的 $\varepsilon$, 即使较小的 $\varepsilon$ 可以减少偏差错误(bias error) (1.4), 因为 $\frac{1}{N^{1/2}\varepsilon^{1+d/4}}$ 关于 $\varepsilon$ 发散. 据 [1] 推测，收敛率(convergence rate) $\frac{1}{N^{1/2}\varepsilon^{1+d/4}}$ 不太可能被提高, 因为在非参数回归(nonparametric regression)中估计二阶导数的速率是相同的.

本文我们证明方差错误(variance error)为 $O(\frac{1}{N^{1/2}\varepsilon^{1/2+d/4}}$, 这对于 (1.5) 中的收敛率来说提高了一个渐近因子 $\sqrt{\varepsilon}$. 此改进源于观察到噪声项 $Wf$ 和 $Df$ 高度相关, 相关系数(correlation coefficient)是 $1-O(\varepsilon)$. 在 $\nabla_M f=0$ 的点收敛率会更快. 并且, 我们通过一个对称参数将偏差细化为 $O(\varepsilon)$ 而非 $O(\varepsilon^{1/2})$. 最后, 将两个误差项作平衡导出 $\varepsilon$ 的一个自然的优化选择,

\[
\varepsilon=\frac{C(\mathcal{M})}{N^{1/(3+d/2)}},
\]

这给出了连续 Laplace-Beltrami 和离散的图 Laplacian 逼近之间的一个最小误差.

常数 $C(\mathcal{M})$ 是关于流形 $\mathcal{M}$ 的一个函数, 其取决于该流形的几何、维数、曲率和体积, 但独立于样本点的数量 $N$. 方差估计(variance error)取决于流形的内在维数及其体积, 它可以用 $O(1/N)$ 的精度作估计[11], 而不是 $O(1/\sqrt{N})$. 然而, 偏差估计(bias error)依赖于流形的局部曲率, 而这并不知道(which is unknown in advance). 因此, 在处理现实生活中的应用程序时, 仍然需要使用一些试错实验, 以便为给定的 $N$ 找到最佳的 $\varepsilon$, 因为 $C(\mathcal{M})$ 是未知的. 当引入一个不同的 $N$, 参数 $\varepsilon$ 可根据 (1.6) 来选取.

我们的结论总结如下:

\[
\frac{1}{\varepsilon}\sum_{j=1}^{N}L_{ij}f(x_j)=\frac{1}{2}\Delta_M f(x_i)+O\biggl(\frac{1}{N^{1/2}\varepsilon^{1/2+d/4}},\varepsilon\biggr).\tag{1.7}
\] 

 

 

 

 

 

关键字:  图拉普拉斯(Graph Laplacians)

 

---

相关阅读

 

(11条消息) 一文读懂 Bias（偏差）、Error（误差）、Variance（方差）_bias偏差_Suprit的博客-CSDN博客

 

本源量子  

本源量子云-本源量子云平台 (originqc.com.cn)  [https://qcloud.originqc.com.cn/]

QPanda 2 — QPanda 文档 (qpanda-tutorial.readthedocs.io)  [https://qpanda-tutorial.readthedocs.io/zh/latest/]

GitHub - OriginQ/QPanda-2: QPanda 2 is an open source quantum computing framework developed by OriginQC that can be used to build, run, and optimize quantum algorithms.    [https://github.com/OriginQ/QPanda-2]

---

七维高科

http://www.7d-soft.com/

1. “领先世界、当今最强大、最易于使用的数值优化分析计算软件平台：最强大的全局优化算法令优化拟合不再成为难题，独特的ODE求解器轻松应对任意边值问题。”
2. “非线性拟合、非线性方程组、参数估算反演、方程求解、微分方程求解、微分方程拟合 、非线性规划、混合整数规划...，令你忘却Matlab、Origin、Lingo、Gams等世界品牌使用时的繁琐、低效与不足。”
3. “遍布世界数十万名科技工作者的选择！”

 

---

元计算

http://www.yuanjisuan.cn/

 

号称中国唯一原创工业仿真解决方案提供商

“FELAC.IDE采用元件化思想自主定制有限元计算的基本工序，使用有限元语言来书写偏微分方程及多种算法，高效、高质量地生成有限元求解器的通用有限元软件集成开发平台。”

 

---

中科大ABACUS

第一性原理量子力学计算软件ABACUS（Atomic-orbital Based Ab-initio Computation at UStc）
http://abacus.ustc.edu.cn/
http://ooe.ustc.edu.cn/news8.html

 

---

MMP

中科院

http://www.mmrc.iss.ac.cn/mmp/index.htm

(已停止研发)

---

maTHμ

清华大学

https://github.com/maTHmU

 

---

baltamatica

北太振寰（重庆）科技有限公司 (baltamatica.com)   [http://www.baltamatica.com/]

北太天元数值计算通用软件Numerical Computation Software

“北太天元数值计算通用软件提供科学计算、可视化、交互式程序设计，具备强大的底层数学函数库，支持数值计算、数据分析、数据可视化、数据优化、算法开发等工作，并通过SDK与API接口，扩展支持各类学科与行业场景，为各领域科学家与工程师提供优质、可靠的科学计算环境。”

 

---

CAE 论坛

全国工程计算软件发展论坛 (caeforum.cn)    [http://caeforum.cn/]

用来计算 $\pi$ 的楚德诺夫斯基(Chudnovsky)公式:

\[
\frac{1}{\pi}=12\sum_{k=0}^{+\infty}\frac{(-1)^k\cdot(6k)!(13591409+545140134k)}{(3k)!(k!)^3\cdot 640320^{3k+\frac{3}{2}}}.
\]

等价于

\[
\frac{(640320)^{\frac{3}{2}}}{12\pi}=\frac{426880\sqrt{10005}}{\pi}=\sum_{k=0}^{+\infty}\frac{(6k)!(545140134k+13591409)}{(3k)!(k!)^3(-262537412640768000)^k}
\]

 

\[
\frac{640320^{3/2}}{12\pi}=\frac{426880\sqrt{10005}}{\pi}=\sum^\infty_{k=0}\frac{(6k)!(545140134k+13591409)}{(3k)!(k!)^3\left(-640320\right)^{3k}}
\]

 

 

---

References:

https://en.wikipedia.org/wiki/Chudnovsky_algorithm

https://bbs.emath.ac.cn/thread-17586-1-1.html

http://numbers.computation.free.fr/Constants/PiProgram/pifast.html

http://www.numberworld.org/y-cruncher/internals/binary-splitting-library.html#pi_chudnovsky

 

编程计算 Gamma 函数的值

\[
\Gamma(s)=\int_0^{+\infty}e^{-t}t^s\cdot\frac{1}{t}dt
\]

GNU GSL

http://www.gnu.org/software/gsl/manual/html_node/Gamma-Functions.html#index-gsl_005fsf_005flngamma-583

 

 

References:

http://stackoverflow.com/questions/15472803/gamma-or-log-gamma-function-in-c-or-c

http://www.nr.com/

http://www.cplusplus.com/forum/general/93834/

Lagrange 插值多项式 (或称 Lagrange 多项式)

Lagrange polynomial (interpolation polynomial in the Lagrange form )

 

是指经过 $n$ 个点的阶小于 $n$ 的如下多项式 $P(x)$. （这 $n$ 个点是 $\{(x_i,y_i)\mid y_i=f(x_i), i=1,2\ldots, n\}$）

\[
P(x):=\sum_{j=1}^{n}P_j(x),
\]

其中

\[
P_j(x)=y_j\prod_{k=1,k\neq j}^{n}\frac{x-x_k}{x_j-x_k},
\]

确切地,

\[
\begin{split}
P(x)=&\frac{(x-x_2)(x-x_3)\cdots (x-x_n)}{(x_1-x_2)(x_1-x_3)\cdots (x_1-x_n)}y_1\\
&+\frac{(x-x_1)(x-x_3)\cdots (x-x_n)}{(x_2-x_1)(x_2-x_3)\cdots (x_2-x_n)}y_2\\
&+\cdots +\frac{(x-x_1)(x-x_2)\cdots (x-x_{n-1})}{(x_n-x_1)(x_n-x_2)\cdots (x_n-x_{n-1})}y_n.
\end{split}
\]

---

矩阵为参数的 Lagrange 多项式

设 $A$ 是可对角化矩阵, 有 $k$ 个不同的特征值 $\lambda_1,\lambda_2,\ldots,\lambda_k$. $A$ 的 $k$ 个 Frobenius covariant(协变量) 是指下面的 $k$ 个矩阵:

\[
A_i:=\prod_{j=1,j\neq i}^{k}\frac{1}{\lambda_i-\lambda_j}(A-\lambda_j I),\quad i=1,2,\ldots, k.
\]

此本质上就是 Lagrange 插值多项式, 只不过 $x$ 换成了矩阵 $A$.

 

考虑 3n+1循环中变换 $T(n)$ 的逆变换

\[ n=S(m)= \begin{cases} 2m&n\text{ is even}\\ (2m-1)/3&n\text{ is odd and }n>1. \end{cases} \] 问是否可由数字1在变换S下生成所有自然数. 一个很简单的事实是, 在S下由1生成了一棵树. 现在的问题是这棵树是否包含了所有自然数?