All-at-once optimization for CP tensor decomposition (with Poblano)

We explain how to use cp_opt_legacy with the Poblano toolbox. The default is to use L-BFGS-B (not Poblano), which is described here.

Contents

Newer version available

This documentation is for the legacy implementation of CP-OPT. We recommend using the newer version, described here.

Poblano Optimization Toolbox

Check that you have Poblano 1.1 (or later) installed. The output of your 'ver' command should look something like the following.

ver
-----------------------------------------------------------------------------------------------------
MATLAB Version: 9.11.0.1873467 (R2021b) Update 3
MATLAB License Number: 40941043
Operating System: Microsoft Windows 11 Pro Version 10.0 (Build 22000)
Java Version: Java 1.8.0_202-b08 with Oracle Corporation Java HotSpot(TM) 64-Bit Server VM mixed mode
-----------------------------------------------------------------------------------------------------
MATLAB                                                Version 9.11        (R2021b) 
Deep Learning Toolbox                                 Version 14.3        (R2021b) 
Image Processing Toolbox                              Version 11.4        (R2021b) 
Mapping Toolbox                                       Version 5.2         (R2021b) 
Optimization Toolbox                                  Version 9.2         (R2021b) 
Parallel Computing Toolbox                            Version 7.5         (R2021b) 
Poblano Toolbox (Sandia National Labs)                Version 1.2                  
Statistics and Machine Learning Toolbox               Version 12.2        (R2021b) 
Symbolic Math Toolbox                                 Version 9.0         (R2021b) 
Tensor Toolbox (Sandia National Labs)                 Version 3.3.a-dev    (R2022a)
Text Analytics Toolbox                                Version 1.8         (R2021b) 
Wavelet Toolbox                                       Version 6.0         (R2021b) 

Create an example problem.

Create an example 50 x 40 x 30 tensor with rank 5 and add 10% noise.

rng('default'); %<- Setting random seed for reproducibility of this script
R = 5;
info = create_problem('Size', [50 40 30], 'Num_Factors', R, 'Noise', 0.10);
X = info.Data;
M_true = info.Soln;

Create initial guess using 'nvecs'

M_init = create_guess('Data', X, 'Num_Factors', R, ...
    'Factor_Generator', 'nvecs');

Set up the optimization parameters

It's genearlly a good idea to consider the parameters of the optimization method. The default options may be either too stringent or not stringent enough. The most important options to consider are detailed here.

% Get the defaults
ncg_opts = ncg('defaults');
% Tighten the stop tolerance (norm of gradient). This is often too large.
ncg_opts.StopTol = 1.0e-6;
% Tighten relative change in function value tolearnce. This is often too large.
ncg_opts.RelFuncTol = 1.0e-20;
% Increase the number of iterations.
ncg_opts.MaxIters = 10^4;
% Only display every 10th iteration
ncg_opts.DisplayIters = 10;
% Display the final set of options
ncg_opts
ncg_opts = 

  struct with fields:

                   Display: 'iter'
              DisplayIters: 10
           LineSearch_ftol: 1.0000e-04
           LineSearch_gtol: 0.0100
    LineSearch_initialstep: 1
         LineSearch_maxfev: 20
         LineSearch_method: 'more-thuente'
         LineSearch_stpmax: 1.0000e+15
         LineSearch_stpmin: 1.0000e-15
           LineSearch_xtol: 1.0000e-15
              MaxFuncEvals: 10000
                  MaxIters: 10000
                RelFuncTol: 1.0000e-20
              RestartIters: 20
                 RestartNW: 0
              RestartNWTol: 0.1000
                   StopTol: 1.0000e-06
                 TraceFunc: 0
            TraceFuncEvals: 0
                 TraceGrad: 0
             TraceGradNorm: 0
              TraceRelFunc: 0
                    TraceX: 0
                    Update: 'PR'

Call the cp_opt_legacy method

Here is an example call to the cp_opt_legacy method. By default, each iteration prints the least squares fit function value (being minimized) and the norm of the gradient. The meaning of any line search warnings can be checked via doc cvsrch.

[M,~,output] = cp_opt_legacy(X, R, 'init', M_init, ...
    'opt', 'ncg', 'opt_options', ncg_opts);
 Iter  FuncEvals       F(X)          ||G(X)||/N        
------ --------- ---------------- ----------------
     0         1   56854.34228466       0.73889204
    10        74     560.29237101       0.16300955
    20       112     555.75709148       0.01918674
    30       151     555.35904396       0.03192486
    40       197     554.86989841       0.02453038
    50       232     554.60386950       0.02480723
    60       278     554.20899944       0.01806078
    70       314     554.02850123       0.01538631
    80       361     553.85012049       0.01281683
    90       393     553.77782494       0.00940403
   100       430     553.70149389       0.01093787
   110       460     553.65835815       0.00742205
   120       502     553.60426584       0.00672556
   130       535     553.57456837       0.00716477
   140       573     553.54912769       0.00722112
   150       602     553.53638723       0.00491183
   160       643     553.50901424       0.00771484
   170       674     553.49543150       0.00529711
   180       711     553.47556496       0.00555256
   190       741     553.46331427       0.00427601
   200       782     553.42704100       0.01005031
   210       810     553.40560402       0.00514693
   220       845     553.38994007       0.00366819
   230       875     553.38210100       0.00463339
   240       915     553.35647341       0.00945372
   250       943     553.33916747       0.00490671
   260       976     553.32642617       0.00409114
   270      1004     553.32141282       0.00365606
   280      1042     553.30667908       0.00489423
   290      1069     553.29806451       0.00440921
   300      1098     553.29113580       0.00322205
   310      1124     553.28734544       0.00291091
   320      1159     553.27683196       0.00513674
   330      1185     553.27171789       0.00190590
   340      1213     553.26887365       0.00168944
   350      1235     553.26750075       0.00202764
   360      1264     553.26520481       0.00205382
   370      1288     553.26396789       0.00118306
   380      1313     553.26303978       0.00103313
   390      1336     553.26245272       0.00090600
   400      1365     553.26095190       0.00251028
   410      1387     553.26027289       0.00056015
   420      1410     553.26000008       0.00024910
   430      1434     553.25982687       0.00052903
   440      1461     553.25942739       0.00108626
   450      1484     553.25916245       0.00036871
   460      1508     553.25905558       0.00056165
   470      1529     553.25898884       0.00032671
   480      1553     553.25889059       0.00060876
   490      1574     553.25882766       0.00032151
   500      1596     553.25876444       0.00023671
   510      1618     553.25872303       0.00033837
   520      1644     553.25862132       0.00054118
   530      1664     553.25854497       0.00028158
   540      1686     553.25847388       0.00037617
   550      1706     553.25843771       0.00022447
   560      1729     553.25837133       0.00021679
   570      1750     553.25829606       0.00012975
   580      1785     553.25685112       0.00119009
   590      1811     553.25583438       0.00062378
   600      1833     553.25570368       0.00034460
   610      1854     553.25565513       0.00030809
   620      1881     553.25553904       0.00041974
   630      1902     553.25548356       0.00026437
   640      1926     553.25544371       0.00023281
   650      1946     553.25542085       0.00020827
   660      1972     553.25536368       0.00042064
   670      1993     553.25532712       0.00015715
   680      2013     553.25529804       0.00017142
   690      2034     553.25527498       0.00022332
   700      2058     553.25520793       0.00065962
   710      2079     553.25515381       0.00014658
   720      2103     553.25512827       0.00019213
   730      2124     553.25511761       0.00014752
   740      2148     553.25509165       0.00019916
   750      2169     553.25507415       0.00017411
   760      2191     553.25505347       0.00019845
   770      2211     553.25503991       0.00016042
   780      2234     553.25501109       0.00027521
   790      2254     553.25499230       0.00013643
   800      2274     553.25497366       0.00015664
   810      2294     553.25496352       0.00014402
   820      2315     553.25494449       0.00020764
   830      2335     553.25493315       0.00013404
   840      2355     553.25492028       0.00009366
   850      2375     553.25491087       0.00014055
   860      2397     553.25487921       0.00038256
   870      2419     553.25485605       0.00008978
   880      2439     553.25484348       0.00004687
   890      2459     553.25483328       0.00010907
   900      2482     553.25479714       0.00049832
   910      2502     553.25476720       0.00010257
   920      2525     553.25473442       0.00038630
   930      2545     553.25471159       0.00008553
   940      2568     553.25470406       0.00010964
   950      2590     553.25469722       0.00008719
   960      2615     553.25467640       0.00020325
   970      2637     553.25465924       0.00010910
   980      2659     553.25465415       0.00007850
   990      2681     553.25464979       0.00007725
  1000      2707     553.25463371       0.00015462
  1010      2729     553.25462078       0.00008769
  1020      2750     553.25461714       0.00006284
  1030      2772     553.25461466       0.00006613
  1040      2798     553.25459825       0.00011030
  1050      2821     553.25458631       0.00007644
  1060      2841     553.25458426       0.00004516
  1070      2861     553.25458259       0.00007473
  1080      2885     553.25457317       0.00014742
  1090      2905     553.25456670       0.00008952
  1100      2925     553.25456407       0.00004139
  1110      2945     553.25456311       0.00005391
  1120      2970     553.25455730       0.00015353
  1130      2990     553.25455358       0.00006328
  1140      3010     553.25455235       0.00003214
  1150      3030     553.25455207       0.00003521
  1160      3055     553.25454938       0.00010005
  1170      3075     553.25454778       0.00003822
  1180      3095     553.25454697       0.00001848
  1190      3117     553.25454661       0.00002347
  1200      3141     553.25454050       0.00017036
  1210      3161     553.25453748       0.00001725
  1220      3181     553.25453722       0.00001960
  1230      3201     553.25453716       0.00001126
  1240      3223     553.25453690       0.00002290
  1250      3243     553.25453672       0.00001412
  1260      3263     553.25453662       0.00001600
  1270      3283     553.25453655       0.00000855
  1280      3305     553.25453640       0.00001281
  1290      3325     553.25453625       0.00000304
  1300      3346     553.25453584       0.00000693
  1310      3366     553.25453550       0.00000215
  1320      3388     553.25453523       0.00000437
  1330      3408     553.25453501       0.00000186
  1340      3428     553.25453497       0.00000191
  1350      3448     553.25453492       0.00000141
  1360      3469     553.25453477       0.00000349
  1370      3489     553.25453464       0.00000227
  1380      3509     553.25453460       0.00000174
  1390      3529     553.25453457       0.00000172
  1400      3550     553.25453447       0.00000432
  1410      3570     553.25453438       0.00000233
  1420      3590     553.25453435       0.00000194
  1430      3610     553.25453433       0.00000125
  1440      3630     553.25453426       0.00000277
  1450      3650     553.25453419       0.00000122
  1460      3670     553.25453417       0.00000133
  1461      3672     553.25453417       0.00000079

Check the output

It's important to check the output of the optimization method. In particular, it's worthwhile to check the exit flag. A zero (0) indicates successful termination with the gradient smaller than the specified StopTol, and a three (3) indicates a successful termination where the change in function value is less than RelFuncTol. The meaning of any other flags can be checked via doc poblano_params.

exitflag = output.ExitFlag
exitflag =

     0

The fit is the percentage of the data that is explained by the model. Because we have noise, we do not expect the fit to be perfect.

fit = output.Fit
fit =

   99.0220

Evaluate the output

We can "score" the similarity of the model computed by CP and compare that with the truth. The score function on ktensor's gives a score in [0,1] with 1 indicating a perfect match. Because we have noise, we do not expect the fit to be perfect. See doc score for more details.

scr = score(M,M_true)
scr =

    0.7998

Overfitting example

Consider the case where we don't know R in advance. We might guess too high. Here we show a case where we guess R+1 factors rather than R.

% Generate initial guess of the corret size
M_plus_init = create_guess('Data', X, 'Num_Factors', R+1, ...
    'Factor_Generator', 'nvecs');
% Loosen the stop tolerance (norm of gradient).
ncg_opts.StopTol = 1.0e-2;
% Run the algorithm
[M_plus,~,output] = cp_opt_legacy(X, R+1, 'init', M_plus_init, ...
    'opt', 'ncg', 'opt_options', ncg_opts);
exitflag = output.ExitFlag
fit = output.Fit
 Iter  FuncEvals       F(X)          ||G(X)||/N        
------ --------- ---------------- ----------------
     0         1   56854.95767322       0.61575601
    10        74     560.53901958       0.13589655
    15        95     556.07779574       0.00781867

exitflag =

     0


fit =

   99.0170

% Check the answer (1 is perfect)
scr = score(M_plus, M_true)
scr =

    0.7998