問題背景：

在求解MTSP問題的時候，因為已知的為各個巡檢點之間路徑耗時長度，而這個具體描述採用無向圖結構可以很好的描述，在matlab中通過函式（graphallshortestpaths）可以得到任意兩個巡檢點之間的距離矩陣

 1 %%得到任意兩個巡檢點之間的路徑時間長度
 2 %W表示從一個巡檢點到另一個巡檢點的路上花費時間
 3 W=[2    1    3    5    1    1    1    4    2    2    1    5    1    2    3    2    6    2    7    2    2    2    6    3    1    2
 
    2    2    1    1    3];
 4 %x表示出發的巡檢點的編號
 5 x=[1    2    2    2    3    3    4    4    5    6    6    6    8    9    9    10    10    11    11    12    13    15    15    16    17    19    20    21    22    23    25];
 6 %y表示到達的巡檢點的編號
 7 y=[2    3    4    19    5    6    21    23    7    8    14    10
 
    17    24    25    11    12    13    15    15    16    18    26    18    25    20    22    22    23    24    26];
 8 A=26;B=31;
 9 m=A;n=A;nzmax=B;
10 DG=sparse(x,y,W,m,n);
11 UG=tril(DG+DG');%得到無向連通圖矩陣
12 D2D=graphallshortestpaths(UG,'directed',false);

。得到距離矩陣以後我們還想根據距離矩陣求出每個巡檢點的座標位置矩陣。也就是說現在問題轉換為這樣的一類問題：我們只知道空間中的點（樣本）的距離，那麼怎麼來重構這些點的相對位置呢？

顯然歐式距離是最直觀的距離，那麼我們就會想使用歐式距離來進行計算重構，我們還希望能夠在不同維度上進行重構，比如2維或者3維。

怎麼做？

有這麼個解決方法叫做MDS 全稱為 Multidimensional Scaling。

MDS解決的問題包括：

當n個物件（object）中各對物件之間的相似性（或距離）給定時，確定這些物件在低維空間中的表示，並使其儘可能與原先的相似性（或距離）“大體匹配”，使得由降維所引起的任何變形達到最小。多維空間中排列的每一個點代表一個物件，因此點間的距離與物件間的相似性高度相關。也就是說，兩個相似的物件由多維空間中兩個距離相近的點表示，而兩個不相似的物件則由多維空間兩個距離較遠的點表示。多維空間通常為二維或三維的歐氏空間，但也可以是非歐氏三維以上空間。

MDS演算法應用的場景包括：

將高維的原始資料投影到低維矩陣上，使得低維矩陣各點之間的相對距離最大程度的接近於原來的高維矩陣傳達的資訊。

在模式識別中，已知一個樣本和另一個樣本在空間的距離，根據空間距離大小進行分類，但是同樣我們只是知道空間上這些點之間的距離，那麼如何根據距離來重構這些點之間的相對位置呢？

演算法的核心思想以及實現思路：

step1：問題重述：

我們有這麼一個距離矩陣，我們通過這個矩陣計算出點的相對位置矩陣X，使得通過X反過來計算距離矩陣與原距離矩陣D差距最小。所以這是一個最優化問題。

step2：

假設X={x1, x2, ..., xn}是一個n×q的矩陣，n為樣本數，q是原始的維度，其中每個xi是矩陣X的一列，xi∈Rq。我們並不知道xi在空間中的具體位置，也就是說對於每個xi，其座標(xi1, xi2, ... , xiq) 都是未知的。我們所知道的僅僅是the pair-wise Euclidean distances for X，我們用一個矩陣DX來表示。因此，對於DX中的每一個元素，可以寫成

或者可以寫成

對於矩陣DX，則有

其中，

這裡的z為

現在讓我們來做平移，從而使得矩陣DX中的點具有zero mean，注意平移操作並不會改變X中各個點的相對關係。為了便於理解，我們先來考察一下AeeT/n和eeTA/n的意義，其中A是一個n×n的方陣。

不難發現AeeT/n中第i行的每個元素都是A中第i行的均值，類似的，我們還可以知道，eeTA/n中第i列的每個元素都是A中第i列的均值。因此，我可以定義centering matrix H如下

所以DXH的作用就是從DX中的每個元素裡減去列均值，HDXH的作用就是在此基礎上再從DX每個元素裡又減去了行均值，因此centering matrix的作用就是把元素分佈的中心平移到座標原點，從而實現zero mean的效果。更重要的是，Let D be a distance matrix, one can transform it to an inner product matrix （Kernel Matrix） by K=-HDH/2, 即

上一步之所以成立，因為

因為BX是一個內積矩陣（Gram Matrix/Kernel Matrix），所以它是對稱的，這樣一來，它就可以被對角化，即BX=U∑UT。

 而我們的最終目的是find a concrete set of n points (Y) in k dimensions so that the pairwise Euclidean distances between all the pairs in the concrete set Y is a close approximation to the pair-wise distances given to us in the matrix DX i.e. we want to find DY such that

Note that after applying the ”double centering” operation to both X and Y，上式服從

最終這個問題的解就是Y=U∑1/2。

下面我們首先來演示在MATLAB中計算上述MDS問題的過程，為了展示其原理，下面的程式碼並不會直接使用MATLAB中內建的用於求解MDS的現成函式。

clear;clc;
%%
%任意兩個城市之間的距離矩陣
D=[[0,2,3,5,4,4,6,6,11,9,11,15,13,5,17,15,7,18,7,9,6,8,9,10,8,11],
    [2,0,1,3,2,2,4,4,9,7,9,13,11,3,15,13,5,16,5,7,4,6,7,8,6,9],
    [3,1,0,4,1,1,3,3,8,6,8,12,10,2,14,12,4,15,6,8,5,7,8,9,5,8],
    [5,3,4,0,5,5,7,7,7,10,12,16,14,6,18,16,8,19,7,5,1,3,4,5,9,12],
    [4,2,1,5,0,2,2,4,9,7,9,13,11,3,15,13,5,16,7,9,6,8,9,10,6,9],
    [4,2,1,5,2,0,4,2,7,5,7,11,9,1,13,11,3,14,7,9,6,8,9,9,4,7],
    [6,4,3,7,2,4,0,6,11,9,11,15,13,5,17,15,7,18,9,11,8,10,11,12,8,11],
    [6,4,3,7,4,2,6,0,5,7,9,13,11,3,11,13,1,13,9,11,8,9,8,7,2,5],
    [11,9,8,7,9,7,11,5,0,12,14,14,16,8,12,17,4,14,8,6,6,4,3,2,3,6],
    [9,7,6,10,7,5,9,7,12,0,2,6,4,6,8,6,8,9,12,14,11,13,14,14,9,12],
    [11,9,8,12,9,7,11,9,14,2,0,8,2,8,7,4,10,7,14,16,13,15,16,16,11,13],
    [15,13,12,16,13,11,15,13,14,6,8,0,9,12,2,7,12,4,18,20,17,18,17,16,11,8],
    [13,11,10,14,11,9,13,11,16,4,2,9,0,10,7,2,12,5,16,18,15,17,18,18,13,13],
    [5,3,2,6,3,1,5,3,8,6,8,12,10,0,14,12,4,15,8,10,7,9,10,10,5,8],
    [17,15,14,18,15,13,17,11,12,8,7,2,7,14,0,5,10,2,20,18,18,16,15,14,9,6],
    [15,13,12,16,13,11,15,13,17,6,4,7,2,12,5,0,14,3,18,20,17,19,20,19,14,11],
    [7,5,4,8,5,3,7,1,4,8,10,12,12,4,10,14,0,12,10,10,9,8,7,6,1,4],
    [18,16,15,19,16,14,18,13,14,9,7,4,5,15,2,3,12,0,21,20,20,18,17,16,11,8],
    [7,5,6,7,7,7,9,9,8,12,14,18,16,8,20,18,10,21,0,2,6,4,5,6,11,14],
    [9,7,8,5,9,9,11,11,6,14,16,20,18,10,18,20,10,20,2,0,4,2,3,4,9,12],
    [6,4,5,1,6,6,8,8,6,11,13,17,15,7,18,17,9,20,6,4,0,2,3,4,9,12],
    [8,6,7,3,8,8,10,9,4,13,15,18,17,9,16,19,8,18,4,2,2,0,1,2,7,10],
    [9,7,8,4,9,9,11,8,3,14,16,17,18,10,15,20,7,17,5,3,3,1,0,1,6,9],
    [10,8,9,5,10,9,12,7,2,14,16,16,18,10,14,19,6,16,6,4,4,2,1,0,5,8],
    [8,6,5,9,6,4,8,2,3,9,11,11,13,5,9,14,1,11,11,9,9,7,6,5,0,3],
    [11,9,8,12,9,7,11,5,6,12,13,8,13,8,6,11,4,8,14,12,12,10,9,8,3,0]];
D([1,22],:)=D([22,1],:);
D(:,[1,22])=D(:,[22,1]);
%%
%求各個城市的座標
DSquare = D.^2;
[N,~]=size(D);
H = eye(N)-ones(N)/N;
K = -0.5*H*DSquare*H;
 
[eigVec, eigVal] = eig(K);
 
Y = eigVec(:,1:2)*sqrt(eigVal(1:2,1:2));%%求解的巡檢點的座標分佈情況。