Arnoldi Iteration 思考

1. 投影平面

假设我们有一个向量q,我们需要关于向量q，构建一个投影平面P，使得给定任何向量v,可以通过公式$p=Pv$，快速得到向量v在投影平面P上的投影向量p.

• 计算向量内积,向量v在向量q上的投影长度|p|
  $$\begin{array}{c}{v}^{T}q=|v||q|\cos \theta =|p||q|\to |p|=\frac{v^{T}q}{|q|}\end{array}$$
• 我们知道，q方向上的单位向量为$\frac{q}{|q|}$,那么投影向量p可得,$v^{T}q$为标量，随便放位置
  $$\begin{array}{c}p=|p|\cdot \frac{q}{|q|}=\frac{v^{T}q}{|q|}\cdot \frac{q}{|q|}=\frac{v^{T}q}{q^{T}q}q\end{array}$$
• 重点！内积可以随便转换，并且标量位置可以随便放！
  $$\begin{array}{c}{v}^{T}q=q^{T}v\end{array}$$
• 整理可得：
  $$\begin{array}{c}p=\frac{q^{T}v}{q^{T}q}q=\frac{q^{T}vq}{q^{T}q}\end{array}$$
• 标量位置随意可得: $q^{T}vq\to q{q}^{T}v$
  $$\begin{array}{c}p=\frac{q^{T}vq}{q^{T}q}=\frac{q{q}^T}{q^{T}q}v\end{array}$$
• 第一个是投影矩阵P
  $$\begin{array}{c}P=\frac{q{q}^T}{q^{T}q},p=Pv\end{array}$$
• 第二，快速计算一个向量v在向量q上的投影p
  $$\begin{array}{c}p=\frac{q^{T}vq}{q^{T}q}\end{array}$$
• 第三，当q为单位向量的时候，$q^{T}q=|q{|}^2=1$,像不像二次型形式，就是这么神奇！
  $$\begin{array}{c}p=q^{T}vq\end{array}$$
• 第四 ，一般情况下计算垂直向量e,向量几何关系可得v=p+e，
  $$\begin{array}{c}e=v-p=v-\frac{q^{T}vq}{q^{T}q}\end{array}$$
  第五，特殊情况下,|q|=1,整理可得：
  $$\begin{array}{c}e=v-q^{T}vq\end{array}$$
  

2. Arnoldi Iteration

arnoldi Iteration的作用是想在原来的krylov 子空间中增加一个向量$A{q}_k$，具体思路如下图所示：

• 小结：arnoldi Iteration 本质上就是新建一个向量$v$，为了让v向量和以前已知的向量$q_1,q_2,\cdots ,q_k$垂直，通过不断迭代，将v向量减去掉所有在$q_1,q_2,\cdots ,q_k$上的投影向量$e_k$，这样最后得到的向量$q_k$就一定是垂直于$q_1,q_2,\cdots ,q_k$

3. python 代码

后续提供详细的，现在直接粘贴吧。

import numpy as npdef arnoldi_iteration(A, b, k):"""Perform Arnoldi iteration to generate an orthonormal basis for the Krylov subspace.Parameters:A : numpy.ndarrayThe input matrix (n x n).b : numpy.ndarrayThe initial vector (n, ).k : intThe number of iterations, which defines the size of the Krylov subspace.Returns:Q : numpy.ndarrayThe orthonormal basis for the Krylov subspace (n x (k+1)).H : numpy.ndarrayThe Hessenberg matrix (k+1 x k)."""n = A.shape[0]Q = np.zeros((n, k + 1))  # Orthonormal basisH = np.zeros((k + 1, k))  # Hessenberg matrix# Normalize the initial vectorQ[:, 0] = b / np.linalg.norm(b)for j in range(k):v = A @ Q[:, j]  # Matrix-vector multiplicationfor i in range(j + 1):H[i, j] = np.dot(Q[:, i].conj(), v)  # Project v onto the current basis vectorsv = v - H[i, j] * Q[:, i]  # Make v orthogonal to Q[:, i]H[j + 1, j] = np.linalg.norm(v)  # Normalize v to get the next basis vectorif H[j + 1, j] != 0 and j + 1 < k:Q[:, j + 1] = v / H[j + 1, j]return Q, H# Example usage
if __name__ == "__main__":# Define a random matrix A and a random vector bA = np.random.rand(5, 5)b = np.random.rand(5)k = 4Q, H = arnoldi_iteration(A, b, k)print("Orthonormal basis Q:\n", Q)print("Hessenberg matrix H:\n", H)