MATLAB实现低秩矩阵填充

MATLAB实现低秩矩阵填充，整合了三种经典的方法

方法名称	核心原理	适用场景	优点	缺点
基于SVD的迭代硬阈值	迭代过程中保留前k个奇异值，强制矩阵低秩	已知矩阵真实秩的情况；数据量不大时	原理简单，收敛性强	需要知道矩阵的秩，计算量大
基于梯度下降的矩阵分解	将矩阵分解为两个小矩阵，优化其乘积与已知元素的误差	大规模矩阵；推荐系统	不需要知道矩阵的秩，速度较快	容易陷入局部最优
基于凸优化的奇异值阈值	用核范数（奇异值之和）近似矩阵的秩，进行优化	矩阵秩未知；理论保证强的恢复	恢复性能好，理论完善	计算速度慢，迭代次数多

基于SVD的迭代硬阈值法

这种方法假设你已知矩阵的真实秩。它通过迭代的方式，在每一步都对矩阵进行SVD分解，只保留前k个最大的奇异值，从而强制矩阵保持低秩状态。

function X_hat = svd_iterative_impute(X, k, max_iters, tol)
    % X: 包含NaN值的观测矩阵
    % k: 假设的矩阵秩
    % max_iters: 最大迭代次数
    % tol: 收敛阈值
    % 初始化：用0填充缺失值
    X_hat = X;
    missing_mask = isnan(X);
    X_hat(missing_mask) = 0;
    
    for iter = 1:max_iters
        X_prev = X_hat;
        % SVD分解
        [U, S, V] = svd(X_hat, 'econ');
        % 硬阈值：只保留前k个奇异值
        S_k = S;
        S_k(k+1:end, k+1:end) = 0;
        % 重构矩阵
        X_hat = U * S_k * V';
        % 确保已知元素不变
        X_hat(~missing_mask) = X(~missing_mask);
        
        % 检查收敛
        if norm(X_hat - X_prev, 'fro') / norm(X_prev, 'fro') < tol
            fprintf('SVD迭代在第%d次收敛。\n', iter);
            break;
        end
    end
end

基于梯度下降的矩阵分解法

该方法将原始矩阵分解为两个小矩阵（用户矩阵和物品矩阵）的乘积，通过梯度下降最小化已知位置上的预测误差。这种方法在推荐系统中非常常用。

function [P, Q] = matrix_factorization_impute(X, k, max_iters, learning_rate, lambda)
    % X: 包含NaN的观测矩阵
    % k: 隐特征维度
    % max_iters: 最大迭代次数
    % learning_rate: 学习率
    % lambda: 正则化参数
    % 初始化P和Q
    [m, n] = size(X);
    P = rand(m, k);
    Q = rand(n, k);
    
    % 找到非NaN元素的位置和值
    [row, col, values] = find(~isnan(X));
    obs_idx = [row, col];
    
    fprintf('开始梯度下降...\n');
    for iter = 1:max_iters
        total_error = 0;
        % 随机打乱观测顺序
        idx = randperm(length(values));
        for i = 1:length(idx)
            r = row(idx(i));
            c = col(idx(i));
            true_val = X(r, c);
            
            % 预测值
            pred_val = P(r, :) * Q(c, :)';
            error = pred_val - true_val;
            total_error = total_error + error^2;
            
            % 梯度更新
            P_temp = P(r, :);
            P(r, :) = P(r, :) - learning_rate * (error * Q(c, :) + lambda * P(r, :));
            Q(c, :) = Q(c, :) - learning_rate * (error * P_temp + lambda * Q(c, :));
        end
        
        % 计算总损失
        total_loss = total_error / length(values) + (lambda/2) * (norm(P, 'fro') + norm(Q, 'fro'));
        if mod(iter, 100) == 0
            fprintf('迭代 %d, 损失: %.6f\n', iter, total_loss);
        end
    end
end

基于凸优化的奇异值阈值法

该方法不直接优化秩（这是个非凸问题），而是优化矩阵的核范数（所有奇异值之和），这是秩的凸近似。通过软阈值方法迭代地更新矩阵。

function X_hat = svt_impute(X, tau, delta, max_iters, tol)
    % X: 包含NaN的观测矩阵
    % tau: 阈值参数
    % delta: 步长
    % max_iters: 最大迭代次数
    % tol: 收敛阈值
    % 初始化
    X_hat = zeros(size(X));
    missing_mask = isnan(X);
    Y = zeros(size(X));
    
    % 将缺失值位置设为0
    X_known = X;
    X_known(missing_mask) = 0;
    
    for iter = 1:max_iters
        X_prev = X_hat;
        % SVD分解
        [U, S, V] = svd(Y, 'econ');
        % 软阈值操作
        S_soft = sign(S) .* max(abs(S) - tau, 0);
        % 更新矩阵估计
        X_hat = U * S_soft * V';
        % 更新Y（只修正已知元素）
        Y = Y + delta * (X_known - X_hat .* ~missing_mask);
        
        % 检查收敛
        if norm(X_hat - X_prev, 'fro') / (norm(X_prev, 'fro') + eps) < tol
            fprintf('SVT在%d次迭代后收敛。\n', iter);
            break;
        end
    end
end

完整演示与对比

下面是一个完整的演示代码，展示如何使用这三种方法并比较它们的效果：

% 生成模拟数据
m = 100; n = 100; % 矩阵维度
true_rank = 5; % 真实秩
A_true = randn(m, true_rank) * randn(true_rank, n); % 生成低秩矩阵

% 随机隐藏部分元素作为缺失值
obs_ratio = 0.3; % 观测比例
mask = rand(m, n) < obs_ratio;
X_obs = A_true;
X_obs(~mask) = NaN;

fprintf('矩阵维度: %d x %d, 真实秩: %d, 观测比例: %.2f\n', ...
        m, n, true_rank, obs_ratio);

% 方法1: 基于SVD的迭代硬阈值
fprintf('\n=== 方法1: SVD迭代硬阈值 ===\n');
tic;
X_svd = svd_iterative_impute(X_obs, true_rank, 500, 1e-6);
time_svd = toc;
error_svd = norm(X_svd - A_true, 'fro') / norm(A_true, 'fro');
fprintf('相对误差: %.6f, 耗时: %.2f秒\n', error_svd, time_svd);

% 方法2: 基于矩阵分解
fprintf('\n=== 方法2: 矩阵分解 ===\n');
tic;
[P, Q] = matrix_factorization_impute(X_obs, true_rank, 1000, 0.01, 0.01);
X_mf = P * Q';
time_mf = toc;
error_mf = norm(X_mf - A_true, 'fro') / norm(A_true, 'fro');
fprintf('相对误差: %.6f, 耗时: %.2f秒\n', error_mf, time_mf);

% 方法3: 奇异值阈值(SVT)
fprintf('\n=== 方法3: 奇异值阈值(SVT) ===\n');
tic;
X_svt = svt_impute(X_obs, 1, 1.2, 500, 1e-6);
time_svt = toc;
error_svt = norm(X_svt - A_true, 'fro') / norm(A_true, 'fro');
fprintf('相对误差: %.6f, 耗时: %.2f秒\n', error_svt, time_svt);

% 结果显示
fprintf('\n=== 结果对比 ===\n');
fprintf('方法                 相对误差          耗时(秒)\n');
fprintf('SVD迭代硬阈值       %.6f          %.2f\n', error_svd, time_svd);
fprintf('矩阵分解            %.6f          %.2f\n', error_mf, time_mf);
fprintf('SVT                 %.6f          %.2f\n', error_svt, time_svt);

参考代码 各种低秩约束矩阵填充方法 www.youwenfan.com/contentcnk/78960.html

关键参数选择建议

1. 目标秩k的选择：

   • 如果知道真实秩，直接使用
   • 否则可以使用特征值衰减观察法：plot(svd(X_known))，寻找拐点

2. SVT方法参数：

   • tau：通常取矩阵大小的倍数，如 tau = 5 * (m+n)/2
   • delta：通常在1.0-1.5之间

3. 矩阵分解参数：

   • 学习率：0.001-0.1，需要尝试
   • 正则化参数：0.01-0.1防止过拟合

实际应用技巧

1. 数据预处理：对已知元素进行归一化可以显著提高性能
2. 缺失模式：对于非均匀缺失的情况，可能需要更多迭代
3. 收敛判断：可以观察重构误差的变化情况来调整迭代次数

这三种方法各有优劣，你可以根据具体问题选择最适合的方案。