Implicit Symmetric CP Decomposition for Symmetric K-Tensors

The function cp_isym computes the symmetric CP decomposition of a symmetric tensor that is in symmetric k-tensor format, meaning that it is the sum of symmetric outer products and stored as a symktensor object. The decomposition is described in the following reference:

Contents

Requirements

This code requires an optimization solver. Our examples use the L-BFGS-B optimization method; see Optimization Options for details of installation and usage as well as other options for optimization solvers.

Create a sample problem based on a Gaussian mixture model.

Here we used the third-order moment (d=3), for a problem with random variables of dimension n=100. We take a mixture of r=5 Gaussians and collect p=750 observations.

d = 3;
n = 100;
r = 5;
p = 250*r;

We construct the means so that they are collinear, i.e., the cosine of the angle between any two vectors is 0.5. This is a more difficult problem than purely random vectors which are close to orthogonal.

rng('default');
collinearity = 0.5;
Atrue = matrandcong(n,r,collinearity);

Each Gaussian gets equal weight in the mixture.

wghts = 1/r * ones(r,1);

The `true` solution is based on the means and their prevalence in the mixture.

Mtrue = symktensor(wghts,Atrue,d);
Mtrue = normalize(Mtrue);

The data tensor is based on noisy observations from the mixture model. The idx determines which Gaussian is sampled. And the samples have noise as specified by stddev. The result is assembled into an observation tensor that in symktensor format.

idx = randsample(r,p,true,wghts);
stddev = 1e-3;
noise = stddev * randn(n,p);
V = Atrue(:,idx) + noise;
X = symktensor(1/p*ones(p,1),V,d);
X = normalize(X,0);

Recommend multiple runs of CP-ISYM

Since this is a non-convex optimization problem, it's critically important to do multiple runs of the optimzation problem. In the code below, we run 10 times. Not every run will converge to the global minimum. We take solution with the lowest function value, which should ideally yield the best match to the true solution. For our artificial problem, we can actually verify the score of how well the computed solution matches with the ideal solution.

rng(7)
nruns = 10;
M = cell(nruns,1);
info = cell(nruns,1);
for i = 1:nruns
    fprintf('------------- Run %d of %d --------------\n', i, nruns);
    [M{i},info{i}] = cp_isym(X,r,'Xnormsqr','exact','printitn',0);
    fprintf('Final F: %.2e, Score: %.2f\n', info{i}.f, score(M{i},Mtrue));
    if i == 1
        ifbest = i;
        fbest = info{i}.f;
        sfbest = score(M{i},Mtrue);
    elseif info{i}.f < fbest
        ifbest = i;
        fbest = info{i}.f;
        sfbest = score(M{i},Mtrue);
    end
end
------------- Run 1 of 10 --------------
Final F: 1.94e-02, Score: 0.38
------------- Run 2 of 10 --------------
Final F: 1.94e-02, Score: 0.41
------------- Run 3 of 10 --------------
Final F: 1.20e-09, Score: 0.95
------------- Run 4 of 10 --------------
Final F: 7.76e-09, Score: 0.95
------------- Run 5 of 10 --------------
Final F: 1.94e-02, Score: 0.41
------------- Run 6 of 10 --------------
Final F: 5.01e-09, Score: 0.95
------------- Run 7 of 10 --------------
Final F: 3.65e-09, Score: 0.95
------------- Run 8 of 10 --------------
Final F: 7.29e-09, Score: 0.95
------------- Run 9 of 10 --------------
Final F: 1.94e-02, Score: 0.41
------------- Run 10 of 10 --------------
Final F: 6.14e-02, Score: 0.00
fprintf('------------- Best of %d runs --------------\n', nruns);
fprintf('Run %d -> Final F: %.2e, Score: %.2f\n', ifbest, fbest, sfbest);
------------- Best of 10 runs --------------
Run 3 -> Final F: 1.20e-09, Score: 0.95

Options for CP-ISYM with L-BFGS-B (default) optimization solver

The CP solution should produce something close to Mtrue, but it will not be exact because of the randomness in building the observation tensor. Here we show what happens when the true solution is provided as an initial guess. The function should converge to zero, and the similarlity score should be close to one. The method defaults to 'lfbgsb'.

[M,info] = cp_isym(X,r,'Xnormsqr','exact','init',Mtrue);
fprintf('\n\t\t*** Final F: %.2e, Score: %.2f ***\n\n', info.f, score(M,Mtrue));
CP-ISYM Implict Symmetric CP Optimzation
L-BFGS-B Optimization (https://github.com/stephenbeckr/L-BFGS-B-C)
Order, Size: 3, 100
Lower bound: -Inf, Upper bound: Inf
Parameters: m=5, ftol=1e-10, gtol=1e-05, maxiters = 1000, subiters = 10

Begin Main Loop
Iter     1, f(x) = 1.798804e-04, ||grad||_infty = 1.54e-02, 4.50e-03
Iter     2, f(x) = 2.009885e-06, ||grad||_infty = 1.15e-03, 5.64e-03
Iter     3, f(x) = 1.320009e-06, ||grad||_infty = 7.90e-04, 6.92e-03
Iter     4, f(x) = 9.160013e-07, ||grad||_infty = 2.02e-04, 8.24e-03
Iter     5, f(x) = 6.929470e-07, ||grad||_infty = 1.16e-04, 9.59e-03
Iter     6, f(x) = 2.975243e-07, ||grad||_infty = 2.70e-04, 1.09e-02
Iter     7, f(x) = 1.647469e-07, ||grad||_infty = 1.26e-04, 1.23e-02
Iter     8, f(x) = 1.401774e-07, ||grad||_infty = 6.75e-05, 1.37e-02
Iter     9, f(x) = 1.275847e-07, ||grad||_infty = 7.95e-05, 1.50e-02
Iter    10, f(x) = 7.883570e-08, ||grad||_infty = 1.39e-04, 1.70e-02
Iter    11, f(x) = 5.170653e-08, ||grad||_infty = 3.19e-05, 1.83e-02
Iter    12, f(x) = 3.582501e-08, ||grad||_infty = 1.63e-05, 1.95e-02
Iter    13, f(x) = 1.921946e-08, ||grad||_infty = 9.09e-05, 2.10e-02
Iter    14, f(x) = 6.336888e-09, ||grad||_infty = 4.28e-05, 2.22e-02
Iter    15, f(x) = 3.691645e-09, ||grad||_infty = 3.04e-05, 2.37e-02
Iter    16, f(x) = 2.252858e-09, ||grad||_infty = 9.17e-06, 2.51e-02
End Main Loop

Final f: 2.2529e-09
Setup time: 0.0024 seconds
Optimization time: 0.026 seconds
Iterations: 16
Total iterations: 35
Exit condition: CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL.

		*** Final F: 2.25e-09, Score: 0.95 ***

It is also possible to specify the initial guess as factor matrix

[M,info] = cp_isym(X,r,'Xnormsqr','exact','init',Mtrue);
fprintf('\n\t\t*** Final F: %.2e, Score: %.2f ***\n\n', info.f, score(M,Mtrue));
CP-ISYM Implict Symmetric CP Optimzation
L-BFGS-B Optimization (https://github.com/stephenbeckr/L-BFGS-B-C)
Order, Size: 3, 100
Lower bound: -Inf, Upper bound: Inf
Parameters: m=5, ftol=1e-10, gtol=1e-05, maxiters = 1000, subiters = 10

Begin Main Loop
Iter     1, f(x) = 1.798804e-04, ||grad||_infty = 1.54e-02, 4.42e-03
Iter     2, f(x) = 2.009885e-06, ||grad||_infty = 1.15e-03, 5.64e-03
Iter     3, f(x) = 1.320009e-06, ||grad||_infty = 7.90e-04, 7.02e-03
Iter     4, f(x) = 9.160013e-07, ||grad||_infty = 2.02e-04, 8.23e-03
Iter     5, f(x) = 6.929470e-07, ||grad||_infty = 1.16e-04, 1.00e-02
Iter     6, f(x) = 2.975243e-07, ||grad||_infty = 2.70e-04, 1.14e-02
Iter     7, f(x) = 1.647469e-07, ||grad||_infty = 1.26e-04, 1.26e-02
Iter     8, f(x) = 1.401774e-07, ||grad||_infty = 6.75e-05, 1.39e-02
Iter     9, f(x) = 1.275847e-07, ||grad||_infty = 7.95e-05, 1.53e-02
Iter    10, f(x) = 7.883570e-08, ||grad||_infty = 1.39e-04, 1.67e-02
Iter    11, f(x) = 5.170653e-08, ||grad||_infty = 3.19e-05, 1.80e-02
Iter    12, f(x) = 3.582501e-08, ||grad||_infty = 1.63e-05, 1.93e-02
Iter    13, f(x) = 1.921946e-08, ||grad||_infty = 9.09e-05, 2.06e-02
Iter    14, f(x) = 6.336888e-09, ||grad||_infty = 4.28e-05, 2.19e-02
Iter    15, f(x) = 3.691645e-09, ||grad||_infty = 3.04e-05, 2.31e-02
Iter    16, f(x) = 2.252858e-09, ||grad||_infty = 9.17e-06, 2.46e-02
End Main Loop

Final f: 2.2529e-09
Setup time: 0.0016 seconds
Optimization time: 0.025 seconds
Iterations: 16
Total iterations: 35
Exit condition: CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL.

		*** Final F: 2.25e-09, Score: 0.95 ***

So far we've been using the exact objective function which is close to zero for the optimal zero. We can, however, ignore the constant |X|^2 term in objective function, which saves some preprocessing cost. The main disadvantage of ignoring the constant term is that the final function value will no longer be near zero.

[M,info] = cp_isym(X,r,'Xnormsqr',0,'init',Atrue,'printitn',1);
fprintf('\n\t\t*** Final F (adjusted): %.2e, Score: %.2f ***\n\n', ...
    info.f + norm(X)^2, score(M,Mtrue));
CP-ISYM Implict Symmetric CP Optimzation
L-BFGS-B Optimization (https://github.com/stephenbeckr/L-BFGS-B-C)
Order, Size: 3, 100
Lower bound: -Inf, Upper bound: Inf
Parameters: m=5, ftol=1e-10, gtol=1e-05, maxiters = 1000, subiters = 10

Begin Main Loop
Iter     1, f(x) = -2.108869e-01, ||grad||_infty = 5.67e-02, 3.00e-03
Iter     2, f(x) = -2.551309e-01, ||grad||_infty = 4.68e-02, 5.40e-03
Iter     3, f(x) = -2.954389e-01, ||grad||_infty = 1.44e-02, 7.86e-03
Iter     4, f(x) = -2.983326e-01, ||grad||_infty = 1.31e-02, 8.97e-03
Iter     5, f(x) = -2.998836e-01, ||grad||_infty = 5.63e-03, 1.03e-02
Iter     6, f(x) = -3.003519e-01, ||grad||_infty = 5.88e-04, 1.20e-02
Iter     7, f(x) = -3.003573e-01, ||grad||_infty = 7.43e-05, 1.36e-02
Iter     8, f(x) = -3.003574e-01, ||grad||_infty = 5.83e-05, 1.49e-02
Iter     9, f(x) = -3.003576e-01, ||grad||_infty = 6.35e-05, 1.62e-02
Iter    10, f(x) = -3.003578e-01, ||grad||_infty = 7.69e-05, 1.76e-02
Iter    11, f(x) = -3.003579e-01, ||grad||_infty = 6.13e-05, 2.02e-02
Iter    12, f(x) = -3.003580e-01, ||grad||_infty = 1.96e-05, 2.14e-02
Iter    13, f(x) = -3.003580e-01, ||grad||_infty = 2.51e-05, 2.26e-02
Iter    14, f(x) = -3.003580e-01, ||grad||_infty = 3.02e-06, 2.38e-02
End Main Loop

Final f: -3.0036e-01
Setup time: 0.0013 seconds
Optimization time: 0.024 seconds
Iterations: 14
Total iterations: 33
Exit condition: CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL.

		*** Final F (adjusted): 4.60e-10, Score: 0.95 ***

By default, the method uses the randomized range finder (rrf) for the initial guess. This is the default initial guess and works well.

rng(1)
[M,info] = cp_isym(X,r,'Xnormsqr','exact','printitn',5,'init','rrf');
fprintf('\n\t\t*** Final F: %.2e, Score: %.2f ***\n\n', info.f, score(M,Mtrue));
CP-ISYM Implict Symmetric CP Optimzation
L-BFGS-B Optimization (https://github.com/stephenbeckr/L-BFGS-B-C)
Order, Size: 3, 100
Lower bound: -Inf, Upper bound: Inf
Parameters: m=5, ftol=1e-10, gtol=1e-05, maxiters = 1000, subiters = 10

Begin Main Loop
Iter     5, f(x) = 4.265680e-02, ||grad||_infty = 1.81e-02, 6.42e-03
Iter    10, f(x) = 2.260594e-02, ||grad||_infty = 1.28e-02, 1.04e-02
Iter    15, f(x) = 1.363733e-02, ||grad||_infty = 1.88e-02, 1.48e-02
Iter    20, f(x) = 1.381487e-03, ||grad||_infty = 1.22e-02, 1.95e-02
Iter    25, f(x) = 5.283198e-05, ||grad||_infty = 1.45e-03, 2.32e-02
Iter    30, f(x) = 6.142064e-07, ||grad||_infty = 1.11e-04, 2.67e-02
Iter    35, f(x) = 1.823079e-08, ||grad||_infty = 1.78e-05, 3.09e-02
Iter    39, f(x) = 8.114700e-10, ||grad||_infty = 3.67e-06, 3.46e-02
End Main Loop

Final f: 8.1147e-10
Setup time: 0.00097 seconds
Optimization time: 0.035 seconds
Iterations: 39
Total iterations: 86
Exit condition: CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL.

		*** Final F: 8.11e-10, Score: 0.95 ***

If we run instead with a random Gaussian initial guess, we are unlikey to converge.

rng(77)
[M,info,M0] = cp_isym(X,r,'Xnormsqr','exact','printitn',5,'init',@randn);
fprintf('\n\t\t*** Final F: %.2e, Score: %.2f ***\n\n', info.f, score(M,Mtrue));
CP-ISYM Implict Symmetric CP Optimzation
L-BFGS-B Optimization (https://github.com/stephenbeckr/L-BFGS-B-C)
Order, Size: 3, 100
Lower bound: -Inf, Upper bound: Inf
Parameters: m=5, ftol=1e-10, gtol=1e-05, maxiters = 1000, subiters = 10

Begin Main Loop
Iter     5, f(x) = 1.568859e+00, ||grad||_infty = 2.43e+03, 4.85e-03
Iter    10, f(x) = 3.002914e-01, ||grad||_infty = 1.21e-02, 9.12e-03
Iter    11, f(x) = 3.002914e-01, ||grad||_infty = 1.68e-03, 9.79e-03
End Main Loop

Final f: 3.0029e-01
Setup time: 0.0032 seconds
Optimization time: 0.01 seconds
Iterations: 11
Total iterations: 24
Exit condition: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH.

		*** Final F: 3.00e-01, Score: 0.00 ***

Sometimes a random initialization does a little better if the convergence tolerances are tightened and the maximum iterations is increased. We also increase the memory in Limited-memory BFGS to see if that helps. It gets a little better, but still no where near the solution with the randomized range finder.

[M,info] = cp_isym(X,r,'Xnormsqr','exact','printitn',5,'init',M0,...
    'ftol',1e-14,'gtol',1e-8,'maxiters',10000,'m',25);
fprintf('\n\t\t*** Final F: %.2e, Score: %.2f ***\n\n', info.f, score(M,Mtrue));
CP-ISYM Implict Symmetric CP Optimzation
L-BFGS-B Optimization (https://github.com/stephenbeckr/L-BFGS-B-C)
Order, Size: 3, 100
Lower bound: -Inf, Upper bound: Inf
Parameters: m=25, ftol=1e-14, gtol=1e-08, maxiters = 10000, subiters = 10

Begin Main Loop
Iter     5, f(x) = 1.568859e+00, ||grad||_infty = 2.43e+03, 4.86e-03
Iter    10, f(x) = 3.002914e-01, ||grad||_infty = 8.83e-03, 8.24e-03
Iter    15, f(x) = 3.002914e-01, ||grad||_infty = 9.08e-04, 1.36e-02
Iter    20, f(x) = 3.002914e-01, ||grad||_infty = 1.39e-02, 1.75e-02
Iter    25, f(x) = 3.002914e-01, ||grad||_infty = 1.59e-01, 2.18e-02
Iter    30, f(x) = 3.002899e-01, ||grad||_infty = 1.81e+00, 2.69e-02
Iter    35, f(x) = 2.191333e-01, ||grad||_infty = 3.14e+02, 3.63e-02
Iter    40, f(x) = 1.203681e-01, ||grad||_infty = 4.07e+03, 4.29e-02
Iter    45, f(x) = 1.116589e-01, ||grad||_infty = 7.46e+02, 5.05e-02
Iter    50, f(x) = 1.041222e-01, ||grad||_infty = 4.91e+01, 6.26e-02
Iter    55, f(x) = 1.034506e-01, ||grad||_infty = 1.35e-03, 6.81e-02
Iter    60, f(x) = 1.034506e-01, ||grad||_infty = 4.54e-02, 7.33e-02
Iter    65, f(x) = 1.034504e-01, ||grad||_infty = 5.71e-01, 7.93e-02
Iter    70, f(x) = 1.034220e-01, ||grad||_infty = 6.39e+00, 8.47e-02
Iter    75, f(x) = 1.021516e-01, ||grad||_infty = 1.85e+02, 9.17e-02
Iter    80, f(x) = 8.415741e-02, ||grad||_infty = 5.40e+02, 9.89e-02
Iter    85, f(x) = 8.404990e-02, ||grad||_infty = 6.41e+00, 1.05e-01
Iter    90, f(x) = 8.404952e-02, ||grad||_infty = 1.19e-01, 1.12e-01
Iter    95, f(x) = 8.404951e-02, ||grad||_infty = 9.01e-03, 1.18e-01
Iter    97, f(x) = 8.404951e-02, ||grad||_infty = 5.18e-05, 1.21e-01
End Main Loop

Final f: 8.4050e-02
Setup time: 0.00064 seconds
Optimization time: 0.12 seconds
Iterations: 97
Total iterations: 214
Exit condition: CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH.

		*** Final F: 8.40e-02, Score: 0.07 ***

Options for CP-ISYM with L-BFGS optimization solver from Poblano

rng(1)
[M,info] = cp_isym(X,r,'Xnormsqr','exact','printitn',5,'method','lbfgs');
fprintf('\n\t\t*** Final F: %.2e, Score: %.2f ***\n\n', info.f, score(M,Mtrue));
CP-ISYM Implict Symmetric CP Optimzation
Poblano L-BFGS Optimization
Order, Size: 3, 100
Parameters: m=5, ftol=1e-10, gtol=1e-05, maxiters = 1000, maxsubiters=10000

Begin Main Loop
 Iter  FuncEvals       F(X)          ||G(X)||/N        
------ --------- ---------------- ----------------
     0         1       3.92027347       0.02425307
     5        18       0.04889771       0.00026180
    10        30       0.02227461       0.00010159
    15        45       0.01590084       0.00016813
    20        60       0.00156181       0.00008811
    25        70       0.00005713       0.00001458
    26        73       0.00001608       0.00000982
End Main Loop

Final f: 1.6075e-05
Setup time: 0.0018 seconds
Optimization time: 0.087 seconds
Iterations: 26
Total iterations: 73
Exit condition: Successful termination based on StopTol

		*** Final F: 1.61e-05, Score: 0.95 ***

Options for CP-ISYM with Quasi-Newton optimization solver from Optimization Toolbox

rng(1)
[M,info] = cp_isym(X,r,'Xnormsqr','exact','printitn',5,'method','fminunc');
fprintf('\n\t\t*** Final F: %.2e, Score: %.2f ***\n\n', info.f, score(M,Mtrue));
CP-ISYM Implict Symmetric CP Optimzation
Quasi-Newton Optimization (via Optimization Toolbox)
Order, Size: 3, 100
Parameters: gtol=1e-05, maxiters = 1000, maxsubiters=10000

Begin Main Loop
                                                        First-order 
 Iteration  Func-count       f(x)        Step-size       optimality
     0           1          3.92027                          1.95
     1           3         0.275394       0.164109         0.0628  
     2           7          0.14028        12.0063           0.11  
     3           8         0.128657              1         0.0684  
     4          10         0.085529            0.5         0.0396  
     5          11        0.0815813              1         0.0462  
     6          12        0.0798902              1        0.00739  
     7          14        0.0776973             10         0.0267  
     8          15        0.0750789              1         0.0245  
     9          17        0.0704473       0.480552          0.042  
    10          18        0.0670865              1         0.0368  
    11          19        0.0616269              1         0.0325  
    12          20         0.055824              1         0.0328  
    13          21        0.0544128              1         0.0339  
    14          22        0.0505937              1         0.0154  
    15          23        0.0481738              1         0.0168  
    16          24        0.0416669              1          0.054  
    17          25        0.0381286              1         0.0336  
    18          26        0.0335275              1         0.0159  
    19          27        0.0317987              1         0.0146  
                                                        First-order 
 Iteration  Func-count       f(x)        Step-size       optimality
    20          28        0.0301083              1         0.0111  
    21          29        0.0293649              1         0.0161  
    22          30         0.027827              1        0.00893  
    23          31        0.0265521              1        0.00975  
    24          32        0.0261249              1        0.00918  
    25          33        0.0251787              1        0.00625  
    26          34          0.02455              1        0.00806  
    27          35         0.023779              1        0.00911  
    28          36        0.0230934              1        0.00845  
    29          37        0.0224019              1        0.00844  
    30          38        0.0220964              1        0.00687  
    31          39        0.0218681              1         0.0049  
    32          40         0.021706              1        0.00375  
    33          41        0.0215136              1         0.0033  
    34          42        0.0213676              1        0.00304  
    35          43        0.0212558              1        0.00323  
    36          44        0.0211222              1        0.00409  
    37          45        0.0208738              1        0.00475  
    38          46        0.0205572              1        0.00656  
    39          47        0.0201352              1         0.0065  
                                                        First-order 
 Iteration  Func-count       f(x)        Step-size       optimality
    40          48        0.0196422              1        0.00387  
    41          49        0.0193932              1        0.00333  
    42          50        0.0191319              1        0.00872  
    43          51        0.0186628              1        0.00878  
    44          52         0.017685              1        0.00907  
    45          53        0.0169716              1        0.00861  
    46          54          0.01577              1        0.00994  
    47          56        0.0131656            0.5         0.0106  
    48          57        0.0117538              1         0.0153  
    49          58       0.00879083              1         0.0134  
    50          59       0.00839261              1          0.012  
    51          60       0.00394004              1        0.00878  
    52          61       0.00304266              1          0.013  
    53          62       0.00244168              1         0.0069  
    54          63       0.00225538              1        0.00627  
    55          64        0.0014968              1        0.00464  
    56          65       0.00105029              1        0.00366  
    57          66      0.000674101              1        0.00456  
    58          67      0.000456871              1         0.0042  
    59          68      0.000245825              1        0.00299  
                                                        First-order 
 Iteration  Func-count       f(x)        Step-size       optimality
    60          69      9.62095e-05              1         0.0024  
    61          70      4.02262e-05              1        0.00114  
    62          71      2.70099e-05              1       0.000664  
    63          72      2.21994e-05              1       0.000556  
    64          73      1.73964e-05              1       0.000481  
    65          74      1.28057e-05              1       0.000433  
    66          75      1.05719e-05              1       0.000388  
    67          76       9.4027e-06              1        0.00034  
    68          77       7.7716e-06              1       0.000308  
    69          78      5.32426e-06              1       0.000286  
    70          79      2.70776e-06              1       0.000259  
    71          80      1.17346e-06              1       0.000219  
    72          81      6.75895e-07              1       0.000149  
    73          82      5.27851e-07              1       0.000111  
    74          83       4.0086e-07              1       9.15e-05  
    75          84      2.47388e-07              1        6.9e-05  
    76          85      1.30409e-07              1       5.92e-05  
    77          86      7.85393e-08              1       4.16e-05  
    78          87      6.17683e-08              1       3.58e-05  
    79          88      5.12695e-08              1       2.56e-05  

Local minimum found.

Optimization completed because the size of the gradient is less than
the selected value of the optimality tolerance.



End Main Loop

Final f: 5.1269e-08
Setup time: 1.1 seconds
Optimization time: 0.63 seconds
Iterations: 79
Total iterations: 88
Exit condition: Gradient norm smaller than tolerance

		*** Final F: 5.13e-08, Score: 0.95 ***

Options for CP-ISYM with Adam stochastic optimization solver

When the number of observations is very large, a stochastic method may be faster. Our example is very small (p=1250), but we use it nonetheless just to demonstrate the capabilities of the stochastic method. The method is called by specifying the 'method' to be 'adam'.

rng(5)
[M,info] = cp_isym(X,r,'method','adam');
fprintf('\n\t\t*** Final F: %.2e, Score: %.2f ***\n\n', f_implicit(M,X,norm(X)^2), score(M,Mtrue));
CP-ISYM Implict Symmetric CP Optimzation
Adam Stochastic Optimization
Order, Size: 3, 100
Function: 1000 samples out of 1250 observations
Gradient: using 100 samples  out of 1250 observations
Max iterations (epochs): 100
Iterations per epoch: 100
Learning rate / decay / maxfails: 0.01 0.1 1
Backup to best known solution after epoch fail? true
beta1 / beta2 / epsilon: 0.9 0.999 1e-08

Begin Main Loop
Epoch  0: f~ = 4.828625e+00
Epoch  1: f~ = -2.657838e-01, step = 0.01
Epoch  2: f~ = -2.832830e-01, step = 0.01
Epoch  3: f~ = -2.982487e-01, step = 0.01
Epoch  4: f~ = -2.985366e-01, step = 0.01
Epoch  5: f~ = -2.984525e-01, step = 0.01, nfails = 1 (resetting to solution from last epoch)
Epoch  6: f~ = -3.001280e-01, step = 0.001
Epoch  7: f~ = -2.996104e-01, step = 0.001, nfails = 2 (resetting to solution from last epoch)
End Main Loop

Final f~: -3.0013e-01
Setup time: 0.01 seconds
Optimization time: 0.40 seconds
Number of epochs: 7
Total iterations: 700
Exit condition: Number of failed epochs > 1

		*** Final F: 1.51e-04, Score: 0.95 ***

Since this problem is small, we can compute the function value exactly.

rng(5)
[M,info] = cp_isym(X,r,'method','adam','fsamp','exact','Xnormsqr','exact');
fprintf('\n\t\t*** Final F: %.2e, Score: %.2f ***\n\n', f_implicit(M,X,norm(X)^2), score(M,Mtrue));
CP-ISYM Implict Symmetric CP Optimzation
Adam Stochastic Optimization
Order, Size: 3, 100
Function: exact
Gradient: using 100 samples  out of 1250 observations
Max iterations (epochs): 100
Iterations per epoch: 100
Learning rate / decay / maxfails: 0.01 0.1 1
Backup to best known solution after epoch fail? true
beta1 / beta2 / epsilon: 0.9 0.999 1e-08

Begin Main Loop
Epoch  0: f = 5.123527e+00
Epoch  1: f = 3.595388e-02, step = 0.01
Epoch  2: f = 1.618031e-02, step = 0.01
Epoch  3: f = 1.507292e-03, step = 0.01
Epoch  4: f = 2.553210e-03, step = 0.01, nfails = 1 (resetting to solution from last epoch)
Epoch  5: f = 1.624327e-04, step = 0.001
Epoch  6: f = 1.601540e-04, step = 0.001
Epoch  7: f = 3.981757e-04, step = 0.001, nfails = 2 (resetting to solution from last epoch)
End Main Loop

Final f: 1.6015e-04
Setup time: 0.01 seconds
Optimization time: 0.41 seconds
Number of epochs: 7
Total iterations: 700
Exit condition: Number of failed epochs > 1

		*** Final F: 1.60e-04, Score: 0.96 ***

We can also change the number of samples in each stochastic gradient ('gsamp') and correspondingly increase the number of iterations per epoch ('subiters').

rng(5)
[M,info] = cp_isym(X,r,'method','adam','fsamp','exact','Xnormsqr','exact',...
    'gsamp',10,'subiters',250);
fprintf('\n\t\t*** Final F: %.2e, Score: %.2f ***\n\n', f_implicit(M,X,norm(X)^2), score(M,Mtrue));
CP-ISYM Implict Symmetric CP Optimzation
Adam Stochastic Optimization
Order, Size: 3, 100
Function: exact
Gradient: using 10 samples  out of 1250 observations
Max iterations (epochs): 100
Iterations per epoch: 250
Learning rate / decay / maxfails: 0.01 0.1 1
Backup to best known solution after epoch fail? true
beta1 / beta2 / epsilon: 0.9 0.999 1e-08

Begin Main Loop
Epoch  0: f = 5.123527e+00
Epoch  1: f = 2.239609e-02, step = 0.01
Epoch  2: f = 6.190717e-03, step = 0.01
Epoch  3: f = 1.955220e-02, step = 0.01, nfails = 1 (resetting to solution from last epoch)
Epoch  4: f = 6.693318e-03, step = 0.001, nfails = 2 (resetting to solution from last epoch)
End Main Loop

Final f: 6.1907e-03
Setup time: 0.00 seconds
Optimization time: 0.26 seconds
Number of epochs: 4
Total iterations: 1000
Exit condition: Number of failed epochs > 1

		*** Final F: 6.19e-03, Score: 0.89 ***

We can further improve by increasing the number of times we decrease the step length ('maxfails').

rng(5)
[M,info] = cp_isym(X,r,'method','adam','fsamp','exact','Xnormsqr','exact',...
    'gsamp',10,'subiters',250,'maxfails',4);
fprintf('\n\t\t*** Final F: %.2e, Score: %.2f ***\n\n', f_implicit(M,X,norm(X)^2), score(M,Mtrue));
CP-ISYM Implict Symmetric CP Optimzation
Adam Stochastic Optimization
Order, Size: 3, 100
Function: exact
Gradient: using 10 samples  out of 1250 observations
Max iterations (epochs): 100
Iterations per epoch: 250
Learning rate / decay / maxfails: 0.01 0.1 4
Backup to best known solution after epoch fail? true
beta1 / beta2 / epsilon: 0.9 0.999 1e-08

Begin Main Loop
Epoch  0: f = 5.123527e+00
Epoch  1: f = 2.239609e-02, step = 0.01
Epoch  2: f = 6.190717e-03, step = 0.01
Epoch  3: f = 1.955220e-02, step = 0.01, nfails = 1 (resetting to solution from last epoch)
Epoch  4: f = 6.693318e-03, step = 0.001, nfails = 2 (resetting to solution from last epoch)
Epoch  5: f = 1.111507e-03, step = 0.0001
Epoch  6: f = 3.658467e-04, step = 0.0001
Epoch  7: f = 1.958052e-04, step = 0.0001
Epoch  8: f = 2.450883e-04, step = 0.0001, nfails = 3 (resetting to solution from last epoch)
Epoch  9: f = 1.802417e-04, step = 1e-05
Epoch 10: f = 1.692018e-04, step = 1e-05
Epoch 11: f = 1.604512e-04, step = 1e-05
Epoch 12: f = 1.384002e-04, step = 1e-05
Epoch 13: f = 1.418440e-04, step = 1e-05, nfails = 4 (resetting to solution from last epoch)
Epoch 14: f = 1.356883e-04, step = 1e-06
Epoch 15: f = 1.352304e-04, step = 1e-06
Epoch 16: f = 1.333213e-04, step = 1e-06
Epoch 17: f = 1.320493e-04, step = 1e-06
Epoch 18: f = 1.304712e-04, step = 1e-06
Epoch 19: f = 1.304975e-04, step = 1e-06, nfails = 5 (resetting to solution from last epoch)
End Main Loop

Final f: 1.3047e-04
Setup time: 0.01 seconds
Optimization time: 1.41 seconds
Number of epochs: 19
Total iterations: 4750
Exit condition: Number of failed epochs > 4

		*** Final F: 1.30e-04, Score: 0.97 ***

We also change the initial learning rate via 'rate'.

rng(5)
[M,info] = cp_isym(X,r,'method','adam','fsamp','exact','Xnormsqr','exact','rate',1e-3);
fprintf('\n\t\t*** Final F: %.2e, Score: %.2f ***\n\n', f_implicit(M,X,norm(X)^2), score(M,Mtrue));
CP-ISYM Implict Symmetric CP Optimzation
Adam Stochastic Optimization
Order, Size: 3, 100
Function: exact
Gradient: using 100 samples  out of 1250 observations
Max iterations (epochs): 100
Iterations per epoch: 100
Learning rate / decay / maxfails: 0.001 0.1 1
Backup to best known solution after epoch fail? true
beta1 / beta2 / epsilon: 0.9 0.999 1e-08

Begin Main Loop
Epoch  0: f = 5.123527e+00
Epoch  1: f = 1.636029e-01, step = 0.001
Epoch  2: f = 6.960158e-02, step = 0.001
Epoch  3: f = 4.769470e-02, step = 0.001
Epoch  4: f = 3.871023e-02, step = 0.001
Epoch  5: f = 3.369702e-02, step = 0.001
Epoch  6: f = 3.012290e-02, step = 0.001
Epoch  7: f = 2.722983e-02, step = 0.001
Epoch  8: f = 2.466563e-02, step = 0.001
Epoch  9: f = 2.200866e-02, step = 0.001
Epoch 10: f = 1.999008e-02, step = 0.001
Epoch 11: f = 1.800170e-02, step = 0.001
Epoch 12: f = 1.660821e-02, step = 0.001
Epoch 13: f = 1.525271e-02, step = 0.001
Epoch 14: f = 1.399967e-02, step = 0.001
Epoch 15: f = 1.185572e-02, step = 0.001
Epoch 16: f = 9.517354e-03, step = 0.001
Epoch 17: f = 6.242648e-03, step = 0.001
Epoch 18: f = 3.505938e-03, step = 0.001
Epoch 19: f = 1.700016e-03, step = 0.001
Epoch 20: f = 7.406634e-04, step = 0.001
Epoch 21: f = 4.796099e-04, step = 0.001
Epoch 22: f = 4.453608e-04, step = 0.001
Epoch 23: f = 1.749952e-04, step = 0.001
Epoch 24: f = 1.949576e-04, step = 0.001, nfails = 1 (resetting to solution from last epoch)
Epoch 25: f = 6.149062e-05, step = 0.0001
Epoch 26: f = 7.074476e-05, step = 0.0001, nfails = 2 (resetting to solution from last epoch)
End Main Loop

Final f: 6.1491e-05
Setup time: 0.01 seconds
Optimization time: 1.13 seconds
Number of epochs: 26
Total iterations: 2600
Exit condition: Number of failed epochs > 1

		*** Final F: 6.15e-05, Score: 0.96 ***