# 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:

E. Acar, D. M. Dunlavy, T. G. Kolda and M. Mørup, Scalable Tensor Factorizations for Incomplete Data, Chemometrics and Intelligent Laboratory Systems 106(1):41-56, March 2011 (doi:10.1016/j.chemolab.2010.08.004)

## Contents

- Important Information
- Create an example problem with missing data.
- Create initial guess using 'nvecs'
- Set up the optimization parameters
- Call the `cp_wopt` method
- Check the output
- Evaluate the output
- Create a SPARSE example problem with missing data.
- Create initial guess using 'nvecs'
- Set up the optimization parameters
- Call the `cp_wopt` method
- Check the output
- Evaluate the output

## 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 7.52e-04 seconds.
    (If zeroing is done in preprocessing, set 'skip_zeroing' to true.)
     Iter  FuncEvals       F(X)          ||G(X)||/N
    ------ --------- ---------------- ----------------
         0         1      42.12686745       0.23070139
        10        37       3.20399740       0.01068598
        20        74       1.86647166       0.01496155
        30       102       1.75795458       0.00124314
        40       122       1.75699489       0.00030707
        50       164       1.59387469       0.03315542
        60       202       0.62497618       0.03131283
        70       229       0.34450247       0.00826425
        80       254       0.31352609       0.00161307
        90       275       0.31288321       0.00018934
       100       295       0.31287112       0.00003616
       110       315       0.31287061       0.00000820
       120       335       0.31287058       0.00000100

## 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.9841

## 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 3.01e-02 seconds.
    (If zeroing is done in preprocessing, set 'skip_zeroing' to true.)
     Iter  FuncEvals       F(X)          ||G(X)||/N
    ------ --------- ---------------- ----------------
         0         1     561.51059983       0.01729762
        10        51      17.54710510       0.00825676
        20        82      16.81335392       0.00164784
        30       112      16.63583752       0.00182422
        40       149      16.09905940       0.00369331
        50       186      12.61496613       0.01236469
        60       223       5.78600003       0.00503093
        70       251       5.42873077       0.00059175
        80       274       5.42075260       0.00006912
        90       294       5.42061561       0.00001694
       100       314       5.42061124       0.00000194
       102       318       5.42061122       0.00000084

## 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.9972