在matlab中使用定点法实现快速优化算法

Implement a fast optimization algorithm using fixed point method in matlab

我正在 matlab 中使用定点法实现快速优化算法。该方法的目标是找到 u 的最优值。表示 u={u_i,i=1..2}。 u 的最优值可以通过以下步骤获得:

对不起我的形象,因为我不能在这里输入数学方程式。 为了完成这项任务,我试图找到你遵循上述步骤。但是,我不知道如何在等式 25 中实现项 \sum_{j!=i} (u_j-1)。这是我的代码。请查看它,您能否给我一些关于我的实施的评论或建议以更正它们。目前,我尝试 运行 该代码,但它给出了错误的答案。

function u = compute_u_TV(Im0, N_class)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
theta=0.001;
gamma=0.01;
tau=0.1;
sigma=0.1;
N_class=2; % only have u1 and u2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Iterative segmentation process
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:N_class
    v(:,:,i) = Im0/max(Im0(:)); % u between 0 and 1.   
    qxv(:,:,i) = zeros(size(Im0));
    qyv(:,:,i) = zeros(size(Im0));
    u(:,:,i) = v(:,:,i);
for iteration=1:10000
    u_temp=u;
    % Update v
    Divqi = ( BackwardX(qxv(:,:,i)) + BackwardY(qyv(:,:,i)) );
    Term = Divqi - u(:,:,i)/ (theta*gamma);
    TermX = ForwardX(Term);
    TermY = ForwardY(Term); 
    Norm = sqrt(TermX.^2 + TermY.^2);
    Denom = 1 + tau*Norm;
    %Equation 24
    qxv(:,:,i) = (qxv(:,:,i) + tau*TermX)./Denom;
    qyv(:,:,i) = (qyv(:,:,i) + tau*TermY)./Denom;  
    v(:,:,i) = u(:,:,i) - theta*gamma* Divqi;  %Equation 23   
    % Update u  
    u(:,:,i) = (v(:,:,i) - theta* gamma* Divqi -theta*gamma*sigma*(sum(u(:))-u(:,:,i)-1))./(1+theta* gamma*sigma);
    u(:,:,i) = max(u(:,:,i),0); 
    u(:,:,i) = min(u(:,:,i),1);
    check=u_temp(:,:,i)-u(:,:,i);
    if(abs(sum(check(:)))<=0.1)
        break;  
    end
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sub-functions- X.Berson
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dx]=BackwardX(u);
[Ny,Nx] = size(u);
dx = u;
dx(2:Ny-1,2:Nx-1)=( u(2:Ny-1,2:Nx-1) - u(2:Ny-1,1:Nx-2) );
dx(:,Nx) = -u(:,Nx-1);

function [dy]=BackwardY(u);
[Ny,Nx] = size(u);
dy = u;
dy(2:Ny-1,2:Nx-1)=( u(2:Ny-1,2:Nx-1) - u(1:Ny-2,2:Nx-1) );
dy(Ny,:) = -u(Ny-1,:);

function [dx]=ForwardX(u);
[Ny,Nx] = size(u);
dx = zeros(Ny,Nx);
dx(1:Ny-1,1:Nx-1)=( u(1:Ny-1,2:Nx) - u(1:Ny-1,1:Nx-1) );

function [dy]=ForwardY(u);
[Ny,Nx] = size(u);
dy = zeros(Ny,Nx);
dy(1:Ny-1,1:Nx-1)=( u(2:Ny,1:Nx-1) - u(1:Ny-1,1:Nx-1) );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End of sub-function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

你应该做

u(:,:,i) = (v(:,:,i) - theta* gamma* Divqi -theta*gamma*sigma* ...
           (sum(u(:,:,1:size(u,3) ~= i),3) -1))./(1+theta* gamma*sigma);

您搜索的部分是

sum(u(:,:,1:size(u,3) ~= i),3)

让我们分解一下:

 1:size(u,3) ~= i

是一个向量,包含从 1 到第三维最大大小 u 的所有值,但 i 除外。

然后

u(:,:,1:size(u,3) ~= i)

u除了j = i

的所有第三维矩阵

最后,

sum(...,3) 

是所有矩阵的三次元之和。

如果有帮助请告诉我!