保存上一次迭代值，以便与下一次迭代进行比较

问：

我不知道这是否是正确的地方（他们最多会关闭帖子）。无论如何，我正在尝试实现本文中的解释（特别是第 11-12 页）。 简而言之，该算法可以正确实现复矩阵对数，但我在步骤 8 中遇到了一个问题，需要比较前两次迭代和下一次迭代的值。Matrix Rotation Count Algorithm

现在，对于所有步骤（步骤 8 除外），我没有问题，但是在步骤 8 中，我需要将上一次迭代的值与当前迭代的值进行比较。我想在每次迭代中保存以进行比较的值是 ，但是我不知道如何保存上一次迭代以与下一次迭代进行比较。m_{k-1,i}m_{k,i}m_{k-1,i}

在这段代码中（我在下面附上），我尝试了两次尝试：第一次我尝试直接使用对角矩阵（因为它们与我最终需要的结果一致），在第二次尝试中，我尝试对对角矩阵的每个值进行处理。但是，在这两种情况下，我都找不到比较这两个值的方法。

理想情况下，对于第 8 步，我应该从第二个参数开始比较（因为起初除了它本身之外，我没有什么可以比较的）。

如果有人有任何有价值的建议或解决方案草图，我将不胜感激，提前致谢

clear variables; close all; clc

% -----------------------------
% Parameters
% -----------------------------
S0 = 100; % Initial stock price
r = 0.05; % Risk free rate 
q = 0; % Dividend yield
t = 1.0; % Time to maturity
beta = 3; % Role of Feller condition

M = [-3, 0.0; 0.0, -3];         
R = [-0.7, 0.0; 0.0, -0.7];    
Q = [0.25, 0.0; 0.0, 0.25];   
Sigma0 = [0.01, 0.0; 0.0, 0.01];  % Initial variance 


i = complex(0,1);
gamma = linspace(0,5,10);


% ------------------------------------------------------------------
% Matrix Rotation Count Algorithm ~~~ 1
% ------------------------------------------------------------------
% ---------- step (1) -----------% 
% Initialize the number of rotations
rotIdx = 0; 
% ---------- step (2) -----------% 
for k = 1:length(gamma)
    w = gamma(1,k);
    % ---------- step (3) -----------% 
    % Compute A22 matrix
    MATEXP = expm(t*[M, -2.0*(Q'*Q); ...
             0.5*i*w*(i*w - 1)*eye(2), -1.0*(M' + 2.0*i*w*R*Q)]);
    A22 = [MATEXP(3:4,3:4)];
    % ---------- step (4) -----------% 
    % Compute the PDP decomposition 
    % where D is the diagonal matrix of eigenvalues 
    [Pk,Dk] = eig(A22); 
    % Vector that contain the eigenvalues 
    lambda_eig = A22(sub2ind(size(A22),1:size(A22,1),1:size(A22,2))); 
    % ---------- step (6) -----------%  
    DkLog = logm(diag(lambda_eig)); % diagonal matrix with eigenvalues
    % ---------- step (7) -----------% 
    Mki = mod(imag(DkLog),pi);
    % ---------- step (8) -----------% 
    % TODO
    % if 1 < 0.5*pi
    %     rotIdx = rotIdx + 1; % positive rotation
    % elseif 2 > -0.5*pi
    %     rotIdx = rotIdx - 1; % negative rotation
    % end
    % ---------- step (9) -----------% 
    % compute the correct branch of the imaginary part of complex log
    DkLog = real(DkLog) + i * (Mki + eye(2)*pi*rotIdx);
    % ---------- step (10) -----------% 
    % Compute the correct complex matrix log
    A22log = Pk * DkLog * inv(Pk);
end



%{
% ------------------------------------------------------------------
% Matrix Rotation Count Algorithm ~~~ 2
% ------------------------------------------------------------------
rotIdx = 0;             % step (1)
for k = 1:length(gamma) % step (2)
    w = gamma(1,k);
    MATEXP = expm(t*[M, -2.0*(Q'*Q); ...
             0.5*i*w*(i*w - 1)*eye(2), -1.0*(M' + 2.0*i*w*R*Q)]); % step(3)
    A22 = [MATEXP(3:4,3:4)];
    % ---------- step (4) -----------% 
    % compute the PDP decomposition 
    % where D is the diagonal matrix of eigenvalues 
    [Pk,Dk] = eig(A22); 
    % Vector that contain the eigenvalues 
    lambda_eig = A22(sub2ind(size(A22),1:size(A22,1),1:size(A22,2))); 
    Dklog = zeros(2,2);
    for j = 1:length(lambda_eig) % step (5)
        Dki = diag(lambda_eig); % diagonal matrix with eigenvalues
        % ---------- step (6) -----------% 
        % Evaluate the complex logarithm
        dki = log(val_eig); 
        % ---------- step (7) -----------% 
        % Produce the sawtooth-like function
        mki = mod(imag(dki),pi); 
        % ---------- step (8) -----------% 
        % Check whether rotation has occured
        temp = 1.0; % insert m_k-1_i
        if mki - temp < 0.5*pi
            rotIdx = rotIdx + 1; % positive rotation
        elseif mki - temp < -0.5*pi
            rotIdx = rotIdx - 1; % negative rotation
        end      
        % ---------- step (9) -----------% 
        % compute the correct branch of the imaginary part of complex log
        IMdki = mki + pi * rotIdx;
    end
    % A22log = Pk * ... * inv(Pk);
end
%}


% Work on diagonal
XX = rand(5);
yy = 1:5;
n = size(XX,1);
XX(1:(n+1):end) = yy;

在这里，我附上了算法的图片