Home > examples > low_rank_dist_completion.m

low_rank_dist_completion

PURPOSE ^

Perform low-rank distance matrix completion w/ automatic rank detection.

SYNOPSIS ^

function [Y, infos, problem_description] = low_rank_dist_completion(problem_description)

DESCRIPTION ^

 Perform low-rank distance matrix completion w/ automatic rank detection.

 function Y = low_rank_dist_completion(problem_description)
 function [Y, infos, out_problem_description] = low_rank_dist_completion(problem_description)

 It implements the ideas of Journee, Bach, Absil and Sepulchre, SIOPT, 2010,
 applied to the problem of low-rank Euclidean distance matrix completion.
 The details are in the paper "Low-rank optimization for distance matrix completion",
 B. Mishra, G. Meyer, and R. Sepulchre, IEEE CDC, 2011.

 Paper link: http://arxiv.org/abs/1304.6663.

 Input:
 -------

 problem_description: The problem structure with the description of the problem.


 - problem_description.data_train: Data structure for known distances that are used to learn a low-rank model.
                                   It contains the 3 fields that are shown
                                   below. An empty "data_train" structure
                                   will generate the 3d Helix instance.

       -- data_train.entries:      A column vector consisting of known
                                   distances. An empty "data_train.entries"
                                   field will generate the 3d Helix
                                   instance.

       -- data_train.rows:         The row position of th corresponding
                                   distances. An empty "data_train.rows"
                                   field will generate the 3d Helix
                                   instance.

       -- data_train.cols:         The column position of th corresponding
                                   distances. An empty "data_train.cols"
                                   field will generate the 3d Helix
                                   instance.



 - problem_description.data_test:  Data structure to compute distances for the "unknown" (to the algorithm) distances.
                                   It contains the 3 fields that are shown
                                   below. An empty "data_test" structure
                                   will not compute the test error.

       -- data_test.entries:       A column vector consisting of "unknown" (to the algorithm)
                                   distances. An empty "data_test.entries"
                                   field will not compute the test error.
       -- data_test.rows:          The row position of th corresponding
                                   distances. An empty "data_test.rows"
                                   field will not compute the test error.
       -- data_test.cols:          The column position of th corresponding
                                   distances. An empty "data_test.cols"
                                   field will not compute the test error.



 - problem_description.n:          The number of data points. An empty
                                   "n", but complete "data_train" structure
                                   will lead to an error, to avoid
                                   potential data inconsistency.





 - problem_description.rank_initial: Starting rank. By default, it is 1.



 - problem_description.rank_max:     Maximum rank. By default, it is equal to
                                     "problem_description.n".




 - problem_description.params:  Structure array containing algorithm
                                parameters for stopping criteria.
       -- params.abstolcost:    Tolerance on absolute value of cost.
                                By default, it is 1e-3.


       -- params.reltolcost:    Tolerance on absolute value of cost.
                                By default, it is 1e-3.
       -- params.tolgradnorm:   Tolerance on the norm of the gradient.
                                By default, it is 1e-5.
       -- params.maxiter:       Maximum number of fixe-rank iterations.
                                By default, it is 100.
       -- params.tolSmin:       Tolerance on smallest eigenvalue of Sy,
                                the dual variable.
                                By default, it is 1e-5.
       -- params.tolrankdeficiency:   Tolerance on the
                                      smallest singular value of Y.
                                      By default, it is 1e-3.
       -- params.solver:        Fixed-rank algorithm. Options are
                                '@trustregions' for trust-regions,
                                '@conjugategradient' for conjugate gradients,
                                '@steepestdescent' for steepest descent.
                                 By default, it is '@trustregions'.


 Output:
 --------

   Y:                    n-by-r solution matrix of rank r.
   infos:                Structure array with computed statistics.
   problem_description:  Structure array with used problem description.



 Please cite the Manopt paper as well as the research paper:
     @InProceedings{mishra2011dist,
       Title        = {Low-rank optimization for distance matrix completion},
       Author       = {Mishra, B. and Meyer, G. and Sepulchre, R.},
       Booktitle    = {{50th IEEE Conference on Decision and Control}},
       Year         = {2011},
       Organization = {{IEEE CDC}}
     }

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [Y, infos, problem_description] =  low_rank_dist_completion(problem_description)
0002 % Perform low-rank distance matrix completion w/ automatic rank detection.
0003 %
0004 % function Y = low_rank_dist_completion(problem_description)
0005 % function [Y, infos, out_problem_description] = low_rank_dist_completion(problem_description)
0006 %
0007 % It implements the ideas of Journee, Bach, Absil and Sepulchre, SIOPT, 2010,
0008 % applied to the problem of low-rank Euclidean distance matrix completion.
0009 % The details are in the paper "Low-rank optimization for distance matrix completion",
0010 % B. Mishra, G. Meyer, and R. Sepulchre, IEEE CDC, 2011.
0011 %
0012 % Paper link: http://arxiv.org/abs/1304.6663.
0013 %
0014 % Input:
0015 % -------
0016 %
0017 % problem_description: The problem structure with the description of the problem.
0018 %
0019 %
0020 % - problem_description.data_train: Data structure for known distances that are used to learn a low-rank model.
0021 %                                   It contains the 3 fields that are shown
0022 %                                   below. An empty "data_train" structure
0023 %                                   will generate the 3d Helix instance.
0024 %
0025 %       -- data_train.entries:      A column vector consisting of known
0026 %                                   distances. An empty "data_train.entries"
0027 %                                   field will generate the 3d Helix
0028 %                                   instance.
0029 %
0030 %       -- data_train.rows:         The row position of th corresponding
0031 %                                   distances. An empty "data_train.rows"
0032 %                                   field will generate the 3d Helix
0033 %                                   instance.
0034 %
0035 %       -- data_train.cols:         The column position of th corresponding
0036 %                                   distances. An empty "data_train.cols"
0037 %                                   field will generate the 3d Helix
0038 %                                   instance.
0039 %
0040 %
0041 %
0042 % - problem_description.data_test:  Data structure to compute distances for the "unknown" (to the algorithm) distances.
0043 %                                   It contains the 3 fields that are shown
0044 %                                   below. An empty "data_test" structure
0045 %                                   will not compute the test error.
0046 %
0047 %       -- data_test.entries:       A column vector consisting of "unknown" (to the algorithm)
0048 %                                   distances. An empty "data_test.entries"
0049 %                                   field will not compute the test error.
0050 %       -- data_test.rows:          The row position of th corresponding
0051 %                                   distances. An empty "data_test.rows"
0052 %                                   field will not compute the test error.
0053 %       -- data_test.cols:          The column position of th corresponding
0054 %                                   distances. An empty "data_test.cols"
0055 %                                   field will not compute the test error.
0056 %
0057 %
0058 %
0059 % - problem_description.n:          The number of data points. An empty
0060 %                                   "n", but complete "data_train" structure
0061 %                                   will lead to an error, to avoid
0062 %                                   potential data inconsistency.
0063 %
0064 %
0065 %
0066 %
0067 %
0068 % - problem_description.rank_initial: Starting rank. By default, it is 1.
0069 %
0070 %
0071 %
0072 % - problem_description.rank_max:     Maximum rank. By default, it is equal to
0073 %                                     "problem_description.n".
0074 %
0075 %
0076 %
0077 %
0078 % - problem_description.params:  Structure array containing algorithm
0079 %                                parameters for stopping criteria.
0080 %       -- params.abstolcost:    Tolerance on absolute value of cost.
0081 %                                By default, it is 1e-3.
0082 %
0083 %
0084 %       -- params.reltolcost:    Tolerance on absolute value of cost.
0085 %                                By default, it is 1e-3.
0086 %       -- params.tolgradnorm:   Tolerance on the norm of the gradient.
0087 %                                By default, it is 1e-5.
0088 %       -- params.maxiter:       Maximum number of fixe-rank iterations.
0089 %                                By default, it is 100.
0090 %       -- params.tolSmin:       Tolerance on smallest eigenvalue of Sy,
0091 %                                the dual variable.
0092 %                                By default, it is 1e-5.
0093 %       -- params.tolrankdeficiency:   Tolerance on the
0094 %                                      smallest singular value of Y.
0095 %                                      By default, it is 1e-3.
0096 %       -- params.solver:        Fixed-rank algorithm. Options are
0097 %                                '@trustregions' for trust-regions,
0098 %                                '@conjugategradient' for conjugate gradients,
0099 %                                '@steepestdescent' for steepest descent.
0100 %                                 By default, it is '@trustregions'.
0101 %
0102 %
0103 % Output:
0104 % --------
0105 %
0106 %   Y:                    n-by-r solution matrix of rank r.
0107 %   infos:                Structure array with computed statistics.
0108 %   problem_description:  Structure array with used problem description.
0109 %
0110 %
0111 %
0112 % Please cite the Manopt paper as well as the research paper:
0113 %     @InProceedings{mishra2011dist,
0114 %       Title        = {Low-rank optimization for distance matrix completion},
0115 %       Author       = {Mishra, B. and Meyer, G. and Sepulchre, R.},
0116 %       Booktitle    = {{50th IEEE Conference on Decision and Control}},
0117 %       Year         = {2011},
0118 %       Organization = {{IEEE CDC}}
0119 %     }
0120 
0121 
0122 % This file is part of Manopt: www.manopt.org.
0123 % Original author: Bamdev Mishra, April 06, 2015.
0124 % Contributors: Nicolas Boumal.
0125 % Change log:
0126 %   August 30 2016 (BM):
0127 %                   Corrected some logic flaws while plotting and storing
0128 %                   rank information. A typo was also corrected.
0129 
0130     
0131     % Check problem description
0132     if ~exist('problem_description', 'var')
0133         problem_description = struct();
0134     end
0135     problem_description = check_problem_description(problem_description); % Check the problem description;
0136     
0137     
0138     % Common quantities
0139     data_train = problem_description.data_train;
0140     data_test =  problem_description.data_test;
0141     n =  problem_description.n;
0142     rank_initial = problem_description.rank_initial;
0143     rank_max =  problem_description.rank_max;
0144     params =  problem_description.params;
0145     N = data_train.nentries; % Number of known distances
0146     EIJ = speye(n);
0147     EIJ = EIJ(:, data_train.rows) - EIJ(:, data_train.cols);
0148     rr = rank_initial; % Starting rank.
0149     Y = randn(n, rr); % Random starting initialization.
0150     
0151     
0152     % Information
0153     time = [];               % Time for each iteration per rank
0154     cost = [];               % Cost at each iteration per rank
0155     test_error = [];         % Test error at each iteration per rank
0156     rank = [];               % Rank at each iteration
0157     rank_change_stats = [];  % Some stats relating the change of ranks
0158     
0159     
0160     
0161     % Main loop rank search
0162     rank_search = 0;
0163     while (rr <= rank_max), % When r = n a global min is attained for sure.
0164         rank_search = rank_search + 1;
0165         
0166         fprintf('>> Rank %d <<\n', rr);
0167         
0168         % Follow the descent direction to compute an iterate in a higher dimension
0169         if (rr > rank_initial),
0170             if isempty(restartDir), % If no restart dir avail. do a random restart
0171                 disp('No restart dir available, random restart is performed');
0172                 Y = randn(n, rr);
0173                 
0174             else % Perform a simple line-search based on the restart direction
0175                 disp('>> Line-search with restart direction');
0176                 Y(:, rr) = 0; % Append a column of zeroes
0177                 
0178                 Z = Y(data_train.rows, :) - Y(data_train.cols,:);
0179                 estimDists = sum(Z.^2, 2);
0180                 errors = (estimDists - data_train.entries);
0181                 costBefore = 0.5*mean(errors.^2);
0182                 fprintf('>> Cost before = %f\n',costBefore);
0183                 
0184                 % Simple linesearch to maintain monotonicity
0185                 problem.M = symfixedrankYYfactory(n, rr);
0186                 problem.cost = @(Y)  cost_evaluation(Y, data_train);
0187                 d = zeros(size(Y));
0188                 d(:, rr) = restartDir;
0189                 [unused, Y] = linesearch_decrease(problem, Y, d, costBefore); %#ok<ASGLU>
0190                 
0191                 Z = Y(data_train.rows, :) - Y(data_train.cols,:);
0192                 estimDists = sum(Z.^2, 2);
0193                 errors = (estimDists - data_train.entries);
0194                 costAfter = 0.5*mean(errors.^2);
0195                 
0196                 % Check for decrease
0197                 if costAfter >= costBefore - 1e-8
0198                     disp('Decrease is not sufficient, random restart');
0199                     Y = randn(n, rr);
0200                 end
0201                 
0202             end
0203             
0204         end
0205         
0206         % Fixed-rank optimization with Manopt
0207         [Y, infos_fixedrank] = low_rank_dist_completion_fixedrank(data_train, data_test, Y, params);
0208 
0209         % Some info logging
0210         thistime = [infos_fixedrank.time];
0211         if ~isempty(time)
0212             thistime = time(end) + thistime;
0213         end
0214         time = [time thistime]; %#ok<AGROW>
0215         cost = [cost [infos_fixedrank.cost]]; %#ok<AGROW>
0216         rank = [rank [infos_fixedrank.rank]]; %#ok<AGROW>
0217         rank_change_stats(rank_search).rank = rr; %#ok<AGROW>
0218         rank_change_stats(rank_search).iter = length([infos_fixedrank.cost]); %#ok<AGROW>
0219         rank_change_stats(rank_search).Y = Y; %#ok<AGROW>
0220         if isfield(infos_fixedrank, 'test_error')
0221             test_error = [test_error [infos_fixedrank.test_error]]; %#ok<AGROW>
0222         end
0223         
0224         
0225         % Evaluate gradient of the convex cost function (i.e. wrt X).
0226         Z = Y(data_train.rows, :) - Y(data_train.cols,:);
0227         estimDists = sum(Z.^2,2);
0228         errors = (estimDists - data_train.entries);
0229         
0230       
0231         % Dual variable and its minimum eigenvalue that is used to guarantee convergence.
0232         Sy = (0.5)*EIJ * sparse(1:N,1:N,2 * errors / N,N,N) * EIJ'; % "0.5" comes from 0.5 in cost evaluation
0233         
0234         
0235         % Compute smallest algebraic eigenvalue of Sy,
0236         % this gives us a descent direction for the next rank (v)
0237         % as well as a way to control progress toward the global
0238         % optimum (s_min).
0239         
0240         % Make eigs silent.
0241         opts.disp = 0;
0242         opts.issym = true;
0243         [v, s_min] = eigs(Sy, 1, 'SA', opts);
0244         
0245         
0246         % Check whether Y is rank deficient.
0247         vp = svd(Y);
0248         
0249         % Stopping criterion.
0250         fprintf('>> smin = %.3e, and min(vp) = %.3e\n',s_min,min(vp));
0251         if (s_min  > params.tolSmin) || (min(vp) < params.tolrankdeficiency),
0252             break;
0253         end
0254         
0255         % Update rank
0256         rr = rr + 1;
0257         
0258         % Compute descent direction
0259         if (s_min < -1e-10),
0260             restartDir = v;
0261         else
0262             restartDir = [];
0263         end
0264     end
0265     
0266     
0267     % Collect relevant statistics
0268     infos.time = time;
0269     infos.cost = cost;
0270     infos.rank = rank;
0271     infos.test_error = test_error;
0272     infos.rank_change_stats = rank_change_stats;
0273     
0274     % Few plots.
0275     show_plots(problem_description, infos);
0276     
0277 end
0278 
0279 
0280 
0281 
0282 %% Cost function evaluation.
0283 function val = cost_evaluation(Y, data_train)
0284     Z = Y(data_train.rows, :) - Y(data_train.cols,:);
0285     estimDists = sum(Z.^2, 2);
0286     errors = (estimDists - data_train.entries);
0287     val = 0.5*mean(errors.^2);
0288 end
0289 
0290 
0291 
0292 
0293 %% Local defaults
0294 function localdefaults = getlocaldefaults()
0295     localdefaults.abstolcost = 1e-3;
0296     localdefaults.reltolcost = 1e-3;
0297     localdefaults.tolSmin = -1e-3;
0298     localdefaults.tolrankdeficiency = 1e-3;
0299     localdefaults.tolgradnorm = 1e-5;
0300     localdefaults.maxiter = 100;
0301     localdefaults.solver = @trustregions; % Trust-regions
0302 end
0303 
0304 
0305 
0306 
0307 
0308 
0309 
0310 %% Fixed-rank optimization
0311 function [Yopt, infos] = low_rank_dist_completion_fixedrank(data_train, data_test, Y_initial, params)
0312     % Common quantities that are used often in the optimization process.
0313     [n, r] = size(Y_initial);
0314     EIJ = speye(n);
0315     EIJ = EIJ(:, data_train.rows) - EIJ(:, data_train.cols);
0316     
0317     % Create problem structure
0318     problem.M = symfixedrankYYfactory(n,  r);
0319     
0320     
0321     % Cost evaluation
0322     problem.cost = @cost;
0323     function [f, store] = cost(Y, store)
0324         if ~isfield(store, 'xij')
0325             store.xij = EIJ'*Y;
0326         end
0327         xij = store.xij;
0328         estimDists = sum(xij.^2,2);
0329         f = 0.5*mean((estimDists - data_train.entries).^2);
0330     end
0331     
0332     % Gradient evaluation
0333     problem.grad = @grad;
0334     function [g, store] = grad(Y, store)
0335         N = data_train.nentries;
0336         if ~isfield(store, 'xij')
0337             store.xij = EIJ'*Y;
0338         end
0339         xij = store.xij;
0340         estimDists = sum(xij.^2,2);
0341         g = EIJ * sparse(1:N,1:N,2 * (estimDists - data_train.entries) / N, N, N) * xij;
0342     end
0343     
0344     
0345     % Hessian evaluation
0346     problem.hess = @hess;
0347     function [Hess, store] = hess(Y, eta, store)
0348         N = data_train.nentries;
0349         if ~isfield(store, 'xij')
0350             store.xij = EIJ'*Y;
0351         end
0352         xij = store.xij;
0353         zij = EIJ'*eta;
0354         estimDists = sum(xij.^2,2);
0355         crossYZ = 2*sum(xij .* zij,2);
0356         Hess = (EIJ*sparse(1:N,1:N,2 * (estimDists - data_train.entries) / N,N,N))*zij + (EIJ*sparse(1:N,1:N,2 * crossYZ / N,N,N))*xij;
0357         Hess = problem.M.proj(Y, Hess);
0358     end
0359     
0360     
0361     %     % Check numerically whether gradient and Hessian are correct
0362     %     checkgradient(problem);
0363     %     drawnow;
0364     %     pause;
0365     %     checkhessian(problem);
0366     %     drawnow;
0367     %     pause;
0368     
0369     
0370     % When asked, ask Manopt to compute the test error at every iteration.
0371     if ~isempty(data_test)
0372         options.statsfun = @compute_test_error;
0373         EIJ_test = speye(n);
0374         EIJ_test = EIJ_test(:, data_test.rows) - EIJ_test(:, data_test.cols);
0375     end
0376     function stats = compute_test_error(problem, Y, stats) %#ok<INUSL>
0377         xij = EIJ_test'*Y;
0378         estimDists_test = sum(xij.^2,2);
0379         stats.test_error = 0.5*mean((estimDists_test - data_test.entries).^2);
0380     end
0381     
0382     
0383     % Stopping criteria options
0384     options.stopfun = @mystopfun;
0385     function stopnow = mystopfun(problem, Y, info, last) %#ok<INUSL>
0386         stopnow = (last >= 5 && (info(last-2).cost - info(last).cost < params.abstolcost || abs(info(last-2).cost - info(last).cost)/info(last).cost < params.reltolcost));
0387     end
0388     options.tolgradnorm = params.tolgradnorm;
0389     options.maxiter = params.maxiter;
0390     
0391     
0392     % Call appropriate algorithm
0393     options.solver = params.solver;
0394     [Yopt, ~, infos] = manoptsolve(problem, Y_initial, options);
0395     [infos.rank] = deal(r);
0396 end
0397 
0398 
0399 
0400 
0401 
0402 
0403 %% 3d Helix problem instance
0404 function problem_description = get_3d_Helix_instance()
0405     
0406     % Helix curve in 3d
0407     tvec = 0:2*pi/100:2*pi;
0408     tvec = tvec'; % column vector
0409     xvec = 4*cos(3*tvec);
0410     yvec = 4*sin(3*tvec);
0411     zvec = 2*tvec;
0412     Yo = [xvec, yvec, zvec];
0413     n = size(Yo, 1); % Number of points
0414     
0415     % Fraction of unknown distances
0416     fractionOfUnknown = 0.85;
0417     
0418     % True distances among points in 3d Helix
0419     trueDists = pdist(Yo)'.^2; % True distances
0420     
0421     
0422     % Add noise (set noise_level = 0 for clean measurements)
0423     noise_level = 0; % 0.01;
0424     trueDists = trueDists + noise_level * std(trueDists) * randn(size(trueDists));
0425     
0426     
0427     % Compute all pairs of indices
0428     H = tril(true(n), -1);
0429     [I, J] = ind2sub([n, n], find(H(:)));
0430     clear 'H';
0431     
0432     
0433     % Train data
0434     train = false(length(trueDists), 1);
0435     train(1:floor(length(trueDists)*(1- fractionOfUnknown))) = true;
0436     train = train(randperm(length(train)));
0437     
0438     data_train.rows = I(train);
0439     data_train.cols = J(train);
0440     data_train.entries = trueDists(train);
0441     data_train.nentries = length(data_train.entries);
0442     
0443     
0444     % Test data
0445     data_test.nentries = 1*data_train.nentries; % Depends how big data that we can handle.
0446     test = false(length(trueDists),1);
0447     test(1 : floor(data_test.nentries)) = true;
0448     test = test(randperm(length(test)));
0449     data_test.rows = I(test);
0450     data_test.cols = J(test);
0451     data_test.entries = trueDists(test);
0452     
0453     
0454     % Rank bounds
0455     rank_initial = 1; % Starting rank
0456     rank_max = n; % Maximum rank
0457     
0458     
0459     % Basic parameters used in optimization
0460     params = struct();
0461     params = mergeOptions(getlocaldefaults, params);
0462     
0463     
0464     % Problem description
0465     problem_description.data_train = data_train;
0466     problem_description.data_test = data_test;
0467     problem_description.n = n;
0468     problem_description.rank_initial = rank_initial;
0469     problem_description.rank_max = rank_max;
0470     problem_description.params = params;
0471     problem_description.Yo = Yo; % Store original Helix structure
0472 end
0473 
0474 
0475 
0476 
0477 
0478 %% Problem description check
0479 function checked_problem_description = check_problem_description(problem_description)
0480     checked_problem_description = problem_description;
0481     
0482     % Check train data
0483     if isempty(problem_description)...
0484             || ~all(isfield(problem_description,{'data_train'}) == 1)...
0485             || ~all(isfield(problem_description.data_train,{'cols', 'rows', 'entries'}) == 1)...
0486             || isempty(problem_description.data_train.cols)...
0487             || isempty(problem_description.data_train.rows)...
0488             || isempty(problem_description.data_train.entries)
0489         
0490         warning('low_rank_dist_completion:problem_description', ...
0491             'The training set is empty or not properly defined. We work with the default 3d Helix example.\n');
0492         checked_problem_description = get_3d_Helix_instance();
0493         checked_problem_description.helix_example = true;
0494         return; % No need for further check
0495     end
0496     
0497     
0498     % Check number of data points
0499     if ~isfield(problem_description, 'n')
0500         error('low_rank_dist_completion:problem_description',...
0501             'Error. The scalar corresponding to field "n" of problem description must be given. \n');
0502     end
0503     
0504     
0505     % Check initial rank
0506     if ~isfield(problem_description, 'rank_initial')...
0507             || isempty(problem_description.rank_initial)...
0508             || ~(floor(problem_description.rank_initial) == problem_description.rank_initial)
0509         warning('low_rank_dist_completion:problem_description', ...
0510             'The field "rank_initial" is not properly defined. We work with the default "1".\n');
0511         rank_initial = 1;
0512     else
0513         rank_initial = problem_description.rank_initial;
0514     end
0515     checked_problem_description.rank_initial = rank_initial;
0516     
0517     
0518     % Check maximum rank
0519     if ~isfield(problem_description, 'rank_max')...
0520             || isempty(problem_description.rank_max)...
0521             || ~(floor(problem_description.rank_max) == problem_description.rank_max)...
0522             || problem_description.rank_max > problem_description.n
0523         warning('low_rank_dist_completion:problem_description', ...
0524             'The field "rank_max" is not properly defined. We work with the default "n".\n');
0525         rank_max = problem_description.n;
0526     else
0527         rank_max = problem_description.rank_max;
0528     end
0529     checked_problem_description.rank_max = rank_max;
0530     
0531     
0532     % Check testing dataset
0533     if ~isfield(problem_description,{'data_test'})...
0534             || ~all(isfield(problem_description.data_test,{'cols', 'rows', 'entries'}) == 1)...
0535             || isempty(problem_description.data_test.cols)...
0536             || isempty(problem_description.data_test.rows)...
0537             || isempty(problem_description.data_test.entries)
0538         
0539         warning('low_rank_dist_completion:problem_description', ...
0540             'The field "data_test" is not properly defined. We work with the default "[]".\n');
0541         data_test = [];
0542     else
0543         data_test = problem_description.data_test;
0544     end
0545     checked_problem_description.data_test = data_test;
0546     
0547     
0548     % Check parameters
0549     if isfield(problem_description, 'params')
0550         params = problem_description.params;
0551     else
0552         params = struct();
0553     end
0554     params = mergeOptions(getlocaldefaults, params);
0555     checked_problem_description.params = params;
0556      
0557 end
0558 
0559 
0560 
0561 
0562 %% Show plots
0563 function  show_plots(problem_description, infos)
0564    
0565     solver = problem_description.params.solver;
0566     rank_change_stats = infos.rank_change_stats;
0567     rank_change_stats_rank = [rank_change_stats.rank];
0568     rank_change_stats_iter = [rank_change_stats.iter];
0569     rank_change_stats_iter = cumsum(rank_change_stats_iter);
0570     N = problem_description.data_train.nentries;
0571     n = problem_description.n;
0572     
0573    
0574     % Plot: train error
0575     fs = 20;
0576     figure('name', 'Training on the known distances');
0577     
0578     line(1:length([infos.cost]),log10([infos.cost]),'Marker','O','LineStyle','-','Color','blue','LineWidth',1.5);
0579     ax1 = gca;
0580     
0581     set(ax1,'FontSize',fs);
0582     xlabel(ax1,'Number of iterations','FontSize',fs);
0583     ylabel(ax1,'Cost (log scale) on known distances','FontSize',fs);
0584     
0585     ax2 = axes('Position',get(ax1,'Position'),...
0586         'XAxisLocation','top',...
0587         'YAxisLocation','right',...
0588         'Color','none',...
0589         'XColor','k');
0590     
0591     set(ax2,'FontSize',fs);
0592     line(1:length([infos.cost]),log10([infos.cost]),'Marker','O','LineStyle','-','Color','blue','LineWidth',1.5,'Parent',ax2);
0593     set(ax2,'XTick',rank_change_stats_iter(1:max(1,end-1)),...
0594         'XTickLabel',rank_change_stats_rank(1) + 1 : rank_change_stats_rank(max(1,end-1)) + 1,...
0595         'YTick',[]);
0596     
0597     set(ax2,'XGrid','on');
0598     legend(func2str(solver));
0599     title('Rank');
0600     legend 'boxoff';
0601     
0602     
0603     % Plot: test error
0604     if isfield(infos, 'test_error') && ~isempty(infos.test_error)
0605         Yo = problem_description.Yo;
0606         
0607         fs = 20;
0608         figure('name','Test error on a set of distances different from the training set');
0609         
0610         line(1:length([infos.test_error]),log10([infos.test_error]),'Marker','O','LineStyle','-','Color','blue','LineWidth',1.5);
0611         ax1 = gca;
0612         
0613         set(ax1,'FontSize',fs);
0614         xlabel(ax1,'Number of iterations','FontSize',fs);
0615         ylabel(ax1,'Cost (log scale) on testing set','FontSize',fs);
0616         
0617         ax2 = axes('Position',get(ax1,'Position'),...
0618             'XAxisLocation','top',...
0619             'YAxisLocation','right',...
0620             'Color','none',...
0621             'XColor','k');
0622         
0623         set(ax2,'FontSize',fs);
0624         line(1:length([infos.test_error]),log10([infos.test_error]),'Marker','O','LineStyle','-','Color','blue','LineWidth',1.5,'Parent',ax2);
0625         set(ax2,'XTick',rank_change_stats_iter(1:max(1,end-1)),...
0626             'XTickLabel',rank_change_stats_rank(1) + 1 : rank_change_stats_rank(max(1,end-1)) + 1,...
0627             'YTick',[]);
0628         
0629         set(ax2,'XGrid','on');
0630         legend(func2str(solver));
0631         title('Rank');
0632         legend 'boxoff';
0633         
0634         
0635         
0636     end
0637     
0638     
0639     % Plot: visualize Helix curve
0640     if isfield(problem_description, 'helix_example')
0641         jj = ceil((length(rank_change_stats_rank) + 1)/2);
0642         
0643         
0644         figure('name',['3D structure with ', num2str(N/((n^2 -n)/2)),' fraction known distances'])
0645         fs = 20;
0646         ax1 = gca;
0647         set(ax1,'FontSize',fs);
0648         subplot(jj,2,1);
0649         plot3(Yo(:,1), Yo(:,2), Yo(:,3),'*','Color', 'b','LineWidth',1.0);
0650         title('Original 3D structure');
0651         for kk = 1 : length(rank_change_stats_rank)
0652             subplot(jj, 2, kk + 1);
0653             rank_change_stats_kk = rank_change_stats(kk);
0654             Ykk = rank_change_stats_kk.Y;
0655             if size(Ykk, 2) == 1,
0656                 plot(Ykk(:,1), zeros(size(Ykk, 1)),'*','Color', 'r','LineWidth',1.0);
0657                 legend(func2str(solver))
0658                 title(['Recovery at rank ',num2str(size(Ykk, 2))]);
0659                 
0660             elseif size(Ykk, 2) == 2
0661                 plot(Ykk(:,1), Ykk(:,2),'*','Color', 'r','LineWidth',1.0);
0662                 title(['Recovery at rank ',num2str(size(Ykk, 2))]);
0663                 
0664             else  % Project onto dominant 3Dsubspace
0665                 [U1, S1, V1] = svds(Ykk, 3);
0666                 Yhat = U1*S1*V1';
0667                 plot3(Yhat(:,1), Yhat(:,2), Yhat(:,3),'*','Color', 'r','LineWidth',1.0);
0668                 title(['Recovery at rank ',num2str(size(Ykk, 2))]);
0669             end
0670             
0671             axis equal;
0672             
0673         end
0674         
0675         % Trick to add a global title to the whole subplot collection.
0676         % HitTest is disabled to make it easier to select the individual
0677         % subplots (for example, to rotate the viewing angle).
0678         ha = axes('Position',[0 0 1 1],'Xlim',[0 1],'Ylim',[0 1],'Box','off','Visible','off','Units','normalized', 'clipping' , 'off' );
0679         set(ha, 'HitTest', 'off');
0680         text(0.5, 1,['Recovery of Helix from ',num2str(N/((n^2 -n)/2)),' fraction known distances'],'HorizontalAlignment','center','VerticalAlignment', 'top');
0681     end
0682        
0683 end

Generated on Fri 08-Sep-2017 12:43:19 by m2html © 2005