This solves the Lagrangean dual of the TSP (as in the Held-Karp formulation) by solving the dual of the dual (an LP problem) by column generation

Routines

Off[General::spell1]

dualfunction[dual_,dir_,al_]:=Module[{nc,dv},
       onetreefind[
          Table[cost[[i,j]]
          -dual[[i]]-al dir[[i]]
          -dual[[j]]-al dir[[j]],{i,nn},{j,nn}]
                     ];
       dv=2-onetreedegree[onetree];
       costcompute[onetree]+(dual+al dir).dv
]

takeheurtour[
        node_]:=Module[
{marked, currnode, 
          minvalue, nextnode, tour, tourlength, 
          val}, 
marked = 
          Table[False, {nn}]; 
currnode = 
          node; 
tour = {currnode}; 
tourlength = 0; 
Do[
          marked[[currnode]] = True; 
          minvalue = Infinity; 
  Do[ 
            val = costs[[currnode + 1, 
                  i + 1]] ; 
     If[! marked[[i]]  &&  (val <  
                    minvalue)
             ,  minvalue = val ; 
              nextnode = 
                i;
               ]  
          , {i, 
              nn}
          ]; 
       currnode = 
            nextnode; 
       AppendTo[tour, 
            currnode]; 
     tourlength += 
            minvalue;  
  , {nn - 
              1}
   ]; 
     {tour, tourlength +  
            costs[[currnode + 1, 
                node + 1]]}
]

takebestheurtour := 
      Module[{tour, value, besttour, 
          bestvalue}, 
   {besttour, bestvalue} =  
          takeheurtour[1]; 
  Do[{tour, value} = 
            takeheurtour[i]; 
        If[
            value  <  bestvalue
        ,  
            bestvalue = 
              value; 
             besttour = 
              tour
             ]
     , {i, 2, 
            nn}
    ]; 
{besttour, 
          bestvalue}
]

onetreedegree[onetree_]:=Block[{values },
values=Table[0,{nn}];
ft=Flatten[onetree];
Do[values[[ft[[i]]]]++
 ,{i,Length[ft]}];
 values];

costcompute[tree_]:=Sum[cost[[tree[[i,1]],tree[[i,2]]]],
             {i,Length[tree]}]

show[a_,b_,c_,d_]:=Block[{b1,dd1,dd2,dd3},
           Which[b==1,b1=0.005,b==2,b1=0.008,b==3,b1=0.015];
           Which[d=="R",dd1=1;dd2=0;dd3=0;,
                 d=="B",dd1=0;dd2=0;dd3=1;,
                 d=="G",dd1=0;dd2=1;dd3=0;,
                 d=="Y",dd1=1;dd2=1;dd3=0.5;,
                 d=="Bk",dd1=0;dd2=0;dd3=0;
                 ];
           If[c==0,
          Graphics[{Thickness[b1],RGBColor[dd1,dd2,dd3],
          Table[Line[{nodes[[a[[i,1]]]],
                         nodes[[a[[i,2]]]]}
                       ]
                   ,{i,Length[a]}]}],
          Graphics[{Thickness[b1],Dashing[{0.01}],RGBColor[dd1,dd2,dd3],
          Table[Line[{nodes[[a[[i,1]]]],
                         nodes[[a[[i,2]]]]}
                       ]
                   ,{i,Length[a]}]}]
                   ]]

minspantree[arcs_]:=Block[{p,r,kk,lt,ki,kj,ri,rj},
root[k_]:=Block[{kp},kp=k;While[kp!=p[[kp]],kp=p[[kp]]];kp];
p=Table[i,{i,nn}];
r=Table[1,{nn}];
For[kk=1;lt=0;tree={},
    (lt<nn-2)&&(kk<=Length[arcs]),
    kk++,
    ki=arcs[[kk,1]];
    kj=arcs[[kk,2]];
    ri=root[ki];
    rj=root[kj];
    If[ri!=rj,
       Which[r[[ri]]>r[[rj]]
            ,p[[rj]]=ri
            ,r[[rj]]>r[[ri]]
            ,p[[ri]]=rj
            ,r[[ri]]==r[[rj]]
            ,p[[rj]]=ri;
             r[[ri]]++
             ];
        lt++;
        AppendTo[tree,arcs[[kk]]];
       ];
     ];
   If[lt<nn-2,Print["disconnected graph !"];
             Do[pp=ii;
                While[!p[[pp]]==pp,pp=p[[pp]]];
                Print[ii," -> ",pp]
                ,{ii,n}]
     ];

]

onetreefind[costmatrix_]:=Block[
{specarcs,arcs,first,second,val},
specarcs={};
Do[AppendTo[specarcs,{costmatrix[[i,j]],{i,j}}],
{i,2,nn-1},{j,i+1,nn}];
arcs=Transpose[Sort[specarcs]][[2]];
minspantree[arcs];
first={10000,{0,0}};
second={10000,{0,0}};
Do[val=costmatrix[[1,i]];
   If[val < second[[1]] 
      ,If[val >= first[[1]],
          second={val,{1,i}},
          second=first;
          first={val,{1,i}}]
      ]
   ,{i,2,nn}];
onetree=Append[tree,first[[2]]];
AppendTo[onetree,second[[2]]];
   ];

How to use this notebook

First run the section Routines. 

Then generate the data as described in the section Data. You have alternative ways of generating the data.

Data

Input mode n.1 
The arc costs are uniformly randomly  generated
between 1 and 99.
The graph nodes are represented on a circle.

nn=6; seed=10 ->  strong duality
nn=6; seed=11 -> duality gap, dual optimum with two onetrees

nn=4;na=nn(nn-1)/2;

SeedRandom[10];

nodes=Table[5 {Sin[2 N[Pi]/nn i],Cos[2 N[Pi]/nn i]},{i,0,nn-1}];

bb=Table[Random[Integer,{1,99}],{na}]

costs=Table[0,{nn+1},{nn+1}];
k=1;
Do[costs[[1,i+1]]=i;costs[[i+1,1]]=i,{i,nn}];
Do[costs[[i+1,j+1]]=bb[[k]];
   costs[[j+1,i+1]]=bb[[k++]],
{i,nn-1},{j,i+1,nn}];
MatrixForm[costs]

cost=Transpose[Drop[Transpose[Drop[costs,1]],1]];

rangeplot={{-7,5.7},{-7,7}};

aaa=Graphics[{PointSize[0.04],Map[Point ,nodes],
              MapThread[Text,{Range[nn],1.13 nodes}]} ];

bbb=Graphics[Table[{Thickness[0.0001(100-costs[[i+1,j+1]])],
                    Line[{nodes[[i]],nodes[[j]]}]},
                   {i,nn-1},{j,i+1,nn}]];

Show[aaa,bbb,AspectRatio->Automatic,
            PlotRange->rangeplot];

the thicker an arc is, the better it is

User provided initial tour for optimum estimate (list of n nodes)

starttour={1,2,3,4}

tourarcs=Table[{starttour[[i]],starttour[[i+1]]},
    {i,Length[starttour]-1}];
    AppendTo[ tourarcs,{starttour[[-1]],starttour[[1]]}]

est=costcompute[tourarcs]

ccc=show[tourarcs,2,0,"Bk"];
Show[ccc,aaa,AspectRatio->Automatic,PlotRange->rangeplot];

Input mode n.2
The arc costs are typed by the user.
The graph nodes are represented on a circle.

Set #1

nn=6;
bb={21,15,16,12,21,
       16,19,17,22,
          16,14,13,
             16,20, 
                15};

Set #2  (comb inequality)

nn=6;
bb={2, 2,10,10, 1,
       2,10, 1,10,
          1,10,10,
             2, 2,
                2};

Set #3

nn=12;
bb={1,2,9,2,9,3,9,9,9,9,9,
      9,2,9,2,9,3,9,9,9,9,
        1,2,9,9,9,3,9,9,9,
          9,2,9,9,9,3,9,9,
            1,9,9,9,9,3,9,
              9,9,9,9,9,3,
                1,2,9,2,9,
                  9,2,9,2,
                    1,2,9,
                      9,2,
                        1};

Setting-up data     YOU    MUST   DO   THIS   FOR  ANY   DATA  SET!

na=nn(nn-1)/2;

nodes=Table[5 {Sin[2 N[Pi]/nn i],Cos[2 N[Pi]/nn i]},{i,0,nn-1}];

bmax=Max[bb];bmin=Min[ bb];a1=99/(bmax-bmin);a2=1-bmin a1;

costs=Table[0,{nn+1},{nn+1}];
k=1;
Do[costs[[1,i+1]]=i;costs[[i+1,1]]=i,{i,nn}];
Do[costs[[i+1,j+1]]=bb[[k]];
   costs[[j+1,i+1]]=bb[[k++]],
{i,nn-1},{j,i+1,nn}];
MatrixForm[costs]

cost=Transpose[Drop[Transpose[Drop[costs,1]],1]];

rangeplot={{-7,7},{-7,7}};

(*aaa=Graphics[{RGBColor[1,1,1],Map[Disk[#,.4] &,nodes],
              RGBColor[0,0,0],AbsoluteThickness[1],
              Map[Circle[#,.4] &,nodes],
              MapThread[Text,{Range[nn],nodes}]} ];*)

aaa=Graphics[{PointSize[0.04],Map[Point ,nodes],
              MapThread[Text,{Range[nn],1.1 nodes}]} ];

(*aaa=Graphics[{Table[{Point[nodes[[i]]],
      Text[ToString[i],nodes[[i]],0.8{-nodes[[i,1]],-nodes[[i,2]]}]},
              {i,Length[nodes]}]}];*)

bbb=Graphics[Table[{RGBColor[0,1,0],Thickness[0.0001(100-(costs[[i+1,j+1]]a1+a2))],
                    Line[{nodes[[i]],nodes[[j]]}]},
                   {i,nn-1},{j,i+1,nn}]
                   ];

Show[bbb,aaa,AspectRatio->Automatic,
            PlotRange->rangeplot];

the thicker an arc is, the better it is

User provided initial tour for optimum estimate (list of n nodes)

starttour={1,3,2,5,6,4}

tourarcs=Table[{starttour[[i]],starttour[[i+1]]},
    {i,Length[starttour]-1}];
    AppendTo[ tourarcs,{starttour[[-1]],starttour[[1]]}]

est=costcompute[tourarcs]

ccc=show[tourarcs,2,0,"B"];
Show[ccc,aaa,AspectRatio->Automatic,PlotRange->rangeplot];

Input mode n.3
First the nodes are randomly generated within a square.
Then the arc costs are taken as Euclidean distances.

nice seeds=40 (nn=10) -> many iterations
              10 (nn=10) -> conv. comb. of two 1-trees,
                            many iterations
              50 (nn=10) -> solution at 48th iteration,
                            nice graphs.
              40 (nn= 8) -> optimal tour !
              40 (nn= 7) -> optimal tour !
              25 (nn= 5) -> very fast with duality gap
              35 (nn= 5) -> fast with no gap
              25 (nn= 6) -> no few iterations, no gap
              25 (nn=50) -> gap less than 1 ! optimum found by manually turning the onetree into a tour

SeedRandom[25]

nn=50;na=nn(nn-1)/2;

nodes={};
Do[AppendTo[nodes,
 {Random[Integer,{1,99}],Random[Integer,{1,99}]}],{nn}]

nodes  (* coordinates *)

bb=Flatten[Table[Floor[N[Sqrt[(nodes[[i,1]]-nodes[[j,1]])^2+
                    (nodes[[i,2]]-nodes[[j,2]])^2]]],
           {i,nn-1},{j,i+1,nn}]];

costs=Table[0,{nn+1},{nn+1}];
k=1;
Do[costs[[1,i+1]]=i;costs[[i+1,1]]=i,{i,nn}];
Do[costs[[i+1,j+1]]=bb[[k]];
   costs[[j+1,i+1]]=bb[[k++]],
{i,nn-1},{j,i+1,nn}];
MatrixForm[costs];

cost=Transpose[Drop[Transpose[Drop[costs,1]],1]];

rangeplot={{-5,105},{-5,105}};

aaa=Graphics[Table[{Point[nodes[[i]]],
      Text[ToString[i],nodes[[i]],{-2,0}]},{i,Length[nodes]}]];

bbb=Graphics[{Dashing[{0.01}],Table[Line[{nodes[[i]],nodes[[j]]}],
    {i,nn-1},{j,i+1,nn}]}];

Show[aaa,AspectRatio->Automatic,PlotRange->rangeplot];

User provided initial tour for optimum estimate (list of n nodes)

starttour={1,18,4,7,10,13,2,5,3,15,9,17,20,14,11,19,6,12,8,16}

tourarcs=Table[{starttour[[i]],starttour[[i+1]]},
    {i,Length[starttour]-1}];
    AppendTo[ tourarcs,{starttour[[-1]],starttour[[1]]}]

length of tour

est=costcompute[tourarcs]

ccc=show[tourarcs,2,0,"B"];
Show[ccc,aaa,AspectRatio->Automatic,PlotRange->rangeplot];

Heuristic initial tour for manual initialization (list of n nodes). The LP execution is faster if you initialize manually

{starttour,value}=takebestheurtour

tourarcs=Table[{starttour[[i]],starttour[[i+1]]},
    {i,Length[starttour]-1}];
    AppendTo[ tourarcs,{starttour[[-1]],starttour[[1]]}];

length of tour

est=costcompute[tourarcs]

ccc=show[tourarcs,2,0,"B"];
Show[ccc,aaa,AspectRatio->Automatic,PlotRange->rangeplot];

Input mode n.4
Euclidean instances generated by the user

provide nn and  coordinates (between 1 and 100)

nn=12;na=nn(nn-1)/2;

nodes= 30 {{0,0.5},{0.5,0},{1,1},{0,2.5},{1,2},{0.5,3},
         {3,0.5},{2,1},{2.5,0},{2.5,3},{3,2.5},{2,2}}+5

bb=Flatten[Table[Floor[N[Sqrt[(nodes[[i,1]]-nodes[[j,1]])^2+
                    (nodes[[i,2]]-nodes[[j,2]])^2]]],
           {i,nn-1},{j,i+1,nn}]];

costs=Table[0,{nn+1},{nn+1}];
k=1;
Do[costs[[1,i+1]]=i;costs[[i+1,1]]=i,{i,nn}];
Do[costs[[i+1,j+1]]=bb[[k]];
   costs[[j+1,i+1]]=bb[[k++]],
{i,nn-1},{j,i+1,nn}];
MatrixForm[costs];

cost=Transpose[Drop[Transpose[Drop[costs,1]],1]];

rangeplot={{-5,105},{-5,105}};

aaa=Graphics[Table[{Point[nodes[[i]]],
      Text[ToString[i],nodes[[i]],{-2,0}]},{i,Length[nodes]}]];

bbb=Graphics[{Dashing[{0.01}],Table[Line[{nodes[[i]],nodes[[j]]}],
    {i,nn-1},{j,i+1,nn}]}];

Show[aaa,AspectRatio->Automatic,PlotRange->rangeplot];

User provided initial tour for manual initialization (list of n nodes)

starttour={1,3,2,5,6,4}

tourarcs=Table[{starttour[[i]],starttour[[i+1]]},
    {i,Length[starttour]-1}];
    AppendTo[ tourarcs,{starttour[[-1]],starttour[[1]]}]

a={{1,2},{2,9},{9,5},{5,6},{6,3},{3,10},{10,4},{4,7},{7,8},{8,1}};

ccc=show[tourarcs,2,0,"B"];
Show[ccc,aaa,AspectRatio->Automatic,PlotRange->rangeplot];

Input mode n.5
Instances derived by Hamiltonian circuit problems

Seeds : 	(nn=10)
            23

nn=10;na=nn(nn-1)/2;

SeedRandom[23];
nodes=Table[5 {Sin[2 N[Pi]/nn i],Cos[2 N[Pi]/nn i]},{i,0,nn-1}];

d=3;
prob=d/(nn-1);
arcs={};
bb={};
Do[If[Random[]<prob,AppendTo[arcs,{i,j}];AppendTo[bb,0]
                   ,AppendTo[bb,1]]
                   ,{i,1,nn-1},{j,i+1,nn}]

arcs  (*list of arcs*)

bb

bbb=Graphics[{RGBColor[0,0,1],Table[
      Line[{nodes[[arcs[[i,1]]]],nodes[[arcs[[i,2]]]]}],
    {i,Length[arcs]}]}];

rangeplot={{-7,7},{-6,6}};

aaa=Graphics[{PointSize[0.02],Table[{Point[nodes[[i]]],
      Text[ToString[i],nodes[[i]],0.5{-nodes[[i,1]],-nodes[[i,2]]}]},{i,Length[nodes]}]}];
      Show[aaa,bbb,AspectRatio->Automatic,
      PlotRange->rangeplot];

Input mode n.5
Other instances

Italian regional headtowns

nn=20;na=nn(nn-1)/2;

nodes=
{{48, 61}, {1, 88}, {47, 48}, {76, 35}, {31, 73}, {58, 39},
 {31, 64},  {11, 72}, {15, 84}, {55, 31}, {47, 1}, {39, 56},
  {67, 29}, {68, 1}, {40, 42}, {3, 80}, {31, 90}, {50, 85},
   {47, 89}, {39, 82}};

bb =  {610,195,465,215,325,320,510,425,390,1110,165,460,875,285,545,440,445,415,305,
       795,1065,395,925,470,245,180,955,1675,625,1100,1440,745,110,400,585,550,445,
               400,400,190,375,600,610,245,965,175,390,730,110,730,625,630,600,490,
                  670,220,720,945,880,260,690,610,145,455,450,1000,895,900,870,760,
                     530,105,295,210,590,1310,260,735,1075,380,330,225,295,260,155,
                           490,720,740,135,825,385,200,590,225,860,755,760,725,620,
                              225,300,490,1205,155,630,970,280,395,315,395,365,255,
                                140,715,1430,380,860,1195,505,170,370,540,505,400,
                                     785,1505,455,930,1270,575,140,245,410,380,270,
                                           735,380,155,500,220,885,800,880,850,740,
                                         1100,680,235,935,1600,1520,1600,1570,1460,
                                                   525,865,170,550,470,480,450,340,
                                                      445,360,1030,945,895,865,755,
                                                      700,1365,1285,1365,1335,1225,                                     
                                                               675,590,670,640,530,
                                                                   355,540,510,400,
                                                                       285,225,195,
                                                                            75,155,
                                                                               125};

costs=Table[0,{nn+1},{nn+1}];
k=1;
Do[costs[[1,i+1]]=i;costs[[i+1,1]]=i,{i,nn}];
Do[costs[[i+1,j+1]]=bb[[k]];
   costs[[j+1,i+1]]=bb[[k++]],
{i,nn-1},{j,i+1,nn}];
MatrixForm[costs];

cost=Transpose[Drop[Transpose[Drop[costs,1]],1]];

rangeplot={{-5,81},{-5,95}};

aaa=Graphics[Table[{Point[nodes[[i]]],
      Text[ToString[i],nodes[[i]],{-2,0}]},{i,Length[nodes]}]];

bbb=Graphics[{Dashing[{0.01}],Table[Line[{nodes[[i]],nodes[[j]]}],
    {i,nn-1},{j,i+1,nn}]}];

Show[aaa,AspectRatio->Automatic,PlotRange->rangeplot];

User provided initial tour for manual initialization (list of n nodes)

starttour={2,9,17,19,18,20,5,1,12,3,6,4,13,14,11,10,15,7,8,16}

tourarcs=Table[{starttour[[i]],starttour[[i+1]]},
    {i,Length[starttour]-1}];
    AppendTo[ tourarcs,{starttour[[-1]],starttour[[1]]}]

ccc=show[tourarcs,2,0,"B"];
Show[ccc,aaa,AspectRatio->Automatic,PlotRange->rangeplot];

Manual computation

First a minimal spanning 1-tree is found for the original cost. onetree is the list of
arcs of the 1-tree. tree is the minimal spamming tree over the nodes
{2,..,nn}

Initialization

dual=Table[0,{nn}];

onetreefind[cost]

onetree

tree

treeshow=show[onetree,2,0,"Bk"];

dualvalue=costcompute[onetree]

Show[treeshow,aaa,AspectRatio->Automatic,PlotRange->rangeplot];

subgradient computation

subgr=2-onetreedegree[onetree]

dualvalue=costcompute[onetree]+dual.subgr

Iteration

dualfunctionplot = False;

If[dualfunctionplot,Plot[dualfunction[dual,subgr,t],{t,0,2},PlotRange->All]]

lambda = 0.5;

coeff=lambda (est-dualvalue)/(subgr.subgr)

dual=dual+coeff subgr

newcost=Table[cost[[i,j]]-dual[[i]]-dual[[j]],{i,nn},{j,nn}];

MatrixForm[newcost];

onetreefind[newcost]

subgr=2-onetreedegree[onetree]

dualvalue=costcompute[onetree]+dual.subgr //N

treeshow=show[onetree,2,0,"Bk"];

Show[aaa,treeshow,AspectRatio->Automatic,PlotRange->rangeplot];

Manual adjustment of onetrees into tours

see the arcs of the onetree

tour = Sort[onetree]

identify "wrong" arcs of the onetree and change them

Position[tour, {16, 35}][[1, 1]]

tour[[28]] = {16, 23}

cost

costcompute[tour]

have a look at the modified solution

tourshow=show[tour,2,0,"Bk"];

Show[aaa,tourshow,AspectRatio->Automatic,PlotRange->rangeplot];