MPCA（Multilinear Principal Component Analysis）算法，这是一种专门用于张量数据降维的强大工具。
MPCA算法核心思想
与传统PCA的区别
| 特性 | 传统PCA | MPCA |
|---|---|---|
| 数据形式 | 向量（2阶张量） | 高阶张量 |
| 数据结构 | 忽略多维结构 | 保留张量结构 |
| 降维方式 | 单一投影矩阵 | 多个模式投影矩阵 |
MPCA数学原理
1. 张量基础概念
% 张量表示：X ∈ ℝ^{I₁ × I₂ × ... × Iₙ}
% 模式-n展开：将张量展开为矩阵
2. MPCA目标函数
MPCA的目标是找到一组投影矩阵 {U⁽¹⁾, U⁽²⁾, …, U⁽ⁿ⁾}，使得投影后的张量在降维的同时最大程度保留原始变异信息。
目标函数：
max{U(n)}∑i=1N∥Yi−Yˉ∥F2\max_{\{U^{(n)}\}} \sum_{i=1}^{N} \| \mathcal{Y}_i - \bar{\mathcal{Y}} \|_F^2max{U(n)}∑i=1N∥Yi−Yˉ∥F2
其中 Yi=Xi×1U(1)T×2U(2)T⋯×NU(N)T\mathcal{Y}_i = \mathcal{X}_i \times_1 U^{(1)T} \times_2 U^{(2)T} \cdots \times_N U^{(N)T}Yi=Xi×1U(1)T×2U(2)T⋯×NU(N)T
MATLAB实现
1. 核心MPCA类实现
classdef MPCA < handle
properties
% 投影矩阵
U_proj % 元胞数组，存储各模式投影矩阵
mean_tensor % 训练数据的均值张量
explained_variance % 各模式解释方差
original_dims % 原始维度
target_dims % 目标维度
end
methods
function obj = MPCA()
% MPCA构造函数
obj.U_proj = {};
obj.mean_tensor = [];
obj.explained_variance = [];
end
2. 张量工具函数
function X_mat = tensor_unfold(obj, X, mode)
% 张量模式展开
dims = size(X);
N = ndims(X);
% 重新排列维度
new_order = [mode, setdiff(1:N, mode)];
X_permuted = permute(X, new_order);
% 重塑为矩阵
X_mat = reshape(X_permuted, dims(mode), prod(dims)/dims(mode));
end
function X_tensor = tensor_fold(obj, X_mat, mode, target_dims)
% 矩阵折叠回张量
N = length(target_dims);
% 计算折叠后的维度
folded_dims = target_dims;
folded_dims(mode) = size(X_mat, 1);
% 重塑为张量
X_reshaped = reshape(X_mat, folded_dims);
% 逆排列恢复原始维度顺序
[~, inverse_order] = sort([mode, setdiff(1:N, mode)]);
X_tensor = permute(X_reshaped, inverse_order);
end
3. MPCA训练算法
function fit(obj, X_train, target_dims, varargin)
% MPCA训练函数
% 输入:
% X_train: 训练张量样本（维度: I₁ × I₂ × ... × Iₙ × N_samples）
% target_dims: 目标维度数组
% varargin: 可选参数（最大迭代次数、容忍度等）
p = inputParser;
addParameter(p, 'max_iter', 100, @isnumeric);
addParameter(p, 'tolerance', 1e-6, @isnumeric);
addParameter(p, 'verbose', true, @islogical);
parse(p, varargin{:});
max_iter = p.Results.max_iter;
tolerance = p.Results.tolerance;
verbose = p.Results.verbose;
% 获取张量维度信息
tensor_dims = size(X_train);
N_modes = length(tensor_dims) - 1; % 最后一个维度是样本数
num_samples = tensor_dims(end);
obj.original_dims = tensor_dims(1:end-1);
obj.target_dims = target_dims;
if verbose
fprintf('开始MPCA训练...\n');
fprintf('原始维度: %s\n', mat2str(obj.original_dims));
fprintf('目标维度: %s\n', mat2str(target_dims));
fprintf('样本数量: %d\n', num_samples);
end
% 步骤1: 中心化数据
obj.mean_tensor = mean(X_train, ndims(X_train));
X_centered = X_train - repmat(obj.mean_tensor, [ones(1, N_modes), num_samples]);
% 步骤2: 初始化投影矩阵
obj.U_proj = cell(1, N_modes);
for n = 1:N_modes
obj.U_proj{n} = randn(tensor_dims(n), target_dims(n));
% 正交化
[obj.U_proj{n}, ~] = qr(obj.U_proj{n}, 0);
end
% 步骤3: 交替优化
prev_total_variance = 0;
for iter = 1:max_iter
if verbose && mod(iter, 10) == 0
fprintf('迭代 %d/%d...\n', iter, max_iter);
end
% 对每个模式交替优化
for n = 1:N_modes
% 构造辅助张量（固定其他模式）
X_projected = obj.partial_project(X_centered, n);
% 模式-n展开
Xn = obj.tensor_unfold(X_projected, n);
% 计算协方差矩阵
cov_matrix = Xn * Xn';
% 特征分解
[U, S, ~] = svd(cov_matrix, 'econ');
% 选择前target_dims(n)个特征向量
obj.U_proj{n} = U(:, 1:target_dims(n));
end
% 计算当前总方差
current_variance = obj.calculate_total_variance(X_centered);
% 检查收敛
if iter > 1 && abs(current_variance - prev_total_variance) < tolerance
if verbose
fprintf('在迭代 %d 收敛\n', iter);
end
break;
end
prev_total_variance = current_variance;
end
% 计算解释方差
obj.calculate_explained_variance(X_centered);
if verbose
fprintf('MPCA训练完成\n');
obj.display_results();
end
end
4. 投影与重构函数
function Y = transform(obj, X)
% 将新张量投影到低维空间
% X: 输入张量（可以是单个张量或张量样本集）
if ndims(X) == length(obj.original_dims)
% 单个张量
X_centered = X - obj.mean_tensor;
Y = X_centered;
for n = 1:length(obj.U_proj)
Y = obj.mode_n_product(Y, obj.U_proj{n}', n);
end
else
% 张量样本集
num_samples = size(X, ndims(X));
Y_dims = [obj.target_dims, num_samples];
Y = zeros(Y_dims);
for i = 1:num_samples
% 提取单个样本
subs = repmat({':'}, 1, ndims(X)-1);
X_i = X(subs{:}, i);
% 中心化并投影
X_centered = X_i - obj.mean_tensor;
Y_i = X_centered;
for n = 1:length(obj.U_proj)
Y_i = obj.mode_n_product(Y_i, obj.U_proj{n}', n);
end
% 存储结果
Y(subs{:}, i) = Y_i;
end
end
end
function X_recon = inverse_transform(obj, Y)
% 从低维表示重构原始张量
if ndims(Y) == length(obj.target_dims)
% 单个低维张量
X_recon = Y;
for n = length(obj.U_proj):-1:1
X_recon = obj.mode_n_product(X_recon, obj.U_proj{n}, n);
end
X_recon = X_recon + obj.mean_tensor;
else
% 低维张量样本集
num_samples = size(Y, ndims(Y));
X_dims = [obj.original_dims, num_samples];
X_recon = zeros(X_dims);
for i = 1:num_samples
subs = repmat({':'}, 1, ndims(Y)-1);
Y_i = Y(subs{:}, i);
X_recon_i = Y_i;
for n = length(obj.U_proj):-1:1
X_recon_i = obj.mode_n_product(X_recon_i, obj.U_proj{n}, n);
end
X_recon_i = X_recon_i + obj.mean_tensor;
X_recon(subs{:}, i) = X_recon_i;
end
end
end
参考代码 MPCA算法，用于张量的降维 www.youwenfan.com/contentcsl/60159.html
应用示例：人脸图像数据集
1. 数据准备与预处理
% 模拟人脸图像数据（3阶张量：高度×宽度×样本数）
function [tensor_data, labels] = generate_face_data(num_samples, image_size)
% 生成模拟人脸数据
tensor_data = zeros([image_size, num_samples]);
labels = randi([1, 5], num_samples, 1); % 5个不同的人
for i = 1:num_samples
% 生成基础人脸模板（高斯分布模拟面部特征）
base_face = generate_base_face(image_size, labels(i));
% 添加个体变异和噪声
individual_variation = 0.1 * randn(image_size);
noise = 0.05 * randn(image_size);
tensor_data(:, :, i) = base_face + individual_variation + noise;
end
end
function base_face = generate_base_face(image_size, person_id)
[X, Y] = meshgrid(1:image_size(2), 1:image_size(1));
center = image_size / 2;
% 根据不同的人ID生成不同的面部特征
switch person_id
case 1
% 圆脸
base_face = exp(-((X-center(2)).^2 + (Y-center(1)).^2) / (0.3*min(image_size))^2);
case 2
% 椭圆脸
base_face = exp(-((X-center(2)).^2/0.4 + (Y-center(1)).^2/0.6) / (0.3*min(image_size))^2);
otherwise
base_face = exp(-((X-center(2)).^2 + (Y-center(1)).^2) / (0.35*min(image_size))^2);
end
end
2. MPCA降维与可视化
% 主演示程序
function demonstrate_MPCA()
clear; close all; clc;
fprintf('=== MPCA算法演示 ===\n');
% 生成模拟数据
image_size = [64, 64]; % 图像尺寸
num_samples = 200;
num_test = 50;
[train_data, train_labels] = generate_face_data(num_samples, image_size);
[test_data, test_labels] = generate_face_data(num_test, image_size);
fprintf('数据维度: %s\n', mat2str(size(train_data)));
% MPCA降维
target_dims = [16, 16]; % 目标维度
mpca = MPCA();
mpca.fit(train_data, target_dims, 'max_iter', 50, 'verbose', true);
% 训练数据投影
train_lowdim = mpca.transform(train_data);
fprintf('降维后维度: %s\n', mat2str(size(train_lowdim)));
% 测试数据投影
test_lowdim = mpca.transform(test_data);
% 重构测试样本
test_recon = mpca.inverse_transform(test_lowdim);
% 可视化结果
visualize_mpca_results(train_data, test_data, test_recon, train_lowdim, mpca);
% 性能评估
evaluate_mpca_performance(train_data, test_data, test_recon, mpca);
end
3. 结果可视化函数
function visualize_mpca_results(original_train, original_test, reconstructed, lowdim, mpca)
figure('Position', [100, 100, 1400, 1000]);
% 1. 原始图像 vs 重构图像
subplot(2,3,1);
sample_idx = 1;
original_img = original_test(:,:,sample_idx);
recon_img = reconstructed(:,:,sample_idx);
montage({original_img, recon_img}, 'Size', [1, 2]);
title('原始图像 vs 重构图像');
% 2. 重构误差分布
subplot(2,3,2);
reconstruction_errors = zeros(size(original_test, 3), 1);
for i = 1:size(original_test, 3)
diff = original_test(:,:,i) - reconstructed(:,:,i);
reconstruction_errors(i) = norm(diff(:)) / norm(original_test(:,:,i)(:));
end
histogram(reconstruction_errors, 20);
xlabel('相对重构误差');
ylabel('频数');
title('重构误差分布');
grid on;
% 3. 低维空间可视化（前2个维度）
subplot(2,3,3);
if ndims(lowdim) == 3
lowdim_2d = reshape(lowdim, [], size(lowdim, 3));
[coeff, score] = pca(lowdim_2d');
scatter(score(:,1), score(:,2), 30, 'filled');
xlabel('第一主成分');
ylabel('第二主成分');
title('低维空间中的样本分布');
grid on;
end
% 4. 解释方差
subplot(2,3,4);
if ~isempty(mpca.explained_variance)
bar(mpca.explained_variance);
xlabel('模式');
ylabel('解释方差比例');
title('各模式解释方差');
grid on;
end
% 5. 投影矩阵可视化
subplot(2,3,5);
U1 = mpca.U_proj{1};
imagesc(U1(:,1:min(10, size(U1,2))));
colorbar;
title('模式-1投影矩阵（前10列）');
xlabel('投影方向');
ylabel('原始维度');
% 6. 压缩比计算
subplot(2,3,6);
original_size = prod(mpca.original_dims);
compressed_size = prod(mpca.target_dims);
compression_ratio = original_size / compressed_size;
labels = {'原始数据', '压缩后数据'};
sizes = [original_size, compressed_size] / 1e3; % 转换为千单位
bar(sizes);
set(gca, 'XTickLabel', labels);
ylabel('数据量 (千元素)');
title(sprintf('压缩比: %.2f:1', compression_ratio));
grid on;
end
性能评估与比较
1. MPCA性能评估
function evaluate_mpca_performance(train_data, test_data, test_recon, mpca)
fprintf('\n=== MPCA性能评估 ===\n');
% 计算重构误差
mse_errors = zeros(size(test_data, 3), 1);
for i = 1:size(test_data, 3)
mse_errors(i) = mean((test_data(:,:,i) - test_recon(:,:,i)).^2, 'all');
end
fprintf('平均重构MSE: %.4f\n', mean(mse_errors));
fprintf('重构MSE标准差: %.4f\n', std(mse_errors));
% 计算压缩比
original_size = prod(mpca.original_dims);
compressed_size = prod(mpca.target_dims);
compression_ratio = original_size / compressed_size;
fprintf('压缩比: %.2f:1\n', compression_ratio);
% 计算信息保留率
total_variance_original = sum(var(reshape(train_data, [], size(train_data,3)), 2));
train_lowdim = mpca.transform(train_data);
total_variance_compressed = sum(var(reshape(train_lowdim, [], size(train_lowdim,3)), 2));
variance_retained = total_variance_compressed / total_variance_original;
fprintf('方差保留率: %.2f%%\n', variance_retained * 100);
end
2. 与传统PCA比较
function compare_with_PCA(tensor_data, target_dims)
% 将张量数据展平为向量
[d1, d2, num_samples] = size(tensor_data);
X_flat = reshape(tensor_data, d1*d2, num_samples)';
% 传统PCA
[coeff, score, ~, ~, explained] = pca(X_flat);
% 选择主成分数量（达到类似压缩比）
target_pc = round(prod(target_dims) * 1.2); % 稍多维度以公平比较
% PCA重构
X_recon_pca = score(:,1:target_pc) * coeff(:,1:target_pc)';
X_recon_pca = reshape(X_recon_pca', d1, d2, num_samples);
% MPCA重构
mpca = MPCA();
mpca.fit(tensor_data, target_dims, 'max_iter', 50, 'verbose', false);
X_recon_mpca = mpca.inverse_transform(mpca.transform(tensor_data));
% 比较重构质量
mse_pca = mean((tensor_data - X_recon_pca).^2, 'all');
mse_mpca = mean((tensor_data - X_recon_mpca).^2, 'all');
fprintf('\n=== MPCA vs PCA 比较 ===\n');
fprintf('PCA重构MSE: %.4f\n', mse_pca);
fprintf('MPCA重构MSE: %.4f\n', mse_mpca);
fprintf('MPCA相对改进: %.2f%%\n', (mse_pca - mse_mpca)/mse_pca*100);
end
MPCA算法优势与应用场景
主要优势
- 结构保持：保留张量的多维结构信息
- 维度灵活性：可独立控制每个模式的降维程度
- 计算效率：相比将张量展平后使用PCA，计算更高效
- 解释性强：每个模式的投影矩阵具有明确的物理意义
典型应用场景
| 应用领域 | 张量结构 | MPCA作用 |
|---|---|---|
| 人脸识别 | 图像×光照×姿态 | 特征提取，降维 |
| 视频分析 | 帧×行×列×时间 | 时空特征提取 |
| 神经科学 | 脑区×时间×被试 | 多被试脑信号分析 |
| 推荐系统 | 用户×物品×时间 | 多维偏好建模 |
这个MPCA实现提供了张量降维的强大工具，特别适合处理具有多维结构的数据。
扫一扫
1559
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
