✅作者简介：热爱科研的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范数最小化算法的步骤如下：
-
计算图像的拉普拉斯梯度。
-
计算L1范数正则化项。
-
计算目标函数。
-
使用梯度下降法更新图像的灰度值。
-
重复步骤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 isNonnegativet0 = 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.');endcase '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;otherwiseerror(['The parameter ''' parameterName ''' is not recognized by the function ''' mfilename '''.']);endendclear varargintimeSteps = 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));endepsilon = c;nz_x = zeros(N,1);if isempty(xk_1)xk_1 = zeros(N,1);gamma_xk = i;elsegamma_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 < maxiteriter = iter+1;% warning('off','MATLAB:divideByZero')gamma_x = gamma_xk;z_x = z_xk;x_k = xk_1;%%%%%%%%%%%%%%%%%%%%%%%%% update on x %%%%%%%%%%%%%%%%%%%%%%%%%% Update directiondel_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*epscount_delta_stop = count_delta_stop + 1;if count_delta_stop >= 500if verbosedisp('stuck in some corner');endbreak;endelsecount_delta_stop = 0;endold_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 measureerrorSteps(iter) = abs(norm(xk_1, 1) - norm(xG, 1)) ;epsSteps(iter) = epsilon;break;endtimeSteps(iter) = toc(t0) ;% MODIFIED -- measure the difference of the L1 norms, instead of the Euclidean distance%errorSteps(iter) = norm(xk_1-xG) ; % original measureerrorSteps(iter) = abs(norm(xk_1, 1) - norm(xG, 1)) ;epsSteps(iter) = epsilon;% compute stopping criteria and test for terminationkeep_going = true;switch stoppingCriterioncase STOPPING_GROUND_TRUTHkeep_going = norm(xk_1-xG)>tolerance;case STOPPING_SPARSE_SUPPORTif delta~=0nz_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 >= 1criterionActiveSet = num_changes_active / num_nz_x;keep_going = (criterionActiveSet > tolerance);endendcase STOPPING_DUALITY_GAPerror('Duality gap is not a valid stopping criterion for Homotopy.');case STOPPING_OBJECTIVE_VALUEif delta~=0% continue if not yeat reached target value tolAprev_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);endcase STOPPING_SUBGRADIENTkeep_going = norm(delta*del_x_vec)>tolerance;case STOPPING_TIMEkeep_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;% endif ~keep_goingbreak;endif ~isempty(out_x)% If an element is removed from gamma_xlen_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 columnsQ11 = 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_xgamma_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 columnsiA11 = 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;endz_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.edufunction [i_delta, delta, out_x] = update_primal(out_x)global N gamma_x z_x xk_temp del_x_vec pk_temp dk epsilon isNonnegativegamma_lc = setdiff(1:N, union(gamma_x, out_x));if isNonnegativedelta1 = inf;elsedelta1_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;endenddelta2_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;endif delta1>delta2delta = delta2;i_delta = gamma_lc(delta2_pos_ind(i_delta2));elsedelta = delta1;i_delta = gamma_lc(delta1_pos_ind(i_delta1));enddelta3_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 isNonnegativewrong_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径向基神经网络时序、回归预测和分类
844
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
