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 np

def arnoldi_iteration(A, b, k):
    """
    Perform Arnoldi iteration to generate an orthonormal basis for the Krylov subspace.

    Parameters:
    A : numpy.ndarray
        The input matrix (n x n).
    b : numpy.ndarray
        The initial vector (n, ).
    k : int
        The number of iterations, which defines the size of the Krylov subspace.

    Returns:
    Q : numpy.ndarray
        The orthonormal basis for the Krylov subspace (n x (k+1)).
    H : numpy.ndarray
        The Hessenberg matrix (k+1 x k).
    """
    
    n = A.shape[0]
    Q = np.zeros((n, k + 1))  # Orthonormal basis
    H = np.zeros((k + 1, k))  # Hessenberg matrix

    # Normalize the initial vector
    Q[:, 0] = b / np.linalg.norm(b)

    for j in range(k):
        v = A @ Q[:, j]  # Matrix-vector multiplication
        for i in range(j + 1):
            H[i, j] = np.dot(Q[:, i].conj(), v)  # Project v onto the current basis vectors
            v = 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 vector
        if 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 b
    A = np.random.rand(5, 5)
    b = np.random.rand(5)
    k = 4

    Q, H = arnoldi_iteration(A, b, k)
    print("Orthonormal basis Q:\n", Q)
    print("Hessenberg matrix H:\n", H)