function [x, ...
SquaredNormOfResidualVector, ...
ResidualVectorbMinusMx, ...
exitflag] = ...
BoundedVariablesLeastSquares0to1GJK(M,b) ;
% Purpose Use the GJK algorithm to solve the bounded variables least squares problem
% Mx=b, where each entry of x is greater than or equal to 0 and less than or equal to 1.
%
% Description The bounded variables least squares (BVLS) problem is to find a vector x that
% comes closest to satisfying Mx=b, where M is a known matrix of m rows and
% n columns, and b is a known column vector with n entries. BVLS constrains
% x so that each entry of x has an upper and lower bound. By a suitable scaling
% and translation, one can always find an equivalent BVLS problem such that
% each entry of x is greater than or equal to 0 and less than or equal to 1.
%
% The algorithm behind this routine is described in [Centore2015]. Rather than
% more common approaches involving quadratic programming (QP), the routine
% implements the Gilbert-Johnson-Keerthi (GJK) algorithm ([Gilbert1988],
% [Rabbitz1994]), which finds the minimal
% distance between two convex polytopes. The routine uses a special case of
% the algorithm, in which the minimal distance is found between a point and a
% convex polytope. The approach in [Centore2015] considers Mx as a subset of
% the vector space R^m. That subset is determined by both the columns of M
% and the constraints on x. With the given constraints, Mx is the zonotope
% generated by the columns of M, where each column is treated as a vector in
% R^m. The target point b is also in R^m. Since a zonotope is a convex
% polytope, the GJK algorithm can find the minimum distance between b and the
% zonotope, GJK also finds the point b_Z on the zonotope that is closest to b,
% and an expression for b_Z as a convex combination of the columns of M. This
% information is sufficient to work backwards to a minimizing vector x. See
% [Centore2015] for further details.
%
% The GJK approach is proposed here as a BVLS algorithm that has the potential
% to improve on current algorithms, but whose effectiveness has not been
% tested. The GJK approach has the advantage of being easy to understand and
% implement, but has the disadvantage of using 2^n points as a generating set
% for the zonotope. The GJK algorithm speed depends approximately linearly on
% the size of the generating set, so the GJK approach could be exponential in
% the number of entries in the vector of x. On the other hand, there is no
% need to solve the zonotope (i.e. determine its vertices), not to express it
% in terms of half-space inequalities. If there are not too many variables, then
% the GJK approach might be preferable. In any event, the algorithm is presented
% as open source so that other researchers may make tests and comparisons as they
% wish, and also add any desired improvements.
%
% References [Centore2015] Paul Centore, A GJK-Based Algorithm for Constrained Least Squares,
% available at wwww.99main.com/~centore, 2015
% [Gilbert1988] Elmer G. Gilbert, Daniel W. Johnson, & Sathiya Keerthi, A Fast
% Procedure for Computing the Distance Between Complex Objects in
% Three-Dimensional Space, IEEE Journal of Robotics and Automation,
% Vol. 4, No. 2, April 1988, pp. 193-203.
% [Rabbitz1994] Rich Rabbitz, Fast Collision Detection of Moving Convex Polyhedra,
% in Section I.8 of Graphics Gems IV (IBM Version), ed. Paul Heckbert,
% Academic Press, 1994.
%
% M,b Entries in the BVLS problem Mx=b. M is a matrix with
% m rows and n columns, and b is a column vector with m entries.
%
% x A column vector of n entries, that minimizes the norm of Mx-b.
%
% ResidualVectorbMinusMx The residual vector, given by b-Mx, where x is the
% vector which minimizes the norm of Mx-b.
%
% SquaredNormOfResidualVector The squared norm of the residual vector b-Mx.
%
% ErrorF Each vertex is a linear combination of the generating vectors.
% The ith row of the matrix VertexCoefficients gives the
% coefficients of the generating vectors for the ith vertex, which
% is given in the ith row of the matrix Vertices.
%
% Author Paul Centore (January 25, 2015)
%
% Copyright 2014 Paul Centore
%
% This file is part of MunsellAndKubelkaMunkToolbox.
%
% MunsellAndKubelkaMunkToolbox is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% MunsellAndKubelkaMunkToolbox is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with MunsellAndKubelkaMunkToolbox. If not, see .
% Initialize variables
x = [] ;
SquaredNormOfResidualVector = NaN ;
ResidualVectorbMinusMx = [] ;
exitflag = 1 ;
% Extract information about the number of rows and columns of M
[m,n] = size(M) ;
% The index matrix IndexJ expresses a generating set for the zonotope P (a generating set
% for P is a set V of vectors whose convex hull is P) as a set of linear combinations of
% the columns of M. A node of a zonotope is a linear combination of the vectors that
% generate the zonotope, where every coefficient in the combination is either 0 or 1. (The
% generators of the zonotope should not be confused with a generating set for the zonotope.
% Every point in the zonotope is a convex combination of vectors in the generating set, but
% many points are not convex combinations of the generators.) The index matrix IndexJ will
% therefore consist of all 2^n possible sequences of 0's and 1's. Each row gives one such
% sequence.
[IndexJ,ind] = combn([0 1],n) ;
% The generators V of P are linear combinations of the columns of M, with coefficients
% given by the index matrix J. Calculate these combinations. The result will be a
% matrix PolytopeGeneratorsV. Each row of PolytopeGeneratorsV is one of the generators
% of P.
NumberOfGenerators = rows(IndexJ) ;
PolytopeGeneratorsV = NaN * ones(NumberOfGenerators,m) ;
for ctr = 1:NumberOfGenerators
% Calculate a linear combination of columns of M, term by term
TempGenerator = zeros(1,m) ;
for ColumnCtr = 1:n
TempGenerator = TempGenerator + ...
IndexJ(ctr,ColumnCtr) .* transpose(M(:,ColumnCtr)) ;
end
PolytopeGeneratorsV(ctr,:) = TempGenerator ;
end
% The constrained least squares solution to Mx=b corresponds to the point on the
% zonotope P that is nearest to the origin. Call the GJK algorithm to calculate
% this point. The return variable alphaCoefficients gives the coefficients in a convex
% combination of the generators of P; this convex combination is the closest point b_Z
% on P to b.
[ b_Z, ...
Distance, ...
PointInPolytope, ...
alphaCoefficients, ...
ErrorFlag] = ...
ClosestPointInConvexPolytopeGJK( ...
PolytopeGeneratorsV, ...
reshape(b,1,length(b))) ;
% If the GJK algorithm did not produce a solution, then relay this failure by setting the
% exit flag of the current routine to 0, and return
if ErrorFlag
exitflag = 0 ;
return
end
% [Centore2015] shows that the entries of x given by dot products of the coefficients alpha
% and the columns of the index matrix J give a BVLS solution that satisfies the constraints.
x = NaN * ones(n,1) ;
for ctr = 1:n
x(ctr,1) = dot(alphaCoefficients, IndexJ(:,ctr)) ;
end
% The residual vector is the vector between the target point b and the point b_Z on the
% zonotope which is closest to b. To be consistent with the returned solution, b_Z is
% calculated as Mx rather than using the GJK output. Following Octave/MATLAB conventions
% for lsqnonneg, this routine also returns the squared norm of the residual vector
ResidualVectorbMinusMx = reshape(b,length(b),1) - (M * reshape(x,length(x),1)) ;
SquaredNormOfResidualVector = (norm(ResidualVectorbMinusMx))^2 ;