基于拉普拉斯梯度法的L1范数最小化的MATLAB实现

 ✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，代码获取、论文复现及科研仿真合作可私信。

🍎个人主页：Matlab科研工作室

🍊个人信条：格物致知。

更多Matlab完整代码及仿真定制内容点击👇

智能优化算法       神经网络预测       雷达通信       无线传感器        电力系统

信号处理              图像处理               路径规划       元胞自动机        无人机 

🔥 内容介绍

摘要

本文提出了一种基于拉普拉斯梯度法的L1范数最小化算法，该算法可以有效地解决图像去噪、图像复原和图像分割等问题。该算法的核心思想是利用拉普拉斯梯度算子来计算图像的梯度信息，然后利用L1范数正则化项来惩罚梯度值较大的像素，从而达到图像去噪和复原的目的。

1. 算法原理

1.1 拉普拉斯梯度算子

拉普拉斯梯度算子是一种二阶微分算子，它可以用来计算图像的梯度信息。拉普拉斯梯度算子的定义如下：

\nabla^2 f(x, y) = \frac{\partial^2 f(x, y)}{\partial x^2} + \frac{\partial^2 f(x, y)}{\partial y^2}

其中，f(x, y)是图像的灰度值，∇^2 f(x, y)是图像的拉普拉斯梯度。

1.2 L1范数正则化项

L1范数正则化项是一种惩罚函数，它可以用来惩罚梯度值较大的像素。L1范数正则化项的定义如下：

\|f(x, y)\|_1 = \sum_{x, y} |f(x, y)|

其中，f(x, y)是图像的灰度值。

1.3 目标函数

基于拉普拉斯梯度法的L1范数最小化算法的目标函数如下：

E(f(x, y)) = \frac{1}{2}\|f(x, y) - g(x, y)\|^2 + \lambda\|f(x, y)\|_1

其中，f(x, y)是图像的灰度值，g(x, y)是原始图像的灰度值，λ是正则化参数。

2. 算法步骤

基于拉普拉斯梯度法的L1范数最小化算法的步骤如下：

1. 计算图像的拉普拉斯梯度。

2. 计算L1范数正则化项。

3. 计算目标函数。

4. 使用梯度下降法更新图像的灰度值。

5. 重复步骤3和步骤4，直到目标函数收敛.

📣 部分代码

%% This function is modified from Matlab Package: L1-Homotopy​% BPDN_homotopy_function.m% % Solves the following basis pursuit denoising (BPDN) problem% min_x  \lambda ||x||_1 + 1/2*||b-Ax||_2^2%% Inputs: % A - m x n measurement matrix% b - measurement vector% lambda - final value of regularization parameter% maxiter - maximum number of homotopy iterations%% Outputs:% x_out - output for BPDN% total_iter - number of homotopy iterations taken by the solver% % Written by: Salman Asif, Georgia Tech% Email: sasif@ece.gatech.edu%%-------------------------------------------+% Copyright (c) 2007.  Muhammad Salman Asif %-------------------------------------------+​function [x_out, total_iter, timeSteps, errorSteps, epsSteps] = SolveHomotopy(A, b, varargin)​global N gamma_x z_x  xk_temp del_x_vec pk_temp dk epsilon isNonnegative​t0 = tic ;​lambda = 1e-6;maxiter = 100;isNonnegative = false;verbose = false;xk_1 = [];​STOPPING_TIME = -2;STOPPING_GROUND_TRUTH = -1;STOPPING_DUALITY_GAP = 1;STOPPING_SPARSE_SUPPORT = 2;STOPPING_OBJECTIVE_VALUE = 3;STOPPING_SUBGRADIENT = 4;STOPPING_DEFAULT = STOPPING_OBJECTIVE_VALUE;​stoppingCriterion = STOPPING_DEFAULT;​% Parse the optional inputs.if (mod(length(varargin), 2) ~= 0 ),    error(['Extra Parameters passed to the function ''' mfilename ''' must be passed in pairs.']);endparameterCount = length(varargin)/2;​for parameterIndex = 1:parameterCount,    parameterName = varargin{parameterIndex*2 - 1};    parameterValue = varargin{parameterIndex*2};    switch lower(parameterName)        case 'stoppingcriterion'            stoppingCriterion = parameterValue;        case 'initialization'            xk_1 = parameterValue;            if ~all(size(xk_1)==[n,1])                error('The dimension of the initial x0 does not match.');            end        case 'groundtruth'            xG = parameterValue;        case 'lambda'            lambda = parameterValue;        case 'maxiteration'            maxiter = parameterValue;        case 'isnonnegative'            isNonnegative = parameterValue;        case 'tolerance'            tolerance = parameterValue;        case 'verbose'            verbose = parameterValue;        case 'maxtime'            maxTime = parameterValue;        otherwise            error(['The parameter ''' parameterName ''' is not recognized by the function ''' mfilename '''.']);    endendclear varargin​timeSteps = nan(1,maxiter) ;errorSteps = nan(1,maxiter) ;epsSteps = nan(1,maxiter);​[K, N] = size(A);​% Initialization of primal and dual sign and supportz_x = zeros(N,1);gamma_x = [];       % Primal support​% Initial stepPrimal_constrk = -A'*b;​if isNonnegative    [c i]  = min(Primal_constrk);    c = max(-c, 0);else    [c i] = max(abs(Primal_constrk));end​epsilon = c;nz_x = zeros(N,1);if isempty(xk_1)    xk_1 = zeros(N,1);    gamma_xk = i;else    gamma_xk = find(abs(xk_1)>eps*10);    nz_x(gamma_xk) = 1;endf = epsilon*norm(xk_1,1) + 1/2*norm(b-A*xk_1)^2;z_x(gamma_xk) = -sign(Primal_constrk(gamma_xk));%Primal_constrk(gamma_xk) = sign(Primal_constrk(gamma_xk))*epsilon;​z_xk = z_x;​% loop parametersiter = 0;out_x = [];old_delta = 0;count_delta_stop = 0;​AtgxAgx = A(:,gamma_xk)'*A(:,gamma_xk);iAtgxAgx = inv(A(:,gamma_xk)'*A(:,gamma_xk));while iter < maxiter    iter = iter+1;    % warning('off','MATLAB:divideByZero')​    gamma_x = gamma_xk;    z_x = z_xk;    x_k = xk_1;​    %%%%%%%%%%%%%%%%%%%%%    %%%% update on x %%%%    %%%%%%%%%%%%%%%%%%%%%        % Update direction    del_x = iAtgxAgx*z_x(gamma_x);    del_x_vec = zeros(N,1);    del_x_vec(gamma_x) = del_x;​    %dk = A'*(A*del_x_vec);    Asupported = A(:,gamma_x);    Agdelx = Asupported*del_x;    dk = A'*Agdelx;        %%% CONTROL THE MACHINE PRECISION ERROR AT EVERY OPERATION: LIKE BELOW.     pk_temp = Primal_constrk;    gammaL_temp = find(abs(abs(Primal_constrk)-epsilon)<min(epsilon,2*eps));    pk_temp(gammaL_temp) = sign(Primal_constrk(gammaL_temp))*epsilon;        xk_temp = x_k;    xk_temp(abs(x_k)<2*eps) = 0;    %%%---        % Compute the step size    [i_delta, delta, out_x] = update_primal(out_x);        if old_delta < 4*eps && delta < 4*eps        count_delta_stop = count_delta_stop + 1;                if count_delta_stop >= 500            if verbose                disp('stuck in some corner');            end            break;        end    else        count_delta_stop = 0;    end    old_delta = delta;        xk_1 = x_k+delta*del_x_vec;    Primal_constrk = Primal_constrk+delta*dk;    epsilon_old = epsilon;    epsilon = epsilon-delta;​    if epsilon <= lambda;%         xk_1 = x_k + (epsilon_old-lambda)*del_x_vec;%         disp('Reach prescribed lambda in SolveHomotopy.');​        gamma_x0 = find(abs(xk_1) > 1e-9);        AtgxAgx0 = A(:,gamma_x0)'*A(:,gamma_x0);        x_temp = AtgxAgx0 \ (A(:,gamma_x0)' * b);        xk_1 = zeros(N,1);        xk_1(gamma_x0) = x_temp;​        timeSteps(iter) = toc(t0) ;        % MODIFIED -- measure the difference of the L1 norms, instead of the Euclidean distance        %errorSteps(iter) = norm(xk_1-xG) ; % original measure        errorSteps(iter) = abs(norm(xk_1, 1) - norm(xG, 1)) ;        epsSteps(iter) = epsilon;               break;    end        timeSteps(iter) = toc(t0) ;    % MODIFIED -- measure the difference of the L1 norms, instead of the Euclidean distance    %errorSteps(iter) = norm(xk_1-xG) ; % original measure  errorSteps(iter) = abs(norm(xk_1, 1) - norm(xG, 1)) ;    epsSteps(iter) = epsilon;        % compute stopping criteria and test for termination    keep_going = true;        switch stoppingCriterion        case STOPPING_GROUND_TRUTH            keep_going = norm(xk_1-xG)>tolerance;        case STOPPING_SPARSE_SUPPORT            if delta~=0                nz_x_prev = nz_x;                nz_x = (abs(xk_1)>eps*10);                num_nz_x = sum(nz_x(:));                num_changes_active = (sum(nz_x(:)~=nz_x_prev(:)));                if num_nz_x >= 1                    criterionActiveSet = num_changes_active / num_nz_x;                    keep_going = (criterionActiveSet > tolerance);                end            end        case STOPPING_DUALITY_GAP            error('Duality gap is not a valid stopping criterion for Homotopy.');        case STOPPING_OBJECTIVE_VALUE            if delta~=0                % continue if not yeat reached target value tolA                prev_f = f;                f = lambda*norm(xk_1,1) + 1/2*norm(b-Asupported*xk_1(gamma_x))^2;                keep_going = (abs((prev_f-f)/prev_f) > tolerance);            end        case STOPPING_SUBGRADIENT            keep_going = norm(delta*del_x_vec)>tolerance;        case STOPPING_TIME            keep_going = timeSteps(iter) < maxTime ;        otherwise,            error('Undefined stopping criterion');    end % end of the stopping criteria switch​    %     if keep_going && norm(xk_1 - x_k)<100*eps%         if verbose%             disp('The iteration is stuck.');%         end%         keep_going = false;%     end        if ~keep_going        break;    end​    if ~isempty(out_x)        % If an element is removed from gamma_x        len_gamma = length(gamma_x);        outx_index = find(gamma_x==out_x(1));        gamma_x(outx_index) = gamma_x(len_gamma);        gamma_x(len_gamma) = out_x(1);        gamma_x = gamma_x(1:len_gamma-1);        gamma_xk = gamma_x;                rowi = outx_index; % ith row of A is swapped with last row (out_x)        colj = outx_index; % jth column of A is swapped with last column (out_lambda)        AtgxAgx_ij = AtgxAgx;        temp_row = AtgxAgx_ij(rowi,:);        AtgxAgx_ij(rowi,:) = AtgxAgx_ij(len_gamma,:);        AtgxAgx_ij(len_gamma,:) = temp_row;        temp_col = AtgxAgx_ij(:,colj);        AtgxAgx_ij(:,colj) = AtgxAgx_ij(:,len_gamma);        AtgxAgx_ij(:,len_gamma) = temp_col;        iAtgxAgx_ij = iAtgxAgx;        temp_row = iAtgxAgx_ij(colj,:);        iAtgxAgx_ij(colj,:) = iAtgxAgx_ij(len_gamma,:);        iAtgxAgx_ij(len_gamma,:) = temp_row;        temp_col = iAtgxAgx_ij(:,rowi);        iAtgxAgx_ij(:,rowi) = iAtgxAgx_ij(:,len_gamma);        iAtgxAgx_ij(:,len_gamma) = temp_col;                AtgxAgx = AtgxAgx_ij(1:len_gamma-1,1:len_gamma-1);                %iAtgxAgx = update_inverse(AtgxAgx_ij, iAtgxAgx_ij,2);        n = size(AtgxAgx_ij,1);        %delete columns        Q11 = iAtgxAgx_ij(1:n-1,1:n-1);        Q12 = iAtgxAgx_ij(1:n-1,n);        Q21 = iAtgxAgx_ij(n,1:n-1);        Q22 = iAtgxAgx_ij(n,n);        Q12Q21_Q22 = Q12*(Q21/Q22);        iAtgxAgx = Q11 - Q12Q21_Q22;                xk_1(out_x(1)) = 0;    else        % If an element is added to gamma_x        gamma_xk = [gamma_x; i_delta];        new_x = i_delta;​        AtgxAnx = A(:,gamma_x)'*A(:,new_x);        AtgxAgx_mod = [AtgxAgx AtgxAnx; AtgxAnx' A(:,new_x)'*A(:,i_delta)];                AtgxAgx = AtgxAgx_mod;                %iAtgxAgx = update_inverse(AtgxAgx, iAtgxAgx,1);        n = size(AtgxAgx,1);        % add columns        iA11 = iAtgxAgx;        iA11A12 = iA11*AtgxAgx(1:n-1,n);        A21iA11 = AtgxAgx(n,1:n-1)*iA11;        S = AtgxAgx(n,n)-AtgxAgx(n,1:n-1)*iA11A12;        Q11_right = iA11A12*(A21iA11/S);        %     Q11 = iA11+ Q11_right;        %     Q12 = -iA11A12/S;        %     Q21 = -A21iA11/S;        %     Q22 = 1/S;        iAtgxAgx = zeros(n);        %iAtB = [Q11 Q12; Q21 Q22];        iAtgxAgx(1:n-1,1:n-1) = iA11+ Q11_right;        iAtgxAgx(1:n-1,n) = -iA11A12/S;        iAtgxAgx(n,1:n-1) =  -A21iA11/S;        iAtgxAgx(n,n) = 1/S;                xk_1(i_delta) = 0;    end​    z_xk = zeros(N,1);    z_xk(gamma_xk) = -sign(Primal_constrk(gamma_xk));    Primal_constrk(gamma_x) = sign(Primal_constrk(gamma_x))*epsilon;endtotal_iter = iter;​x_out = xk_1;​timeSteps = timeSteps(1:total_iter) ;errorSteps = errorSteps(1:total_iter) ;epsSteps = epsSteps(1:total_iter);​% Debiasing%x_out = zeros(N,1);%x_out(gamma_x) = A(:,gamma_x)\b;​% update_primal.m%% This function computes the minimum step size in the primal update direction and% finds change in the primal or dual support with that step.%% Inputs:% gamma_x - current support of x% gamma_lambda - current support of lambda% z_x - sign sequence of x% z_lambda - sign sequence of lambda% del_x_vec - primal update direction% pk_temp% dk% epsilon - current value of epsilon% out_lambda - element removed from support of lambda in previous step (if any)%% Outputs:% i_delta - index corresponding to newly active primal constraint (new_lambda)% out_x - element in x shrunk to zero% delta - primal step size%% Written by: Salman Asif, Georgia Tech% Email: sasif@ece.gatech.edu​function [i_delta, delta, out_x] = update_primal(out_x)​global N gamma_x z_x  xk_temp del_x_vec pk_temp dk epsilon isNonnegative​gamma_lc = setdiff(1:N, union(gamma_x, out_x));​if isNonnegative    delta1 = inf;else    delta1_constr = (epsilon-pk_temp(gamma_lc))./(1+dk(gamma_lc));    delta1_pos_ind = find(delta1_constr>0);    delta1_pos = delta1_constr(delta1_pos_ind);    [delta1 i_delta1] = min(delta1_pos);    if isempty(delta1)        delta1 = inf;    endend​delta2_constr = (epsilon+pk_temp(gamma_lc))./(1-dk(gamma_lc));delta2_pos_ind = find(delta2_constr>0);delta2_pos = delta2_constr(delta2_pos_ind);[delta2 i_delta2] = min(delta2_pos);if isempty(delta2)    delta2 = inf;end​if delta1>delta2    delta = delta2;    i_delta = gamma_lc(delta2_pos_ind(i_delta2));else    delta = delta1;    i_delta = gamma_lc(delta1_pos_ind(i_delta1));end​delta3_constr = (-xk_temp(gamma_x)./del_x_vec(gamma_x));delta3_pos_index = find(delta3_constr>0);[delta3 i_delta3] = min(delta3_constr(delta3_pos_index));out_x_index = gamma_x(delta3_pos_index(i_delta3));​out_x = [];if ~isempty(delta3) && (delta3 > 0) && (delta3 <= delta)    delta = delta3;    out_x = out_x_index;end​%%% THESE ARE PROBABLY UNNECESSARY %%% NEED TO REMOVE THEM. ​% The following checks are just to deal with degenerate cases when more% than one elements want to enter or leave the support at any step% (e.g., Bernoulli matrix with small number of measurements)​% This one is ONLY for those indices which are zero. And we don't know where% will its dx point in next steps, so after we calculate dx and its in opposite% direction to z_x, we will have to remove that index from the support.xk_1 = xk_temp+delta*del_x_vec;xk_1(out_x) = 0;​wrong_sign = find(sign(xk_1(gamma_x)).*z_x(gamma_x)==-1);if isNonnegative    wrong_sign = union(wrong_sign, find(xk_1(gamma_x)<0));endif ~isempty(gamma_x(wrong_sign))    delta = 0;    % can also choose specific element which became non-zero first but all    % that matters here is AtA(gx,gl) doesn't become singular.    [val_wrong_x ind_wrong_x] =  sort(abs(del_x_vec(gamma_x(wrong_sign))),'descend');    out_x = gamma_x(wrong_sign(ind_wrong_x));end​% If more than one primal constraints became active in previous iteration i.e.,% more than one elements wanted to enter the support and we added only one.% So here we need to check if those remaining elements are still active.i_delta_temp = gamma_lc(abs(pk_temp(gamma_lc)+delta*dk(gamma_lc))-(epsilon-delta) >= 10*eps);if ~isempty(i_delta_temp)​    i_delta_more = i_delta_temp;    if (length(i_delta_more)>=1) && (~any((i_delta_temp==i_delta)))        % ideal way would be to check that incoming element doesn't make AtA        % singular!        [v_temp i_temp] = max(-pk_temp(i_delta_more)./dk(i_delta_more));        i_delta = i_delta_more(i_temp);        delta = 0;        out_x = [];    endend

⛳️ 运行结果

🔗 参考文献

[YGZ+10] A. Yang, A. Ganesh, Z. Zhou, S. Sastry, and Y. Ma. 

Fast L1-Minimization Algorithms for Robust Face Recognition. 

arXiv:1007.3753 [cs.CV] -- https://arxiv.org/abs/1007.3753

🎈 部分理论引用网络文献，若有侵权联系博主删除

🎁  关注我领取海量matlab电子书和数学建模资料

👇  私信完整代码和数据获取及论文数模仿真定制

1 各类智能优化算法改进及应用

生产调度、经济调度、装配线调度、充电优化、车间调度、发车优化、水库调度、三维装箱、物流选址、货位优化、公交排班优化、充电桩布局优化、车间布局优化、集装箱船配载优化、水泵组合优化、解医疗资源分配优化、设施布局优化、可视域基站和无人机选址优化、背包问题、 风电场布局、时隙分配优化、 最佳分布式发电单元分配、多阶段管道维修、 工厂-中心-需求点三级选址问题、 应急生活物质配送中心选址、 基站选址、 道路灯柱布置、 枢纽节点部署、 输电线路台风监测装置、 集装箱船配载优化、 机组优化、 投资优化组合、云服务器组合优化、 天线线性阵列分布优化

2 机器学习和深度学习方面

2.1 bp时序、回归预测和分类

2.2 ENS声神经网络时序、回归预测和分类

2.3 SVM/CNN-SVM/LSSVM/RVM支持向量机系列时序、回归预测和分类

2.4 CNN/TCN卷积神经网络系列时序、回归预测和分类

2.5 ELM/KELM/RELM/DELM极限学习机系列时序、回归预测和分类

2.6 GRU/Bi-GRU/CNN-GRU/CNN-BiGRU门控神经网络时序、回归预测和分类

2.7 ELMAN递归神经网络时序、回归\预测和分类

2.8 LSTM/BiLSTM/CNN-LSTM/CNN-BiLSTM/长短记忆神经网络系列时序、回归预测和分类

2.9 RBF径向基神经网络时序、回归预测和分类

2.10 DBN深度置信网络时序、回归预测和分类

2.11 FNN模糊神经网络时序、回归预测

2.12 RF随机森林时序、回归预测和分类

2.13 BLS宽度学习时序、回归预测和分类

2.14 PNN脉冲神经网络分类

2.15 模糊小波神经网络预测和分类

2.16 时序、回归预测和分类

2.17 时序、回归预测预测和分类

2.18 XGBOOST集成学习时序、回归预测预测和分类

方向涵盖风电预测、光伏预测、电池寿命预测、辐射源识别、交通流预测、负荷预测、股价预测、PM2.5浓度预测、电池健康状态预测、用电量预测、水体光学参数反演、NLOS信号识别、地铁停车精准预测、变压器故障诊断

2.图像处理方面

图像识别、图像分割、图像检测、图像隐藏、图像配准、图像拼接、图像融合、图像增强、图像压缩感知

3 路径规划方面

旅行商问题（TSP）、车辆路径问题（VRP、MVRP、CVRP、VRPTW等）、无人机三维路径规划、无人机协同、无人机编队、机器人路径规划、栅格地图路径规划、多式联运运输问题、 充电车辆路径规划（EVRP）、 双层车辆路径规划（2E-VRP）、 油电混合车辆路径规划、 船舶航迹规划、 全路径规划规划、 仓储巡逻

4 无人机应用方面

无人机路径规划、无人机控制、无人机编队、无人机协同、无人机任务分配、无人机安全通信轨迹在线优化、车辆协同无人机路径规划

5 无线传感器定位及布局方面

传感器部署优化、通信协议优化、路由优化、目标定位优化、Dv-Hop定位优化、Leach协议优化、WSN覆盖优化、组播优化、RSSI定位优化

6 信号处理方面

信号识别、信号加密、信号去噪、信号增强、雷达信号处理、信号水印嵌入提取、肌电信号、脑电信号、信号配时优化

7 电力系统方面

微电网优化、无功优化、配电网重构、储能配置、有序充电

8 元胞自动机方面

交通流 人群疏散 病毒扩散 晶体生长 金属腐蚀

9 雷达方面

卡尔曼滤波跟踪、航迹关联、航迹融合