How can I optimize Clarke-Wright heuristic for TSP in Matlab?

Question

I have implemented the Clarke-Wright huristic to solve TSP (based on the pseudo-code here). I have attached my implementation in Matlab. However it is not fast enough for me and takes O(n²) space (because of pairwise distances). I wonder if there is any theoretic or practical optimization I can apply to reduce the complexity (specially, space complexity). It would be grateful if anyone can help me.

function [tour, length] = clarke_wright (data)

n=size(data,1); % number of records

center = mean(data,1); % mean of data

hubIdx = knnsearch(data,center,'k',1);  % nearest record to center

distances = dist(data,data');   % this requires O(n^2) space :(

savings = zeros(n);  % place to store the saving after adding an edge %

% Can be more vectorized? %
for i=1:n    
    if i==hubIdx
        continue;
    end
        savings(i,(i+1):n)=distances(i,hubIdx)+distances(hubIdx,(i+1):n)-distances(i,(i+1):n);
end

minParent = 1:n;

[~,si] = sort(savings(:),'descend');
si=si(1:(end/2));

Vh = zeros(1,n);
Vh(hubIdx) = 1;
VhCount = n-1;
degrees = zeros(1,n);

selectedIdx = 1;  % edge to try for insertion

tour = zeros(n,2);
curEdgeCount = 1;

while VhCount>2
    i = mod(si(selectedIdx)-1,n)+1;
    j = floor((si(selectedIdx)-1)/n)+1;

    if Vh(i)==0 && Vh(j)==0 && (minParent(i)~=minParent(j)) && i~=j && i~=hubIdx && j~=hubIdx     % always all degrees are <= 2, so it is not required to check them
%     if (minParent(i)~=minParent(j)) && isempty(find(degrees>2, 1)) && i~=j && i~=hubIdx && j~=hubIdx && Vh(i)==0 && Vh(j)==0
        degrees(i)=degrees(i)+1;
        degrees(j)=degrees(j)+1;
        tour(curEdgeCount,:) = [i,j];

        if minParent(i)<minParent(j)
            minParent(minParent==minParent(j))=minParent(i);
        else
            minParent(minParent==minParent(i))=minParent(j);            
        end


        curEdgeCount = curEdgeCount + 1;

        if degrees(i)==2
            Vh(i) = 1;
            VhCount = VhCount - 1;
        end

        if degrees(j)==2
            Vh(j) = 1;
            VhCount = VhCount - 1;
        end
    end
    selectedIdx = selectedIdx + 1;

end

remain = find(Vh==0);
n1=remain(1);
n2=remain(2);

tour(curEdgeCount,:) = [hubIdx n1];
curEdgeCount = curEdgeCount + 1;

tour(curEdgeCount,:) = [hubIdx n2];

tour = stitchTour(tour);
tour=tour(:,1)';
length=distances(tour(end),tour(1));
for i=1:n-1  % how can I vectorize these lines?
    length=length+distances(tour(i),tour(i+1));
end
end


function tour = stitchTour(t) % uniforms the tour [a b; b c; c d; d e;.... ]

n=size(t,1);

[~,nIdx] = sort(t(:,1));
t=t(nIdx,:);

tour(1,:) = t(1,:);
t(1,:) = -t(1,:);
lastNodeIdx = tour(1,2);

for i=2:n
    nextEdgeIdx = find(t(:,1)==lastNodeIdx,1);
    if ~isempty(nextEdgeIdx)
        tour(i,:) = t(nextEdgeIdx,:);
        t(nextEdgeIdx,:)=-t(nextEdgeIdx,:);
    else
        nextEdgeIdx = find(t(:,2)==lastNodeIdx,1);
        tour(i,:) = t(nextEdgeIdx,[2 1]);
        t(nextEdgeIdx,:)=-t(nextEdgeIdx,:);
    end
    lastNodeIdx = tour(i,2);
end


end

Answer 1

This is what you can do if space is an issue (will probably reduce calculation speed a bit).

I have not really looked into your code, but judging from the pseudo code this should do the trick:

For each pair or points, calculate the savings created by connecting them.

If this is better than the best savings found so far, update the best savings, and remember the two points.

After checking all pairs just implement the best savings.

This way you will barely require extra space at all.