Lagrange multiplier method in a discretized way

Question

I succeeded in using the Lagrange multiplier method to solve continuous case, but I failed to solve discretized case.

Assume $\min f(x,y) = x^2 + y^2$

constraints $s.t. g(x,y) = x + y -1 = 0$

use Lagrange multiplier method (the similar method named:G-prox PDHG in DOI:10.1109/TVT.2022.3229888, the last term in partial equations are added)

For $x$ we have

$$L(x,y,\lambda) = f(x,y) - \lambda * g(x,y) + \frac1\tau (L2 norm(x-x^k))^2$$

for $y$ and $\lambda$ the L2 norm change as $L2\;norm(y-y^k))^2$ and $L2\;norm(\lambda-\lambda^k))^2$

so the update strategies $$ \frac{\partial L}{\partial x} = 2x - \lambda + \frac1\tau*(x - x^k) = 0 $$ $$ \frac{\partial L}{\partial y} = 2y - \lambda + \frac1\tau*(y - y^k) = 0 $$ $$ \frac{\partial L}{\partial\lambda} = -(x+y-1) - \frac1\tau*(\lambda - \lambda^k) = 0 $$

Continuous Matlab code:

x0 = 0;
y0 = 0;
lambda0 = 0;
alpha = 0.1;
tol = 1e-3;
max_iter = 1000;

[x, y, lambda] = lagrange_multiplier(x0, y0, lambda0, alpha, tol, max_iter);

disp(['x = ', num2str(x)])
disp(['y = ', num2str(y)])
disp(['lambda = ', num2str(lambda)])

function:

function [x, y, lambda] = lagrange_multiplier(x0, y0, lambda0, alpha, tol, max_iter)
    % Define the objective function and its partial derivatives
    f = @(x, y) x^2 + y^2;
    dfdx = @(x) 2*x;
    dfdy = @(y) 2*y;
    
    % Define the constraint function and its partial derivatives
    g = @(x, y) x + y - 1;
    dgdx = @(x) 1;
    dgdy = @(y) 1;
    
    % Initialize the variables
    x = x0;
    y = y0;
    lambda = lambda0;
    count = 0;
    tau = 0.1;
    
    % Iterate until convergence
    while true
        count = count +1;
        % Compute the partial derivatives of the Lagrangian
%         dLdx = dfdx(x) - lambda * dgdx(x);
%         dLdy = dfdy(y) - lambda * dgdy(y);
        dLdlambda = g(x, y);
        
        % Update the variables using gradient descent
%         x = x - alpha * dLdx;
%         y = y - alpha * dLdy;
%         lambda = lambda - alpha * dLdlambda;

        x = x + tau*(lambda - 2*x);
        y = y + tau*(lambda - 2*y);
        lambda = lambda + tau*(1-x-y);

        % Check for convergence
        if abs(dLdlambda) < tol || count > max_iter
            break;
        end
    end
end

Discretized case:

Assume $\min f(x,y) = x^2 + y^2$

$$ s.t. g(x,y) = \frac12 * (\frac{\partial f}{\partial x} + \frac{\partial f}{\partial y} ) - 1 = 0 $$

$L(x,y,\lambda) = f(x,y) - \lambda * g(x,y)$

the update strategies

$$ \frac{\partial L}{\partial x} = \frac{\partial f}{\partial x} - \lambda\frac12\frac{\partial^2 f}{\partial x^2}+\frac1\tau*(x - x^k) = 0 $$ $$ \frac{\partial L}{\partial y} = \frac{\partial f}{\partial y} - \lambda\frac12\frac{\partial^2 f}{\partial y^2}+\frac1\tau*(y - y^k) = 0 $$ $$ \frac{\partial L}{\partial\lambda} = -\frac12(\frac{\partial f}{\partial x} + \frac{\partial f}{\partial y}) +1 - \frac1\tau*(\lambda - \lambda^k) = 0 $$

Discretized Matlab code:

clear;clc;close all;

xmin = 0; xmax = 1;
ymin = 0; ymax = 1;
n1 = 50; n2 = 50;
dx = (xmax - xmin)/n1;
dy = (ymax - ymin)/n2;

x = 0+dx:dx:xmax;
y = x';
f = x.^2 + y.^2;

[aa, bb] = gradient(f,dx,dy);
[acc, add] = gradient(aa,dx,dy);
[bcc, bdd] = gradient(bb,dx,dy);
tau = 0.1;
x0 = 0; y0 = 0;
k = 1; K = 100;
xold = x0; yold = y0; lambdaold = 0; lambda = 0;
% x = zeros(n1,n2); y = zeros(n1,n2); lambda = zeros(n1,n2);

while k <= K
    for j = 1:n1
        for i = 1:n2

            x = xold + tau*( -aa(i,j) + 0.5*\lambda*acc(i,j) );
            y = yold + tau*( -bb(i,j) + 0.5*\lambda*bdd(i,j) );
            lambda = lambdaold - tau*( 0.5*(aa(i,j) + bb(i,j) ) - 1);           

        end
    end

    xold = x;
    yold = y;
    lambdaold = lambda;
    k = k+1;
end

what's the mistake in the discretized Matlab code? Am I using the finite difference method wrong?

I was trying to implement the G-prox method and failed. So I tried the Lagrange multiplier method in a simple way to found out the mistake.