"K-means" is the most simple and famous clustering algorithm, which has numerous applications. For a given as an input number of clusters it segments set of points in R^n to that given number of clusters. It minimizes the so-called "inertia" i.e. sum distances^2 to clusters centers = $\sum_{i ~ - ~ cluster~ number} \sum_{X - points~ in ~i-th ~ cluster} |X_{in ~ i-th ~ cluster} - center_{i-th~ cluster} |^2 $

By some reasons let me ask, what happens for the plane i.e. there is no any natural clusters, but still we can pose minimization task and it will produce something. Let us look on the example:

So, most clusters look like hexagons. Especially the most central one which is colored in red. Well, boundary spoils things, also may be not enough sample size/iteration number - simulation is not a perfect thing - but I made many and pictures are similar...
Hexagonal lattice appears in many somewhat related topics, so it might be that some reason exists.

Question 0 What is known on "inertia" minimization on the plane/torus ? (torus - to avoid boundary effects.) (Any references/ideas are welcome). Do hexagons arise as generic clusters ?

Question 1 Consider a torus of sizes R1,R2 , consider the number of clusters to be mn , is it true that hexagonal lattice will provide the global minima for "inertia" ? (At least for consistent values of R1, R2, m,n (R1=am, R2 = a*n) ).

(Instead of finite number of points we can consider the continuous case and substitute summation over points by the integral. Or we can sample large enough uniform datacloud - as done in simulation).

Let me mention beautiful survey by Henry Cohn at ICM2010, where lots of optimization problems of somewhat related spirit are discussed and which sound simple, but remain unsolved for years. That question is not discussed there unfortunately.

The Python code for the simulation above. One can use colab.research.google.com - to run it - no need to install anything - can use google's powers for free.

from sklearn.cluster import  KMeans
import numpy as np
import matplotlib.pyplot as plt
from scipy.spatial.distance import cdist
import time
#import datetime

t0 = time.time()
N = 10**5 # Number of uniformly scattered point 
d = 2 # dimension of space 
X = np.random.rand(N,d) # Generate random uniform N poins on [0,1]^d
n_clusters = 225 # Number of clusters for Kmeans
clustering = KMeans(n_clusters=n_clusters,  
      init='k-means++', n_init=10, max_iter=600, # max_iter increased twice from default  
      tol=0.0001,  random_state=None,  algorithm= 'full' ).fit(X) # Run K-means with default params 
      # https://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans.html#sklearn.cluster.KMeans

print(time.time() - t0, ' secs passed' ) # Print time passed 
 
cluster_centers_ = clustering.cluster_centers_ # 
predicted_clusters = clustering.labels_ #

####################################################################
# Choose the most central classter - hope boundary effect on it would be negligble 
central_point = 0.5 * np.ones(d)  # Choose central pint  
idx_most_central_cluster  = np.argmin( cdist( central_point.reshape(1,-1), cluster_centers_ ) ) # Find cluster most close to central point 
coords_most_central_cluster_center = cluster_centers_[idx_most_central_cluster,: ] 
mask = idx_most_central_cluster  == predicted_clusters 
Xm = X[mask,: ] # Select only points from the most central cluster

#######################################################################
# Plotting 
fig = plt.figure( figsize= (20,10 ) ) # 20 - horizontal size, 6 - vertical size 
plt.scatter( X[:,0], X[:,1], c = predicted_clusters )  # scatter plot all the points  colored according to different clusters
plt.scatter( cluster_centers_[:,0],  cluster_centers_[:,1], c = 'blue' ) # Mark centers of the clusters
plt.scatter( Xm[:,0], Xm[:,1], c = 'red' ) # Color by red most central cluster 
plt.title('n_sample = '+str(N) + ' n_cluster = ' + str(n_clusters))
plt.show()