浅谈：函数逼近~曲线拟合~最小二乘法~DFT
nekko
·
·
个人记录
最小二乘法
最小化：f_0(\vec{x})=\Vert A\vec{x}-\vec{b} \Vert_2
即最小化：f(\vec{x})=\Vert A\vec{x}-\vec{b} \Vert_2^2=(A\vec{x}-\vec{b})^T(A\vec{x}-\vec{b})
即：f(\vec{x})=\vec{x}^TA^TA\vec{x}-\vec{b}^TA\vec{x}-\vec{x}^TA^T\vec{b}+\vec{b}^T\vec{b}
令偏导为 0，得到（A^TA\hat{\bold{x}}=A^T\bold{b}）：
\frac{\partial f}{\partial \vec{x}}=2A^TA\vec{x}-2A^T\vec{b}=0 \Rightarrow A^TA\vec{x}=A^T\vec{b}
矩阵求导：
\begin{cases}
\frac{\partial}{\partial \vec{x}}\vec{b}^T\vec{x}=\frac{\partial}{\partial \vec{x}}\vec{x}^T\vec{b}=\vec{b} \\ \frac{\partial}{\partial \vec{x}} \vec{x}^TA^TA\vec{x}=A^TA\vec{x}+(A^TA)^T\vec{x}=2A^TA\vec{x}
\end{cases}
是否是极小值，即海森矩阵 H(f)=2A^TA 是否正定？
考虑弱化条件：A 列线性无关
则 \forall \vec{x} \ne \vec{0},\vec{x}^T(A^TA)\vec{x}=\Vert A\vec{x}\Vert_2^2>0
所以 \hat{\bold{x}} 为极小值点
设 m 个线性无关向量 a_1,\cdots,a_m \in \mathbb{R}^n，其基矩阵 A=\begin{bmatrix} a_1 & \cdots & a_m \end{bmatrix}。
对于 f \in \mathbb{R}^n，设 w=Ax \in \text{span}(A) 满足 \Vert f-w \Vert 最小，则：
A^TAx=A^T f
特别的，当 f \in \text{span}(A) 时，线性方程组 Ax=f 有唯一解。
函数逼近
考虑最佳平方逼近。
设 f,g \in C[a,b]，定义内积 (f,g)=\int_a^b f(x)g(x)dx
考虑基函数集 \Phi=\{\varphi_0,\cdots,\varphi_n\} \subseteq C[a,b]，现在想求 f 在 \varphi=\text{span}(\Phi) 的投影
即求 S^* \in \varphi，满足 \Vert f-S^* \Vert_2^2=\inf_{S \in \varphi} \Vert f-S \Vert_2^2
即求 I(a_0,\cdots,a_n)=\int_a^b \left( \sum_{i=0}^{n}a_i\varphi_i(x)-f(x) \right)^2dx 的最小值
\begin{aligned}
&\frac{\partial I}{\partial a_k}=2\int_a^b\left( \sum_{i=0}^{n}a_i\varphi_i(x)-f(x) \right)\varphi_k(x)dx=0 \\
\Rightarrow &\sum_{i=0}^{n}a_i\int_a^b\varphi_i(x)\varphi_k(x)dx=\int_a^bf(x)\varphi_k(x)dx \\
\Rightarrow &\sum_{i=0}^{n}(\varphi_i,\varphi_k)a_i=(f,\varphi_k) \\
\Rightarrow &\hat\Phi_{nn}\vec{a}=\vec{b}
\end{aligned}
解线性方程组！
特别的，当选取的基函数 \Phi 是两两正交的时候，则有 a_i=\frac{(f,\varphi_k)}{\Vert \varphi_i \Vert}
比如说傅里叶变换（基函数 \{1,\sin x,\cos x,\cdots\} 在 [-\pi,\pi] 两两正交）：
\begin{cases}
f(x)=\frac{a_0}{2}+\sum_{n=1}^{+\infty}(a_n\cos nx+b_n\sin nx) \\
a_n=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)\cos nxdx \\
b_n=\frac{1}{\pi}\int_{-\pi}^{\pi}f(x)\sin nxdx
\end{cases}
设 f \in C[a,b]，定义 (f,g)=\int_a^bf(x)g(x)dx。考虑一组基 \phi_0,\cdots,\phi_m \in C[a,b]，找一组系数 a，使得 \Vert \Phi a-f\Vert 最小。
类似的（多元函数求极值），也可以导出最小二乘解：
\Phi^T\Phi a=\Phi^Tf
特别的，当 \Phi 中函数两两正交时，有：
a_i=\frac{(\phi_i,f)}{(\phi_i,\phi_i)}
此时有广义傅里叶级数：
f(x) \sim \sum_{k=0}^{\infty} \frac{(\phi_i,f)}{(\phi_i,\phi_i)}\phi_i
在一般区间 x \in [a,b] 上的傅里叶展开可以转换到 t \in [-1,1] 上，只需要作变换：
x=\frac{b-a}{2}t+\frac{b+a}{2} \quad (-1 \le t \le 1)
曲线拟合
设一组实验数据 (x_0,y_0),\cdots,(x_m,y_m) \in \mathbb{R}^2，令 f 为该数据的拟合曲线，有误差：
\begin{aligned}
&\delta_i=f(x_i)-y_i \\
& \Vert \delta \Vert_2^2=\min_{f} \sum_{i=0}^{m}(f(x_i)-y_i)^2
\end{aligned}
考虑一组基 \phi_0,\cdots,\phi_m \in C[a,b]，找一组系数 a，令 f(x)=\sum_{k=0}^{m}a_k\phi_k(x)。
定义 (\phi_i,\phi_j)=\sum_{k=0}^{m} \phi_i(x_k)\phi_j(x_k)，类似的可以得到：
\Phi^T\Phi a=\Phi^T y
特别的，当 \Phi 中函数关于点集 \{x_i\} 两两正交时，有：
a_i=\frac{(\phi_i,y)}{(\phi_i,\phi_i)}
傅里叶变换
最佳平方三角逼近
考虑一组正交基：
1,\cos x,\sin x,\cos 2x,\sin 2x,\cdots,\cos kx, \sin kx,\cdots
设函数 f \in C[0,2\pi]，将其在 [0,2\pi] 上展开：
f(x) \sim \frac{a_0}{2}+\sum_{k=1}^{\infty}(a_k\cos kx+ b_k\sin kx)
其中：
\begin{aligned}
&\frac{a_0}{2}=\frac{(f,1)}{(1,1)}=\frac{\int_0^{2\pi}f(x)dx}{\int_0^{2\pi}dx}=\frac{1}{2\pi}\int_0^{2\pi}f(x)dx\\
&a_k=\frac{(f,\cos kx)}{(\cos kx, \cos kx)}=\frac{\int_0^{2\pi}f(x)\cos kx dx}{\frac{1}{2}\int_0^{2\pi} (1+\cos(2kx))dx}=\frac{1}{\pi}\int_0^{2\pi}f(x)\cos kxdx \\
&b_k=\frac{1}{\pi}
\end{aligned}
三角插值
考虑函数族：
1,\cos x,\sin x,\cos 2x,\sin 2x,\cdots,\cos mx,\sin mx
其在点集 \{\frac{2\pi k}{2m+1} \mid k=0,1,\cdots,2m\} 上两两正交。
先考虑：
\sum_{k=0}^{2m}\sin px_k\sin qx_k=-\frac{1}{2}\sum_{k=0}^{2m}\cos(p+q)x_k-\cos (p-q) x_k
考虑到：
\begin{aligned}
\sum_{k=0}^{2m}\sin tx_k&=\sum_{k=0}^{2m}\Im(e^{itx_k})=\Im \left(\sum_{k=0}^{2m}e^{(2it\pi k)/(2m+1)}\right) \\
&=\Im \left(\frac{1-e^{2\pi i t \cdot (2m+1)/(2m+1)}}{1-e^{2\pi i t/(2m+1)}}\right)=\Im \left(\frac{1-1}{1-e^{2\pi it/(2m+1)}}\right)=0 \\
\sum_{k=0}^{2m}\cos(tx_k)&=\sum_{k=0}^{2m}\Re(e^{itx_k})=0
\end{aligned}
所以：
\sum_{k=0}^{2m}\sin px_k\sin qx_k=
\begin{cases}
\frac{2m+1}{2} \quad &p = q \ne 0 \\
0 \quad &\text{otherwise}
\end{cases}
类似的，可以证明：
\begin{aligned}
&\sum_{k=0}^{2m}\sin px_k\cos qx_k=\frac{1}{2}\sum_{k=0}^{2m}\sin(p+q)x_k+\sin(p-q)x_k=0 \\
&\sum_{k=0}^{2m}\sin px_k=\sum_{k=0}^{2m}\cos px_k=0 \\
&\sum_{k=0}^{2m}\cos px_k\cos qx_k=\frac{1}{2}\sum_{k=0}^{2m}\cos(p+q)x_k+\cos(p-q)x_k=
\begin{cases}
\frac{2m+1}{2} \quad &p=q \ne 0 \\
2m+1 \quad &p=q=0 \\
0 \quad & \text{otherwise}
\end{cases}
\end{aligned}
同理可知，偶数个均等分的离散点集上，三角函数族仍满足正交性。
更一般的情况，复函数族：
1,e^{ix},e^{2ix},\cdots,e^{(n-1)ix}
在点集 \{\frac{2\pi k}{n} \mid k=0,1,\cdots,n-1\} 上（广义）两两正交。
一方面：
\begin{aligned}
e^{pix}e^{qix}
&=(\cos px+i\sin px)(\cos qx+i\sin qx) \\
&=\cos px \cos qx-\sin px \sin qx+i\sin px \cos qx+i\sin qx\cos px
\end{aligned}
另一方面：
\begin{aligned}
\sum_{k=0}^{n-1}e^{pix_k}e^{qix_k}
&=\sum_{k=0}^{n-1}e^{i(p+q)2\pi k/n} \\
&=
\begin{cases}
\frac{1-e^{i(p+q)2\pi}}{1-e^{i(p+q)2\pi /n}} \quad &p+q \ne 0 \\
n \quad & p+q =0
\end{cases} \\
&=
\begin{cases}
0 \quad &p+q \ne 0 \\
n \quad & p+q =0
\end{cases}
\end{aligned}
即该函数族在内积 (f,g)=\sum_{k=0}^{n-1}f(x_k)g(-x_k) 上具有正交性。
同时可以看到上式蕴含单位根反演系数：
[n|m]=\frac{1}{n}\sum_{k=0}^{n-1}e^{2\pi ikm/n}
可以得到三角插值（离散傅里叶变换，DFT）：
\begin{aligned}
&f(x) = \sum_{k=0}^{n-1}a_ke^{ikx} \\
&a_k=\frac{(y,e^{ikx})}{(e^{ikx},e^{ikx})}=\frac{1}{n}\sum_{j=0}^{n-1}y_je^{-ikx_j}=\frac{1}{n}\sum_{j=0}^{n-1}y_je^{-2\pi ikj/n} \\
&y_j=\sum_{k=0}^{n-1}a_ke^{ikx_j}=\sum_{k=0}^{n-1}a_ke^{2\pi ikj /n}
\end{aligned}
设有一多项式 A(x)=\sum_{k=0}^{n-1}a_kx^k，则：
- 通过DFT可以得到点值表示：(e^{ikx_j},y_j)
- 通过IDFT可以得到系数表示：A(x)=\sum_{k=0}^{n-1}a_kx^k
因为 kj \in \mathbb{Z}，因此 e^{2\pi ikj/n} 只有 n 个取值。
设 n 次单位根 w_n=e^{2\pi i/n}，则 w_n^{kn+r}=w_n^r，只需要考虑计算 C_j=\sum_{k=0}^{n-1}x_kw^{jk}
FFT有空再写，咕咕咕
例题
51nod 1747 近似多项式
求一个 n 次多项式 f(x)=\sum_{i=0}^{n}a_ix^i，最小化：
\int_{0}^{1}(f(x)-e^x))^2 \mathrm{d} x
其中 f(x) 应可以写为：
f(x)=\sum_{i=0}^{n}a_ix^i=\sum_{i=0}^{n} (\sum_{j=0}^{\infty}p_{i,j}e^j)x^i
其中 p_{i,j} 是有理数，只需要输出 p_{i,1},p_{i,0} 对 998244353 取模的值即可
1 \le n \le 200
最佳平方逼近
\begin{bmatrix}
(\phi_0,\phi_0) & (\phi_0,\phi_1) & \cdots & (\phi_0,\phi_n) \\
(\phi_1,\phi_0) & (\phi_1,\phi_1) & \cdots & (\phi_1,\phi_n) \\
\vdots & \vdots & \ddots & \vdots \\
(\phi_n,\phi_0) & (\phi_n,\phi_1) & \cdots & (\phi_n,\phi_n)
\end{bmatrix}
\begin{bmatrix}
a_0 \\
a_1 \\
\vdots \\
a_n
\end{bmatrix}
=
\begin{bmatrix}
(\phi_0,e^x) \\
(\phi_1,e^x) \\
\vdots \\
(\phi_n,e^x)
\end{bmatrix}
其中
\begin{cases}
\int_0^{1}\phi_i^2dx=\int_0^{1}x^{2i}dx=\frac{1}{2i+1} \\
\int_0^{1}\phi_i\phi_jdx=\int_0^{1}x^{i+j}dx=\frac{1}{i+j+1} \\
\int_0^{1}\phi_ne^xdx=\int_0^{1}x^ne^xdx=\left(e\sum_{i=0}^{n} (-1)^{n-i} \frac{n!}{i!}\right)+(-1)^{n+1}n!
\end{cases}
即 (\varphi_n,e^x)=b_ne+c_n，实际上 p_{i,j} 也只有 j=0,1。
当然也可以不用计算上面那个定积分，由分部积分可得 \int x^ne^xdx=e^xA(x)+C，其中 A(x) 是某 n 次多项式。
求导得：e^x(A+A')=x^ne^x，即 A+A'=x^n
设 A(x)=\sum_{i=0}^{n}a_ix^i，比较系数得：
\begin{cases}
a_n=1 \\
a_{i}+(i+1)a_{i+1}=0
\end{cases}
即 a_i=(-1)^{n-i} n^{\underline {n-i}}
回到 (\varphi_n,e^x)，可以发现：
\begin{cases}
b_n=1-nb_{n-1} \\
c_n=-nc_{n-1}
\end{cases}
因此 \hat\Phi_{nn}\vec{a}=\vec{b} 的右侧就知道了。
考虑到 b=ke+r，不妨设 a=pe+q，那么只需要解两次线性方程组即可。