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_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

关键参数选择建议

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. 收敛判断：可以观察重构误差的变化情况来调整迭代次数

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