Code covered by the BSD License

# Submodular Function Optimization

### Andreas Krause (view profile)

28 Jun 2008 (Updated )

This toolbox provides functions for maximizing and minimizing submodular set functions.

sfo_pspiel_fixed_r(F,V,Q,D,R,dists,pdallok,Vroot)
```% Helper function for sfo_pspiel, by Andreas Krause (krausea@gmail.com)
% Example: See sfo_tutorial.m
% solve pSPIEL for fixed value of R
% last parameter pdallok controls whether it's ok to use the trivial PD or not
function [A, E, result] = sfo_pspiel_fixed_r(F,V,Q,D,R,dists,pdallok,Vroot)

result.failed = 0;

done = 0; trial = 0; maxTrials = 10;
while ~done && trial<maxTrials;
trial = trial+1;
[pd,done,a,successrate] = sfo_pspiel_pd(V,dists,R,1/2);
if ~pdallok && length(pd)==1 && isempty(sfo_setdiff_fast(V,pd{1}))
% we got a trivial partition. Let's avoid that by reducing R
done = 0;
R = R/2;
elseif ~done
% we didn't find a PD that achieves success probability >= 1/2.
% Increase R.
R = 2*R;
end
end

% construct modular approximation graph (MAG)

% for rooted k-MST, we build the MAG on the residual submodular function
% hence the reward of the root node is ignored --> need to reduce quota
if Vroot>0
valRoot = F(Vroot);
else
valRoot = 0;
end
% compute the Quota-MST solution on the MAG

if isempty(allels)
% failed to satisfy quota!
A = []; E = [];
result = prepareResult(1,[],[],[],0,1e99,0,Q,R,pd,a,successrate,gs,gimp);
return
else
% transfering solution from MAG back to original graph
selectednodes = sfo_pspiel_transfer_back(allels,gs,node2cluster);

% now greedily augment to satisfy quota
selectednodes = sfo_pspiel_augment(F,V,selectednodes,Q,dists);

% solve Steiner tree problem to connect selectednodes
[totalCost,selectededges,supportnodes] = sfo_pspiel_get_cost(selectednodes,D,dists);

A = [selectednodes supportnodes];
E = selectededges;

selectedVal = F(selectednodes); totalVal = F(A);

% generate result structure
result = prepareResult(0,selectednodes,supportnodes,selectededges,totalVal,totalCost,selectedVal,Q,R,pd,a,successrate,gs,gimp);
end

%% Now we set up the modular approximation graph
nnodes = 0; nedges = 0; %these count the #nodes and #edges used

reward = []; % modular reward per node (takes indices in MAG)
node2cluster = []; %first column is cluster id (in pd), second column is node # in cluster acc. to greedy

edge = []; % indices from to in MAG
edgeweight = []; % edgeweights in MAG
nodeindex =  {}; % nodeindex{i} contains MAG indicies of nodes in cluster i
roots = []; % set of cluster centers, i.e., roots(i) = nodeindex{i}(1)
gs = {}; %greedy sequences (indices in V)
gimp = {}; %greedy improvements

if Vroot>0
% need to make sure root is part of MAG
% first find cluster that contains root
rootPDind = [];
for i = 1:length(pd)
if sum(pd{i}==Vroot)>0 && length(pd{i})>1
rootPDind = i;
break;
end
end
% make root a separate cluster
if ~isempty(rootPDind)
% remove
pd{rootPDind} = sfo_setdiff_fast(pd{rootPDind},Vroot);
end
pd{end+1} = Vroot;
% now define a new residual submodular function F'(A) = F(A u root)-F(root)
F = sfo_fn_residual(F,Vroot);
end

for i = 1:length(pd);
% compute greedy subsets and modular approximation per cluster
ni = length(pd{i}); %number of clusters

if pd{i}(1) ~= Vroot || length(pd{i})>1
[gs{i},scores] = sfo_greedy_lazy(F,pd{i},ni);
else
% in case we have a root, that will give value 0 and greedy will
% not pick it.
gs{i} = Vroot;
scores = 1e-99;
end
if isempty(gs{i}) % require that the greedy algorithm picks at least 1 element per cluster. Maybe F == 0?
error('greedy algorithm did not pick any elements!');
end

% compute greedy improvements F(Aj)-F(A_{j-1}) for i-th cluster
gimp{i} = scores(1:end)-[0 scores(1:(end-1))];

% compute increase in MST costs per cluster
mstweight = zeros(1,ni);
for j = 1:length(gs{i})
mstweight(j) = sfo_pspiel_mst(dists(gs{i}(1:j),gs{i}(1:j)));

% add a new node to MAG
nnodes = nnodes+1;
reward(nnodes) = gimp{i}(j); % assign greedy improvment score
node2cluster(nnodes,:)=[i,j]; % i.e., MAG node nnodes is j-th node in i-th cluster

if j>1
% not cluster center; need to attach node to center by an edge
nedges = nedges+1;
edge(nedges,:)=[nnodes-1,nnodes];

% edge cost is difference in MST costs
edgeweight(nedges) = max(0,mstweight(j)-mstweight(j-1));
nodeindex{i} = [nodeindex{i},nnodes];
else
% new cluster center. Will be a possible root for the k-MST
% algorithm
nodeindex{i} = nnodes;
roots = [roots,nnodes];
end

end
end
if Vroot>0
roots = roots(end);
end

% now create the "core" of G', which is a complete subgraph (clique) connecting the
% cluster centers.
for i = 1:length(pd)
for j = (i+1):length(pd)
nedges = nedges+1;
% connect centers of clusters i and j
edge(nedges,:) = [nodeindex{i}(1),nodeindex{j}(1)];

% edge weight is shortest path distance between cluster centers
edgeweight(nedges) = dists(gs{i}(1),gs{j}(1));
end
end

% generate adjacency matrix representation for k-MST algorithm
INFINITY = 1e99; %guardian value for non-existent edges
for j=1:nedges
end
for i = 1:nnodes
end

%% Now transfer back from MAG to original graph
% A is the MAG nodes selected
function selectednodes = sfo_pspiel_transfer_back(A,gs,node2cluster)
% now we need to figure out how many nodes are selected in each cluster
numpercl = zeros(1,length(gs));
for y = A
% for indexing
clind = node2cluster(y,1); % the cluster index
clpos = node2cluster(y,2); % position in the cluster
numpercl(clind)=max(numpercl(clind),clpos);
end
selectednodes = []; % the real indices selected in V

% now select nodes in V: for each cluster i pick the first numpercl(i) greedy elements
for i=1:length(gs)
selectednodes = [selectednodes, gs{i}(1:numpercl(i))];
end

%% Greedily augment the network to achieve quota F(A)>=Q
% Due to (R,gamma)-locality, the selected nodes might achieve quota Q on
% the MAG, but less than Q on F. Here, we simply greedily add nodes until
% the quota is met.
function A = sfo_pspiel_augment(F,V,A,Q,dists)
while 1
currentScore = F(A);
if (currentScore>=Q) % satisfied quota
break
end
Vc = sfo_setdiff_fast(V,A); % elements that are still available
Fcur = init(F,A); % optimization for speed
scores  = zeros(1,length(Vc)); % scores for adding each available element
for i = 1:length(Vc)
Fcur = init(Fcur,A);
newCost = min(dists(A,Vc(i))); % cost of adding = cheapest shortest path connection cost
newScore = max(Q,inc(Fcur,A,Vc(i))); % benefit of adding; don't care if > Q
scores(i) = (newScore-currentScore)/max(newCost,1e-10); % score = benefit / cost ratio
end
[tmp best] = max(scores); % find best element
A = [A Vc(best)];
end

%% prepare the result output structure
result.failed = failed;
result.selectedNodes = selectednodes;
result.supportNodes = supportnodes;
result.selectedEdges = selectededges;
result.totalVal = totalVal;
result.totalWeight = totalWeight;
result.selectedVal = selectedVal;
result.parameters.quotaMI = quota;
result.parameters.localityR = R;