All-at-once optimization for CP tensor decomposition

We explain how to use cp_opt function which implements the CP-OPT method that fits the CP model using direct or all-at-once optimization. This is in contrast to the cp_als function which implements the CP-ALS that fits the CP model using alternating optimization. The CP-OPT method is described in the following reference:

Contents

Third-party optimization software

The cp_opt method uses third-party optimization software to do the optimization. You can use either

The remainder of these instructions assume L-BFGS-B is being used. See here for instructions on using cp_opt with Poblano.

Check that the software is installed.

Be sure that lbfgsb is in your path.

help lbfgsb
  x = lbfgsb( fcn, l, u )
    uses the lbfgsb v.3.0 library (fortran files must be installed;
        see compile_mex.m ) which is the L-BFGS-B algorithm.
    The algorithm is similar to the L-BFGS quasi-Newton algorithm,
    but also handles bound constraints via an active-set type iteration.
    This version is based on the modified C code L-BFGS-B-C, and so has 
    a slightly different calling syntax than previous versions.
 
   The minimization problem that it solves is:
        min_x  f(x)     subject to   l <= x <= u
 
  'fcn' is a function handle that accepts an input, 'x',
    and returns two outputs, 'f' (function value), and 'g' (function gradient).
 
  'l' and 'u' are column-vectors of constraints. Set their values to Inf
    if you want to ignore them. (You can set some values to Inf, but keep
    others enforced).
 
  The full format of the function is:
  [x,f,info] = lbfgsb( fcn, l, u, opts )
    where the output 'f' has the value of the function f at the final iterate
    and 'info' is a structure with useful information
        (self-explanatory, except for info.err. The first column of info.err
         is the history of the function values f, and the second column
         is the history of norm( gradient, Inf ).  )
 
    The 'opts' structure allows you to pass further options.
    Possible field name values:
 
        opts.x0     The starting value (default: all zeros)
        opts.m      Number of limited-memory vectors to use in the algorithm
                        Try 3 <= m <= 20. (default: 5 )
        opts.factr  Tolerance setting (see this source code for more info)
                        (default: 1e7 ). This is later multiplied by machine epsilon
        opts.pgtol  Another tolerance setting, relating to norm(gradient,Inf)
                        (default: 1e-5)
        opts.maxIts         How many iterations to allow (default: 100)
        opts.maxTotalIts    How many iterations to allow, including linesearch iterations
                        (default: 5000)
        opts.printEvery     How often to display information (default: 1)
        opts.errFcn         A function handle (or cell array of several function handles)
                        that computes whatever you want. The output will be printed
                        to the screen every 'printEvery' iterations. (default: [] )
                        Results saved in columns 3 and higher of info.err variable
 
  Stephen Becker, srbecker@alumni.caltech.edu
  Feb 14, 2012
  Updated Feb 21 2015, Stephen Becker, stephen.becker@colorado.edu

Create an example problem.

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

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');

Call the cp_opt 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.

[M,M0,output] = cp_opt(X, R, 'init', M_init);
Iter    10, f(x) = 1.239376e+03, ||grad||_infty = 9.52e+01
Iter    20, f(x) = 1.204254e+03, ||grad||_infty = 6.93e-01
Iter    30, f(x) = 1.033479e+03, ||grad||_infty = 9.20e+01
Iter    40, f(x) = 5.932041e+02, ||grad||_infty = 8.66e+00
Iter    50, f(x) = 5.891379e+02, ||grad||_infty = 2.58e-01
Iter    60, f(x) = 5.891340e+02, ||grad||_infty = 7.49e-03
Iter    61, f(x) = 5.891340e+02, ||grad||_infty = 5.58e-03

Check the output

It's important to check the output of the optimization method. In particular, it's worthwhile to check the exit message.

exitmsg = output.ExitMsg
exitmsg =

    'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH.'

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

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

Overfitting example

Re-using the same example as before, 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');
% Run the algorithm
[M_plus,~,output] = cp_opt(X, R+1, 'init', M_plus_init);
exitmsg = output.ExitMsg
fit = output.Fit
Iter    10, f(x) = 1.239651e+03, ||grad||_infty = 9.54e+01
Iter    20, f(x) = 1.204230e+03, ||grad||_infty = 1.30e+00
Iter    30, f(x) = 1.202606e+03, ||grad||_infty = 1.08e+01
Iter    40, f(x) = 7.488881e+02, ||grad||_infty = 1.42e+02
Iter    50, f(x) = 5.919198e+02, ||grad||_infty = 4.48e+00
Iter    60, f(x) = 5.889069e+02, ||grad||_infty = 1.93e+00
Iter    70, f(x) = 5.881156e+02, ||grad||_infty = 8.79e-01
Iter    80, f(x) = 5.877669e+02, ||grad||_infty = 9.43e-01
Iter    90, f(x) = 5.874421e+02, ||grad||_infty = 2.62e+00
Iter   100, f(x) = 5.870874e+02, ||grad||_infty = 2.98e-01
Iter   110, f(x) = 5.869636e+02, ||grad||_infty = 2.99e+00
Iter   120, f(x) = 5.868945e+02, ||grad||_infty = 4.64e-01
Iter   130, f(x) = 5.868383e+02, ||grad||_infty = 4.58e-01
Iter   140, f(x) = 5.867528e+02, ||grad||_infty = 3.39e-01
Iter   150, f(x) = 5.867107e+02, ||grad||_infty = 7.51e-01
Iter   160, f(x) = 5.866676e+02, ||grad||_infty = 2.29e-01
Iter   170, f(x) = 5.866355e+02, ||grad||_infty = 2.30e-01
Iter   180, f(x) = 5.866268e+02, ||grad||_infty = 1.77e-01
Iter   190, f(x) = 5.866142e+02, ||grad||_infty = 1.19e-01
Iter   200, f(x) = 5.865911e+02, ||grad||_infty = 2.12e-01
Iter   210, f(x) = 5.865839e+02, ||grad||_infty = 1.83e-01
Iter   220, f(x) = 5.865713e+02, ||grad||_infty = 3.90e-01
Iter   230, f(x) = 5.865528e+02, ||grad||_infty = 8.10e-01
Iter   240, f(x) = 5.865372e+02, ||grad||_infty = 1.05e-01
Iter   250, f(x) = 5.865300e+02, ||grad||_infty = 2.82e-01
Iter   260, f(x) = 5.865229e+02, ||grad||_infty = 1.66e-01
Iter   270, f(x) = 5.865135e+02, ||grad||_infty = 1.06e-01
Iter   280, f(x) = 5.865011e+02, ||grad||_infty = 8.57e-02
Iter   290, f(x) = 5.864903e+02, ||grad||_infty = 3.23e-01
Iter   300, f(x) = 5.864816e+02, ||grad||_infty = 1.39e-01
Iter   310, f(x) = 5.864728e+02, ||grad||_infty = 1.16e-01
Iter   320, f(x) = 5.864642e+02, ||grad||_infty = 1.15e-01
Iter   330, f(x) = 5.864525e+02, ||grad||_infty = 3.76e-01
Iter   340, f(x) = 5.864428e+02, ||grad||_infty = 2.33e-01
Iter   350, f(x) = 5.864069e+02, ||grad||_infty = 1.92e-01
Iter   360, f(x) = 5.863788e+02, ||grad||_infty = 4.34e-01
Iter   370, f(x) = 5.863443e+02, ||grad||_infty = 1.84e-01
Iter   380, f(x) = 5.863030e+02, ||grad||_infty = 1.94e-01
Iter   390, f(x) = 5.862826e+02, ||grad||_infty = 2.97e-01
Iter   400, f(x) = 5.862708e+02, ||grad||_infty = 1.20e-01
Iter   410, f(x) = 5.862519e+02, ||grad||_infty = 4.23e-01
Iter   420, f(x) = 5.862373e+02, ||grad||_infty = 1.72e-01
Iter   430, f(x) = 5.862314e+02, ||grad||_infty = 4.63e-01
Iter   440, f(x) = 5.862213e+02, ||grad||_infty = 9.90e-02
Iter   450, f(x) = 5.862163e+02, ||grad||_infty = 6.01e-01
Iter   460, f(x) = 5.862098e+02, ||grad||_infty = 1.04e-01
Iter   470, f(x) = 5.862055e+02, ||grad||_infty = 7.14e-02
Iter   480, f(x) = 5.862028e+02, ||grad||_infty = 1.27e-01
Iter   490, f(x) = 5.862019e+02, ||grad||_infty = 8.58e-02
Iter   500, f(x) = 5.861960e+02, ||grad||_infty = 1.39e-01
Iter   510, f(x) = 5.861923e+02, ||grad||_infty = 1.04e-01
Iter   520, f(x) = 5.861899e+02, ||grad||_infty = 5.11e-02
Iter   530, f(x) = 5.861883e+02, ||grad||_infty = 7.26e-02
Iter   540, f(x) = 5.861877e+02, ||grad||_infty = 2.21e-02
Iter   550, f(x) = 5.861871e+02, ||grad||_infty = 4.50e-02
Iter   560, f(x) = 5.861869e+02, ||grad||_infty = 6.50e-02
Iter   570, f(x) = 5.861868e+02, ||grad||_infty = 2.86e-02
Iter   580, f(x) = 5.861860e+02, ||grad||_infty = 2.24e-02
Iter   590, f(x) = 5.861859e+02, ||grad||_infty = 2.18e-02
Iter   600, f(x) = 5.861857e+02, ||grad||_infty = 1.93e-02
Iter   610, f(x) = 5.861856e+02, ||grad||_infty = 3.24e-02
Iter   620, f(x) = 5.861854e+02, ||grad||_infty = 1.61e-02
Iter   630, f(x) = 5.861853e+02, ||grad||_infty = 2.29e-02
Iter   640, f(x) = 5.861852e+02, ||grad||_infty = 8.46e-03
Iter   650, f(x) = 5.861851e+02, ||grad||_infty = 9.60e-03
Iter   660, f(x) = 5.861851e+02, ||grad||_infty = 2.05e-02
Iter   664, f(x) = 5.861851e+02, ||grad||_infty = 3.68e-02

exitmsg =

    'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH.'


fit =

   99.0247

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

    0.9983

Nonnegative factorization

We can employ lower bounds to get a nonnegative factorization.

Create an example problem.

Create an example 50 x 40 x 30 tensor with rank 5 and add 10% noise. We select nonnegative factor matrices and lambdas. The create_problem doesn't really know how to add noise without going negative, so we hack it to make the observed tensor be nonzero.

R = 5;
info = create_problem('Size', [50 40 30], 'Num_Factors', R, 'Noise', 0.10,...
    'Factor_Generator', 'rand', 'Lambda_Generator', 'rand');
X = info.Data .* (info.Data > 0); % Force it to be nonnegative
M_true = info.Soln;

Generate initial guess of the corret size

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

Call the cp_opt method

Here we specify a lower bound of zero with the last two arguments.

[M,M0,output] = cp_opt(X, R, 'init', M_init,'lower',0);
Iter    10, f(x) = 1.373055e+02, ||grad||_infty = 1.52e+01
Iter    20, f(x) = 5.094910e+01, ||grad||_infty = 4.21e+00
Iter    30, f(x) = 4.495606e+01, ||grad||_infty = 2.25e+00
Iter    40, f(x) = 4.154133e+01, ||grad||_infty = 1.92e+00
Iter    50, f(x) = 4.067009e+01, ||grad||_infty = 2.14e+00
Iter    60, f(x) = 3.996698e+01, ||grad||_infty = 1.80e+00
Iter    70, f(x) = 3.964063e+01, ||grad||_infty = 1.45e+00
Iter    80, f(x) = 3.941584e+01, ||grad||_infty = 5.82e-01
Iter    90, f(x) = 3.935186e+01, ||grad||_infty = 2.06e-01
Iter   100, f(x) = 3.933107e+01, ||grad||_infty = 1.76e-01
Iter   110, f(x) = 3.932266e+01, ||grad||_infty = 1.61e-01
Iter   120, f(x) = 3.931633e+01, ||grad||_infty = 2.04e-01
Iter   130, f(x) = 3.931059e+01, ||grad||_infty = 2.72e-01
Iter   140, f(x) = 3.930927e+01, ||grad||_infty = 2.68e-01
Iter   150, f(x) = 3.930873e+01, ||grad||_infty = 2.33e-01
Iter   160, f(x) = 3.930837e+01, ||grad||_infty = 1.91e-01
Iter   170, f(x) = 3.930823e+01, ||grad||_infty = 1.91e-01
Iter   180, f(x) = 3.930818e+01, ||grad||_infty = 1.87e-01
Iter   190, f(x) = 3.930817e+01, ||grad||_infty = 1.84e-01
Iter   200, f(x) = 3.930816e+01, ||grad||_infty = 1.81e-01
Iter   202, f(x) = 3.930816e+01, ||grad||_infty = 1.81e-01

Check the output

exitmsg = output.ExitMsg
exitmsg =

    'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH.'

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

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