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