MInimal Distence index for ellipsoidal clusters

Input: partition Π={(Subscript[π, 1]),⃛,(Subscript[π, k])} of the set 𝒜 ⊂ (ℝ^n), centers (Subscript[c, 1]),⃛,(Subscript[c, k]) ∈ (ℝ^n), 
covariance matrices (∑_j)∈(ℝ^(n×n)), j=1⃛k;
Output: MMD -index;
Submodules:  IntersectTest[],MinR[]

Clear[MMD]
MMD[PI_,z_,Scov_,Ind2_:0]:=Module[{k=Length[PI],min,P1,P2},
DDD=Table[0,{r,k},{s,k}];
PI0=Table[PI[[j,All,2]],{j,k}];
Ell=Table[Ellipsoid[z[[j]],kappa  Scov[[j]]],{j,k}];
	 (* M-udaljenosti (najprije se testira je li presjek prazan) *)
Do[
Do[
If[IntersectTest[Ell[[r,1]],Ell[[r,2]],Ell[[s,1]],Ell[[s,2]]],
min=FindMinimum[{EuclideanDistance[x,y],{x∈Ell[[r]],y∈Ell[[s]]}}, {x,y}];If[Ind2≠ 0,Print[min];{P1,P2}={x/.min[[2]],y/.min[[2]]};
		];
DDD[[r,s]]=min[[1]],
{P1,P2}={Ell[[r,1]],Ell[[r,1]]};DDD[[r,s]]=0
	];
If[Ind2≠ 0, 
If[n==2,
Print[Show[
Graphics[{boja[[r]],EdgeForm[Directive[Thickness[.004],boja[[r]]]],Opacity[0.075],Ell[[r]]}],
Graphics[{boja[[s]],EdgeForm[Directive[Thickness[.004],boja[[s]]]],Opacity[0.075],Ell[[s]]}],
ListPlot[PI0,PlotStyle->Table[{PointSize[.008],boja[[j]]},{j,Length[sol[[2]]]}]], 
ListPlot[{P1,P2}, PlotStyle->{Red,PointSize[.025]}], 
Graphics[{Red,Thin,Line[{P1,P2}]}],
Axes->True, AxesOrigin->{0,0},AxesStyle->{Thickness[.005],Thickness[.005]},GridLines->{tabx,taby},PlotRange->AG,AspectRatio->Automatic,ImageSize->250]],
Print[Show[
Graphics3D[{RGBColor[z[[r]]],Opacity[0.5],Ell[[r]]}],Graphics3D[{RGBColor[z[[s]]],Opacity[0.5],Ell[[s]]}],
ListPointPlot3D[PI0,PlotStyle->Table[{PointSize[.001],Opacity[.5],RGBColor[z[[j]]]},{j,Length[sol[[2]]]}]],
ListPointPlot3D[{P1,P2}, PlotStyle->{Red,PointSize[Large]}],
Graphics3D[{Red,Thin,Line[{P1,P2}]}],
Axes->True,AxesStyle->{Red,Green,Blue},AxesLabel->{"R","G","B"},ImageSize->250]]
		];
			];
,{s,r+1,k}]
,{r,k}];
Print["DDD=",MatrixForm[DDD]," - distance matrix between clusters"];
m=Sum[Length[PI[[j]]],{j,Length[PI]}];
mwf=MWF[PI,z,Scov];
(m-k)MinR[DDD] /mwf
]

Example

Mahalanobis distance-like function and Normalized Mahalanobis distance-like function

Input: data points x,y ∈ (ℝ^n), symmetric positive definite matrix ∑ ∈ (ℝ^(n×n)); 
Output dm: Mahalanobis distance between the points x,y ∈ (ℝ^n);
Output dM: normalized Mahalanobis distance between the points x,y ∈ (ℝ^n);
Submodules: none needed

dm[x_,y_,S_]:=(x-y).Inverse[S].(x-y);
dM[x_,y_,S_]:=(Det[S]^(1/n))(x-y).Inverse[S].(x-y);

Minimal element of the upper triangular matrix

Input: The upper triangular matrix ℝ;
Output: Minimal element;
Submodules: None

MinR[R_]:=Module[{n=Length[R[[1]]],min=Infinity},
Do[
red=Drop[R[[i]],{1,i}];
minn=Min[red];
If[minn<min,min=minn];
,{i,n-1}];
min]

Testing if two ellipsoids is not intersecting

Input: Ell1[c1,cov1],Ell2[c2,cov2];
Output: True (not intersecting)or False (intersecting)
Submodules: None

IntersectTest[c1_,cov1_,c2_,cov2_]:=Module[{FF,max},
FF[lambda_]:=lambda (1-lambda) (c2-c1).Inverse[lambda cov2+(1-lambda) cov1].(c2-c1);
max=First[FindMaximum[{FF[x],0<x<1},x]];
max>1
]

Mahalanobis weighted 𝓀 - means algorithm

Input: data set 𝒜 ⊂ (ℝ^n), data weights w ∈ (ℝ^n), initial centers (Subscript[z, 1]),⃛,(Subscript[z, k]) ∈ (ℝ^n), ;
Output: partition Π={(Subscript[π, 1]),⃛,(Subscript[π, k])} (each partition consists of clus-ters with corresponding weights associated 
to their elements), centers (Subscript[c, 1]),⃛,(Subscript[c, k]) ∈ (ℝ^n), covariance matrices (∑_j)∈(ℝ^(n×n)), j=1⃛k, number of iterations It;
Submodules: dM[],  DIRECT[]

(* It is used dM *)
MWKMeans[w_,A_,z_,Ind_:0]:=Module[{m=Length[w],k=Length[z],W=Total[w],PI0,PIpriv,tab,min,imin,Pod,c,c0,n=Length[z[[1]]],Scov,pom,corr,It=0,epsilon=.000000001},
Pod=Table[{w[[i]],A[[i]]},{i,m}];
c0=z;
						(* Partition calculating *)
PI=Table[{},{j,k}];Scov=Table[0,{j,k}];pom=Table[0,{j,k}];mj=Table[0,{j,k}];
Do[
tab=Table[Norm[Pod[[i]][[2]]-z[[j]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
PI=DeleteCases[PI,{}];
k=Length[PI];
                                                           (* Centers calculating *)
c=Table[Sum[PI[[j,s,1]]*PI[[j,s,2]],{s,Length[PI[[j]]]}]/Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}],{j,k}];
							(* LOOP *)
While[
Chop[Norm[c-c0]]≠0,c0=c;It=It+1;
(* It<1,c0=c;It=It+1;*)
				(* Covariance matrices Scov calculating *)
Do[
mj[[j]]=Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}];
Scov[[j]]=Sum[PI[[j,s,1]]*KroneckerProduct[(PI[[j,All,2]][[s]]-c[[j]]),(PI[[j,All,2]][[s]]-c[[j]])],{s,Length[PI[[j]]]}]/mj[[j]]
,{j,k}];
Do[
If[Min[Eigenvalues[Scov[[j]]]]<epsilon,Print["Cov.Mat. ",j,"-tog centra (",c[[j]],"): ",Scov[[j]]," je singularna"];Print["|PI[j]|: ",Length[PI[[j]]],", Eigen_min=",Min[Eigenvalues[Scov[[j]]]]];Scov[[j]]=IdentityMatrix[n] ]
,{j,k}];
PIpriv=PI;
PI=Table[{},{i,k}];
corr=Table[Log[m^2 Det[Scov[[jj]]]/mj[[jj]]^2],{jj,k}];
Do[
Do[
tab=Table[dM[PIpriv[[jj,s,2]],c[[j]],Scov[[j]]],{j,k}];
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],PIpriv[[jj,s]] ],
{s,Length[PIpriv[[jj]]]}],
{jj,k}];
PI=DeleteCases[PI,{}];
k=Length[PI];
c=Table[Sum[PI[[j,s,1]]*PI[[j,s,2]],{s,Length[PI[[j]]]}]/Sum[PI[[j,s,1]],{s,Length[PI[[j]]]}],{j,k}];
]; 
				(* The case of reducing the number of clusters in a partition *)
If[Length[c]<Length[c0],
F2[u_,v_]:=Sum[w[[i]]Min[Append[Table[Norm[A[[i]]-c[[s]]]^2,{s,Length[c]}],Norm[A[[i]]-{u,v}]^2]],{i,m}];
dom={xG,yG};
								(* The new DIRECT center *)
minp=DIRECT[F2,dom,DMin->10.^(-3),MaxIter->50, Eps->10^(-4)]//Timing; min=minp[[2]];Print["DIRECT: ",min];
c=Append[c,min[[1]]];
Print["Novi optimirani centar: ",c];
k=Length[c]; Print["Korigirani k: ",k];
								(* The new partition *)
PI=Table[{},{j,k}];
Do[
tab=Table[Norm[Pod[[i]][[2]]-c[[j]]],{j,k}]//N;
imin=Ordering[tab,1][[1]];
PI[[imin]]=Append[PI[[imin]],Pod[[i]]],
{i,m}];
];
{PI,c,Scov,It}
]

Execution

SeedRandom[7]; 
tabx=taby=Table[i,{i,20}];
boja={Blue,Green,Orange,Brown,Purple,Magenta}
m0=300; k=4; m=k m0; n=2;      
p=.7;  kappa=Quantile[ChiSquareDistribution[n],p];
Print["p = ",p,";  kappa = ",kappa];
CC={{2.5,2.5},{6,2.5},     {3.,5.5},{7.5,7.5}};
AG={xG,yG}={{0,10},{0,10}};
slC=ListPlot[CC,PlotMarkers->{"✶",15},PlotStyle->{Gray,PointSize[.025],Opacity[.9]}];
cov={.7{{2,1},{1,1}}, .3{{4,1.5},{1.5,1.4}},  .3{{1.6,2},{2,3.6}},1.3{{.8,.4},{.4,.8}}};
PI0=Table[RandomReal[MultinormalDistribution[CC[[j]],cov[[j]]],m0],{j,k}];
A=Flatten[PI0,1];
w=Table[1,{i,Length[A]}]; 
z0={{1,1},{5,2},{3,6},{8,8}}; 
sol=MWKMeans[w,A,z0,0];
				(*  *)
MMD[sol[[1]],sol[[2]],sol[[3]],1]