function [x, ...
SquaredNormOfResidualVector, ...
ResidualVectorbMinusMx, ...
exitflag] = ...
lsqnonnegGJK(M,b) ;
% Purpose Solve the non-negative least squares problem Mx=b, where each entry of x is
% zero or positive. Use the GJK-based algorithm introduced in [Centore2015].
%
% Description The non-negative least squares (NNLS) problem is to find a vector x, such
% that the norm of Mx-b is minimized, where M is a known matrix, and b is a
% known vector; there is a requirement, however, that each entry of x be
% zero or positive. If the minimum norm is 0, then x satisfies Mx=b exactly,
% if not, then x can be thought of as a vector which comes as close as possible
% to satisfying Mx=b.
%
% This routine uses an algorithm introduced in [Centore2015]. That algorithm
% is based on the Gilbert-Johnson-Keerthi (GJK) algorithm that was introduced
% in 1988 ([Gilbert1988]). Readable descriptions appear in [Rabbitz1994] (although
% the term GJK algorithm is not used there), [Ericson2004], and [McCutchan2006].
% In its general form, the algorithm finds the minimum distance between two
% convex polytopes. We actually need only a special case in which one polytope
% is a point.
%
% The algorithm in [Centore2015] considers the term Mx in the NNLS equation as
% a set in a vector space. Since each x is non-negative, Mx is the convex cone
% generated by the columns of M. b can be thought of
% as a target vector. The least squares problem finds the point on the convex
% cone that is closest to b.
%
% Rather than dealing with the cone as a whole, it is truncated at a certain
% distance Gamma. The truncated cone is the convex polytope P_Gamma generated by
% (adjusted) column vectors of M, and the zero vector. The GJK algorithm
% can then be used to find the nearest point on P_Gamma to b, and the shortest
% distance. As Gamma increases, the distance (as a function of Gamma) will
% reach a plateau, and the nearest point will be identical for any larger
% Gammas. When this happens, that point will also be the minimizing point
% for the entire convex cone. This routine starts off with a small Gamma,
% and increases it by a factor of 10 (up to some maximum), until the convergence
% is noted, at which time it is concluded that a solution has been found.
%
% This novel algorithm is presented here because it seems to have the
% potential to improve on current algorithms, both because the GJK algorithm is
% very fast, and because the calculation of the necessary convex polytopes is
% almost immediate from the matrix M. No speed tests or formal comparisons
% have been made with other algorithms, however, so any improvements are, for
% now at least, just promising hypotheses. The new algorithm does have the
% advantage of being easy to understand and interpret. Some issues with
% numerical stability have also been noted, which might cause problems for any
% algorithm. Attempts were made to mitigate instability, but doubtless more
% sophisticated approaches are available. Other researchers with
% more experience in NNLS problems are encourage to use, test, and modify this
% open source code, in the process of which the algorithm s usefulness can be
% assessed.
%
% 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.
% [Ericson2004] Christer Ericson, The Gilbert-Johnson-Keerthi Algorithm,
% http://realtimecollisiondetection.net/pubs/SIGGRAPH04_Ericson_GJK_notes.pdf
% [McCutchan2006] John McCutchan, Introduction to GJK, 9 November 2006,
% http://www.cas.mcmaster.ca/~carette/SE3GB3/2006/notes/gjk1_pres.pdf
%
% Variables M,b Entries in the NNLS 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.
%
% SquaredNormOfResidualVector The squared norm of the residual vector b-Mx.
%
% ResidualVectorbMinusMx The residual vector, given by b-Mx, where x is the
% non-negative vector which minimizes the norm of Mx-b.
%
% exitflag A performance flag that is 0 when the algorithm fails to
% produce a result, and a positive integer otherwise.
%
% Author Paul Centore (January 8, 2015)
%
% Copyright 2015 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 = 0 ;
% Minimum distance for some comparisons, to avoid numerical roundoff error
EpsDist = 1e-12 ;
[Rows, Cols] = size(M) ;
% Loop through a set of truncated cones, whose size increases by a factor of 10 at each
% iteration. Apply the GJK algorithm to find shortest distance from the truncated cone to
% the target vector b. When the distance stops decreasing, the constant distance is the
% distance to the convex cone Mx, and the GJK algorithm will have found the closest point
% on Mx to the target b; this is (or leads to) an NNLS solution.
SolutionFound = false ;
PowerOfTen = 0 ;
MaxPowerOfTen = 15 ;
% The following matrices keep track of information over each iteration. The ith row of
% each matrix corresponds to the ith truncated cone.
ClosestPoints = [] ;
Distances = [] ;
PointInPolytopes = [] ;
CoefficientsList = [] ;
NormsOfMxMinusb = [] ;
while not(SolutionFound) && PowerOfTen <= MaxPowerOfTen
% A truncated cone is the convex hull of the columns of M and the zero vector; the
% truncated cone is scaled by a factor which increases tenfold at each iteration.
PolytopeGenerators = [transpose(M); zeros(1,Rows)] ;
PolytopeGenerators = (10^PowerOfTen) * PolytopeGenerators ;
% Apply the GJK algorithm to a truncated cone, to find the nearest point to b on
% that truncated cone, and related information
[ClosestPoint, ...
Distance, ...
PointInPolytope,...
Coefficients] = ...
ClosestPointInConvexPolytopeGJK(...
PolytopeGenerators,...
reshape(b,1,length(b)) );
% The variable Coefficients give the coefficients in a convex combination of the
% polytope generators (the zero vector and the scaled columns of M) that produces the
% nearest point to b on the truncated cone. They must be rescaled to give similar
% coefficients in terms of the original columns of M,
Coefficients = Coefficients*(10^PowerOfTen) ;
% Record the information for the new truncated cone as a new row in the information matrices
MxMinusb = M *reshape(Coefficients(1:(end-1)),length(Coefficients)-1,1) - reshape(b,length(b),1) ;
NormsOfMxMinusb = [NormsOfMxMinusb; norm(MxMinusb)] ;
ClosestPoints = [ClosestPoints; ClosestPoint] ;
Distances = [Distances; Distance ] ;
PointInPolytopes = [PointInPolytopes; PointInPolytope] ;
CoefficientsList = [CoefficientsList; reshape(Coefficients,1,length(Coefficients))] ;
% Check whether the point b is inside or on a truncated cone, in which case it is also
% in the convex cone Mx. If so, then a solution has been found.
if PointInPolytope
SolutionFound = true ;
% If not, then see whether the distance to the truncated cone has stopped decreasing.
% If it has, then there will be no improvement as the truncated cones increase in
% size, and we have found the nearest point
elseif length(Distances) >= 2
if abs(Distances(end) - Distances(end-1)) < EpsDist
SolutionFound = true ;
end
end
% Make next truncated cone 10 times the size of the current one
PowerOfTen = PowerOfTen + 1 ;
end
% If a solution has been found, then calculate the returned values, such as the residual
% vector and its norm.
if SolutionFound
% The next sorting step was introduced because some numerical instabilities were
% seen. Since each truncated cone properly contains the previous truncated cone, the
% minimum distance to the point b should constantly be decreasing. The fact that it
% is not indicates numerical problems. To try to mitigate instability when it occurs,
% check directly for the cone with the smallest distance, and use information from
% that cone.
[SortedNorms, SortedIndices] = sort(NormsOfMxMinusb) ;
x = CoefficientsList(SortedIndices(1),1:(end-1)) ;
ResidualVectorbMinusMx = reshape(b,length(b),1) - (M * reshape(x,length(x),1)) ;
SquaredNormOfResidualVector = (Distances(SortedIndices(1)))^2 ;
exitflag = 1 ;
end