基于交替方向乘子法（ADMM）的RPCA MATLAB实现

function [L, S] = rpca_admm(X, lambda, rho, max_iter, tol)
    % RPCA的ADMM实现（低秩+稀疏分解）
    % 输入参数：
    % X - 输入矩阵 (m x n)
    % lambda - 稀疏正则化参数 (默认1/sqrt(max(size(X))))
    % rho - ADMM惩罚参数 (默认1.0)
    % max_iter - 最大迭代次数 (默认1000)
    % tol - 收敛容差 (默认1e-7)
    % 输出：
    % L - 低秩矩阵
    % S - 稀疏矩阵
​
    [m, n] = size(X);
    if nargin < 2 || isempty(lambda)
        lambda = 1 / sqrt(max(m, n));  % 经验公式end
    if nargin < 3 || isempty(rho)
        rho = 1.0;  % 默认惩罚参数end
    if nargin < 4 || isempty(max_iter)
        max_iter = 1000;
    end
    if nargin < 5 || isempty(tol)
        tol = 1e-7;
    end
​
    % 初始化变量
    L = zeros(m, n);
    S = zeros(m, n);
    Y = zeros(m, n);  % 拉格朗日乘子
    mu = 0.01 / norm(X, 'fro');  % 自适应步长mu_bar = mu * 1e7;  % 最大步长
    dual_norm = 0;
​
    % 主迭代循环
    for iter = 1:max_iter
        % 更新低秩矩阵L
        [U, SVD, V] = svd(X - S + Y/rho, 'econ');
        diagS = diag(SVD);
        diagS = max(diagS - 1/rho, 0);  % 奇异值收缩svp = length(find(diagS > 1e-10));  % 有效秩估计
        L = U(:, 1:svp) * diag(diagS(1:svp)) * V(:, 1:svp)';
        
        % 更新稀疏矩阵S
        residual = X - L + Y/rho;
        S = sign(residual) .* max(abs(residual) - lambda/rho, 0);
        
        % 更新拉格朗日乘子
        Y = Y + rho * (X - L - S);
        
        % 收敛性检查
        primal_residual = norm(X - L - S, 'fro');
        dual_residual = rho * norm(S - S_prev, 'fro');
        if iter > 1 && primal_residual < tol && dual_residual < tol
            break;
        end
        S_prev = S;
        
        % 自适应调整rhoif mod(iter, 10) == 0
            mu = min(mu * 1.2, mu_bar);
        end
    end
end

% 生成测试数据
[m, n] = (100, 150);
L = randn(m, 10)*0.1;  % 低秩成分
S = 5*ones(m, n);      % 稀疏噪声
X = L + S;
​
% 执行RPCA分解
tic;
[L_est, S_est] = rpca_admm(X, 1/sqrt(n), 1.2, 2000, 1e-6);
toc;
​
% 可视化结果
figure;
subplot(131); imagesc(L); title('真实低秩');
subplot(132); imagesc(S); title('真实稀疏');
subplot(133); imagesc(L_est); title('估计低秩');

% 加载图像
img = imread('lena.png');
img_gray = rgb2gray(img);
img_double = im2double(img_gray);
​
% 添加高斯噪声
noise_level = 0.2;
noisy_img = imnoise(img_double, 'gaussian', 0, noise_level);
​
% RPCA去噪
tic;
[low_rank, sparse_noise] = rpca_admm(noisy_img, 0.1, 1.5, 1500, 1e-7);
denoised_img = low_rank;
toc;
​
% 显示结果
figure;
subplot(131); imshow(noisy_img); title('含噪图像');
subplot(132); imshow(low_rank); title('RPCA去噪');
subplot(133); imshow(img_gray); title('原始图像');

X_gpu = gpuArray(X);
[L_gpu, S_gpu] = rpca_admm(X_gpu, ...);
L = gather(L_gpu);

parfor i = 1:block_num
    block = X((i-1)*block_size+1:i*block_size, :);
    [L_block, S_block] = rpca_admm(block, ...);
    L((i-1)*block_size+1:i*block_size, :) = L_block;
end

X_sparse = sparse(X);
[L_sparse, S_sparse] = rpca_admm(X_sparse, ...);

v = L + (iter/(iter+3))*(L - L_prev);

[U, ~, ~] = svd(X, 'econ'); X_reduced = U(:, 1:100) * X;