矩阵求逆算法

Jack315

有如下 \(n\times n\) 的矩阵：
\[A=\begin{pmatrix}
1+k&-2k&k&0&0&0&\dots&0&0&0&0&0&0\\
-2k&1+5k&-4k&k&0&0&\dots&0&0&0&0&0&0\\
k&-4k&1+6k&-4k&k&0&\dots&0&0&0&0&0&0\\
0&k&-4k&1+6k&-4k&k&\dots&0&0&0&0&0&0\\
0&0&k&-4k&\ddots&\ddots&\ddots&0&0&0&0&0&0\\
0&0&0&k&\ddots&\ddots&\ddots&\ddots&0&0&0&0&0\\
0&0&0&0&\ddots&\ddots&\ddots&\ddots&\ddots&0&0&0&0\\

0&0&0&0&0&\ddots&\ddots&\ddots&\ddots&k&0&0&0\\
0&0&0&0&0&0&\ddots&\ddots&\ddots&-4k&k&0&0\\
0&0&0&0&0&0&\dots&k&-4k&1+6k&-4k&k&0\\
0&0&0&0&0&0&\dots&0&k&-4k&1+6k&-4k&k\\
0&0&0&0&0&0&\dots&0&0&k&-4k&1+5k&-2k\\
0&0&0&0&0&0&\dots&0&0&0&k&-2k&1+k\\
\end{pmatrix}\]
按定义逆阵为：
\[A^{-1}=\frac{1}{|A|}
\begin{pmatrix}
A_{11}&A_{21}&\dots&A_{n1}\\
A_{12}&A_{22}&\dots&A_{n2}\\
\vdots&\vdots&\ddots&A_{n2}\\
A_{12}&A_{22}&\dots&A_{n2}\\
\end{pmatrix}\]
如果按定义求逆阵，效率比较低。
注意到矩阵的对角对称性，并且是一个稀疏矩阵，
有没有更好的算法能降低运算的时间和空间复杂度？

Jack315

【问题的由来】
在时间序列分析中，常需要进行趋势滤波，使得“残差”成为一个平稳随机过程。
有很多种进行趋势滤波的方法，其中比较常用的一种方法是用 Hodrick-Prescott Filter (HP 滤波器) 。

HP 滤波器的基本原理如下：
设时间序列 \(\begin{matrix}y_t,&t=1,2,\dots,n\end{matrix}\) 由趋势分量 \(\tau_t\) 和循环分量 \(c_t\) 组成，即：\(y_t=\tau_t+c_t\) 。
HP滤波器通过最小化下列目标函数而将趋势分量和循环分量分离：
\[h=\lim_{\tau_1,\tau_2,\dots\tau_n}\bigg\{\sum_{t=1}^n(y_t-\tau_t)^2-\lambda\sum_{t=2}^{n-1}(\tau_{t-1}-2\tau_t+\tau_{t+1})\bigg\}\]
其中：\(\lambda\) 为一正平滑参数，用于对趋势中的变异进行惩罚。

最小化目标函数的解为：
\[\hat\tau_t=(I+\lambda F)^{-1}y_t=(I+\lambda D^TD)^{-1}y_t\]
其中：
\(I\) 为 \(n\times n\) 单位阵：
\[I=\begin{pmatrix}
1&0&\dots&0&0\\
0&1&\dots&0&0\\
\vdots&\vdots&\ddots&\vdots&\vdots\\
0&0&\dots&1&0\\
0&0&\dots&0&1
\end{pmatrix}\]
\(D\) 为 \((n-2)\times n\) 矩阵：
\[D=\begin{pmatrix}
1&-2&1&0&0&\dots&0&0&0&0&0\\
0&1&-2&1&0&\dots&0&0&0&0&0\\
0&0&1&-2&1&\dots&0&0&0&0&0\\

0&0&0&1&\ddots&\ddots&0&0&0&0&0\\
0&0&0&0&\ddots&\ddots&\ddots&0&0&0&0\\

0&0&0&0&0&\ddots&\ddots&1&0&0&0\\
0&0&0&0&0&\dots&1&-2&1&0&0\\
0&0&0&0&0&\dots&0&1&-2&1&0\\
0&0&0&0&0&\dots&0&0&1&-2&1\\
\end{pmatrix}\]
\(F=D^TD\) 计算的结果为 \(n\times n\) 矩阵：
\[F=\begin{pmatrix}
1&-2&1&0&0&0&\dots&0&0&0&0&0&0\\
-2&5&-4&1&0&0&\dots&0&0&0&0&0&0\\
1&-4&6&-4&1&0&\dots&0&0&0&0&0&0\\
0&1&-4&6&-4&1&\dots&0&0&0&0&0&0\\
0&0&1&-4&\ddots&\ddots&\ddots&0&0&0&0&0&0\\
0&0&0&1&\ddots&\ddots&\ddots&\ddots&0&0&0&0&0\\
0&0&0&0&\ddots&\ddots&\ddots&\ddots&\ddots&0&0&0&0\\

0&0&0&0&0&\ddots&\ddots&\ddots&\ddots&1&0&0&0\\
0&0&0&0&0&0&\ddots&\ddots&\ddots&-4&1&0&0\\
0&0&0&0&0&0&\dots&1&-4&6&-4&1&0\\
0&0&0&0&0&0&\dots&0&1&-4&6&-4&1\\
0&0&0&0&0&0&\dots&0&0&1&-4&5&-2\\
0&0&0&0&0&0&\dots&0&0&0&1&-2&1\\
\end{pmatrix}\]
想用VB、C或C++自己编程实现 HP 滤波器，这样就有了 1# 中的问题。
对于几百个数据点以下的情况，用定义求逆阵没什么问题，
但当 \(n\) 增长到上千以后算法的时间和空间复杂度就变得至关重要。

这个问题的实质是一个优化（求最小化参数）问题。
如果从这个角度出发有好的算法也希望能推荐。

Jack315

有大佬发文说 HP 滤波器不合适：
《Why You Should Never Use the Hodrick-Prescott Filter》
HPFilter.part1.rar (250 KB, 下载次数: 0)

2023-10-8 07:13 上传

点击文件名下载附件

HPFilter.part2.rar (190.78 KB, 下载次数: 0)

2023-10-8 07:13 上传

点击文件名下载附件

但现在不管这个，想把各类趋势滤波器都玩一下 。

mathe

三对角阵可以通过追赶法来解决。
你这里的五对角阵也可以用类似方法去解决，我们假设
\(\begin{pmatrix}a_{11}&a_{12}&a_{13}&0&0&0&0&\cdots&0&0&0\\
a_{21}&a_{22}&a_{23}&a_{24}&0&0&0&\cdots&0&0&0\\
a_{31}&a_{32}&a_{33}&a_{34}&a_{35}&0&0&\cdots&0&0&0\\
0&a_{42}&a_{43}&a_{44}&a_{45}&a_{46}&0&\cdots&0&0&0\\
\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots\\
0&0&0&0&0&0&0&\cdots&a_{n,n-2}&a_{n,n-1}&a_{n,n}
\end{pmatrix}=LU\)
其中
L=\(\begin{pmatrix}1&0&0&0&0&0&0&\cdots&0&0&0\\
l_{21}&1&0&0&0&0&0&\cdots&0&0&0\\
l_{31}&l_{32}&1&0&0&0&0&\cdots&0&0&0\\
0&l_{42}&a_{43}&1&0&0&0&\cdots&0&0&0\\
\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots\\
0&0&0&0&0&0&0&\cdots&l_{n,n-2}&l_{n,n-1}&1
\end{pmatrix}\)

U=\(\begin{pmatrix}u_{11}&u_{12}&u_{13}&0&0&0&0&\cdots&0&0&0\\
0&u_{22}&u_{23}&u_{24}&0&0&0&\cdots&0&0&0\\
0&0&u_{33}&u_{34}&u_{35}&0&0&\cdots&0&0&0\\
0&0&0&u_{44}&u_{45}&u_{46}&0&\cdots&0&0&0\\
\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots&\cdots\\
0&0&0&0&0&0&0&\cdots&0&0&u_{n,n}
\end{pmatrix}\)
得出连个矩阵L,U后就容易计算了

补充一下，好像大家看不懂
第一步我们需要将U第一列乘上L,结果要求等于A的第一列，得到方程组
\(\begin{cases}u_{11}=a_{11}\\ u_{11}l_{21}=a_{21}\\u_{11}l_{31}=a_{31}\end{cases}\)
由此可以解得U和L的第一列\(u_{11},l_{21},l_{31}\)
第二步将U的第二列乘上L,要求等于A的第二列，得到方程组
\(\begin{cases}u_{12}=a_{22}\\ u_{22}=a_{22}-u_{12}l_{21}\\u_{22}l_{32}=a_{32}-l_{31}u_{12}\\u_{22}l_{42}=a_{42}\end{cases}\)
由此得到U和L的第二列\(u_{12},u_{22},l_{32},l_{42}\)
...
这样经过n步，每步最多5个方程5个变量，解得L,U中所有非零参数。
有了LU分解形式后，后面无论求逆还是解方程等都很简单了

nyy

一切交给软件，这样求解最简单

ShuXueZhenMiHu

A和F也不一样啊。 另外公式显示的很奇怪，根本没法看哪。

ShuXueZhenMiHu

F矩阵不可逆啊，因为对10行10列的F求阶梯型矩阵以后，我们有

Jack315

ShuXueZhenMiHu 发表于 2023-10-13 10:52
A和F也不一样啊。 另外公式显示的很奇怪，根本没法看哪。

运算是关于求矩阵 A 的逆：
\[A=(I+\lambda F)\]

因为矩阵 D （参考 2#）是二阶差分形成的矩阵，因此相关矩阵关于两条对角线都对称。
整个矩阵只有左上角和右下角的两个 2 x 2 的矩阵元素略有差异，其它元素都遵循相同的模式。
如果把这两个对角的 2 x 2 矩阵划掉，剩下的余子式是一个标准的 TOEPLITZ 矩阵。
其快速算法的运算复杂度大致为 \(O(n^2 )\) 。
但比 TOEPLITZ 矩阵更进一步的是矩阵 A 是关于两条对角线都对称的。
因此猜想可能有更好的算法。

A 的逆阵有下列特点：
1) 关于两条对角线对称；
2) 每一行所有元素的和为 1 ；
3) （由于对称性），每一列所有元素的和为 1 。

一般的软件使用 LU 分解，其运算复杂度为 \(O(n^3 )\) 。
由于前述的对称性，希望能找到更佳的算法。

Matlab 中有 hpfilter 函数，运算速度比较快。
但无法用来生成 C++ 代码，也就无法一窥其具体的算法。

如果能找出这个算法，并使用稀疏矩阵，估计空间和运算复杂度都能大幅降低。
对于运算复杂度，应该能优于 \(O(n^2 )\)，非常可能是 \(O(n\log_2 n )\) 甚至更优。

这是自己在 Matlab 里编的 hpfilter 函数：

1. function t = trend(y, k)
   
2. % trend 函数用平滑系数 k 从时间序列 y 中提取趋势分量 t 。
   
3.    
   
4.     % 时间序列长度。
   
5.     n = length(y);
   
6.    
   
7.     switch n
   
8.         case {1, 2}
   
9.             t = y;
   
10.         case 3
    
11.             t = [5 * y(1) + 2 * y(2) - y(3); ...
    
12.                 2 * (y(1) + y(2) + y(3)); ...
    
13.                 -y(1) + 2 * y(2) + 5 * y(3)
    
14.                 ] / 6;
    
15.         case 4
    
16.             t = [7 * y(1) + 4 * y(2) + y(3) - 2 * y(4); ...
    
17.                 4 * y(1) + 3 * y(2) + 2 * y(3) + y(4); ...
    
18.                 y(1) + 2 * y(2) + 3 * y(3) + 4 * y(4); ...
    
19.                 -2 * y(1) + y(2) + 4 * y(3) + 7 * y(4)
    
20.                 ] / 10;
    
21.         otherwise
    
22.             % 生成系数矩阵。
    
23.             x = zeros(n);
    
24.             for i = 1 : n
    
25.                 switch i
    
26.                     case 1
    
27.                         x(i, 1 : 3) = [k + 1; -2 * k; k];
    
28.                     case 2
    
29.                         x(i, 1 : 4) = [-2 * k; 5 * k + 1; -4 * k; k];
    
30.                     case n - 1
    
31.                         x(i, n - 3 : n) = [k; -4 * k; 5 * k + 1; -2 * k];
    
32.                     case n
    
33.                         x(i, n - 2 : n) = [k; -2 * k; k + 1];
    
34.                     otherwise
    
35.                         x(i, i - 2 : i + 2) = [k; -4 * k; 6 * k + 1; -4 * k; k];
    
36.                 end
    
37.             end
    
38. 
    
39.             % 求逆阵并还原为全元素矩阵。
    
40.             x = full(inv(x));
    
41. 
    
42.             % 计算趋势势分量。
    
43.             t = x * y;
    
44.     end
    
45. end

复制代码

其运算速度和占用空间都不如 Matlab 自带的 hpfilter 函数。

Matlab 中 hpfilter 函数的参考：
hpfilter

Jack315

ShuXueZhenMiHu 发表于 2023-10-13 11:21
F矩阵不可逆啊，因为对10行10列的F求阶梯型矩阵以后，我们有

给个 Excel 的例子。
分别为 n=10 (偶数) 和 n=9 (奇数) 两种情况。
两个表格显示了：
1) 运算过程；
2) 非零元素的排列；
3) 矩阵关于两条对角线的对称性。
HP滤波器.rar (17.5 KB, 下载次数: 1)

2023-10-14 07:52 上传

点击文件名下载附件

Jack315

mathe 发表于 2023-10-8 21:56
三对角阵可以通过追赶法来解决。
你这里的五对角阵也可以用类似方法去解决，我们假设
\(\begin{pmatrix}a_{ ...

“追赶法”搞明白了，是个优秀的算法。准备试一下。
再次感谢前辈指导！