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

- E. Acar, D. M. Dunlavy and T. G. Kolda, A Scalable Optimization Approach for Fitting Canonical Tensor Decompositions, J. Chemometrics 25(2):67-86, February 2011 (doi:10.1002/cem.1335)

## Contents

- Third-party optimization software
- Check that the software is installed.
- Create an example problem.
- Create initial guess using 'nvecs'
- Call the `cp_opt` method
- Check the output
- Evaluate the output
- Overfitting example
- Nonnegative factorization
- Create an example problem.
- Generate initial guess of the corret size
- Call the `cp_opt` method
- Check the output
- Evaluate the output

## Third-party optimization software

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

- L-BFGS-B by Stephen Becker (preferred), or
- POBLANO Version 1.1 by Evrim Acar, Daniel Dunlavy, and Tamara Kolda.

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