列更新后的Krylov子空间Arnoldi分解高效更新方法

问题描述 #

已知通过Arnoldi算法得到矩阵A的Krylov子空间分解结果：

[V, H] = Arnoldi(A, v1, m)

在不存储原矩阵A的前提下，当给A的第k列添加向量dAk（即矩阵更新为 (A_{\text{new}} = A_{\text{old}} + dAk \cdot e_k^T)，其中(e_k)是第k个标准单位向量），如何直接更新已有的正交基V和上Hessenberg矩阵H，无需重新运行完整的Arnoldi算法？

原示例代码（带中文注释） #

function arnoldi_demo
    n = 5;
    A = randn(n);

    b = randn(n,1);
    m = 5;
    v1 = randn(n,1);

    % 初始Arnoldi分解
    [V, H] = arnoldi(A, v1, m);

    % 对矩阵A的第k列进行秩1扰动
    k = 2;
    dAk = 0.1*randn(n,1);

    % 直接重新计算Arnoldi（作为对比基准）
    A(:,k) = A(:,k) + dAk;
    [V_new_ref, H_new_ref] = arnoldi(A, v1, m);
end

function [V, H] = arnoldi(A, v1, m)
    if nargin == 0
        n = 10;
        A = randn(n);
        v1 = randn(n,1);
        m = 4;
    end
    % 输入：
    % A  - 待分解的方阵
    % v1 - Krylov子空间的初始向量（维度与A匹配）
    % m  - Arnoldi迭代次数，即子空间的维度
    %
    % 输出：
    % V  - Krylov子空间的正交规范基（n×m矩阵）
    % H  - 对应的上Hessenberg矩阵（m×m矩阵）

    n = size(A, 1);     % 矩阵A的维度
    V = zeros(n, m+1);  % 存储正交基（预留第m+1个向量用于残差）
    H = zeros(m+1, m);  % 存储上Hessenberg矩阵（预留最后一行存储残差范数）

    % 初始化第一个正交基向量
    V(:, 1) = v1 / norm(v1);
    for j = 1:m
        % 计算A与当前基向量的乘积
        w = A * V(:, j);
        % 格拉姆-施密特正交化
        for i = 1:j
            H(i, j) = V(:, i)' * w;
            w = w - H(i, j) * V(:, i);
        end
        % 计算残差的范数
        H(j+1, j) = norm(w);
        if H(j+1, j) ~= 0
            % 生成下一个正交基向量
            V(:, j+1) = w / H(j+1, j);
        else
            % 残差为0，子空间已收敛，提前终止
            break;
        end
    end
    % 裁剪到m维的基和矩阵
    V = V(:, 1:m);
    H = H(1:m, 1:m);
end

增量更新解决方案 #

核心原理 #

矩阵的更新属于秩1修正：(A_{\text{new}} = A_{\text{old}} + dAk \cdot e_k^T)。利用原Arnoldi分解的正交性，通过以下步骤完成增量更新：

1. 原Arnoldi分解满足：(A_{\text{old}} V_{\text{old}} = V_{\text{old}} H_{\text{old}} + f e_m^T)，其中(f)是迭代残差向量。
2. 代入更新后的矩阵，得到：(A_{\text{new}} V_{\text{old}} = V_{\text{old}} H_{\text{old}} + f e_m^T + dAk \cdot (V_{\text{old}}(k,:)))。
3. 将增量项投影到原正交基上，得到对H的秩1修正，再通过正交化调整基向量V，最终得到新的Arnoldi分解。

实现代码 #

function [V_new, H_new] = update_arnoldi_after_col_change(V_old, H_old, v1, m, k, dAk)
    n = size(V_old, 1);
    % 计算dAk在原正交基上的投影与正交残差
    c = V_old' * dAk;  % m维投影系数向量
    proj_dAk = V_old * c;
    s = norm(dAk - proj_dAk);  % 正交残差的范数
    r = V_old(k, :);  % V_old的第k行，即e_k^T V_old

    % 对H进行秩1修正
    H_temp = H_old + c * r;

    % 构造A_new*V_old的近似矩阵
    AV_new = V_old * H_temp;
    % 添加秩1更新的增量项
    for j = 1:m
        AV_new(:,j) = AV_new(:,j) + dAk * r(j);
    end

    % 重新对AV_new执行格拉姆-施密特正交化，生成新的V和H
    V_new = zeros(n, m+1);
    H_new = zeros(m+1, m);
    V_new(:,1) = v1 / norm(v1);  % 初始向量保持不变

    for j = 1:m
        w = AV_new(:,j);
        % 正交化到已生成的基向量
        for i = 1:j
            H_new(i,j) = V_new(:,i)' * w;
            w = w - H_new(i,j)*V_new(:,i);
        end
        H_new(j+1,j) = norm(w);
        if H_new(j+1,j) ~= 0
            V_new(:,j+1) = w / H_new(j+1,j);
        else
            break;
        end
    end

    % 裁剪到m维输出
    V_new = V_new(:,1:m);
    H_new = H_new(1:m,1:m);
end

说明 #

• 该方法无需存储原矩阵A，仅依赖已有的V、H，以及更新的列增量dAk和列索引k即可完成更新。
• 相比重新运行完整Arnoldi算法，增量更新的计算量更小，尤其当m远小于n时优势显著。
• 如果原Arnoldi迭代因残差为0提前终止，可跳过残差相关的额外处理。

内容的提问来源于stack exchange，提问作者bardo