Contents
% Find point in intersection of two polyhedra, given by % { x | A1 x <= b1 } and { x | A2 x <= b2 }.
Generate problem data
randn('state', 0); rand('state', 0); n = 5; % dimension of variable m1 = 10; % number of faces for polyhedra 1 m2 = 12; % number of faces for polyhedra 2 c1 = 10*randn(n,1); % center of polyhedra 1 c2 = -10*randn(n,1); % center of polyhedra 2 % consider the following picture: % % a1 % c ---------> x % % from the center "c", we travel along vector "a1" (not necessarily a unit % vector) until we reach x. at "x", a1'x = b. a point y is to the left of x % if a1'y <= b. % % pick m1 random directions with different magnitudes A1 = diag(1 + rand(m1,1))*randn(m1,n); % the value of b is found by traveling from the center along the normal % vectors in A1 and taking its inner product with A1. b1 = diag(A1*(c1*ones(1,m1) + A1')); % pick m2 random directions with different magnitudes A2 = diag(1 + rand(m2,1))*randn(m2,n); % the value of b is found by traveling from the center along the normal % vectors in A1 and taking its inner product with A1. b2 = diag(A2*(c2*ones(1,m2) + A2')); % find the distance between the two polyhedra--make sure they overlap by % checking if the distance is 0 cvx_begin quiet variables x(n) y(n) minimize sum_square(x - y) subject to A1*x <= b1 A2*y <= b2 cvx_end % if the distance is not 0, expand A1 and A2 by a little more than half the % distance if norm(x-y) > 1e-4, A1 = (1 + 0.5*norm(x-y))*A1; A2 = (1 + 0.5*norm(x-y))*A2; % recompute b's as appropriate b1 = diag(A1*(c1*ones(1,m1) + A1')); b2 = diag(A2*(c2*ones(1,m2) + A2')); end
Solve problem
[x history] = polyhedra_intersection(A1, b1, A2, b2, 1.0, 1.0);
iter r norm eps pri s norm eps dual objective 1 6.8713 0.1385 13.8243 0.0689 0.00 2 4.6972 0.1809 3.3124 0.0374 0.00 3 0.9614 0.1692 1.6155 0.0354 0.00 4 0.4124 0.1753 1.3753 0.0313 0.00 5 0.8436 0.1839 0.9909 0.0229 0.00 6 1.0713 0.1884 0.5223 0.0122 0.00 7 1.0909 0.1888 0.0446 0.0012 0.00 8 0.1018 0.1888 0.0001 0.0002 0.00 Elapsed time is 2.392706 seconds.
Reporting
K = length(history.objval); h = figure; plot(1:K, history.objval, 'k', 'MarkerSize', 10, 'LineWidth', 2); ylabel('f(x^k) + g(z^k)'); xlabel('iter (k)'); g = figure; subplot(2,1,1); semilogy(1:K, max(1e-8, history.r_norm), 'k', ... 1:K, history.eps_pri, 'k--', 'LineWidth', 2); ylabel('||r||_2'); subplot(2,1,2); semilogy(1:K, max(1e-8, history.s_norm), 'k', ... 1:K, history.eps_dual, 'k--', 'LineWidth', 2); ylabel('||s||_2'); xlabel('iter (k)');
Compare to alternating projections
% MAX_ITER = 10; % x = zeros(n,1); % z = zeros(n,1); % for k = 1:MAX_ITER % % % x-update % % use cvx to find point in first polyhedra % cvx_begin quiet % variable x(n) % minimize (sum_square(x - z)) % subject to % A1*x <= b1 % cvx_end % % % z-update with relaxation % zold = z; % % use cvx to find point in second polyhedra % cvx_begin quiet % variable z(n) % minimize (sum_square(x - z)) % subject to % A2*z <= b2 % cvx_end % % history1.r_norm(k) = norm(x - z); % history1.s_norm(k) = norm((z - zold)); % % end % % g = figure % subplot(2,1,1); % semilogy(1:MAX_ITER, max(1e-8, history1.r_norm), 'k', 1:K, max(1e-8, history.r_norm), 'r'); % ylabel('||r||_2'); % % subplot(2,1,2); % semilogy(1:MAX_ITER, max(1e-8, history1.s_norm), 'k', 1:K, max(1e-8, history.s_norm), 'r'); % ylabel('||s||_2'); xlabel('iter (k)');