Home > examples > shapefit_smoothed.m

shapefit_smoothed

PURPOSE ^

ShapeFit formulation for sensor network localization from pair directions

SYNOPSIS ^

function [T_hub, T_lsq, T_cvx] = shapefit_smoothed(V, J)

DESCRIPTION ^

 ShapeFit formulation for sensor network localization from pair directions

 function [T_hub, T_lsq, T_cvx] = shapefit_smoothed(V, J)

 This example in based on the paper http://arxiv.org/abs/1506.01437:
 ShapeFit: Exact location recovery from corrupted pairwise directions, 2015
 by Paul Hand, Choongbum Lee and Vladislav Voroninski.

 The problem is the following: there are n points t_1, ..., t_n in R^d
 which need to be estimated (localized). To this end, we are given
 measurements of some of the pairwise directions,
 v_ij = (t_i - t_j) / norm(t_i - t_j) + noise.
 Assume there are m such pairwise measurements, defining a graph with m
 edges over n nodes. J is the signed incidence matrix of this graph (see
 in code). To build J from lists I, J in R^m of nodes, use:
 J = sparse([I ; J], [(1:m)' ; (1:m)'], [ones(m, 1), -ones(m, 1)], n, m, 2*m);

 The measurements are arranged in the matrix V of size d x m. From V, we
 attempt to estimate t_1, ..., t_n, arranged in T, a matrix of size d x n.
 The estimation can only be done up to translation and scaling. The
 returned T's are centered: the columns sum to zero.

 ShapeFit is a formulation of this estimation problem which is robust to
 outliers. It is a nonsmooth, convex optimization problem over an affine
 space, i.e., a linear manifold. We smooth the cost using the pseudo-Huber
 loss cost and solve the problem using Manopt. This requires two
 ingredients: (1) a factory to describe the affine space, see
 shapefitfactory; (2) defining the cost and its derivative, and minimizing
 it while progressively tightening the smooth approximation (with
 warm-start).

 Simply run the example to see the results on random data. It compares the
 smoothed ShapeFit formulation against a least-squares formulation, when
 the measurements include outliers. See in code to vary the noise
 parameters, dimension d, number of nodes n, number of measurements m, ...

 Note: since the problem is convex, this returns the global optimum.
 This example also illustrates the use of Manopt for optimization under
 linear constraints: admittedly a simple subcase of optimization on
 manifolds.


 See also: shapefitfactory

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [T_hub, T_lsq, T_cvx] = shapefit_smoothed(V, J)
0002 % ShapeFit formulation for sensor network localization from pair directions
0003 %
0004 % function [T_hub, T_lsq, T_cvx] = shapefit_smoothed(V, J)
0005 %
0006 % This example in based on the paper http://arxiv.org/abs/1506.01437:
0007 % ShapeFit: Exact location recovery from corrupted pairwise directions, 2015
0008 % by Paul Hand, Choongbum Lee and Vladislav Voroninski.
0009 %
0010 % The problem is the following: there are n points t_1, ..., t_n in R^d
0011 % which need to be estimated (localized). To this end, we are given
0012 % measurements of some of the pairwise directions,
0013 % v_ij = (t_i - t_j) / norm(t_i - t_j) + noise.
0014 % Assume there are m such pairwise measurements, defining a graph with m
0015 % edges over n nodes. J is the signed incidence matrix of this graph (see
0016 % in code). To build J from lists I, J in R^m of nodes, use:
0017 % J = sparse([I ; J], [(1:m)' ; (1:m)'], [ones(m, 1), -ones(m, 1)], n, m, 2*m);
0018 %
0019 % The measurements are arranged in the matrix V of size d x m. From V, we
0020 % attempt to estimate t_1, ..., t_n, arranged in T, a matrix of size d x n.
0021 % The estimation can only be done up to translation and scaling. The
0022 % returned T's are centered: the columns sum to zero.
0023 %
0024 % ShapeFit is a formulation of this estimation problem which is robust to
0025 % outliers. It is a nonsmooth, convex optimization problem over an affine
0026 % space, i.e., a linear manifold. We smooth the cost using the pseudo-Huber
0027 % loss cost and solve the problem using Manopt. This requires two
0028 % ingredients: (1) a factory to describe the affine space, see
0029 % shapefitfactory; (2) defining the cost and its derivative, and minimizing
0030 % it while progressively tightening the smooth approximation (with
0031 % warm-start).
0032 %
0033 % Simply run the example to see the results on random data. It compares the
0034 % smoothed ShapeFit formulation against a least-squares formulation, when
0035 % the measurements include outliers. See in code to vary the noise
0036 % parameters, dimension d, number of nodes n, number of measurements m, ...
0037 %
0038 % Note: since the problem is convex, this returns the global optimum.
0039 % This example also illustrates the use of Manopt for optimization under
0040 % linear constraints: admittedly a simple subcase of optimization on
0041 % manifolds.
0042 %
0043 %
0044 % See also: shapefitfactory
0045 
0046 % This file is part of Manopt: www.manopt.org.
0047 % Original author: Nicolas Boumal, June 18, 2015.
0048 % Contributors:
0049 % Change log:
0050 
0051 
0052     % Generic useful functions
0053     center_cols = @(A) bsxfun(@minus, A, mean(A, 2));
0054     normalize_cols = @(A) bsxfun(@times, A, 1./sqrt(sum(A.^2, 1)));
0055     sqnorm_cols = @(A) sum(A.^2, 1);
0056 
0057     
0058     % DATA GENERATION
0059     %
0060     % If no inputs are specified, generate some random data for
0061     % illustration purposes.
0062     if nargin == 0
0063 
0064         % We estimate n points in R^d
0065         d =   2;
0066         n = 500;
0067 
0068         % Those points are the columns of T : they are what we need to
0069         % estimate, up to scaling and translation. We center T for
0070         % convenience.
0071         T_tru = center_cols(rand(d, n));
0072 
0073         % We get a measurement of some pairs of relative directions.
0074         % Which pairs is encoded in this graph, with J being the (signed,
0075         % transposed) incidence matrix. J is n x m, sparse.
0076         % There are roughly edge_fraction * n * (n-1) / 2 measurements.
0077         edge_fraction = 0.10;
0078         % [ii, jj] = erdosrenyi(n, edge_fraction);
0079         [ii, jj] = randomgraph(n, edge_fraction*nchoosek(n, 2));
0080         m = length(ii);
0081         J = sparse([ii ; jj], [(1:m)' ; (1:m)'], ...
0082                    [ones(m, 1), -ones(m, 1)], n, m, 2*m);
0083 
0084         % The measurements give the directions from one point to another.
0085         % That is: we get the position difference, normalized. Here, with
0086         % Gaussian noise. Least-squares will be well-suited for this.
0087         sigma = .0;
0088         V = normalize_cols(T_tru*J + sigma*randn(d, m)); % d x m
0089 
0090         % Outliers: we replace some of the direction measurements by
0091         % uniformly random unit-norm vectors.
0092         outlier_fraction = .3;
0093         outliers = rand(1, m) < outlier_fraction;
0094         V(:, outliers) = normalize_cols(randn(d, sum(outliers)));
0095         
0096     end % done generating random data
0097     
0098     
0099     
0100     
0101     
0102     [d, m] = size(V);
0103     n = size(J, 1);
0104     assert(size(J, 2) == m, 'J must be n x m, with V of size d x m.');
0105 
0106     VJt = full(V*J');
0107 
0108     % This "manifold" describes the Euclidean space of matrices T of size
0109     % d x n such that <VJt, T> = 1 and T has centered columns: T1 = 0.
0110     problem.M = shapefitfactory(VJt);
0111 
0112     % This linear operator computes the orthogonal projection of each
0113     % difference ti - tj on the orthogonal space to v_ij.
0114     % If the alignment is compatible with the data, then this is zero.
0115     % A(T) is a d x m matrix.
0116     function AT = A(T)
0117         TJ = T*J;
0118         AT = TJ - bsxfun(@times, V, sum(V .* TJ, 1));
0119     end
0120 
0121     % Need the adjoint of A, too. Input is d x m, output is d x n.
0122     Astar = @(W) (W - bsxfun(@times, V, sum(V.*W, 1)))*J';
0123 
0124     
0125     
0126     % LEAST-SQUARES
0127     %
0128     % First, work with a least-squares formulation of the problem.
0129     % That is, we minimize a (very nice) convex cost over an affine space.
0130     % Since the smooth solvers in Manopt converge to critical points, this
0131     % means they converge to global optimizers.
0132     problem.cost  = @(T) 0.5*norm(A(T), 'fro')^2;
0133     problem.egrad = @(T) Astar(A(T));
0134     problem.ehess = @(T, Tdot) Astar(A(Tdot));
0135 
0136     T_lsq = trustregions(problem);
0137     
0138 
0139     
0140     % PSEUDO-HUBER SMOOTHED SHAPEFIT
0141     %
0142     % Now solve the same, but with a pseudo-Huber loss instead of
0143     % least-squares.
0144     % We iteratively sharpen the Huber function, i.e., reduce delta.
0145     % It is important to warm start in such a fashion: trying to optimize
0146     % with a random initial guess and a very small delta is typically slow.
0147     % How fast one should decrease delta, and how accurately one should
0148     % optimize at each intermediate stage, is open for research.
0149     delta = 1;
0150     T_hub = []; % We could use T_lsq as initial guess, too.
0151     problem = rmfield(problem, 'ehess');
0152     warning('off', 'manopt:getHessian:approx');
0153     for iter = 1 : 12
0154         
0155         delta = delta / 2;
0156         
0157         h = @(x2) sqrt(x2 + delta^2) - delta; % pseudo-Huber loss
0158 
0159         problem.cost  = @(T) sum(h(sqnorm_cols(A(T))));
0160         problem.egrad = @(T) Astar(bsxfun(@times, A(T), ...
0161                                     1./sqrt(sqnorm_cols(A(T)) + delta^2)));
0162 
0163         % Solve, using the previous solution as initial guess.
0164         T_hub = trustregions(problem, T_hub);
0165         
0166     end
0167     
0168     
0169     
0170     % CVX SHAPEFIT
0171     %
0172     % Actual ShapeFit cost (nonsmooth), with CVX.
0173     % You can get CVX from http://cvxr.com/.
0174     use_cvx_if_available = false;
0175     if use_cvx_if_available && exist('cvx_version', 'file')
0176         T_cvx = shapefit_cvx(V, J);
0177     else
0178         T_cvx = NaN(d, n);
0179     end
0180     
0181     
0182     
0183     % VISUALIZATION
0184     %
0185     % If T_true is available, for display, we scale the estimators to match
0186     % the norm of the target. The scaling factor is obtained by minimizing
0187     % the norm of the discrepancy : norm(T_tru - scale*T_xxx, 'fro').
0188     % A plot is produced if d is 2 or 3.
0189     if exist('T_tru', 'var') && (d == 2 || d == 3)
0190         
0191         T_lsq = T_lsq * trace(T_tru'*T_lsq) / norm(T_lsq, 'fro')^2;
0192         T_hub = T_hub * trace(T_tru'*T_hub) / norm(T_hub, 'fro')^2;
0193         T_cvx = T_cvx * trace(T_tru'*T_cvx) / norm(T_cvx, 'fro')^2;
0194 
0195     
0196         switch d
0197             case 2
0198                 plot(T_tru(1, :), T_tru(2, :), 'bo', ...
0199                      T_lsq(1, :), T_lsq(2, :), 'rx', ...
0200                      T_hub(1, :), T_hub(2, :), 'k.', ...
0201                      T_cvx(1, :), T_cvx(2, :), 'g.');
0202             case 3
0203                 plot3(T_tru(1, :), T_tru(2, :), T_tru(3, :), 'bo', ...
0204                       T_lsq(1, :), T_lsq(2, :), T_lsq(3, :), 'rx', ...
0205                       T_hub(1, :), T_hub(2, :), T_hub(3, :), 'k.', ...
0206                       T_cvx(1, :), T_cvx(2, :), T_cvx(3, :), 'g.');
0207         end
0208 
0209         legend('ground truth', 'least squares', ...
0210                sprintf('pseudo-huber, \\delta = %.1e', delta), ...
0211                'CVX ShapeFit');
0212            
0213         title(sprintf(['ShapeFit problem : d = %d, n = %d, edge ' ...
0214                        'fraction = %.2g, sigma = %.2g, outlier ' ...
0215                        'fraction = %.2g'], d, n, edge_fraction, sigma, ...
0216                        outlier_fraction));
0217         axis equal;
0218     
0219     end
0220 
0221 end
0222 
0223 
0224 % If CVX is available, it can be used to solve the nonsmooth problem
0225 % directly, very elegantly.
0226 function T_cvx = shapefit_cvx(V, J)
0227     d = size(V, 1);
0228     n = size(J, 1); %#ok<NASGU>
0229     VJt = full(V*J');
0230     cvx_begin
0231         variable T_cvx(d, n)
0232         % We want to minimize this:
0233         % minimize sum( norms( A(T_cvx), 2, 1 ) )
0234         % But unfortunately, CVX doesn't handle bsxfun. Instead, we use
0235         % repmat, which is slower, and hence hurts the comparison in
0236         % disfavor of CVX.
0237         minimize sum( norms( T_cvx*J - V .* repmat(sum(V .* (T_cvx*J), 1), [d, 1])  , 2, 1 ) )
0238         sum(T_cvx, 2) == zeros(d, 1); %#ok<NODEF,EQEFF>
0239         VJt(:).' * T_cvx(:) == 1; %#ok<EQEFF>
0240     cvx_end
0241 end
0242 
0243 
0244 function [I, J, A] = erdosrenyi(n, p) %#ok<DEFNU>
0245 % Generate a random Erdos-Renyi graph with n nodes and edge probability p.
0246 %
0247 % [I, J, A] = erdosrenyi(n, p)
0248 %
0249 % Returns a list of edges (I(k), J(k)) for a random, undirected Erdos-Renyi
0250 % graph with n nodes and edge probability p. A is the adjacency matrix.
0251 %
0252 % I(k) < J(k) for all k, i.e., all(I<J) is true.
0253 %
0254 % The memory requirements for this simple implementation scale as O(n^2).
0255 
0256     X = rand(n);
0257     mask = X <= p;
0258     X( mask) = 1;
0259     X(~mask) = 0;
0260     X = triu(X, 1);
0261 
0262     % A is the adjacency matrix
0263     A = X + X';
0264     
0265     [I, J] = find(X);
0266 
0267 end
0268 
0269 
0270 function [I, J, A] = randomgraph(n, m)
0271 % Generates a random graph over n nodes with at most m edges.
0272 %
0273 % function [I, J, A] = randomgraph(n, m)
0274 %
0275 % Selects m (undirected) edges from a graph over n nodes, uniformly at
0276 % random, with replacement. The self edges and repeated edges are then
0277 % discarded. The remaining number of edges is at most m, and should be
0278 % close to m if m is much smaller than nchoosek(n, 2).
0279 %
0280 % The output satisfies all(I < J). A is the corresponding adjacency matrix.
0281 %
0282 % Uses O(m) memory (not O(n^2)), making it fit for large, sparse graphs.
0283 
0284     % Generate m edges at random, with replacement, and remove repetitions.
0285     IJ = unique(sort(randi(n, m, 2), 2), 'rows');
0286     
0287     % Remove self-edges if any.
0288     IJ = IJ(IJ(:, 1) ~= IJ(:, 2), :);
0289     
0290     % Actual number of selected edges
0291     k = size(IJ, 1);
0292     
0293     I = IJ(:, 1);
0294     J = IJ(:, 2);
0295     
0296     % Build the adjacency matrix of the graph.
0297     A = sparse([I ; J], [J ; I], ones(2*k, 1), n, n, 2*k);
0298 
0299 end

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