Arnoldi Iteration 思考

1. 投影平面

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

• 计算向量内积,向量v在向量q上的投影长度|p|
  v T q = ∣ v ∣ ∣ q ∣ cos ⁡ θ = ∣ p ∣ ∣ q ∣ → ∣ p ∣ = v T q ∣ q ∣
  vTq=|v||q|cos⁡θ=|p||q|→|p|=vTq|q|" role="presentation" style="position: relative;">vTq=|v||q|cosθ=|p||q|→|p|=vTq|q|$$v^{T}q=|v||q|\cos \theta =|p||q|\to |p|=\frac{v^{T}q}{|q|}$$
  vTq=∣v∣∣q∣cosθ=∣p∣∣q∣→∣p∣=∣q∣vTq​​​
• 我们知道，q方向上的单位向量为$\frac{q}{|q|}$,那么投影向量p可得,$v^{T}q$为标量，随便放位置
  p = ∣ p ∣ ⋅ q ∣ q ∣ = v T q ∣ q ∣ ⋅ q ∣ q ∣ = v T q q T q q
  p=|p|⋅q|q|=vTq|q|⋅q|q|=vTqqTqq" role="presentation" style="position: relative;">p=|p|⋅q|q|=vTq|q|⋅q|q|=vTqqTqq$$p=|p|\cdot \frac{q}{|q|}=\frac{v^{T}q}{|q|}\cdot \frac{q}{|q|}=\frac{v^{T}q}{q^{T}q}q$$
  p=∣p∣⋅∣q∣q​=∣q∣vTq​⋅∣q∣q​=qTqvTq​q​​
• 重点！内积可以随便转换，并且标量位置可以随便放！
  v T q = q T v
  vTq=qTv" role="presentation" style="position: relative;">vTq=qTv$$v^{T}q=q^{T}v$$
  vTq=qTv​​
• 整理可得：
  p = q T v q T q q = q T v q q T q
  p=qTvqTqq=qTvqqTq" role="presentation" style="position: relative;">p=qTvqTqq=qTvqqTq$$p=\frac{q^{T}v}{q^{T}q}q=\frac{q^{T}vq}{q^{T}q}$$
  p=qTqqTv​q=qTqqTvq​​​
• 标量位置随意可得: $q^{T}vq\to q{q}^{T}v$
  p = q T v q q T q = q q T q T q v
  p=qTvqqTq=qqTqTqv" role="presentation" style="position: relative;">p=qTvqqTq=qqTqTqv$$p=\frac{q^{T}vq}{q^{T}q}=\frac{q{q}^T}{q^{T}q}v$$
  p=qTqqTvq​=qTqqqT​v​​
• 第一个是投影矩阵P
  P = q q T q T q , p = P v
  P=qqTqTq,p=Pv" role="presentation" style="position: relative;">P=qqTqTq,p=Pv$$P=\frac{q{q}^T}{q^{T}q},p=Pv$$
  P=qTqqqT​,p=Pv​​
• 第二，快速计算一个向量v在向量q上的投影p
  p = q T v q q T q
  p=qTvqqTq" role="presentation" style="position: relative;">p=qTvqqTq$$p=\frac{q^{T}vq}{q^{T}q}$$
  p=qTqqTvq​​​
• 第三，当q为单位向量的时候，$q^{T}q=|q{|}^2=1$,像不像二次型形式，就是这么神奇！
  p = q T v q
  p=qTvq" role="presentation" style="position: relative;">p=qTvq$$p=q^{T}vq$$
  p=qTvq​​
• 第四 ，一般情况下计算垂直向量e,向量几何关系可得v=p+e，
  e = v − p = v − q T v q q T q
  e=v−p=v−qTvqqTq" role="presentation" style="position: relative;">e=v−p=v−qTvqqTq$$e=v-p=v-\frac{q^{T}vq}{q^{T}q}$$
  e=v−p=v−qTqqTvq​​​
  第五，特殊情况下,|q|=1,整理可得：
  e = v − q T v q
  e=v−qTvq" role="presentation" style="position: relative;">e=v−qTvq$$e=v-q^{T}vq$$
  e=v−qTvq​​
  

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)