function [W,L] = TransformationMatrix(Incidence)
%F. de Meijer and R. Sotirov - 31 July 2020
%Code belonging to Algorithm 1 of the paper 'SDP-based bounds for the
%quadratic cycle cover problem via cutting plane augmented Lagrangian
%methods and reinforcement learning'
%This code computes a sparse expression for the transformation matrix
%needed in the facial reduction of our SDP relaxations for the QCCP. It
%exploits the bipartite representation to do so and is valid for any
%directed (simple) graph, provided that the instance is feasible for the QCCP.
% |-------------+-----------------------------------------------------------|
% | Input | Description |
% |-------------+-----------------------------------------------------------|
% | Incidence | m-by-n incidence matrix of a directed graph, where each |
% | | column represents an arc from the node with entry -1 |
% | | to the node with entry 1. |
% |-------------+-----------------------------------------------------------|
% |-------------+-----------------------------------------------------------|
% | Output | Description |
% |-------------+-----------------------------------------------------------|
% | W | sparse expression of the transformation matrix needed in |
% | | the facial reduction for SDP1, SDP2 and SDP3 |
% |-------------+-----------------------------------------------------------|
n = size(Incidence,1);
m = size(Incidence,2);
U = double(Incidence < 0);
V = double(Incidence > 0);
%% Pre-processing step:
%We first need to make sure that Assumption 1 is satisfied. In other words,
%we need to remove all arcs that are never used in a cycle cover. We do this
%by solving m linear CCPs to check whether there exist a cycle cover containing
%each arc. If such a cycle cover does not exist, we set the corresponding entry
%in jVector equal to 1, otherwise it is set to 0.
jVector = zeros(m,1);
for i = 1:m
f = zeros(m,1); f(i) = 1;
[~, val] = linprog(-f, [], [], [U; V], ones(2*n,1), zeros(m,1));
if val == 0
jVector(i) = 1;
end
end
%We now remove the arcs corresponding to a 1 in jVector:
jVector = jVector > 0;
Incidence(:,jVector) = [];
m = m - sum(jVector);
%% Computation of the sparse transformation matrix:
%Initialize W:
W = [];
%The bipartite representation of the graph is called H. Its incidence
%matrix is as follows:
IncidenceH = [U; V];
%Creating H and its corresponding adjacency matrix:
Adjacency = IncidenceH*IncidenceH';
Adjacency = Adjacency - diag(diag(Adjacency));
H = graph(Adjacency);
%Finding connected components in bipartite graph H, using built-in Matlab
%function conncomp:
components = conncomp(H);
numComponents = max(components);
%We treat each connected component separately:
for i = 1:numComponents
%Find vertices and edges in current component. We do this by
%considering the subgraph of the incidence matrix of H implied by the
%current component:
vertices = find(components == i);
subIncidence = IncidenceH(vertices,:);
edges = find(sum(subIncidence,1)>0);
if length(edges) == length(vertices) - 1
%Component is a spanning tree, no columns of W can be found here.
else
%We want to find a spanning tree in the connected component. We do
%this using the built-in function minspantree, that finds a MST
%w.r.t. weight function 1:
M = minspantree(graph(Adjacency(vertices,vertices)));
mstEdges = table2array(M.Edges);
%T is the array that contains the edges in the spanning tree. We
%consruct T iteratively:
T = [];
for j = 1:size(mstEdges,1)
currentEdge = find(subIncidence(mstEdges(j,1),:) == 1 & subIncidence(mstEdges(j,2),:) == 1);
T = [T, currentEdge];
end
%notUsedEdges contains remaining edges in current component that are
%not in T
notUsed = ~ismember(edges, T);
notUsedEdges = edges(notUsed);
%For every edge in the current component that is not in T, we now
%construct a vector w^e. These vectors are linearly independent and
%form a basis for the flow space of H. We construct these iteratively
%using the following for-loop:
for j = 1:length(notUsedEdges)
%Initialize w:
w = zeros(m,1);
thisEdge = notUsedEdges(j);
%Add the current edge to T to get T_new. By construction, T_new
%contains exactly 1 cycle:
T_new = [T, thisEdge];
%smallIncidence equals the incidence matrix of the subgraph
%implied by the edges in T_new:
smallIncidence = IncidenceH(vertices,T_new);
%To detect the cycle, we perform a breadth-first search,
%starting at one of the vertices incident to the added edge:
firstVertices = find(smallIncidence(:,end));
firstVertex = firstVertices(1);
currentVertex = firstVertices(2);
%Initialize the list of visited vertices:
listVisited = [];
listToBeVisited = currentVertex;
%The vector predecessors saves the last edge (in first column)
%and last vertex (in second column) that led to each vertex in
%BFS algorithm:
predecessors = zeros(length(vertices),2);
predecessors(currentVertex,1) = length(T_new); %Since we put thisEdge on last position of T_new
predecessors(currentVertex,2) = firstVertex;
%Since the current subgraph has a cycle, each vertex in the
%cycle can be reached by two distinct paths. In order to
%properly detect the vertices in the cycle, we artificially
%'remove' the last added edge from the current subgraph:
smallIncidence(:,end) = zeros(length(vertices),1);
%Start the BFS:
while ~ismember(firstVertex,listToBeVisited) %We stop when the initial vertex is discovered.
%currentVertex is the new vertex that is visited. offSpring
%contains the edges that are incident to currentVertex:
currentVertex = listToBeVisited(1);
offSpring = find(smallIncidence(currentVertex,:));
%Consider all neighbours of currentVertex not yet visited:
for k = 1:length(offSpring)
thisEdge = offSpring(k);
thisVertex = find(smallIncidence(:,thisEdge)==1);
thisVertex = thisVertex(thisVertex ~= currentVertex);
%Save the predecessor edge and predecessor vertex that
%led to the discovery of thisVertex:
predecessors(thisVertex,1) = thisEdge;
predecessors(thisVertex,2) = currentVertex;
%Put the new vertex on the list of vertices that needs to be
%visited and artificially 'remove' the edge from the
%subgraph:
listToBeVisited = [listToBeVisited, thisVertex];
smallIncidence(:,thisEdge) = zeros(length(vertices),1);
end
%Remove currentVertex from list and put it on list of
%visited vertices.
listToBeVisited(1) = [];
listVisited = [listVisited, currentVertex];
end
%With backward recursion, we want to find the edges in the
%cycle containing edge T_new(end).
thisEdge = 0;
thisVertex = firstVertex;
count = 1;
while thisEdge ~= length(T_new)
%Find last edge before current vertex:
thisEdge = predecessors(thisVertex,1);
%Put 1 or -1 in vector w. NOTE: We do this on position
%T_new(thisEdge), not on position thisEdge, since above we
%worked with a subset of edges. Therefore, we changed
%indexing.
w(T_new(thisEdge)) = (-1)^count;
%Find last vertex before current vertex:
thisVertex = predecessors(thisVertex,2);
count = count + 1;
end
%Add a vector to W:
W = [W, [0; w]];
end
end
end
%Finally, the first column of W is determined by a feasible cycle cover of
%the graph. This can be determined by solving a linear CCP with a constant arc-weight:
f = ones(m,1);
x = linprog(f, [], [], [U; V], ones(2*n,1), zeros(m,1));
W = [[1; x],W];
end