Shifted Power Method for Generalized Tensor Eigenproblem

>> Tensor Toolbox >> Tensor Eigenproblem >> GEAP

The method is described in the following paper:

* T. G. Kolda, J. R. Mayo, An Adaptive Shifted Power Method for
  Computing Generalized Tensor Eigenpairs, SIAM J. Matrix Analysis and
  Applications, 35:1563-1582, 2014, http://dx.doi.org/0.1137/140951758

Contents

rng('default'); %<- Setup for reproducibility

Data tensor for Example 5.1 in Kolda and Mayo (2014)

Tensor taken from Example 1 in E. Kofidis and P. A. Regalia, On the best rank-1 approximation of higher-order supersymmetric tensors, SIAM J. Matrix Anal. Appl., 23:863-884, 2002, http://dx.doi.org/10.1137/S0895479801387413.

% Unique values in lexiographical order
Avals = [0.2883 -0.0031 0.1973 -0.2485 -0.2939 0.3847 0.2972 0.1862 ...
    0.0919 -0.3619 0.1241 -0.3420 0.2127 0.2727 -0.3054 ];

% Create standard tensor object directly
% because that's what EIG_GEAP expects.
A = full(symtensor(Avals, 4,3));

% Display the tensor as a symmetric tensor
disp(symtensor(A),'A')

% Save order and dimension
m = ndims(Asym);
n = size(Asym,1);

A is a symmetric tensor with 4 modes of dimension 3
	(1,1,1,1)    0.2883
	(1,1,1,2)   -0.0031
	(1,1,1,3)    0.1973
	(1,1,2,2)   -0.2485
	(1,1,2,3)   -0.2939
	(1,1,3,3)    0.3847
	(1,2,2,2)    0.2972
	(1,2,2,3)    0.1862
	(1,2,3,3)    0.0919
	(1,3,3,3)   -0.3619
	(2,2,2,2)    0.1241
	(2,2,2,3)   -0.3420
	(2,2,3,3)    0.2127
	(2,3,3,3)    0.2727
	(3,3,3,3)   -0.3054

Create corresponding "identity" tensor

Create an identity tensor of order m and dimension n. Here we use a somewhat convoluted construction method, but it aligns with the mathematical derivation.

if mod(m,2) ~= 0
    error('Identity tensor on exists for even order');
end
Bsym = symtensor(@zeros,m,n);
uniqidx = indices(Bsym);
Bvals = zeros(size(uniqidx,1),1);
for i = 1:size(uniqidx,1)
    pidx = perms(uniqidx(i,:));
    pidxodd = pidx(:,1:2:m-1);
    pidxeven = pidx(:,2:2:m);
    pidxresult = pidxodd == pidxeven;
    Bvals(i) = sum(all(pidxresult,2))/factorial(m);
end
Bsym = symtensor(Bvals, m,n);
B = full(Bsym);
disp(Bsym,'B')

B is a symmetric tensor with 4 modes of dimension 3
	(1,1,1,1)    1.0000
	(1,1,1,2)         0
	(1,1,1,3)         0
	(1,1,2,2)    0.3333
	(1,1,2,3)         0
	(1,1,3,3)    0.3333
	(1,2,2,2)         0
	(1,2,2,3)         0
	(1,2,3,3)         0
	(1,3,3,3)         0
	(2,2,2,2)    1.0000
	(2,2,2,3)         0
	(2,2,3,3)    0.3333
	(2,3,3,3)         0
	(3,3,3,3)    1.0000

Demonstrate identity property

A tensor identity, E, satisfies the following mathematical property for any n-dimensional vector x.

$$\|x\|^2 x = E x^{\otimes 2}$$

Note that it is not scale invariant. When we test this property, it should be close to machine precision.

for i = 1:10
    x = rand(n,1);
    lhs = norm(x)^(m-2) * x;
    rhs = ttsv(B,x,-1);
    fprintf('Identity property error: %g\n', norm(lhs-rhs));
end

Identity property error: 2.22045e-16
Identity property error: 2.498e-16
Identity property error: 5.08768e-16
Identity property error: 4.996e-16
Identity property error: 0
Identity property error: 1.27192e-16
Identity property error: 3.14018e-16
Identity property error: 3.14018e-16
Identity property error: 0
Identity property error: 6.3596e-17

Call eig_geap to find eigenpair

Default is to find a maxima, i.e., a convex solution.

info = eig_geap(A, B, 'MaxIts', 100, 'Display',1);

Generalized Adaptive Tensor Eigenpair Power Method: Convex  
----  --------- ----- ------------ -----
Iter  Lambda    Diff  |newx-x|     Shift
----  --------- ----- ------------ -----
   1  -0.410701 +e+00 2.009472e-01  3.06
   2   0.115407 +e+00 2.554734e-01  2.62
   3   0.333322 +e-01 2.053365e-01  1.80
   4   0.362461 +e-02 9.743037e-02  1.03
   5   0.363244 +e-03 1.767961e-02  0.78
   6   0.363291 +e-04 3.655241e-03  0.79
   7   0.363302 +e-05 1.808619e-03  0.80
   8   0.363305 +e-06 9.127993e-04  0.81
   9   0.363306 +e-06 4.650736e-04  0.81
  10   0.363306 +e-07 2.380818e-04  0.81
  11   0.363306 +e-07 1.221727e-04  0.81
  12   0.363306 +e-08 6.277042e-05  0.81
  13   0.363306 +e-08 3.227074e-05  0.81
  14   0.363306 +e-09 1.659599e-05  0.81
  15   0.363306 +e-10 8.536291e-06  0.81
  16   0.363306 +e-10 4.391091e-06  0.81
  17   0.363306 +e-11 2.258888e-06  0.81
  18   0.363306 +e-11 1.162055e-06  0.81
  19   0.363306 +e-12 5.978110e-07  0.81
  20   0.363306 +e-12 3.075414e-07  0.81
  21   0.363306 +e-13 1.582139e-07  0.81
  22   0.363306 +e-14 8.139288e-08  0.81
  23   0.363306 +e-14 4.187246e-08  0.81
  24   0.363306 +e-15 2.154124e-08  0.81
  25   0.363306 +e-15 1.108187e-08  0.81
Successful Convergence

Demonstrate generalized eigenpair property

$$A x^{\otimes 2} = \lambda B x^{\otimes 2}$$

x = info.x;
lambda = info.lambda;
lhs = ttsv(A,x,-1);
rhs = lambda * ttsv(B,x,-1);
fprintf('Generalized eigenpair identity norm: %g\n', norm(lhs-rhs));

Generalized eigenpair identity norm: 6.70732e-09

Call eig_geap to find another eigenpair

Find a minima (by specifying 'Concave' = true)

info = eig_geap(A, B, 'MaxIts', 100, 'Concave', true, 'Display',1);

Generalized Adaptive Tensor Eigenpair Power Method: Concave 
----  --------- ----- ------------ -----
Iter  Lambda    Diff  |newx-x|     Shift
----  --------- ----- ------------ -----
   1   0.375422 -e-01 1.845569e-01 -2.27
   2   0.052994 -e+00 2.368751e-01 -1.99
   3  -0.041959 -e-01 1.620719e-01 -1.55
   4  -0.045031 -e-03 3.357753e-02 -1.26
   5  -0.045085 -e-04 4.020527e-03 -1.21
   6  -0.045091 -e-05 1.369335e-03 -1.20
   7  -0.045092 -e-06 4.682150e-04 -1.20
   8  -0.045092 -e-07 1.597911e-04 -1.20
   9  -0.045092 -e-08 5.449789e-05 -1.20
  10  -0.045092 -e-09 1.858280e-05 -1.20
  11  -0.045092 -e-10 6.335930e-06 -1.20
  12  -0.045092 -e-11 2.160222e-06 -1.20
  13  -0.045092 -e-12 7.365167e-07 -1.20
  14  -0.045092 -e-13 2.511109e-07 -1.20
  15  -0.045092 -e-14 8.561464e-08 -1.20
  16  -0.045092 -e-15 2.918975e-08 -1.20
  17  -0.045092 -e-16 9.952054e-09 -1.20
Successful Convergence

Test generalized eigenvector properties

x = info.x;
lambda = info.lambda;
lhs = ttsv(A,x,-1);
rhs = lambda * ttsv(B,x,-1);
fprintf('Generalized eigenpair identity norm: %g\n', norm(lhs-rhs));

Generalized eigenpair identity norm: 4.22637e-09

Reproduce tensors from Example 5.3 from Kolda and Mayo (2014)

Avals = [0.4982 -0.0582 -1.1719 0.2236 ...
    -0.0171 0.4597 0.4880 0.1852 ...
    -0.4087 0.7639 0.0000 -0.6162 ...
    0.1519 0.7631 2.6311];

Bvals = [3.0800 0.0614 0.2317 0.8140 ...
    0.0130 2.3551 0.0486 0.0616 ...
    0.0482 0.5288 1.9321 0.0236 ...
    1.8563 0.0681 16.0480];


A = full(symtensor(Avals, 4,3));
B = full(symtensor(Bvals, 4,3));

Compute an eigenpair and test it

info = eig_geap(A, B, 'MaxIts', 100, 'Display',1);

Generalized Adaptive Tensor Eigenpair Power Method: Convex  
----  --------- ----- ------------ -----
Iter  Lambda    Diff  |newx-x|     Shift
----  --------- ----- ------------ -----
   1   0.202754 +e-02 1.169885e-01  0.09
   2   0.203746 +e-03 2.568496e-02  0.14
   3   0.205219 +e-03 3.063915e-02  0.15
   4   0.207435 +e-03 3.717155e-02  0.15
   5   0.210711 +e-02 4.467892e-02  0.16
   6   0.215404 +e-02 5.282495e-02  0.17
   7   0.221730 +e-02 6.069324e-02  0.18
   8   0.229403 +e-02 6.641258e-02  0.20
   9   0.237235 +e-02 6.704464e-02  0.21
  10   0.243478 +e-02 5.985900e-02  0.23
  11   0.247275 +e-02 4.601160e-02  0.27
  12   0.249217 +e-03 3.184013e-02  0.31
  13   0.250186 +e-03 2.171951e-02  0.34
  14   0.250689 +e-03 1.523781e-02  0.37
  15   0.250963 +e-04 1.102013e-02  0.39
  16   0.251118 +e-04 8.165653e-03  0.41
  17   0.251208 +e-04 6.163198e-03  0.42
  18   0.251262 +e-04 4.717275e-03  0.43
  19   0.251294 +e-04 3.649284e-03  0.43
  20   0.251314 +e-05 2.846352e-03  0.44
  21   0.251327 +e-05 2.234267e-03  0.44
  22   0.251334 +e-05 1.762551e-03  0.45
  23   0.251339 +e-05 1.395871e-03  0.45
  24   0.251342 +e-06 1.108889e-03  0.45
  25   0.251344 +e-06 8.830627e-04  0.45
  26   0.251346 +e-06 7.045917e-04  0.45
  27   0.251346 +e-06 5.630598e-04  0.45
  28   0.251347 +e-06 4.505125e-04  0.45
  29   0.251347 +e-06 3.608169e-04  0.46
  30   0.251347 +e-07 2.892073e-04  0.46
  31   0.251348 +e-07 2.319560e-04  0.46
  32   0.251348 +e-07 1.861323e-04  0.46
  33   0.251348 +e-07 1.494218e-04  0.46
  34   0.251348 +e-07 1.199907e-04  0.46
  35   0.251348 +e-08 9.638168e-05  0.46
  36   0.251348 +e-08 7.743418e-05  0.46
  37   0.251348 +e-08 6.222201e-05  0.46
  38   0.251348 +e-08 5.000509e-05  0.46
  39   0.251348 +e-08 4.019126e-05  0.46
  40   0.251348 +e-09 3.230629e-05  0.46
  41   0.251348 +e-09 2.597006e-05  0.46
  42   0.251348 +e-09 2.087774e-05  0.46
  43   0.251348 +e-09 1.678470e-05  0.46
  44   0.251348 +e-09 1.349459e-05  0.46
  45   0.251348 +e-10 1.084972e-05  0.46
  46   0.251348 +e-10 8.723432e-06  0.46
  47   0.251348 +e-10 7.013982e-06  0.46
  48   0.251348 +e-10 5.639603e-06  0.46
  49   0.251348 +e-10 4.534587e-06  0.46
  50   0.251348 +e-10 3.646122e-06  0.46
  51   0.251348 +e-11 2.931758e-06  0.46
  52   0.251348 +e-11 2.357371e-06  0.46
  53   0.251348 +e-11 1.895526e-06  0.46
  54   0.251348 +e-11 1.524170e-06  0.46
  55   0.251348 +e-11 1.225571e-06  0.46
  56   0.251348 +e-12 9.854735e-07  0.46
  57   0.251348 +e-12 7.924141e-07  0.46
  58   0.251348 +e-12 6.371771e-07  0.46
  59   0.251348 +e-12 5.123524e-07  0.46
  60   0.251348 +e-12 4.119816e-07  0.46
  61   0.251348 +e-13 3.312740e-07  0.46
  62   0.251348 +e-13 2.663772e-07  0.46
  63   0.251348 +e-13 2.141939e-07  0.46
  64   0.251348 +e-13 1.722334e-07  0.46
  65   0.251348 +e-13 1.384930e-07  0.46
  66   0.251348 +e-13 1.113623e-07  0.46
  67   0.251348 +e-14 8.954655e-08  0.46
  68   0.251348 +e-14 7.200449e-08  0.46
  69   0.251348 +e-14 5.789891e-08  0.46
  70   0.251348 +e-14 4.655659e-08  0.46
  71   0.251348 +e-14 3.743623e-08  0.46
  72   0.251348 +e-15 3.010253e-08  0.46
  73   0.251348 +e-15 2.420549e-08  0.46
  74   0.251348 +e-15 1.946368e-08  0.46
Successful Convergence

Test generalized eigenvector properties

x = info.x;
lambda = info.lambda;
lhs = ttsv(A,x,-1);
rhs = lambda * ttsv(B,x,-1);
fprintf('Generalized eigenpair identity norm: %g\n', norm(lhs-rhs));

Generalized eigenpair identity norm: 4.69887e-08

Compute another eigenpair and test it

info = eig_geap(A, B, 'MaxIts', 100, 'Display',1);

Generalized Adaptive Tensor Eigenpair Power Method: Convex  
----  --------- ----- ------------ -----
Iter  Lambda    Diff  |newx-x|     Shift
----  --------- ----- ------------ -----
   1   0.311558 +e+00 3.548427e-01  0.38
   2   0.471060 +e-01 2.202178e-01  1.16
   3   0.488515 +e-02 5.394642e-02  1.07
   4   0.503547 +e-02 5.040663e-02  1.06
   5   0.515362 +e-02 4.543136e-02  1.04
   6   0.523706 +e-02 3.881159e-02  1.01
   7   0.529043 +e-02 3.152998e-02  0.99
   8   0.532173 +e-03 2.448560e-02  0.97
   9   0.533881 +e-03 1.829630e-02  0.95
  10   0.534759 +e-03 1.324809e-02  0.94
  11   0.535192 +e-03 9.359745e-03  0.92
  12   0.535397 +e-04 6.491499e-03  0.92
  13   0.535493 +e-04 4.442031e-03  0.91
  14   0.535537 +e-04 3.010740e-03  0.91
  15   0.535557 +e-05 2.027133e-03  0.90
  16   0.535565 +e-05 1.358672e-03  0.90
  17   0.535569 +e-05 9.078326e-04  0.90
  18   0.535571 +e-06 6.053311e-04  0.90
  19   0.535572 +e-06 4.030639e-04  0.90
  20   0.535572 +e-06 2.681327e-04  0.90
  21   0.535572 +e-07 1.782605e-04  0.90
  22   0.535572 +e-07 1.184625e-04  0.90
  23   0.535572 +e-08 7.870213e-05  0.90
  24   0.535572 +e-08 5.227724e-05  0.90
  25   0.535572 +e-08 3.472049e-05  0.90
  26   0.535572 +e-09 2.305812e-05  0.90
  27   0.535572 +e-09 1.531224e-05  0.90
  28   0.535572 +e-09 1.016806e-05  0.90
  29   0.535572 +e-10 6.751913e-06  0.90
  30   0.535572 +e-10 4.483415e-06  0.90
  31   0.535572 +e-10 2.977052e-06  0.90
  32   0.535572 +e-11 1.976792e-06  0.90
  33   0.535572 +e-11 1.312603e-06  0.90
  34   0.535572 +e-11 8.715749e-07  0.90
  35   0.535572 +e-12 5.787287e-07  0.90
  36   0.535572 +e-12 3.842773e-07  0.90
  37   0.535572 +e-13 2.551609e-07  0.90
  38   0.535572 +e-13 1.694272e-07  0.90
  39   0.535572 +e-13 1.124999e-07  0.90
  40   0.535572 +e-14 7.470004e-08  0.90
  41   0.535572 +e-14 4.960090e-08  0.90
  42   0.535572 +e-14 3.293504e-08  0.90
  43   0.535572 +e-15 2.186889e-08  0.90
  44   0.535572 +e-15 1.452096e-08  0.90
Successful Convergence

Test generalized eigenvector properties

x = info.x;
lambda = info.lambda;
lhs = ttsv(A,x,-1);
rhs = lambda * ttsv(B,x,-1);
fprintf('Generalized eigenpair identity norm: %g\n', norm(lhs-rhs));

Generalized eigenpair identity norm: 5.15774e-08

Tensor Toolbox for MATLAB: www.tensortoolbox.org.