When i use cvx to solve this convex problem, something wrong, can you help me?

In the first two iterations, there is no error. However, there is a error in the third iteration.
i have prove this problem is convex, but i have know why the error occurs. Can you help me?

clc,clear;
close all;

%%
K=2;
Nt = 4;
M = 10;
iter_num = 100;
inner_iter_num = 1000;
h_irs_u = rand(M,K)+1i*rand(M,K);
h_ai = rand(M,Nt)+1i*rand(M,Nt);
W = zeros(Nt,K,iter_num+1);

val = zeros(1,iter_num+1);
mu = zeros(K,iter_num+1);
v = zeros(K,iter_num+1);
Phi1= zeros(M,M,iter_num+1);
Phi1(:,:,1) = sqrt(0.5)*diag(exp(1i*2*pi*rand(M,1)));
Phi2= zeros(M,M,iter_num+1);
Phi2(:,:,1) = sqrt(0.5)*diag(exp(1i*2*pi*rand(M,1)));
Pmax = 3;
sigma = 1;
delta = 1e-3;
% initialization
W(:,:,1) = rand(Nt,K)+1i*rand(Nt,K);
mu(:,1) = rand(K,1);
v(:,1) = rand(K,1);
rho = 1.5;
for iter = 1:iter_num

% optimize mu and v
% mu_k = zeros(K,1);
% v_k  = zeros(K,1);
H = [h_irs_u(:,1)'*Phi1(:,:,iter)*h_ai;
    h_irs_u(:,2)'*Phi2(:,:,iter)*h_ai]';
for kk=1:K
   rho = real(conj(v(kk,iter))*H(:,kk)'*W(:,kk,iter));
   mu(kk,iter+1)= rho/2*(rho+sqrt(rho^2+4));
   v(kk,iter+1) = sqrt(1+mu(kk,iter+1))*H(:,kk)'*W(:,kk,iter)/(norm(H(:,kk)'*W(:,:,iter))^2+sigma);
end
   
% optimize W
cvx_begin quiet
variable W0(Nt,K) complex
maximize (2*sqrt(1+mu(1,iter+1))*real(v(1,iter+1)'*H(:,1)'*W0(:,1))-abs(v(1,iter+1))^2*(square_pos(norm(H(:,1)'*W0))+sigma)...
    +2*sqrt(1+mu(2,iter+1))*real(v(2,iter+1)'*H(:,2)'*W0(:,2))-abs(v(2,iter+1))^2*(square_pos(norm(H(:,2)'*W0))+sigma))

subject to
norm(W0)<=Pmax;
cvx_end
W(:,:,iter+1)=W0;

% optimize phase shift at IRS
    Theta1 = zeros(M,inner_iter_num+1);
    Theta1(:,1) = diag(Phi1(:,:,iter));
    Theta2 = zeros(M,inner_iter_num+1);
    Theta2(:,1) = diag(Phi2(:,:,iter));
    psi1=zeros(M,inner_iter_num+1);
    psi1(:,1) = exp(1i*2*pi*rand(M,1));
    psi2= zeros(M,inner_iter_num+1);
    psi2(:,1) = exp(1i*2*pi*rand(M,1));
    u1 = zeros(M,inner_iter_num+1);
    u1(:,1) = rand(M,1);
    u2 = zeros(M,inner_iter_num+1);
    u2(:,1) = rand(M,1);
    x = zeros(1,inner_iter_num+1);
    x(1) = sqrt(0.5);
    y = zeros(1,inner_iter_num+1);
    y(1)=sqrt(0.5);
    for inner_iter=1:inner_iter_num
        % optimize theta1 and theta2
        cvx_begin quiet
        variable theta1(M,1) complex
        variable theta2(M,1) complex
        H1 = [h_irs_u(:,1)'*diag(theta1)*h_ai;
        h_irs_u(:,2)'*diag(theta2)*h_ai]';
        maximize(2*sqrt(1+mu(1,iter+1))*real(v(1,iter+1)'*H1(:,1)'*W0(:,1))-abs(v(1,iter+1))^2*(square_pos(norm(H1(:,1)'*W0))+sigma)...
        +2*sqrt(1+mu(2,iter+1))*real(v(2,iter+1)'*H1(:,2)'*W0(:,2))-abs(v(2,iter+1))^2*(square_pos(norm(H1(:,2)'*W0))+sigma)...
         -rho/2*(square_pos(norm(theta1-x(inner_iter)*psi1(:,inner_iter)+u1(:,inner_iter)))+square_pos(norm(theta2-y(inner_iter)*psi2(:,inner_iter)+u2(:,inner_iter)))))
        cvx_end
        
        Theta1(:,inner_iter+1) = theta1;
        Theta2(:,inner_iter+1) = theta2;
        % optimze psi1 and psi2
        psi1(:,inner_iter+1) = exp(1i*angle(theta1+u1(:,inner_iter)));
        psi2(:,inner_iter+1) = exp(1i*angle(theta2+u2(:,inner_iter)));
        
        % optimize x and y
        a = 2*real((Theta1(:,inner_iter+1)+u1(:,inner_iter))'*psi1(:,inner_iter+1));
        c = 2*real((Theta2(:,inner_iter+1)+u2(:,inner_iter))'*psi2(:,inner_iter+1));
        d = -(M+norm(Theta1(:,inner_iter+1)+u1(:,inner_iter))^2+norm(Theta2(:,inner_iter+1)+u2(:,inner_iter))^2);
        x(inner_iter+1) = a/sqrt(a^2+c^2);
        y(inner_iter+1) = c/sqrt(a^2+c^2);
        % optimze u1 and u2
        u1(:,inner_iter+1) = u1(:,inner_iter) + Theta1(:,inner_iter+1) -x(inner_iter+1)*psi1(:,inner_iter+1);
        u2(:,inner_iter+1) = u2(:,inner_iter) + Theta2(:,inner_iter+1) -y(inner_iter+1)*psi2(:,inner_iter+1);
%         disp([norm(Theta1(:,inner_iter+1)-x(inner_iter+1)*psi1(:,inner_iter+1)), norm(Theta2(:,inner_iter+1)-y(inner_iter+1)*psi2(:,inner_iter+1))]);
        err1 = norm(Theta1(:,inner_iter+1)-x(inner_iter+1)*psi1(:,inner_iter+1));
        err2 = norm(Theta2(:,inner_iter+1)-y(inner_iter+1)*psi2(:,inner_iter+1));
        if err1<=delta && err2<=delta
            break
        end
    end
Phi1(:,:,iter+1) = diag(Theta1(:,inner_iter+1));
Phi2(:,:,iter+1) = diag(Theta2(:,inner_iter+1));
H = [h_irs_u(:,1)'*Phi1(:,:,iter+1)*h_ai;
    h_irs_u(:,2)'*Phi2(:,:,iter+1)*h_ai]';
val(iter+1) = log2(1+abs(H(:,1)'*W0(:,1))^2/(abs(H(:,1)'*W0(:,2))^2+sigma)) + log2(1+abs(H(:,2)'*W0(:,2))^2/(abs(H(:,2)'*W0(:,1))^2+sigma));
disp([iter,val(iter+1) ])
if abs(val(iter+1) - val(iter))<=delta
    break
end
end

figure
plot(val(2:iter+1),'b-','linewidth',1.5);
xlabel('iteration index');
ylabel('Objective value');
grid on

Check the numerical values of the coefficients used in the objective function on the iteration in which the error occurs.

You can read many posts by me on this forum of the hazards of crude unsafeguarded SCA. What you have experienced is perhaps one such manifestation.

ok, thank you very much