smallestinconvexhull.m

function [u_norm, coeffs, u] = smallestinconvexhull(M, x, U, tol)
% Computes a minimal norm convex combination of given tangent vectors in Manopt.
%
% function [u_norm, coeffs, u] = smallestinconvexhull(M, x, U)
% function [u_norm, coeffs, u] = smallestinconvexhull(M, x, U, tol)
%
% M is a manifold as returned by a Manopt factory.
% x is a point on this manifold.
% U is a cell containing N tangent vectors U{1} to U{N} at x.
% tol (default: 1e-8): tolerance for solving the quadratic program.
% 
% This function computes u, a tangent vector at x contained in the convex
% hull spanned by the N vectors U{i}, with minimal norm (according to the
% Riemannian metric on M). This is obtained by solving a convex quadratic
% program involving the Gram matrix of the given tangent vectors.
% The quadratic program is solved using Matlab's built-in quadprog,
% which requires the optimization toolbox. If this toolbox is not
% available, consider replacing with CVX for example.
%
%
% u_norm is the norm of the smallest vector u.
% coeffs is a vector of length N with entries in [0, 1] summing to 1.
% u is the sought vector: u = coeffs(1)*U{1} + ... + coeffs(N)*U{N}.
%
% Nicolas Boumal, Feb. 19, 2013
% Modified April 6, 2016 to work with Manopt.

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, June 28, 2016.
% Contributors: 
% Change log: 
%
%   June 28, 2016 (NB):
%       Adapted for Manopt from original code by same author (Feb. 19, 2013)

% Example code: pick a manifold, a point, and a collection of tangent
% vectors at that point, then get the smallest vector in the convex hull
% of those:
% 
% M = spherefactory(5);
% x = M.rand();
% N = 3;
% U = cell(N,1);
% for k = 1 : N, U{k} = M.randvec(x); end
% [u_norm, coeffs, u] = smallestinconvexhull(M, x, U)

    % We simply need to solve the following quadratic program:
    % minimize ||u||^2 such that u = sum_i s_i U_i, 0 <= s_i <= 1
    %                            and sum_i s_i = 1
    %
    % This is equivalent to solving:
    %  min s'*G*s s.t. 0 <= s <= 1, s'*ones = 1, with G(i, j) = <U_i, U_j> (Gram matrix)
    % Then our solution is s_1 U_1 + ... + s_N U_N.
    
    
    % Compute the Gram matrix of the given tangent vectors
    N = numel(U);
    G = grammatrix(M, x, U);
    
    % Solve the quadratic program.
    % If the optimization toolbox is not available, consider replacing with
    % CVX.
    
    if ~exist('tol', 'var') || isempty(tol)
        tol = 1e-8;
    end
    
    opts = optimset('Display', 'off', 'TolFun', tol);
    [s_opt, cost_opt] ...
          = quadprog(G, zeros(N, 1),     ...  % objective (squared norm)
                     [], [],             ...  % inequalities (none)
                     ones(1, N), 1,      ...  % equality (sum to 1)
                     zeros(N, 1),        ...  % lower bounds (s_i >= 0)
                     ones(N, 1),         ...  % upper bounds (s_i <= 1)
                     [],                 ...  % we do not specify an initial guess
                     opts);

    % Norm of the smallest tangent vector in the convex hull:
    u_norm = real(sqrt(2*cost_opt));

    % Keep track of optimal coefficients
    coeffs = s_opt;
    
    % If required, construct the vector explicitly.
    if nargout >= 3
        u = lincomb(M, x, U, coeffs);
    end
                 
end