Symmetric Tucker Decomposition

The function tucker_sym computes the best symmetric rank-(R,R,...,R) approximation of a symmetric tensor S, according to the specified rank R. The result returned in T is a ttensor with all factors equal. S must be a tensor (not a symtensor) and must be exactly symmetric.

The method is based on the algorithm described in: * Phillip A. Regalia, Monotonically Convergent Algorithms for Symmetric Tensor Approximation, Linear Algebra and its Applications 438(2):875-890, 2013, http://dx.doi.org/10.1016/j.laa.2011.10.033

Contents

Create a symmetric tensor of size 4x4x4

rng(0); %<-- Set seed for reproducibility
S = tensor(@randn,[4,4,4]);
S = symmetrize(S);

Compute a symmetric Tucker decomposition with rank 2 with default parameters

Except for the rank, all parameters are optional. We specify here that we want to print the fit every 10 iterations. The usual default is every iteration. By default, the initial guess is random. Since this is a nonconvex problem, we will run the algorithm multiple times to find the best solution.

best_err = inf;
best_T = [];
errs = zeros(3,1); % Store errors for the first 3 examples
for i = 1:10
	rng(i); % Set different seed each time
	T_tmp = tucker_sym(S,2,'printitn',10);
	err = norm(S - full(T_tmp)) / norm(S);
	if err < best_err
		best_err = err;
		best_T = T_tmp;
	end
end
errs(1) = best_err;
fprintf('\n\n*** Best relative error over 10 runs: %g ***\n', best_err);
disp(best_T)
Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: random
Return type: ttensor

 Iter 10: rel. change in X = 1.993136e+00
 Iter 20: rel. change in X = 2.367469e-05
 Iter 30: rel. change in X = 2.961379e-07
 Iter 40: rel. change in X = 3.704293e-09
 Iter 49: rel. change in X = 7.181217e-11
Final relative error: 8.334750e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: random
Return type: ttensor

 Iter 10: rel. change in X = 3.968407e-04
 Iter 20: rel. change in X = 8.389887e-07
 Iter 30: rel. change in X = 5.553515e-09
 Iter 37: rel. change in X = 8.519612e-11
Final relative error: 7.961686e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: random
Return type: ttensor

 Iter 10: rel. change in X = 4.967143e-05
 Iter 20: rel. change in X = 2.022365e-09
 Iter 23: rel. change in X = 9.745062e-11
Final relative error: 8.198031e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: random
Return type: ttensor

 Iter 10: rel. change in X = 1.057355e-03
 Iter 20: rel. change in X = 7.013479e-07
 Iter 30: rel. change in X = 2.409872e-10
 Iter 32: rel. change in X = 4.545516e-11
Final relative error: 7.378318e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: random
Return type: ttensor

 Iter 10: rel. change in X = 6.805640e-04
 Iter 20: rel. change in X = 9.664977e-08
 Iter 29: rel. change in X = 9.171225e-11
Final relative error: 6.765288e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: random
Return type: ttensor

 Iter 10: rel. change in X = 1.773809e-03
 Iter 20: rel. change in X = 1.884608e-06
 Iter 30: rel. change in X = 1.624699e-09
 Iter 34: rel. change in X = 8.947762e-11
Final relative error: 6.766387e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: random
Return type: ttensor

 Iter 10: rel. change in X = 3.004369e-01
 Iter 20: rel. change in X = 2.070062e-03
 Iter 30: rel. change in X = 3.471168e-06
 Iter 40: rel. change in X = 5.826140e-09
 Iter 47: rel. change in X = 6.649890e-11
Final relative error: 6.663817e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: random
Return type: ttensor

 Iter 10: rel. change in X = 7.402992e-02
 Iter 20: rel. change in X = 7.122327e-04
 Iter 30: rel. change in X = 4.405379e-06
 Iter 40: rel. change in X = 2.716151e-08
 Iter 50: rel. change in X = 1.674617e-10
 Iter 52: rel. change in X = 6.052126e-11
Final relative error: 7.891046e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: random
Return type: ttensor

 Iter 10: rel. change in X = 2.381543e-02
 Iter 20: rel. change in X = 1.025224e-04
 Iter 30: rel. change in X = 3.346359e-07
 Iter 40: rel. change in X = 1.091047e-09
 Iter 45: rel. change in X = 6.229899e-11
Final relative error: 8.071541e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: random
Return type: ttensor

 Iter 10: rel. change in X = 1.144115e-02
 Iter 20: rel. change in X = 3.687818e-04
 Iter 30: rel. change in X = 1.248129e-05
 Iter 40: rel. change in X = 4.411647e-07
 Iter 50: rel. change in X = 1.606548e-08
 Iter 60: rel. change in X = 5.909997e-10
 Iter 66: rel. change in X = 6.653694e-11
Final relative error: 7.406452e-01


*** Best relative error over 10 runs: 0.666382 ***
ans is a ttensor of size 4 x 4 x 4
	ans.core is a tensor of size 2 x 2 x 2
		ans.core(:,:,1) = 
	   -0.8818    1.4306
	    1.4306   -0.2805
		ans.core(:,:,2) = 
	    1.4306   -0.2805
	   -0.2805   -2.2030
	ans.U{1} = 
		   -0.5433    0.5931
		   -0.5974   -0.2862
		    0.0982    0.7404
		   -0.5817   -0.1350
	ans.U{2} = 
		   -0.5433    0.5931
		   -0.5974   -0.2862
		    0.0982    0.7404
		   -0.5817   -0.1350
	ans.U{3} = 
		   -0.5433    0.5931
		   -0.5974   -0.2862
		    0.0982    0.7404
		   -0.5817   -0.1350

We do the same test again, but provide user-specified initial guesses

best_err = inf;
best_T = [];
for i = 1:10
	rng(i); % Set different seed each time
	X0 = randn(4,2); % Initial guess for factor matrix
	T_tmp = tucker_sym(S,2,'init',X0,'printitn',10);
	err = norm(S - full(T_tmp)) / norm(S);
	if err < best_err
		best_err = err;
		best_T = T_tmp;
	end
end
errs(2) = best_err;
fprintf('\n\n*** Best relative error over 10 runs (user init): %g ***\n', best_err);
disp(best_T)
Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: user-provided
Return type: ttensor

 Iter 10: rel. change in X = 4.119526e-02
 Iter 20: rel. change in X = 3.047757e-02
 Iter 30: rel. change in X = 5.084022e-02
 Iter 40: rel. change in X = 1.702917e-02
 Iter 50: rel. change in X = 1.477118e-04
 Iter 60: rel. change in X = 1.028134e-06
 Iter 70: rel. change in X = 7.212640e-09
 Iter 79: rel. change in X = 8.319420e-11
Final relative error: 8.638623e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: user-provided
Return type: ttensor

 Iter 10: rel. change in X = 1.599719e-03
 Iter 20: rel. change in X = 9.599357e-04
 Iter 30: rel. change in X = 1.058735e-03
 Iter 40: rel. change in X = 2.205713e-03
 Iter 50: rel. change in X = 1.936684e-02
 Iter 60: rel. change in X = 2.020894e-03
 Iter 70: rel. change in X = 4.614973e-10
 Iter 71: rel. change in X = 9.920835e-11
Final relative error: 8.178410e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: user-provided
Return type: ttensor

 Iter 10: rel. change in X = 4.952388e-06
 Iter 19: rel. change in X = 4.212445e-11
Final relative error: 8.930580e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: user-provided
Return type: ttensor

 Iter 10: rel. change in X = 2.748723e-02
 Iter 20: rel. change in X = 1.995712e-03
 Iter 30: rel. change in X = 1.198167e-04
 Iter 40: rel. change in X = 6.574993e-06
 Iter 50: rel. change in X = 5.009613e-07
 Iter 60: rel. change in X = 2.309549e-08
 Iter 70: rel. change in X = 1.895766e-09
 Iter 80: rel. change in X = 9.578569e-11
Final relative error: 7.268388e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: user-provided
Return type: ttensor

 Iter 10: rel. change in X = 9.924103e-03
 Iter 20: rel. change in X = 8.831319e-05
 Iter 30: rel. change in X = 7.842128e-07
 Iter 40: rel. change in X = 6.963572e-09
 Iter 49: rel. change in X = 9.917060e-11
Final relative error: 8.915486e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: user-provided
Return type: ttensor

 Iter 10: rel. change in X = 5.525827e-05
 Iter 20: rel. change in X = 4.083307e-09
 Iter 24: rel. change in X = 9.030115e-11
Final relative error: 8.863422e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: user-provided
Return type: ttensor

 Iter 10: rel. change in X = 3.944911e-03
 Iter 20: rel. change in X = 5.363932e-05
 Iter 30: rel. change in X = 6.962089e-07
 Iter 40: rel. change in X = 9.036423e-09
 Iter 50: rel. change in X = 1.172881e-10
 Iter 51: rel. change in X = 7.595905e-11
Final relative error: 8.566810e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: user-provided
Return type: ttensor

 Iter 10: rel. change in X = 1.292455e-02
 Iter 20: rel. change in X = 8.400218e-05
 Iter 30: rel. change in X = 3.583701e-07
 Iter 40: rel. change in X = 1.528702e-09
 Iter 45: rel. change in X = 9.984268e-11
Final relative error: 8.169955e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: user-provided
Return type: ttensor

 Iter 10: rel. change in X = 1.335792e-02
 Iter 20: rel. change in X = 8.125965e-09
 Iter 24: rel. change in X = 2.657894e-11
Final relative error: 8.221815e-01

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: user-provided
Return type: ttensor

 Iter 10: rel. change in X = 4.458519e-05
 Iter 20: rel. change in X = 1.272609e-09
 Iter 23: rel. change in X = 5.514294e-11
Final relative error: 8.429423e-01


*** Best relative error over 10 runs (user init): 0.726839 ***
ans is a ttensor of size 4 x 4 x 4
	ans.core is a tensor of size 2 x 2 x 2
		ans.core(:,:,1) = 
	    0.2591   -1.8347
	   -1.8347   -0.0258
		ans.core(:,:,2) = 
	   -1.8347   -0.0258
	   -0.0258    0.1413
	ans.U{1} = 
		   -0.6941   -0.0516
		    0.0120   -0.5852
		   -0.6632    0.3509
		   -0.2796   -0.7292
	ans.U{2} = 
		   -0.6941   -0.0516
		    0.0120   -0.5852
		   -0.6632    0.3509
		   -0.2796   -0.7292
	ans.U{3} = 
		   -0.6941   -0.0516
		    0.0120   -0.5852
		   -0.6632    0.3509
		   -0.2796   -0.7292

We do the same test again, but use the n-vecs initialization method

This initialization is more expensive but generally works very well. This initialization method is deterministic, so we don't need to set the random seed.

T = tucker_sym(S,2,'init','nvecs','printitn',10);
errs(3) = norm(S - full(T)) / norm(S);
Computing 2 leading e-vectors.

Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: nvecs
Return type: ttensor

 Iter 10: rel. change in X = 1.498092e-06
 Iter 20: rel. change in X = 3.673933e-11
Final relative error: 6.816876e-01

Comparison of errors with different initialization methods

The first two methods are repeated runs with random initializations, while the third method is a single run with the n-vecs initialization. Generally, the n-vecs initialization is better than a random initialization, but it is not guaranteed to find the best solution. Repeated runs with random initializations may find a better solution than a single run with n-vecs.

fprintf('\nSummary of best relative errors for each initialization method:\n');
fprintf('  Default random init (best of 10): %g\n', errs(1));
fprintf('  User random init (best of 10):    %g\n', errs(2));
fprintf('  N-vecs init (single run):         %g\n', errs(3));
Summary of best relative errors for each initialization method:
  Default random init (best of 10): 0.666382
  User random init (best of 10):    0.726839
  N-vecs init (single run):         0.681688

Example: Symmetric tensor with known rank-2 decomposition

If we know the rank, then tucker_sym should compute an exact decomposition with an error of zero (or very close to zero due to numerical errors).

% Create a symmetric factor matrix
A = rand(4,2);
% Construct a symmetric core tensor (2x2x2)
G = tensor(rand(2,2,2));
% Build the symmetric Tucker tensor
S_known = ttensor(G, {A, A, A});
% Symmetrize to ensure exact symmetry
S_known = symmetrize(full(S_known));
% Compute the symmetric Tucker decomposition
T_est = tucker_sym(S_known, 2);
disp('Relative error between original and estimated tensor:');
disp(norm(S_known - full(T_est)) / norm(S_known));
Symmetric Tucker:
Tensor size: 4 x 4 x 4
Rank (R): 2
Tolerance (matrix rel diff): 1.00e-10
Max iterations: 1000
Initialization: random
Return type: ttensor

 Iter  1: rel. change in X = 1.002407e+00
 Iter  2: rel. change in X = 2.283370e-14
Final relative error: 9.196364e-16
Relative error between original and estimated tensor:
   9.1964e-16