Alternating least squares for CANDECOMP/PARAFAC (CP) Decomposition

函数 cp_als 使用著名的交替最小二乘算法（alternating least-squares algorithm）计算张量 $X$ 的最佳 rank-R CP 模型估计值（更多信息请参见 Kolda and Bader, SIAM Review, 2009）。输入 $\mathbf{X}$ 可以是几乎任何类型的张量，包括 tensor、sptensor、ktensor 或 ttensor。输出 CP 模型是一个 ktensor。

请注意，需要较新版本的 MATLAB。下载或克隆 Tensor Toolbox 后，请使用 addpath 将 tensor_toolbox 目录添加到您的路径中。键入 help tensor_toolbox 可查看工具箱提供的函数和类的列表。键入 doc tensor_toolbox 以访问此文档。更多信息请参见 README.md。

加载一些数据

我们使用了来自 Andersson and Bro 的著名氨基酸数据集。它包含了 5 个样本的荧光测量值，这些样本含有 3 种氨基酸：色氨酸 (Tryptophan)、酪氨酸 (Tyrosine) 和苯丙氨酸 (Phenylalanine)。每种氨基酸对应一个秩一分量 (rank-one component)。该张量的大小为 5 x 51 x 201，来自 5 个样本、51 个激发波长和 201 个发射波长。更多详情可在此处找到：http://www.models.life.ku.dk/Amino_Acid_fluo。请引用以下论文来获取此数据：Rasmus Bro, PARAFAC: Tutorial and applications, Chemometrics and Intelligent Laboratory Systems, 1997, 38, 149-171。该数据集可以在 doc 目录中找到。

load aminoacids

rng('default') %<- 设置随机种子以保证此脚本的可复现性
M1 = cp_als(X,3); %<- 调用该方法

方法的基本调用，指定数据张量及其秩

这将使用一个随机的初始猜测。在每次迭代中，它会报告‘拟合度’ (fit)，其定义为 $1-(\text{norm}(X-M)/\text{norm}(X))$ ，这大致是 CP 模型所描述的数据比例，即拟合度为 $1$ 时表示完美拟合。

CP_ALS (CP Alternating Least Squares):Tensor size: [5 201 61]Tensor type: tensorR = 3, maxiters = 50, tol = 1.000000e-04dimorder = [1 2 3]init = randomIter  1: f = 5.191033e-01 f-delta = 5.2e-01Iter  2: f = 6.617871e-01 f-delta = 1.4e-01Iter  3: f = 7.615352e-01 f-delta = 1.0e-01Iter  4: f = 8.286421e-01 f-delta = 6.7e-02Iter  5: f = 8.813646e-01 f-delta = 5.3e-02Iter  6: f = 9.349259e-01 f-delta = 5.4e-02Iter  7: f = 9.592629e-01 f-delta = 2.4e-02Iter  8: f = 9.671789e-01 f-delta = 7.9e-03Iter  9: f = 9.702331e-01 f-delta = 3.1e-03Iter 10: f = 9.715826e-01 f-delta = 1.3e-03Iter 11: f = 9.722850e-01 f-delta = 7.0e-04Iter 12: f = 9.727265e-01 f-delta = 4.4e-04Iter 13: f = 9.730506e-01 f-delta = 3.2e-04Iter 14: f = 9.733119e-01 f-delta = 2.6e-04Iter 15: f = 9.735326e-01 f-delta = 2.2e-04Iter 16: f = 9.737227e-01 f-delta = 1.9e-04Iter 17: f = 9.738876e-01 f-delta = 1.6e-04Iter 18: f = 9.740309e-01 f-delta = 1.4e-04Iter 19: f = 9.741554e-01 f-delta = 1.2e-04Iter 20: f = 9.742635e-01 f-delta = 1.1e-04Iter 21: f = 9.743572e-01 f-delta = 9.4e-05Final f = 9.743572e-01 

• Iter 1 to Iter 21:
  • 这部分显示了算法的每一次迭代。算法总共运行了21次迭代。
• f = …:
  • 这是拟合度（fit）值。它衡量了分解得到的模型M对原始数据X的解释程度，通常定义为1-norm(X-M)/norm(X)。这个值越接近1，表示模型对数据的拟合越好。
• f-delta = …:
  • 这是拟合度的变化量，即当前迭代的f值与上一次迭代f值的差。例如，在Iter2中，f-delta≈0.6617-0.5191=0.1426。

我们通常可以达到 $f=0.97$ 的最终拟合度。当拟合度的变化小于指定的容差时，该方法停止，默认容差为 $1e-4$。

可视化结果

使用 ktensor/viz 函数来可视化结果。

vizopts = {'PlotCommands',{'bar','line','line'},...'ModeTitles',{'Concentration','Emission','Excitation'},...'BottomSpace',0.10,'HorzSpace',0.04,'Normalize',0};
info1 = viz(M1,'Figure',1,vizopts{:});

• 'ModeTitles', {'Concentration', 'Emission', 'Excitation'} 这为三个模式（维度）的图表设置了自定义标题。
• 'HorzSpace', 0.04，设置子图 (subplots) 之间的水平间距为 4% 的窗口宽度，用于调整布局。

增加最大迭代次数

请注意，前一次运行在达到指定的收敛容差之前，仅迭代 50 次就退出了。让我们增加最大迭代次数，并使用相同的初始猜测再次尝试。

M2 = cp_als(X, 3, 'maxiters', 100, 'init', U2);

CP_ALS (CP Alternating Least Squares):Tensor size: [5 201 61]Tensor type: tensorR = 3, maxiters = 100, tol = 1.000000e-04dimorder = [1 2 3]init = user-specifiedIter  1: f = 4.656532e-01 f-delta = 4.7e-01Iter  2: f = 6.849381e-01 f-delta = 2.2e-01Iter  3: f = 7.873635e-01 f-delta = 1.0e-01Iter  4: f = 9.084899e-01 f-delta = 1.2e-01Iter  5: f = 9.477155e-01 f-delta = 3.9e-02Iter  6: f = 9.633209e-01 f-delta = 1.6e-02Iter  7: f = 9.692590e-01 f-delta = 5.9e-03Iter  8: f = 9.715954e-01 f-delta = 2.3e-03Iter  9: f = 9.726257e-01 f-delta = 1.0e-03Iter 10: f = 9.731706e-01 f-delta = 5.4e-04Iter 11: f = 9.735185e-01 f-delta = 3.5e-04Iter 12: f = 9.737732e-01 f-delta = 2.5e-04Iter 13: f = 9.739743e-01 f-delta = 2.0e-04Iter 14: f = 9.741386e-01 f-delta = 1.6e-04Iter 15: f = 9.742746e-01 f-delta = 1.4e-04Iter 16: f = 9.743880e-01 f-delta = 1.1e-04Iter 17: f = 9.744824e-01 f-delta = 9.4e-05Final f = 9.744824e-01 

更改输出频率

使用 'printitn' 选项来更改输出频率。

M2alt2 = cp_als(X, 3, 'maxiters', 100, 'init', U2, 'printitn', 10);

CP_ALS (CP Alternating Least Squares):Tensor size: [5 201 61]Tensor type: tensorR = 3, maxiters = 100, tol = 1.000000e-04dimorder = [1 2 3]init = user-specifiedIter 10: f = 9.731706e-01 f-delta = 5.4e-04Iter 17: f = 9.744824e-01 f-delta = 9.4e-05Final f = 9.744824e-01 

更改容差

也可以放宽或收紧拟合度变化的容差。您可能需要增加迭代次数以使其收敛。

rng('default') % Reproducibility
M5 = cp_als(X,3,'init','nvecs','tol',1e-6,'maxiters',1000,'printitn',10);
CP_ALS (CP Alternating Least Squares):

 Tensor size: [5 201 61]Tensor type: tensorR = 3, maxiters = 1000, tol = 1.000000e-06dimorder = [1 2 3]init = nvecsIter 10: f = 9.334888e-01 f-delta = 3.5e-03Iter 20: f = 9.604549e-01 f-delta = 1.9e-03Iter 30: f = 9.712518e-01 f-delta = 5.8e-04Iter 40: f = 9.741285e-01 f-delta = 1.4e-04Iter 50: f = 9.747733e-01 f-delta = 2.9e-05Iter 60: f = 9.749128e-01 f-delta = 6.4e-06Iter 70: f = 9.749430e-01 f-delta = 1.4e-06Iter 73: f = 9.749461e-01 f-delta = 8.8e-07Final f = 9.749461e-01 

控制因子矩阵的符号模糊性

cp_als 的默认行为是调用 fixsigns 来固定因子矩阵的符号模糊性。您可以在调用 cp_als 时通过传递 fixsigns 参数值为 false 来关闭此行为。

Y = ktensor({[1; 1], [1; 1], [1; -10]});
X = tensor(Y);
M = cp_als(X, 2, 'printitn', 0, 'init', Y); % 默认行为，调用 fixsigns
M = cp_als(X, 2, 'printitn', 0, 'init', Y, 'fixsigns', false); % -fixsigns 未被调用

• ​​'fixsigns', false​​:这会​​禁用​​自动的符号标准化步骤。算法会返回它找到的"原始"解，而不进行任何符号调整。

时间和收敛轨迹

您还可以通过传递 'trace' 参数值为 true 来对方法计时。这将在输出结构体中返回运行该方法所花费的时间，以及每次迭代的拟合度（相对误差是 $1$ 减去拟合度）。

rng('default') % 可复现性
[M6, ~, output6] = cp_als(X, 3, 'init', 'nvecs', 'trace', true, 'printitn', 10);
figure(3); clf;
plot(output6.time_trace, 1-output6.fit_trace, '-o');
title('相对误差 vs. 时间');
xlabel('时间 (秒)');
ylabel('相对误差');

CP_ALS (CP Alternating Least Squares):Tensor size: [5 201 61]Tensor type: tensorR = 3, maxiters = 50, tol = 1.000000e-04dimorder = [1 2 3]init = nvecsIter 10: f = 9.334888e-01 f-delta = 3.5e-03Iter 20: f = 9.604549e-01 f-delta = 1.9e-03Iter 30: f = 9.712518e-01 f-delta = 5.8e-04Iter 40: f = 9.741285e-01 f-delta = 1.4e-04Iter 43: f = 9.744312e-01 f-delta = 8.6e-05Final f = 9.744312e-01