Weighted optimization for CP tensor decomposition with incomplete data

We explain how to use cp_wopt with the POBLANO toolbox. The method is described in the following article:

Contents

Important Information

It is critical to zero out the values in the missing entries of the data tensor. This can be done by calling cp_wopt(X.*P,P,...). This is a frequent source of errors in using this method.

Create an example problem with missing data.

Here we have 25% missing data and 10% noise.

R = 2;
info = create_problem('Size', [15 10 5], 'Num_Factors', R, ...
    'M', 0.25, 'Noise', 0.10);
X = info.Data;
P = info.Pattern;
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_wopt method

Here is an example call to the cp_opt 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_wopt(X, P, R, 'init', M_init, ...
    'opt', 'ncg', 'opt_options', ncg_opts);
Running CP-WOPT...
Time for zeroing out masked entries of data tensor is 4.20e-04 seconds.
(If zeroing is done in preprocessing, set 'skip_zeroing' to true.)
 Iter  FuncEvals       F(X)          ||G(X)||/N        
------ --------- ---------------- ----------------
     0         1      25.10574260       0.17813973
    10        35       0.66553174       0.02681209
    20        68       0.28614990       0.00668033
    30        91       0.27446548       0.00071779
    40       111       0.27426663       0.00011737
    50       131       0.27425951       0.00003687
    60       151       0.27425900       0.00000478
    70       171       0.27425899       0.00000078

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

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.9427

Create a SPARSE example problem with missing data.

Here we have 95% missing data and 10% noise.

R = 2;
info = create_problem('Size', [150 100 50], 'Num_Factors', R, ...
    'M', 0.95, 'Sparse_M', true, 'Noise', 0.10);
X = info.Data;
P = info.Pattern;
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_wopt method

Here is an example call to the cp_opt 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_wopt(X, P, R, 'init', M_init, ...
    'opt', 'ncg', 'opt_options', ncg_opts);
Running CP-WOPT...
Time for zeroing out masked entries of data tensor is 6.03e-02 seconds.
(If zeroing is done in preprocessing, set 'skip_zeroing' to true.)
 Iter  FuncEvals       F(X)          ||G(X)||/N        
------ --------- ---------------- ----------------
     0         1    6245.25896387       0.04522498
    10        60    1783.65135517       0.03570645
    20       120     622.38293705       0.61597552
    30       163      61.17553615       0.00263828
    40       185      61.16563723       0.00000672
    44       193      61.16563716       0.00000080

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

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.9990