Andreas R. Krause

• About me
• Research
• Publications
• CV

Matlab Toolbox for Submodular Function Optimization (v 0.99)

Tutorial script and all implementations Andreas Krause (krausea@gmail.com). Slides and detailed references available at http://www.submodularity.org.

Tested in MATLAB 7.0.1 (R14), 7.2.0 (R2006a), 7.4.0 (R2007a, MAC).

Download the toolbox.

Welcome to this introduction!

We will discuss how we can use the toolbox to optimize submodular set functions, i.e., functions that take a subset A of a finite ground set V to the real numbers, satisfying

Some information on conventions:

All algorithms will use function pointers @F, where F is a function that takes a set (an array) A and returns a real number (see examples below). They will also take an index set V, and A must be a subset of V.

Implemented algorithms:

1) Minimization:

2) Maximization:

3) Miscellaneous

4) Submodular functions:

Here is an overview reference:

A. Krause, C. Guestrin. Near-optimal Observation Selection Using Submodular Functions. Survey paper, Proc. of 22nd Conference on Artificial Intelligence (AAAI) 2007 -- Nectar Track

Contents

clear
randn('state',0); rand('state',0);

PART 1) MINIMIZATION OF SUBMODULAR FUNCTIONS

Queyranne's algorithm for minimizing symmetric submodular functions

We will first explore Queyranne's algorithm for minimizing a symmetric submodular function.

Example: F is a cut function on an undirected graph with adjacency matrix G_un on vertices V_G:

G_un=[0 1 1 0 0 0; 1 0 1 0 0 0; 1 1 0 1 0 0; 0 0 1 0 1 1; 0 0 0 1 0 1; 0 0 0 1 1 0]
V_G = 1:6;

G_un =

 0 1 1 0 0 0
 1 0 1 0 0 0
 1 1 0 1 0 0
 0 0 1 0 1 1
 0 0 0 1 0 1
 0 0 0 1 1 0

Here's the cut function:

F_cut_un = @(A) sfo_fn_cutfun(G_un,A);

Now we can run Queyranne's algorithm, which will return the minimum cut in graph G_un:

A = sfo_queyranne(F_cut_un,V_G)

A =

 4 5 6

Minimizing general submodular functions

We can also minimize general (not necessarily symmetric) submodular functions. This toolbox implements Fujishige's min-norm-point algorithm. We try it out on a test function by Iwata, described, e.g., in Fujishige et al. 06, on a set of size 100

V_iw = 1:100;
F_iw = @(A) sfo_fn_iwata(length(V_iw),A);

Now let's run the algorithm:

A = sfo_min_norm_point(F_iw,V_iw)

suboptimality bound: -11571.000000 <= min_A F(A) <= F(A_best) = -6669.000000; delta<=4902.000000
suboptimality bound: -6834.000000 <= min_A F(A) <= F(A_best) = -6834.000000; delta<=0.000000
suboptimality bound: -6834.000000 <= min_A F(A) <= F(A_best) = -6834.000000; delta<=0.000000

A =

 Columns 1 through 15 

 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

 Columns 16 through 30 

 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63

 Columns 31 through 45 

 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78

 Columns 46 through 60 

 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93

 Columns 61 through 67 

 94 95 96 97 98 99 100

Finding a minimum s-t-cut using general submodular function minimization

Instead of unconstrained minimization

we can also solve

Example: We want to find minimum s-t-cuts in a directed graph with adjacency matrix:

G_dir=[0 1 1.2 0 0 0; 1.3 0 1.4 0 0 0; 1.5 1.6 0 1.7 0 0; 0 0 1.8 0 1.9 2; 0 0 0 2.1 0 2.2; 0 0 0 2.3 2.4 0]

G_dir =

 0 1.0000 1.2000 0 0 0
 1.3000 0 1.4000 0 0 0
 1.5000 1.6000 0 1.7000 0 0
 0 0 1.8000 0 1.9000 2.0000
 0 0 0 2.1000 0 2.2000
 0 0 0 2.3000 2.4000 0

We define the directed cut function:

F_cut_dir = @(A) sfo_fn_cutfun(G_dir,A);

Now let's try to find the minimum cut separating node 1 and 6 in graph G_dir:

A16 = sfo_s_t_mincut(F_cut_dir,V_G,1,6)

suboptimality bound: -0.900000 <= min_A F(A) <= F(A_best) = 4.300000; delta<=5.200000
suboptimality bound: -0.500000 <= min_A F(A) <= F(A_best) = 1.700000; delta<=2.200000
suboptimality bound: -0.500000 <= min_A F(A) <= F(A_best) = 1.700000; delta<=2.200000
suboptimality bound: 1.700000 <= min_A F(A) <= F(A_best) = 1.700000; delta<=0.000000

A16 =

 2 3 1

Now let's separate node 4 and 6:

A46 = sfo_s_t_mincut(F_cut_dir,V_G,4,6)

suboptimality bound: -4.400000 <= min_A F(A) <= F(A_best) = 7.000000; delta<=11.400000
suboptimality bound: -1.800000 <= min_A F(A) <= F(A_best) = 3.900000; delta<=5.700000
suboptimality bound: -1.800000 <= min_A F(A) <= F(A_best) = 3.900000; delta<=5.700000
suboptimality bound: 3.900000 <= min_A F(A) <= F(A_best) = 3.900000; delta<=0.000000

A46 =

 1 2 3 4

Clustering using mutual information

Now we use the Greedy Splitting algorithm for clustering. We first generate 30 data points at random:

n = 10; C1 = randn(2,n)+5; C2 = randn(2,n); C3 = randn(2,n)-5; X = [C1';C2';C3'];

Here's the data points:

figure
plot(X(:,1),X(:,2),'b.'); hold on

Now we use a Gaussian kernel function to measure similarity of the data points:

D = dist(X')/2; sigma_cl = exp(-D.^2)+eye(3*n)*.01;

These points will make up our ground set V:

V = 1:(3*n);

We use entropy as a measure of cluster inhomogeneity:

E = @(A) sfo_fn_entropy(sigma_cl,A);

Now do greedy splitting with 2 clusters:

P = sfo_greedy_splitting(E,V,2)
figure
plot(X(P{1},1),X(P{1},2),'bx'); hold on
plot(X(P{2},1),X(P{2},2),'ro');

iteration 1, cluster 1

P = 

 [1x10 double] [1x20 double]

Now do greedy splitting with 3 clusters:

P = sfo_greedy_splitting(E,V,3)

figure
plot(X(P{1},1),X(P{1},2),'bx'); hold on
plot(X(P{2},1),X(P{2},2),'ro');
plot(X(P{3},1),X(P{3},2),'ks'); hold off

iteration 1, cluster 1
iteration 2, cluster 1
iteration 2, cluster 2

P = 

 [1x10 double] [1x10 double] [1x10 double]

Image denoising using submodular function minimization

We will now do inference in a Markov Random Field using submodular function minimization.

Let's create a binary image of size 40x40, with a white square in front: of a black square

D = 40; S = 10;
img = zeros(D); img(S:(D-S),S:(D-S))=1;

now we perturb it with 20% noise:

noiseProb = 0.2;
mask = rand(size(img))>(1-noiseProb);
imgN = img; imgN(mask)=1-imgN(mask);

Our ground set will be the set of all pixels:

V = 1:numel(img);

Here's the original image and the noisy copy:

figure;
subplot(121); imshow(img); title('original image')
subplot(122); imshow(imgN); title('noisy image')

We define potential functions for the Markov Random Field (ising model). Here, coeffPix is the coefficient for the pixel potentials, the other coefficients enforce smoothness.

coeffPix = 1; coeffH = 1; coeffV = 1; coeffD = .0;

Now we define a submodular function measuring the "energy" of smooth approximations to the noisy image imgN:

val0 = sfo_fn_ising(imgN,[], coeffPix, coeffH, coeffV, coeffD);
F = @(A) sfo_fn_ising(imgN,A, coeffPix, coeffH, coeffV, coeffD)-val0;

Initialize with noisy image:

Ainit = find(imgN(:));

Now we minimize F using the min-norm-point algorithm. We allow ourselves to be at most 1 unit of energy (1 pixel) off the optimal solution:

tic; AD = sfo_min_norm_point(F,V,Ainit,1.999); toc

suboptimality bound: -843.000000 <= min_A F(A) <= F(A_best) = 1725.000000; delta<=2568.000000
suboptimality bound: -327.513820 <= min_A F(A) <= F(A_best) = 789.000000; delta<=1116.513820
suboptimality bound: -250.470243 <= min_A F(A) <= F(A_best) = -19.000000; delta<=231.470243
suboptimality bound: -238.981194 <= min_A F(A) <= F(A_best) = -119.000000; delta<=119.981194
suboptimality bound: -232.609701 <= min_A F(A) <= F(A_best) = -126.000000; delta<=106.609701
suboptimality bound: -226.994977 <= min_A F(A) <= F(A_best) = -155.000000; delta<=71.994977
suboptimality bound: -220.836858 <= min_A F(A) <= F(A_best) = -167.000000; delta<=53.836858
suboptimality bound: -213.784843 <= min_A F(A) <= F(A_best) = -167.000000; delta<=46.784843
suboptimality bound: -206.959990 <= min_A F(A) <= F(A_best) = -175.000000; delta<=31.959990
suboptimality bound: -206.361698 <= min_A F(A) <= F(A_best) = -185.000000; delta<=21.361698
suboptimality bound: -205.753533 <= min_A F(A) <= F(A_best) = -187.000000; delta<=18.753533
suboptimality bound: -204.835702 <= min_A F(A) <= F(A_best) = -189.000000; delta<=15.835702
suboptimality bound: -203.340390 <= min_A F(A) <= F(A_best) = -191.000000; delta<=12.340390
suboptimality bound: -201.515647 <= min_A F(A) <= F(A_best) = -191.000000; delta<=10.515647
suboptimality bound: -200.872117 <= min_A F(A) <= F(A_best) = -191.000000; delta<=9.872117
suboptimality bound: -199.959495 <= min_A F(A) <= F(A_best) = -191.000000; delta<=8.959495
suboptimality bound: -198.246040 <= min_A F(A) <= F(A_best) = -191.000000; delta<=7.246040
suboptimality bound: -196.203150 <= min_A F(A) <= F(A_best) = -191.000000; delta<=5.203150
suboptimality bound: -195.752745 <= min_A F(A) <= F(A_best) = -191.000000; delta<=4.752745
suboptimality bound: -195.352582 <= min_A F(A) <= F(A_best) = -191.000000; delta<=4.352582
suboptimality bound: -194.957251 <= min_A F(A) <= F(A_best) = -191.000000; delta<=3.957251
suboptimality bound: -194.847216 <= min_A F(A) <= F(A_best) = -191.000000; delta<=3.847216
suboptimality bound: -194.720262 <= min_A F(A) <= F(A_best) = -191.000000; delta<=3.720262
suboptimality bound: -194.534848 <= min_A F(A) <= F(A_best) = -192.000000; delta<=2.534848
suboptimality bound: -194.284503 <= min_A F(A) <= F(A_best) = -192.000000; delta<=2.284503
suboptimality bound: -194.043847 <= min_A F(A) <= F(A_best) = -192.000000; delta<=2.043847
suboptimality bound: -193.978602 <= min_A F(A) <= F(A_best) = -192.000000; delta<=1.978602
suboptimality bound: -193.978602 <= min_A F(A) <= F(A_best) = -192.000000; delta<=1.978602
Elapsed time is 80.151602 seconds.

Now let's plot the reconstruction:

imgD = zeros(size(img)); imgD(AD)=1;
figure;
subplot(131); imshow(img); title('original image')
subplot(132); imshow(imgN); title('noisy image')
subplot(133); imshow(imgD); title('reconstructed image')

PART 2) MAXIMIZATION OF SUBMODULAR FUNCTIONS

We can also (near-optimally) solve some maximization problems using submodular functions. We first study constrained optimization:

The lazy greedy algorithm

We will first explore the lazy greedy algorithm on an experimental design problem. More specifically, we want to choose sensing locations to optimally predict pH values on a lake in Merced, California.

As objective function, we use the mutual information criterion of Caselton & Zidek '84 based on a Gaussian Process that we trained on the pH data.

Let's load the data, which contains covariance matrix of the GP, merced_data.sigma, and locations (coordinates) merced_data.coords:

merced_data = load('merced_data');

Define the ground set for the submodular functions:

V_sigma = 1:size(merced_data.sigma,1);

These correspond to the following locations that we can select for sensing:

figure
plot(merced_data.coords(:,1),merced_data.coords(:,2),'k.');
xlabel('Horizontal location along transect'); ylabel('Vertical location (depth)'); title('Possible sensing locations on the lake');

Here's the mutual information criterion:

F_mi = @(A) sfo_fn_mi(merced_data.sigma,A);

Note: mutual information is submodular, but not monotonic! But for k << n it's approximately monotonic, and that's enough (see Krause et al., JMLR '08).

We greedily pick 15 sensor locations for maximizing mutual information:

k = 15;
[A,scores,evals] = sfo_greedy_lazy(F_mi,V_sigma,k); A

A =

 4 15 3 37 82 34 63 35 81 83 60 40 17 19 76

We can find out how many evaluations the naive greedy algorithm would have taken:

nevals = sum(length(V_sigma):-1:(length(V_sigma)-k+1)); evals = sum(evals);
disp(sprintf('Lazy evaluations: %d, naive evaluations: %d, savings: %f%%',evals,nevals,100*(1-evals/nevals)));

Lazy evaluations: 337, naive evaluations: 1185, savings: 71.561181%

Now let's display the chosen locations:

figure
plot(merced_data.coords(:,1),merced_data.coords(:,2),'k.'); hold on
plot(merced_data.coords(A,1),merced_data.coords(A,2),'bs','markerfacecolor','blue');
xlabel('Horizontal location along transect'); ylabel('Vertical location (depth)'); title('Chosen sensing locations (blue squares)');

We can also compute bounds on how far away the greedy solution is from the optimal solution:

bound = sfo_maxbound(F_mi,V_sigma,A,k);
disp(sprintf('Greedy score F(A) = %f; \nNemhauser (1-1/e) bound: %f; \nOnline bound: %f',scores(end),scores(end)/(1-1/exp(1)),bound));

Greedy score F(A) = 17.510471; 
Nemhauser (1-1/e) bound: 27.701158; 
Online bound: 25.974139

Here's how the mutual information increases as we greedily pick more and more sensors:

figure
plot(0:k,[0 scores]); xlabel('Number of elements'); ylabel('submodular utility');

Optimizing greedily over matroids

Sidenote: We also can do greedy optimization over a matroid. To do that, we define a function that takes a set A and outputs 1 if A is independent, 0 otherwise.

Example: The uniform matroid: A independent if and only if length(A)<=k'

checkIndep = @(A) (length(A)<=k);

Now let's run the greedy algorithm, given infinite budget:

C = ones(1,length(V_sigma)); % unit cost
[A,scores,evals] = sfo_greedy_lazy(F_mi,V_sigma,inf,C,0,checkIndep); A

A =

 4 15 3 37 82 34 63 35 81 83 60 40 17 19 76

The lazy greedy coverage algorithm

We can also use the greedy algorithm for solving problems of the form

Example: Picking best sensor locations to achieve a specified amount of mutual information.

Here we want to achieve quota Q:

Q = 5;

Let's run the covering algorithm:

[A,stat] = sfo_cover(F_mi,V_sigma,Q);
disp(sprintf('Coverage possible: %d, Cost: %f',stat,length(A))); A

Coverage possible: 1, Cost: 3.000000

A =

 4 15 1

Now let's try to shoot for a higher quota:

Q = 30;
[A,stat] = sfo_cover(F_mi,V_sigma,Q);
disp(sprintf('Coverage possible: %d, Cost: %f',stat,length(A))); A

Coverage possible: 0, Cost: 38.000000

A =

 Columns 1 through 15 

 4 15 3 37 82 34 63 35 81 83 60 40 17 19 76

 Columns 16 through 30 

 57 50 14 16 42 32 74 20 43 73 29 66 8 75 49

 Columns 31 through 38 

 51 28 72 65 58 52 25 55

The CELF algorithm for budgeted maximization

We can also near-optimally solve problems of the form

where C(A) is an additive cost function. Here, we just make up some cost function:

C = 1:length(V_sigma);

First use a small budget, B = 2. Here, the unit cost greedy solution is better.

A = sfo_celf(F_mi,V_sigma,2,C)

Unit cost solution = 1.721224, Cost-benefit solution = 1.707856

A =

 2

Now use a large budget, B = 15. Here, the cost/benefit greedy solution is better.

A = sfo_celf(F_mi,V_sigma,15,C)

Unit cost solution = 3.815122, Cost-benefit solution = 4.386112

A =

 1 2 8 4

CELF will always choose the better of unit cost or cost/benefit greedy solutions.

Submodular-supermodular procedure of Narasimhan & Bilmes

The submodular-supermodular procedure is a general purpose algortihm for minimizing the difference between two submodular functions. The algorithm is guaranteed to converge to a locally optimal solution.

Here, we will use the submodular-supermodular procedure to do experimental design.

We want to choose subset A of sensor locations in the Lake Merced ph-estimation data set') that maximizes F_mi(A)-|A|, where F_mi is the mutual information criterion. We set F(A) = A, and G(A) = F_mi(A), and minimize F(A)-G(A)

G = F_mi;
F = @(A) length(A);

Now let's run the submodular-supermodular procedure:

A = sfo_sssp(F,G,V_sigma)
disp(sprintf('\n\nMutual information MI(A) = %f, Cost |A| = %d',G(A),F(A)));

suboptimality bound: -2.328092 <= min_A F(A) <= F(A_best) = -2.328092; delta<=-0.000000
suboptimality bound: -2.328092 <= min_A F(A) <= F(A_best) = -2.328092; delta<=0.000000
suboptimality bound: -2.504810 <= min_A F(A) <= F(A_best) = -2.504810; delta<=-0.000000
suboptimality bound: -2.504810 <= min_A F(A) <= F(A_best) = -2.504810; delta<=0.000000
suboptimality bound: -2.572455 <= min_A F(A) <= F(A_best) = -2.572455; delta<=0.000000
suboptimality bound: -2.572455 <= min_A F(A) <= F(A_best) = -2.572455; delta<=0.000000
suboptimality bound: -2.902415 <= min_A F(A) <= F(A_best) = -2.902415; delta<=-0.000000
suboptimality bound: -2.902415 <= min_A F(A) <= F(A_best) = -2.902415; delta<=0.000000
suboptimality bound: -2.942741 <= min_A F(A) <= F(A_best) = -2.942741; delta<=0.000000
suboptimality bound: -2.942741 <= min_A F(A) <= F(A_best) = -2.942741; delta<=0.000000
suboptimality bound: -2.942741 <= min_A F(A) <= F(A_best) = -2.942741; delta<=0.000000
suboptimality bound: -2.942741 <= min_A F(A) <= F(A_best) = -2.942741; delta<=0.000000

A =

 2 5 11 34 37 76 84



Mutual information MI(A) = 9.942741, Cost |A| = 7

Let's again plot the chosen sensing locations:

figure
plot(merced_data.coords(:,1),merced_data.coords(:,2),'k.'); hold on
plot(merced_data.coords(A,1),merced_data.coords(A,2),'b*','markerfacecolor','blue','markersize',10);
xlabel('Horizontal location along transect'); ylabel('Vertical location (depth)'); title('Chosen locations by SSSP (blue stars)');

The Data-Correcting algorithm for maximizing general submodular functions

So far, we've seen algorithms for approximate maximization of submodular functions. We can also try to find the optimal solution (even though that could take exponential time).

The data-correcting algorithm is a branch and bound algorithm for solving

up to a desired accuracy \delta.

Here, we will use the Data Correcting algorithm for finding a maximum directed cut in a graph:

A = sfo_max_dca_lazy(F_cut_dir,V_G,0);
disp(sprintf('Cut value = %f',F(A)));

new best solution: F(A) = 8.800000.

A =

 4 2 6

new best solution: F(A) = 9.000000.

A =

 1 3 6

new best solution: F(A) = 9.200000.

A =

 2 3 6

new best solution: F(A) = 9.500000.

A =

 6 3

Cut value = 2.000000

We can also do constrained maximization, e.g., to find the best cut of size k:

To do that, we define a new submodular function that's -inf if length(A) > k:

FT = @(A) F_cut_dir(A) - 1e10*max(0,length(A)-1);

Now let's run the data correcting algorithm:

A = sfo_max_dca_lazy(FT,V_G,0);
disp(sprintf('Cut value = %f',FT(A)));

new best solution: F(A) = 5.700000.

A =

 4

Cut value = 5.700000

The SATURATE algorithm for robust optimization

SATURATE is a general purpose algorithm for approximately solving problems of the the form

for a collection of monotonic submodular utility functions F_i.

Example: Here let's use the SATURATE algorithm to minimize worst-case variance in GP regression on Merced Lake.

First pick 10 locations greedily using the greedy algorithm:

k = 10;

As objective function, we use variance reduction:

F_var = @(A) sfo_fn_sat_varred(merced_data.sigma,A,0);

Note: Variance reduction is not always submodular, but under certain conditions on the covariance matrix (see Das & Kempe, STOC '08).

We now use the greedy algorithm to minimize the average variance:

AG = sfo_greedy_lazy(F_var,V_sigma,k)

AG =

 24 66 41 12 82 51 15 63 73 34

Now let's use the SATURATE algorithm to minimize the worst case variance.

We need to define the thresholded version of variance reduction function:

F_sat_var = @(A,thr) sfo_fn_sat_varred(merced_data.sigma,A,thr);

Now we can run SATURATE:

AS = sfo_saturate(F_sat_var,V_sigma,k)

interval [0.000000,3.684634], size(ssetMax,Min)=0,0
interval [0.000000,1.842317], size(ssetMax,Min)=0,0
interval [0.000000,0.921159], size(ssetMax,Min)=0,0
interval [0.000000,0.460579], size(ssetMax,Min)=0,0
interval [0.000000,0.230290], size(ssetMax,Min)=3,0
interval [0.115145,0.230290], size(ssetMax,Min)=3,11
interval [0.115145,0.172717], size(ssetMax,Min)=7,11
interval [0.115145,0.143931], size(ssetMax,Min)=8,11
interval [0.115145,0.129538], size(ssetMax,Min)=9,11
interval [0.115145,0.122341], size(ssetMax,Min)=10,11
interval [0.115145,0.118743], size(ssetMax,Min)=10,11
interval [0.115145,0.116944], size(ssetMax,Min)=10,11
interval [0.116044,0.116944], size(ssetMax,Min)=10,11
interval [0.116044,0.116494], size(ssetMax,Min)=10,11
interval [0.116269,0.116494], size(ssetMax,Min)=10,11
interval [0.116269,0.116382], size(ssetMax,Min)=10,11
interval [0.116269,0.116326], size(ssetMax,Min)=10,11
interval [0.116269,0.116297], size(ssetMax,Min)=10,11
interval [0.116283,0.116297], size(ssetMax,Min)=10,11
interval [0.116283,0.116290], size(ssetMax,Min)=10,11
interval [0.116283,0.116287], size(ssetMax,Min)=10,11
interval [0.116285,0.116287], size(ssetMax,Min)=10,11

AS =

 31 70 49 16 81 72 55 8 83 58

Here's the locations selected by both algorithms:

figure
plot(merced_data.coords(:,1),merced_data.coords(:,2),'k.'); hold on
plot(merced_data.coords(AG,1),merced_data.coords(AG,2),'bs','markerfacecolor','blue');
xlabel('Horizontal location along transect'); ylabel('Vertical location (depth)'); title('Greedy locations (blue squares) and SATURATE locations (green diamonds)');
plot(merced_data.coords(AS,1),merced_data.coords(AS,2),'gd','markerfacecolor','green');hold off

Here's the maximum remaining variance under both solutions:

disp(sprintf('max remaining variance: Greedy = %f, SATURATE = %f',sfo_eval_maxvar(merced_data.sigma,AG),sfo_eval_maxvar(merced_data.sigma,AS)));

max remaining variance: Greedy = 0.157456, SATURATE = 0.116368

PART 3) MISCELLANEOUS OTHER FUNCTIONS

The toolbox also provides some other functions that are useful when working with submodular functions.

The Lovasz extension

Every submodular function F induces a convex function, the Lovasz extension, which generalizes F to the unit cube (and the positive orthant).

We will explore the Lovasz extension on the 2 element example function from the tutorial slides at http://www.submodularity.org:

F_ex = @sfo_fn_example;
V_ex = 1:2;

We can now define the Lovasz extension:

g = @(w) sfo_lovaszext(F_ex,V_ex,w);

Now we want to compate F_ex({b}) which should equal g([0,1]). Let's first compute the characteristic vector for set A={b}:

A = [2];
w = sfo_charvector(V_ex,A)
% Now we comparing the function F_ex(A) with Lovasz extension g(w):
disp(sprintf('A = {b}; F(A) = %f; g(wA) = %f',F_ex(A),g(w)));

w =

 0 1

A = {b}; F(A) = 2.000000; g(wA) = 2.000000

Now we plot the Lovasz extension:

[X,Y] = meshgrid(0:.05:1,0:.05:1);
Z = zeros(size(X));
for i = 1:size(X,1)
 for j = 1:size(X,2)
 Z(i,j) = g([X(i,j),Y(i,j)]);
 end
end
figure
surf(X,Y,Z)
xlabel('w_{a}'); ylabel('w_{b}'); zlabel('g(w)')
title('Lovasz extension for example from tutorial: F(\{\})=0,F(\{a\})=-1,F(\{b\})=2,F(\{a,b\})=0 ');

Bounds on optimal solution for minimization

We can also get bounds on the optimal solution for minimization. Let's suppose we guessed solution A={a}:

A = [1];

Let's compute the bound:

bound = sfo_minbound(F_ex,V_ex,A);
disp(sprintf('bound = %f <= min F(A) <= F(B) = %f; ==> A is optimal!',bound,F_ex(A)));

bound = -1.000000 <= min F(A) <= F(B) = -1.000000; ==> A is optimal!

The polyhedron greedy algorithm

Here's an example of the polyhedron greedy algorithm that solves the LP

Let's start with set A = {a}:

A = [1];

Get the characteristic vector:

w = sfo_charvector(V_ex,A);

Run the polyhedron greedy algorithm:

xw = sfo_polyhedrongreedy(F_ex,V_ex,w)

xw =

 -1 1

And now it should hold that

w*xw'

ans =

 -1

Published with MATLAB® 7.2

© 2007 Andreas Krause (krausea at cs dot cmu dot edu)