利用 Lanczos 方法实现张量的 HOSVD 分解

1. 特征值分解（EVD）

如果 $A$ 是一个 $m\times m$ 的 实对称矩阵（$A=A^T$） ，如果存在 $m$ 维列向量 $q$ 和实数 $\lambda$ 满足 $Aq=\lambda q$，则称 $q$ 为 $A$ 的 特征向量、$\lambda$ 为 $A$ 的 特征值。

如果我们求出 $A$ 的 $m$ 个特征值与特征向量，并将所有特征向量标准化，则可以得到标准正交矩阵 $Q=[q_1,q_2,\cdots ,q_m],Q\in {\mathbb{R}}^{m\times m}$，与对角矩阵 $\Sigma =diag[{\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_m]$，他们满足如下形式：
$$A=Q\Sigma Q^T=Q\left[\begin{array}{cccc}{\lambda }_1& \cdots & \cdots & \cdots \\ \cdots & {\lambda }_2& \cdots & \cdots \\ \cdots & \cdots & \ddots & \cdots \\ \cdots & \cdots & \cdots & {\lambda }_m\end{array}\right]Q^T$$

numpy 中对一个矩阵进行特征值分解使用 np.linalg.eig 函数，从下例可以看出，该函数求解得到的特征值大小是随机排列的，如果需要前 $n$ 个大的特征值与特征向量，需要进一步进行处理：

#创建一个对称矩阵
sym_matrix = np.array([[1,2,3,4],
                      [2,5,6,7],
                      [3,6,8,9],
                      [4,7,9,10]])
#对矩阵求解特征值(se)与特征向量(q)
se, q = np.linalg.eig(sym_matrix)    
print('se:')
print(se)
print('q:')
print(q)
"""
se:
[24.06253512 -0.80548492  0.5580362   0.18491359]
q:
[[-0.22593827 -0.72531367 -0.56047033 -0.32976505]
 [-0.44322186 -0.31846973  0.77633831 -0.3153256 ]
 [-0.57278788 -0.14246073 -0.05844028  0.805111  ]
 [-0.6514755   0.59346614 -0.28241204 -0.37897369]]
"""

在特征值分解中，对要分解的矩阵有着较高的要求，它需要被分解的矩阵 $A$ 为实对称矩阵，但是现实中我们遇到的矩阵并不都是实对称的，对于一个 $m\times n$ 的矩阵，如果我们想分解为如上的形式，我们需要使用 SVD分解。

2. 奇异值分解（SVD）

对于一个 $m\times n$ 的实数矩阵 $A$，我们定义其奇异值分解（SVD）的形式为：
$$A=U\Sigma V^T$$

其中 $U\in {\mathbb{R}}^{m\times m}$ 和 $V\in {\mathbb{R}}^{n\times n}$ 为单位正交矩阵，即有 $U{U}^T=I$ 和 $V{V}^T=I$， $U$ 为左奇异矩阵，$V$ 为右奇异矩阵，$\Sigma \in {\mathbb{R}}^{m\times n}$ 仅在主对角线上有值，我们称它们为奇异值。分解示意图为：

奇异值分解的求解通过下式：
$$A{A}^T=U\Sigma V^{T}V{\Sigma }^T{U}^T=U\Sigma {\Sigma }^T{U}^T$$

$$A^{T}A=V{\Sigma }^T{U}^{T}U\Sigma V^T=V{\Sigma }^T\Sigma V^T$$

因为 $A{A}^T$ 和 $A{A}^T$ 是对称矩阵，则通过对上面两个等式求它们的特征值和特征向量，就可以求出 $U,V,\Sigma$ 。

我们可以使用 numpy 中的 np.linalg.svd 函数对矩阵求解 SVD ：

#创建一个5×3的矩阵
T = np.array([[16,25,35],
              [45,23,6],
              [7,12,9],
              [10,18,12],
              [15,14,15]])
#进行SVD分解      
u,ss,vh = np.linalg.svd(T, full_matrices=False)

print('u:')
print(u)

print('ss:')
print(ss)

print('vh')
print(vh)

#输出 u * ss * vh
print('u*ss*vh:')
print(u.dot(np.diag(ss)).dot(vh))
"""
u:
[[-0.57781685 -0.62450998  0.36719022]
 [-0.63112018  0.74837707  0.09682015]
 [-0.21931497 -0.12213601 -0.43502282]
 [-0.31447256 -0.15727956 -0.7906354 ]
 [-0.34759596 -0.10131625  0.20358785]]
ss:
[72.50580717 29.74834802  6.47639687]
vh
[[-0.65566262 -0.58091643 -0.48233041]
 [ 0.66347606 -0.13833268 -0.73529829]
 [ 0.3604248  -0.80212229  0.47612372]]
u*ss*vh:
[[16. 25. 35.]
 [45. 23.  6.]
 [ 7. 12.  9.]
 [10. 18. 12.]
 [15. 14. 15.]]
"""

3. QR 分解

一个矩阵 $A\in {\mathbb{R}}^{m\times n}(m\le n)$ 可以被分解为 $A=QR$ 的形式，其中：

• $Q\in {\mathbb{R}}^{m\times m}$ 是一个正交矩阵，即 $Q^{T}Q=Q{Q}^T=I$
• $R=\left[\begin{array}{c}\hat{R}\\ 0\end{array}\right]\in {\mathbb{R}}^{m\times n}$，$\hat{R}\in {\mathbb{R}}^{m\times m}$是上三角矩阵

numpy 中，进行 QR 分解使用 np.linalg.qr 函数：

A = np.array([[16,25,35],
              [45,23,6],
              [7,12,9],
              [10,18,12],
              [15,14,15]])
q,r = np.linalg.qr(A)
print('q:')
print(q)

print('r:')
print(r)

print('q*r')
print(q.dot(r)) 
"""
q:
[[-0.31051867  0.63947628  0.59444254]
 [-0.87333376 -0.44331977 -0.09225387]
 [-0.13585192  0.33011618 -0.35331103]
 [-0.19407417  0.51220924 -0.66985902]
 [-0.29111125  0.15232425  0.25414072]]
r:
[[-51.52669211 -37.04875904 -24.02638224]
 [  0.          21.10425203  31.12417123]
 [  0.           0.          12.84596908]]
q*r
[[16. 25. 35.]
 [45. 23.  6.]
 [ 7. 12.  9.]
 [10. 18. 12.]
 [15. 14. 15.]]
"""

4. Lanczos 方法

Lanczos 方法是求解大规模稀疏矩阵特征值问题最常用的方法之一，其主要思想是将对称矩阵三对角化，通过对三对角矩阵进行特征分解，从而得到原对称矩阵的特征值与特征向量。

由于三对角矩阵是稀疏矩阵，因此对其进行特征值求解的时间复杂度相对是很低的。

4.1 Lanczos 方法思路

设矩阵 $A\in {\mathbb{R}}^{n\times n}$，如果 $A=A^T$，那么一定存在 正交矩阵 $Q_n$ ，使得 $\mathbf{T}={\mathbf{Q}}_{\mathbf{n}}^{\mathbf{T}}\mathbf{A}{\mathbf{Q}}_{\mathbf{n}}$，其中 $T$ 是对称三对角矩阵。

可以将以上分解中的 $Q$ 写为：
$$Q=[q_1,q_2,\cdots ,q_n]$$

其中 $q_i$ 为 $Q$ 的列向量。

将 $T$ 记为：
$$T=\left[\begin{array}{ccccccccc}{\alpha }_1& {\beta }_1& 0& 0& 0& 0& \cdots & 0& 0\\ {\beta }_1& {\alpha }_2& {\beta }_2& 0& 0& 0& \cdots & 0& 0\\ 0& {\beta }_2& {\alpha }_3& {\beta }_3& 0& 0& \cdots & 0& 0\\ 0& 0& {\beta }_3& {\alpha }_4& {\beta }_4& 0& \cdots & 0& 0\\ 0& 0& 0& {\beta }_4& {\alpha }_5& {\beta }_5& \cdots & 0& 0\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 0& 0& 0& 0& 0& 0& \cdots & {\alpha }_{n-1}& {\beta }_{n-1}\\ 0& 0& 0& 0& 0& 0& \cdots & {\beta }_{n-1}& {\alpha }_n\end{array}\right]$$

从上面我们可以知道，对于一个矩阵 $A$ ，我们如果想要求其三对角矩阵，要求的是 $T$ 中的 $\alpha 、\beta$ 参数，以及 $Q$ 矩阵，因此通过如下的变形，我们求其中的变量。

由于 $Q$ 是正交的，我们可以将求解三对角矩阵的式子写为：
$$AQ=QT$$

对左边的式子，我们取 $Q$ 的第一列与 $A$ 作乘积，得到左边结果的第一列为 $A{q}_1$。同时我们取右式 $T$ 的第一列与 $Q$ 作乘积，得到 ${\alpha }_1{q}_1+{\beta }_1{q}_2$ ，即
$$A{q}_1={\alpha }_1{q}_1+{\beta }_1{q}_2$$

同理，有：
$$A{q}_2={\beta }_1{q}_1+{\alpha }_2{q}_2+{\beta }_2{q}_3$$

$$A{q}_3={\beta }_2{q}_2+{\alpha }_3{q}_3+{\beta }_3{q}_4$$

由此可以得出：
$$A{q}_i={\beta }_{i-1}{q}_{i-1}+{\alpha }_i{q}_i+{\beta }_i{q}_{i+1}，i=1,2,3,\cdots ,n$$

由于 ${\beta }_0,{\beta }_n,q_0,q_{n+1}$ 都是不存在的，所以我们令 ${\beta }_0{q}_0={\beta }_n{q}_{n+1}=0$。

我们将上式两边同时左乘以 $q_i^T$，由于 $Q$ 是正交的，所以其列向量互相正交，且 $q_i^T{q}_i=1$，我们可以得到：
$$q_i^{T}A{q}_i={\alpha }_i$$

通过上式可以看出，如果我们知道了 $q_i$ ，就能求出 ${\alpha }_i$，这里就是求解 $\mathbf{\alpha }$ 的方法。

我们再对等式两边移项可以得到：
$${\beta }_i{q}_{i+1}=A{q}_i-{\beta }_{i-1}{q}_{i-1}-{\alpha }_i{q}_i$$

上式左侧是一个向量，由于 $Q$ 是正交矩阵，所以 $q_{i+1}$ 是一个归一化的单位向量，我们令
$$r_i={\beta }_i{q}_{i+1}=A{q}_i-{\beta }_{i-1}{q}_{i-1}-{\alpha }_i{q}_i$$

有 ${\beta }_i=\Vert {\beta }_i{q}_{i+1}{\Vert }_2$ ，即

$${\beta }_i=\Vert r_i{\Vert }_2=\Vert {\beta }_i{q}_{i+1}{\Vert }_2=\Vert A{q}_i-{\beta }_{i-1}{q}_{i-1}-{\alpha }_i{q}_i{\Vert }_2$$

从上式可以看出，如果我们知道了 ${\alpha }_i、q_i、{\beta }_{i-1}、q_{i-1}$ 就能求出 ${\beta }_i$，又因为 $q$ 都是单位向量，所以我们还能通过 $\Vert {\beta }_i{q}_{i+1}{\Vert }_2$求出 $q_{i+1}$，这样就能循环求解了。

我们从初始条件看，如果我们知道了 $q_1$ 的值，就能根据上面求 $\mathbf{\alpha }$ 的方法求出 ${\alpha }_1$，因为我们假设 ${\beta }_0{q}_0=0$，所以我们就能通过 ${\beta }_0{q}_0、{\alpha }_1、q_1$ 求出 ${\beta }_1$，再通过令 $q_{i+1}=\frac{r_i}{\Vert r_i{\Vert }_2}$，这时 $q_{i+1}$ 就是一个单位向量，我们就得到了下一轮循环的 $q$。因此所有的重点就在 $q_1$ ，在实际求解中，我们可以随机初始化一个单位向量作为 $q_1$，由此我们可以得到计算 $A$ 的三对角矩阵的流程：

1. 随机初始化一个单位向量 $q_1$
2. 计算 ${\alpha }_1=q_1^{T}A{q}_1$
3. $r_1=A{q}_1-{\alpha }_1{q}_1-{\beta }_0{q}_0$
4. ${\beta }_1=\Vert r_1{\Vert }_2$
5. $q_2=\frac{r_1}{{\beta }_1}$
6. 重复以上步骤，依次求出 { q 1 , q 2 , ⋯   , q n } \{q_1,q_2,\cdots,q_n\} {q1​,q2​,⋯,qn​}、${\alpha }_1,{\beta }_1,{\alpha }_2,{\beta }_2,\cdots ,{\alpha }_n$

此即 Lanczos 迭代，其中 $q_i$ 为 Lanczos 向量，若其中某一步 ${\beta }_i=0$ ，则此时得到的 $T_j$ 的特征值都将是 $A$ 的特征值。

到这里，我们就得到了 $T=Q^{T}AQ$，如果此时对 $T$ 进行特征分解得到 $T=WS{W}^T$，所以有
$$A=(QW)S(W^T{Q}^T)$$

我们得到了通过将矩阵 $A$ 变换为三对角矩阵，从而求解特征值与特征向量的方法。将矩阵转换为三对角矩阵的 python 实现为：

#生成一个对称矩阵
sym_matrix = np.array([[1,2,3,4],
                      [2,5,6,7],
                      [3,6,8,9],
                      [4,7,9,10]])
Q = []
A = []
B = []
sum_dim = sym_matrix.shape[0]
    
#随机初始化q1
q1 = np.random.random((sum_dim,1))
q1 = q1/np.linalg.norm(q1)
Q.append(q1)

#计算a0
a0 = (q1.transpose()[0]).dot(sym_matrix).dot(q1)[0]
A.append(a0)

#计算r0
r = sym_matrix.dot(q1) - a0*q1

#计算b0
b0 = np.linalg.norm(r)
B.append(b0)

#循环迭代，得到A、B和Q
for i in range(1,sum_dim):
    qi = r/B[i-1]
    Q.append(qi)

    ai = qi.transpose()[0].dot(sym_matrix).dot(qi)[0]
    A.append(ai)

    r = sym_matrix.dot(qi) - ai*qi - B[i-1]*Q[i-1]
    if i < sum_dim - 1:
        b_next = np.linalg.norm(r)
        B.append(b_next)

#组成三对角矩阵
L = np.diag(A) + np.diag(B,1) + np.diag(B,-1)

print(L)
"""
[[10.99996429  6.7727395   0.          0.          0.        ]
 [ 6.7727395  11.4584147   5.03926679  0.          0.        ]
 [ 0.          5.03926679 -6.36314172  4.98605503  0.        ]
 [ 0.          0.          4.98605503  1.00512957  1.8902924 ]
 [ 0.          0.          0.          1.8902924   0.89963316]]
"""

4.2 三对角矩阵的特征值分解

三对角化后的矩阵 $T$ 有三种常用的方法进行特征分解：

1. QR 法，即对 $T$ 进行 QR 分解，令 $T=QR$，再令 $T^{\prime }=RQ$，如此反复，直至收敛。其中产生的 $Q$ 矩阵的连乘得到的结果就是 $T$ 的特征矩阵。
2. 分治法。这个方法将 $T$ 划均匀分成两块三对角阵和一个秩为 1 的 $2\times 2$ 的块状方阵的和。分别对两个三对角阵分解，然后合并。合并时对特征向量作秩1的修正。
3. MRRR 方法。这个方法理论上是最快的方式，复杂度为 $O(l^2)$，主要采用逆迭代来加快收敛。

4.3 利用 QR 算法求解特征值

给定一个实对称矩阵 $X$，假设其通过特征值分解得到 $X=PS{P}^T$，其中 $P$ 为正交阵，$S$ 是对角阵，利用 QR 算法求解其特征值分解的算法为：

$X_1=X$
$P=E$
$fork=1,2,\cdots$
$X_k=Q_k\times R_k$
$X_{k+1}=R_k\times Q_k$
$P=P\times Q_k$

这个算法很简单，每次迭代只需要做一次 QR 分解和一个矩阵乘法即可，用 python 实现为：

#创建对称矩阵
sym_matrix = np.array([[1,2,3,4],
                      [2,5,6,7],
                      [3,6,8,9],
                      [4,7,9,10]])
dim = sym_matrix.shape
#初始化P为单位矩阵
P = np.eye(dim[0])
X_n = sym_matrix
i = 1
while i<20:   #迭代次数可以根据需要改变
    i += 1
    (Q, R) = np.linalg.qr(X_n)  #QR分解
    X_n = R.dot(Q)
    P = P.dot(Q)

print("QR分解求解特征矩阵")
print(P)

se, q = np.linalg.eig(sym_matrix)
print("特征值函数求解特征矩阵：")
print(q)

"""
QR分解求解特征矩阵
[[-0.22593827 -0.72442125 -0.56162332 -0.32976505]
 [-0.44322186 -0.31970419  0.77583077 -0.3153256 ]
 [-0.57278788 -0.1423676  -0.05866681  0.805111  ]
 [-0.6514755   0.5939146  -0.28146771 -0.37897369]]
特征值函数求解特征矩阵：
[[-0.22593827 -0.72531367 -0.56047033 -0.32976505]
 [-0.44322186 -0.31846973  0.77633831 -0.3153256 ]
 [-0.57278788 -0.14246073 -0.05844028  0.805111  ]
 [-0.6514755   0.59346614 -0.28241204 -0.37897369]]
"""