Mercurial > hg > camir-aes2014

annotate toolboxes/distance_learning/mlr/util/lcprog.m @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 function [z, w, status, z0] = lcprog(M, q, varargin)
wolffd@0 2 %LCPROG - Linear complementary programming with Lemke's algorithm.
wolffd@0 3 %
wolffd@0 4 % [z, w] = lcprog(M, q);
wolffd@0 5 % [z, w, status, z0] = lcprog(M, q, varargin);
wolffd@0 6 %
wolffd@0 7 % Inputs:
wolffd@0 8 % M - Matrix (n x n matrix)
wolffd@0 9 % q - Vector (n vector)
wolffd@0 10 %
wolffd@0 11 % Optional Inputs:
wolffd@0 12 % Optional name-value pairs are:
wolffd@0 13 % 'maxiter' Maximum number of iterations before termination of algorithm (default is (1+n)^3).
wolffd@0 14 % 'tiebreak' Method to break ties for pivot selection (default is 'first'). Options are:
wolffd@0 15 % 'first' Select pivot with lowest index.
wolffd@0 16 % 'last' Select pivot with highest index.
wolffd@0 17 % 'random' Select a pivot at random.
wolffd@0 18 % 'tolpiv' Pivot tolerance below which a number is considered to be zero (default is 1.0e-9).
wolffd@0 19 %
wolffd@0 20 % Outputs:
wolffd@0 21 % z - Solution vector (n vector)
wolffd@0 22 % w - Complementary solution vector (n vector)
wolffd@0 23 % status - Status of optimization with values (integer)
wolffd@0 24 % 6 Ray termination
wolffd@0 25 % 1 Convergence to a solution (normal termination)
wolffd@0 26 % 0 Termination prior to convergence
wolffd@0 27 % -2 Infeasible problem
wolffd@0 28 % -3 Bad pivot to an infeasible solution
wolffd@0 29 % z0 - Artificial variable which is 0 upon normal termination (scalar)
wolffd@0 30 %
wolffd@0 31 % Comments:
wolffd@0 32 % Given a matrix M and a vector q, the linear complementary programming (LCP) problem seeks an
wolffd@0 33 % orthogonal pair of non-negative vectors w and z such that
wolffd@0 34 % w = M * z + q
wolffd@0 35 % with
wolffd@0 36 % w' * z = 0
wolffd@0 37 % w >= 0
wolffd@0 38 % z >= 0 .
wolffd@0 39 %
wolffd@0 40 % The algorithm starts with the problem
wolffd@0 41 % w = M * z + q + z0 * 1
wolffd@0 42 % with an artificial variable z0 > 0 such that q + z0*1 >= 0. The algorithm goes through a series
wolffd@0 43 % of pivot steps that seeks to drive z0 to 0. The algorithm can stop in two ways with either
wolffd@0 44 % normal termination or ray termination.
wolffd@0 45 %
wolffd@0 46 % With normal termination, the artificial variable z0 is zero so that the original problem is
wolffd@0 47 % solved.
wolffd@0 48 %
wolffd@0 49 % With ray termination, the artificial variable z0 is non-negative with "ray" solution in the form
wolffd@0 50 % w = wR + lambda * dw
wolffd@0 51 % z = zR + lambda * dz
wolffd@0 52 % z0 = z0R + lambda * dz0
wolffd@0 53 % that satisfies
wolffd@0 54 % w = M * z + q + z0 * 1
wolffd@0 55 % w' * z = 0
wolffd@0 56 % for any lambda >= 0. If ray termination occurs, the returned variables w, z, and z0 are the
wolffd@0 57 % points on the ray with lambda = 0.
wolffd@0 58 %
wolffd@0 59 % If tolpiv is too small, the algorithm may identify a pivot that is effectively zero. If this
wolffd@0 60 % situation occurs, the status code -3 is triggered. Try a larger pivot tolerance 'tolpiv' although
wolffd@0 61 % a pivot that is too large can trigger false convergence due to failure to find a pivot.
wolffd@0 62 %
wolffd@0 63
wolffd@0 64 %
wolffd@0 65 % Copyright 2007 The MathWorks, Inc.
wolffd@0 66 % $Revision: 1.1.6.4 $ $ Date: $
wolffd@0 67
wolffd@0 68 % initialization
wolffd@0 69
wolffd@0 70 if nargin < 2 || isempty(M) || isempty(q)
wolffd@0 71 error('finance:lcprog:MissingInputArgument', ...
wolffd@0 72 'One or more missing required input arguments M, q.');
wolffd@0 73 end
wolffd@0 74
wolffd@0 75 q = q(:);
wolffd@0 76 n = size(q,1);
wolffd@0 77
wolffd@0 78 % process name-value pairs (if any)
wolffd@0 79
wolffd@0 80 if mod(nargin-2,2) ~= 0
wolffd@0 81 error('finance:lcprog:InvalidNameValuePair', ...
wolffd@0 82 'Invalid name-value pairs with either a missing name or value.');
wolffd@0 83 end
wolffd@0 84
wolffd@0 85 names = { 'maxiter', 'tiebreak', 'tolpiv' }; % names
wolffd@0 86 values = { (1 + n)^3, 'first', 1.0e-9 }; % default values
wolffd@0 87 try
wolffd@0 88 [maxiter, tiebreak, tolpiv] = parsepvpairs(names, values, varargin{:});
wolffd@0 89
wolffd@0 90 catch E
wolffd@0 91 E.throw
wolffd@0 92 end
wolffd@0 93
wolffd@0 94 % check optional arguments
wolffd@0 95
wolffd@0 96 if maxiter <= 0
wolffd@0 97 error('finance:lcprog:NonPositiveInteger', ...
wolffd@0 98 'Maximum number of iterations (maxiter) must be a positive integer.');
wolffd@0 99 end
wolffd@0 100
wolffd@0 101 if ~isempty(strmatch(tiebreak,{'default','first'}))
wolffd@0 102 tb = 1;
wolffd@0 103 elseif ~isempty(strmatch(tiebreak,'last'))
wolffd@0 104 tb = 2;
wolffd@0 105 elseif ~isempty(strmatch(tiebreak,'random'))
wolffd@0 106 tb = 3;
wolffd@0 107 else
wolffd@0 108 tb = 1;
wolffd@0 109 end
wolffd@0 110
wolffd@0 111 if tolpiv < 2*eps
wolffd@0 112 error('finance:lcprog:InvalidTolerance', ...
wolffd@0 113 'Unrealistically small tolerance (tolpiv) specified.');
wolffd@0 114 end
wolffd@0 115
wolffd@0 116 % trivial solution
wolffd@0 117
wolffd@0 118 if all(q >= 0)
wolffd@0 119 z = zeros(n,1);
wolffd@0 120 w = q;
wolffd@0 121 status = 1;
wolffd@0 122 z0 = 0;
wolffd@0 123 return
wolffd@0 124 end
wolffd@0 125
wolffd@0 126 % set up general problem
wolffd@0 127
wolffd@0 128 w = zeros(n,1);
wolffd@0 129 z = zeros(n,1);
wolffd@0 130 u = ones(n,1);
wolffd@0 131 z0 = max(-q);
wolffd@0 132
wolffd@0 133 tableau = [ eye(n), -M, -u, q];
wolffd@0 134 basic = (1:n)';
wolffd@0 135
wolffd@0 136 % initial pivot and initial basis
wolffd@0 137
wolffd@0 138 r = find(q == min(q),1); % r are pivot rows
wolffd@0 139 cx = r; % w(r) leaves the basis
wolffd@0 140 c = 2*n + 1; % initial pivot column
wolffd@0 141 basic(r) = c; % z0 enters the basis
wolffd@0 142
wolffd@0 143 % pivot step
wolffd@0 144
wolffd@0 145 pr = (1/tableau(r,c)) * tableau(r,:);
wolffd@0 146 pc = tableau(:,c);
wolffd@0 147 tableau = tableau - pc * pr;
wolffd@0 148 tableau(r,:) = pr;
wolffd@0 149
wolffd@0 150 % main loop
wolffd@0 151
wolffd@0 152 for iter = 2:maxiter
wolffd@0 153
wolffd@0 154 % get new basic variable
wolffd@0 155
wolffd@0 156 if cx <= n
wolffd@0 157 c = n + cx;
wolffd@0 158 elseif cx <= 2*n
wolffd@0 159 c = cx - n;
wolffd@0 160 end
wolffd@0 161
wolffd@0 162 % do ratio test
wolffd@0 163
wolffd@0 164 d = tableau(:,c);
wolffd@0 165 mind = NaN(n,1);
wolffd@0 166 for i = 1:n
wolffd@0 167 if d(i) > tolpiv
wolffd@0 168 mind(i) = tableau(i,end)/d(i);
wolffd@0 169 end
wolffd@0 170 end
wolffd@0 171 iind = find(mind == nanmin(mind));
wolffd@0 172
wolffd@0 173 % ray termination
wolffd@0 174
wolffd@0 175 if isempty(iind) % cannot find a new basic variable
wolffd@0 176 for i = 1:n
wolffd@0 177 if basic(i) <= n
wolffd@0 178 ii = basic(i);
wolffd@0 179 w(ii) = tableau(i,end);
wolffd@0 180 z(ii) = 0;
wolffd@0 181 elseif basic(i) <= 2*n
wolffd@0 182 ii = basic(i) - n;
wolffd@0 183 w(ii) = 0;
wolffd@0 184 z(ii) = tableau(i,end);
wolffd@0 185 else
wolffd@0 186 z0 = tableau(i,end);
wolffd@0 187 end
wolffd@0 188 end
wolffd@0 189 if lcprealitycheck(M, q, w, z)
wolffd@0 190 status = -3;
wolffd@0 191 z = NaN(n,1);
wolffd@0 192 w = NaN(n,1);
wolffd@0 193 z0 = NaN;
wolffd@0 194 else
wolffd@0 195 status = 6;
wolffd@0 196 end
wolffd@0 197 return
wolffd@0 198 end
wolffd@0 199
wolffd@0 200 % identify next basic variable to be replaced using ratio test
wolffd@0 201
wolffd@0 202 if tb == 3
wolffd@0 203 ii = ceil(rand()*numel(iind));
wolffd@0 204 if ii < 1
wolffd@0 205 ii = 1;
wolffd@0 206 end
wolffd@0 207 r = iind(ii);
wolffd@0 208 elseif tb == 2
wolffd@0 209 r = iind(end);
wolffd@0 210 elseif tb == 1
wolffd@0 211 r = iind(1);
wolffd@0 212 end
wolffd@0 213 for i = 1:numel(iind) % always select z0 if it is in the index set
wolffd@0 214 if basic(iind(i)) == (2*n + 1)
wolffd@0 215 r = iind(i);
wolffd@0 216 end
wolffd@0 217 end
wolffd@0 218
wolffd@0 219 % new basis
wolffd@0 220
wolffd@0 221 cx = basic(r); % get variable leaving the basis
wolffd@0 222 basic(r) = c; % move new basic variable into the basis
wolffd@0 223
wolffd@0 224 % pivot step
wolffd@0 225
wolffd@0 226 pr = (1/tableau(r,c)) * tableau(r,:);
wolffd@0 227 pc = tableau(:,c);
wolffd@0 228 tableau = tableau - pc * pr;
wolffd@0 229 tableau(r,:) = pr;
wolffd@0 230
wolffd@0 231 if ~isfinite(tableau(1,end)) % trap overflow/underflow
wolffd@0 232 warning('finance:lcprog:InfeasibleProblem', ...
wolffd@0 233 'Tableau overflow/underflow due to bad pivot.');
wolffd@0 234 status = -3;
wolffd@0 235 z = NaN(n,1);
wolffd@0 236 w = NaN(n,1);
wolffd@0 237 z0 = NaN;
wolffd@0 238 return
wolffd@0 239 end
wolffd@0 240
wolffd@0 241 % convergence test
wolffd@0 242
wolffd@0 243 if all(basic <= 2*n) % test if z0 left the basis -> convergence
wolffd@0 244 for i = 1:n
wolffd@0 245 if basic(i) <= n
wolffd@0 246 ii = basic(i);
wolffd@0 247 w(ii) = tableau(i,end);
wolffd@0 248 z(ii) = 0;
wolffd@0 249 else
wolffd@0 250 ii = basic(i) - n;
wolffd@0 251 w(ii) = 0;
wolffd@0 252 z(ii) = tableau(i,end);
wolffd@0 253 end
wolffd@0 254 end
wolffd@0 255 z0 = 0;
wolffd@0 256 if lcprealitycheck(M, q, w, z)
wolffd@0 257 status = -3;
wolffd@0 258 z = NaN(n,1);
wolffd@0 259 w = NaN(n,1);
wolffd@0 260 z0 = NaN;
wolffd@0 261 else
wolffd@0 262 status = 1;
wolffd@0 263 end
wolffd@0 264 return
wolffd@0 265 end
wolffd@0 266 end
wolffd@0 267
wolffd@0 268 % too many iterations
wolffd@0 269
wolffd@0 270 for i = 1:n
wolffd@0 271 if basic(i) <= n
wolffd@0 272 ii = basic(i);
wolffd@0 273 w(ii) = tableau(i,end);
wolffd@0 274 z(ii) = 0;
wolffd@0 275 elseif basic(i) <= 2*n
wolffd@0 276 ii = basic(i) - n;
wolffd@0 277 w(ii) = 0;
wolffd@0 278 z(ii) = tableau(i,end);
wolffd@0 279 else
wolffd@0 280 z0 = tableau(i,end);
wolffd@0 281 end
wolffd@0 282 end
wolffd@0 283 if lcprealitycheck(M, q, w, z)
wolffd@0 284 status = -3;
wolffd@0 285 z = NaN(n,1);
wolffd@0 286 w = NaN(n,1);
wolffd@0 287 z0 = NaN;
wolffd@0 288 else
wolffd@0 289 status = 0;
wolffd@0 290 end
wolffd@0 291
wolffd@0 292
wolffd@0 293 function notok = lcprealitycheck(M, q, w, z)
wolffd@0 294 %LCPREALITYCHECK - Make sure candidate solution z and w is complementary basic feasible.
wolffd@0 295
wolffd@0 296 u = w - M*z - q;
wolffd@0 297
wolffd@0 298 if (max(u) - min(u)) > 1
wolffd@0 299 warning('finance:lcprog:InfeasibleProblem', ...
wolffd@0 300 'Candidate solution is infeasible due to a bad pivot.');
wolffd@0 301 notok = true;
wolffd@0 302 else
wolffd@0 303 notok = false;
wolffd@0 304 end
wolffd@0 305
wolffd@0 306
wolffd@0 307 % [EOF]