% 生成 3阶测试张量
clear;
start_t = datetime('now');
X = create_symmetric_tensor(5,3);
% 设置 TR 秩 (注意首尾秩需相等)
ranks = [2,3,3]; % ranks(1) = ranks(end)
% 正确参数设置
params = struct(...
'maxiters', 200, ... % 最大迭代次数
'tol', 1e-8, ... % 收敛阈值
'verbose', true, ... % 显示迭代信息
'conv_crit', 'relative error' ... % 收敛准则
);
% 执行分解
[cores, conv_history] = tr_als(X, ranks, params);
cores1=cat(3,cores{1,:});
cores2=cat(3,cores{2,:});
cores3=cat(3,cores{3,:});
% 验证结果
reconstructed = cores_2_tensor(cores);
error = norm(X(:) - reconstructed(:)) / norm(X(:));
disp(['相对重构误差: ', num2str(error)]);
% 绘制图形
duration = milliseconds(datetime('now') - start_t)function [cores, varargout] = tr_als(X, ranks, varargin)
%tr_als Compute tensor ring decomposition via alternating least squares
%
%cores = tr_als(X, ranks) computes a tensor ring (TR) decomposition of the
%input N-dimensional array X. ranks is a length-N vector containing the
%target ranks. The output cores is a cell containing the N cores tensors,
%each represented as a 3-way array.
%
%cores = tr_als(___, 'tol', tol) is an optional argument that controls the
%termination tolerance. If the change in the relative error is less than
%tol at the conclusion of a main loop iteration, the algorithm terminates.
%Default is 1e-3.
%
%cores = tr_als(___, 'maxiters', maxiters) is an optional argument that
%controls the maximum number of main loop iterations. The default is 50.
%
%cores = tr_als(___, 'verbose', verbose) is an optional argument that
%controls the amount of information printed to the terminal during
%execution. Setting verbose to true will result in more print out. Default
%is false.
%
%This algorithm is based on Alg. 2 in paper by
%Q. Zhao, G. Zhou, S. Xie, L. Zhang, A. Cichocki. "Tensor ring
%decomposition". arXiv:1606.05535, 2016.
%% Handle optional inputs
params = inputParser;
addParameter(params, 'conv_crit', 'none');
addParameter(params, 'tol', 1e-3, @isscalar);
addParameter(params, 'maxiters', 50, @(x) isscalar(x) & x > 0);
addParameter(params, 'verbose', false, @isscalar);
parse(params, varargin{:});
conv_crit = params.Results.conv_crit;
tol = params.Results.tol;
maxiters = params.Results.maxiters;
verbose = params.Results.verbose;
%% Initialize cores
sz = size(X);
cores = initialize_cores(sz, ranks);
if nargout > 1 && tol > 0 && (strcmp(conv_crit, 'relative error') || strcmp(conv_crit, 'norm'))
conv_vec = zeros(1, maxiters);
end
%% Main loop
% Iterate until convergence, for a maximum of maxiters iterations
N = length(sz);
er_old = Inf;
for it = 1:maxiters
for n = 1:N
% Compute G_{[2]}^{\neq n} from cores
G = subchain_matrix(cores, n);
% Compute relevant unfolding of data tensor X
XnT = mode_unfolding(X, n).';
% Solve LS problem and update core
Z = (G \ XnT).';
cores{n} = classical_mode_folding(Z, 2, size(cores{n}));
end
% Check convergence: Relative error
if tol > 0 && strcmp(conv_crit, 'relative error')
% Compute full tensor corresponding to cores
Y = cores_2_tensor(cores);
% Compute current relative error
er = norm(X(:)-Y(:))/norm(X(:));
if verbose
fprintf('\tRelative error after iteration %d: %.8f\n', it, er);
end
% Save current error to conv_vec if required
if nargout > 1
conv_vec(it) = er;
end
% Break if change in relative error below threshold
if abs(er - er_old) < tol
if verbose
fprintf('\tRelative error change below tol; terminating...\n');
end
break
end
% Update old error
er_old = er;
% Check convergence: Norm change
elseif tol > 0 && strcmp(conv_crit, 'norm')
% Compute norm of TR tensor and change if it > 1
norm_new = normTR(cores);
if it == 1
norm_change = Inf;
else
norm_change = abs(norm_new - norm_old);
if verbose
fprintf('\tNorm change after iteration %d: %.8f\n', it, norm_change);
end
end
% Save current norm_change to conv_vec
if nargout > 1
conv_vec(it) = norm_change;
end
% Break if change in relative error below threshold
if norm_change < tol
if verbose
fprintf('\tNorm change below tol; terminating...\n');
end
break
end
% Update old norm
norm_old = norm_new;
% Just print iteration count
else
if verbose
fprintf('\tIteration %d complete\n', it);
end
end
end
if nargout > 1 && exist('conv_vec', 'var')
varargout{1} = conv_vec(1:it);
else
varargout{1} = nan;
end
end
function cores = initialize_cores(sz, ranks, varargin)
%initialize_cores Initializes cores using Gaussian distribution
%
%cores = initialize_cores(sz, ranks) returns a length-N cell with N cores,
%each with entires drawn iid from the standard Gaussian distribution. sz is
%a length-N vector with the sizes, and ranks is a length-N vector with the
%outgoing ranks.
%
%cores = initialize_cores(___, 'init_zero', init_zero) will just initialize
%all cores to zero if init_zero is true. Default is false.
% Handle optional inputs
params = inputParser;
addParameter(params, 'init_zero', false, @isscalar);
parse(params, varargin{:});
init_zero = params.Results.init_zero;
% Main code
N = length(sz);
cores = cell(N,1);
for n = 1:N
R1 = ranks(n);
if n == 1
R0 = ranks(end);
else
R0 = ranks(n-1);
end
if init_zero
cores{n} = zeros(R0, sz(n), R1);
else
cores{n} = randn(R0, sz(n), R1);
end
end
end,修改代码,使得初始点是固定的一个随机种子生成,而不是高斯分布生成的
最新发布