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_itersX_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') < tolfprintf('SVD迭代在第%d次收敛。\n', iter);break;endend
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_iterstotal_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) == 0fprintf('迭代 %d, 损失: %.6f\n', iter, total_loss);endend
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));% 将缺失值位置设为0X_known = X;X_known(missing_mask) = 0;for iter = 1:max_itersX_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) < tolfprintf('SVT在%d次迭代后收敛。\n', iter);break;endend
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
关键参数选择建议
-
目标秩k的选择：
- 如果知道真实秩，直接使用
- 否则可以使用特征值衰减观察法：
plot(svd(X_known))，寻找拐点
-
SVT方法参数：
tau：通常取矩阵大小的倍数，如tau = 5 * (m+n)/2delta：通常在1.0-1.5之间
-
矩阵分解参数：
- 学习率：0.001-0.1，需要尝试
- 正则化参数：0.01-0.1防止过拟合
实际应用技巧
- 数据预处理：对已知元素进行归一化可以显著提高性能
- 缺失模式：对于非均匀缺失的情况，可能需要更多迭代
- 收敛判断：可以观察重构误差的变化情况来调整迭代次数
这三种方法各有优劣，你可以根据具体问题选择最适合的方案。
