Two-Point Local Optimality in $k$-Means via Boundary-Point Screening

Wenlong Lyu; Xujie Xiao; Yuheng Jia

cs.LGOct 5, 2026

Two-Point Local Optimality in $k$ -Means via Boundary-Point Screening

Authors: Wenlong Lyu, Xujie Xiao, Yuheng Jia

Organizations: School of Computer Science and Engineering Southeast University Nanjing, Jiangsu 211189, China

Abstract

Lloyd's algorithm and the discrete local (D-local) optimization method (Li et al., 2025) for $k$ -means provide only weak local-optimality guarantees, and their solution quality remains sensitive to initialization. In this paper, we introduce $r$ -point local optimality, under which no reassignment of at most $r$ samples decreases the objective function, and focus on $r=2$ . The main computational obstacle is the $\mathcal{O}(n^2(k^2+d))$ cost of exhaustive two-point certification for $n$ samples in $d$ dimensions and $k$ clusters. To address this challenge, we prove that (i) every improving two-point move of a D-local optimum must involve a cluster shared by both reassignments, and (ii) only certificate-defined boundary points can participate in an improving pair. Exploiting this structure, we propose Boundary-Point-Screened Two-Point Local Search (BPS-2PLS), which terminates at a two-point local optimum. For fixed $k,d$ and nonvanishing cluster occupancy, the number $m$ of retained candidates satisfies $m=\mathcal{O}_{\mathbb{P}}(\log n)$ under i.i.d. sampling from a bounded-support distribution with bounded density or from a Gaussian mixture. Across twelve benchmarks, BPS-2PLS attains the lowest available mean WCSS on ten. In a subsampling study, screening retains 0.10% to 2.81% of samples on average at the largest tested sizes. The code is available at https://github.com/lwl-learning/BPS-2PLS.

Figures & tables

Configuration $\tau$	Partner constraint	Destination and feasibility	Exact gain $G(u,v)$
Same source, same destination	$j\in{\mathbb{I}}_{s}\setminus\{i\}$	$b=a$ , $n_{s}>2$	$R_{s}(i,j)-A_{a}(i,j)$
Same source, different destinations	$j\in{\mathbb{I}}_{s}\setminus\{i\}$	$b\notin\{s,a\}$ , $n_{s}>2$	$R_{s}(i,j)-A_{a}(i)-A_{b}(j)$
Forward chain	$j\in{\mathbb{I}}_{a}$	$b\notin\{s,a\}$ , $n_{s}>1$	$R_{s}(i)+T_{a}(j,i)-A_{b}(j)$
Swap	$j\in{\mathbb{I}}_{a}$	$b=s$	$T_{s}(i,j)+T_{a}(j,i)$
Different sources, same destination	$j\notin{\mathbb{I}}_{s}\cup{\mathbb{I}}_{a}$	$b=a$ , $n_{s}>1$ , $n_{z_{j}}>1$	$R_{s}(i)+R_{z_{j}}(j)-A_{a}(i,j)$
Backward chain	$j\notin{\mathbb{I}}_{s}\cup{\mathbb{I}}_{a}$	$b=s$ , $n_{z_{j}}>1$	$R_{z_{j}}(j)+T_{s}(i,j)-A_{a}(i)$

Table 1: Search spaces for the partner

v=(j:z_{j}\to b)

and exact gains for fixed

u=(i:s\to a)

Configuration	Pivot ${\bm{o}}_{e}$	Key $e$	Partner superset $\widetilde{{\mathbb{J}}}_{e}$
Same source (both configurations)	${\bm{\mu}}_{s}$	$(\tau,s,a)$	${\mathbb{I}}_{s}$
Forward chain / swap	$({\bm{\mu}}_{s}+{\bm{\mu}}_{a})/2$	$(\tau,s,a)$	${\mathbb{I}}_{a}$
Different sources, same destination	${\bm{\mu}}_{a}$	$(\tau,a)$	$\{j:z_{j}\neq a,\ n_{z_{j}}>1\}$
Backward chain	${\bm{\mu}}_{s}$	$(\tau,s)$	$\{j:z_{j}\neq s,\ n_{z_{j}}>1\}$

Table 2: Reusable envelope families for

u=(i:s\to a)

, with

s\neq a

Panel of Figure 1 — Figure 1: Schematic pivot choices for the six connected configurations: (a) both same-source configurations; (b) forward chains and swaps; (c) different sources with the same destination; and (d) backward chains. Shaded regions denote clusters.

Dataset	Metric	Lloyd	lite-KM	Hartigan	D-LO	Min-D-LO	BPS-2PLS
Yale	WCSS	3328.242	3328.242	3197.414	3214.216	3214.216	3185.016
	Time	0.011	0.009	0.356	0.219	0.246	0.327
ORL	WCSS	2902.673	2902.673	2723.878	2735.73	2735.73	2691.304
	Time	0.025	0.018	1.481	0.935	1.055	2.403
USPS-1000	WCSS	22400.15	22400.15	22234.94	22254.42	22254.42	22229.35
	Time	0.057	0.020	4.910	1.035	1.193	1.280

Table 3: Mean WCSS and runtime (seconds) over ten shared initializations. Bold: lowest available WCSS. Dashes mark Hartigan runs omitted because of excessive runtime.

Figure 2: BPS-2PLS descent trajectories on five datasets. The

x

-axis shows the iteration count after Lloyd convergence, and the

y

-axis shows WCSS as a percentage of the WCSS at Lloyd convergence.

Method	Variant	Yale		ORL		USPS-1000		COIL20		MNIST-2000
		WCSS	Time	WCSS	Time	WCSS	Time	WCSS	Time	WCSS	Time
Min-D-LO	Baseline	3214.216	0.246	2735.73	1.055	22254.42	1.193	32783.05	1.385	76501.68	2.345
	Inner	3185.016	0.327	2691.304	2.403	22229.35	1.280	32772.27	1.253	76499.59	2.521
LocalSearch	Baseline	5204.958	0.074	4594.401	0.196	35930.89	0.106	51511.93	0.439	123605.5	0.352
	Post	3180.593	0.422	2670.928	2.439	22235.22	1.003	31799.29	2.329	76382.11	3.156
MultiSwap	Baseline	3347.979	0.045	2839.742	0.178	22363.17	0.090	33014.32	0.356	76669.54	0.488

Table 4: Heuristic enhancement: mean WCSS and runtime (seconds). Bold: lowest WCSS within each method and dataset, including ties; dashes: omitted runs.

Figure 3: Retained candidate count

m

versus sample size

n

at the first D-local state. Markers show means over five random subsets, each averaged over three initializations; error bars indicate one standard deviation across subsets.

Appendix figures & tables9 assets

Supplementary material from the paper’s appendix.

Appendix

	Full sample-pair cache	Lazy candidate-pair blocks
Pair-distance preprocessing	$\mathcal{O}(n^{2}d)$	None
Screening and search per round	$\mathcal{O}(k^{2}d+nk\log n+m^{2}k)$	$\mathcal{O}(k^{2}d+nk\log n+m^{2}(d+k))$
Peak memory	$\mathcal{O}(nd+nk+k^{2}+n^{2})$	$\mathcal{O}(nd+nk+k^{2}+wm)$

Appendix

Table 5: Distance-policy costs for BPS-2PLS. Screening-and-search time excludes one-point descent and accepted-move state updates; peak memory includes the input data and clustering workspace. The block row limit is

w

Configuration $\tau$	$C_{\tau}(u)$	$Q_{e}(j)$	$\eta_{e}$
Same source same destination	$\frac{n_{s}}{n_{s}-2}D_{is}-\frac{n_{a}}{n_{a}+2}D_{ia}$	$\frac{n_{s}}{n_{s}-2}D_{js}-\frac{n_{a}}{n_{a}+2}D_{ja}$	$-\frac{1}{n_{s}-2}-\frac{1}{n_{a}+2}$
Same source, different destinations	$\frac{n_{s}}{n_{s}-2}D_{is}-A_{a}(i)$	$\frac{n_{s}}{n_{s}-2}D_{js}-A_{b^{\star}(j;s,a)}(j)$	$-\frac{1}{n_{s}-2}$
Forward chain	$R_{s}(i)-D_{ia}$	$D_{ja}-A_{b^{\star}(j;s,a)}(j)$	$\frac{1}{n_{a}}$
Swap	$D_{is}-D_{ia}$	$D_{ja}-D_{js}$	$\frac{1}{n_{s}}+\frac{1}{n_{a}}$
Different sources, same destination	$R_{s}(i)-\frac{n_{a}}{n_{a}+2}D_{ia}$	$R_{z_{j}}(j)-\frac{n_{a}}{n_{a}+2}D_{ja}$	$-\frac{1}{n_{a}+2}$
Backward chain	$D_{is}-A_{a}(i)$	$R_{z_{j}}(j)-D_{js}$	$\frac{1}{n_{s}}$

Appendix

Table 6: Separable gains for

u=(i:s\to a)

, with

s=z_{i}

and envelope key

e=e(\tau,u)

Configuration	Interaction $I(u,v)$	Excess upper bound
Same source, same destination	$I_{s}^{-}(i,j)+I_{a}^{+}(i,j)$	$(16/n_{s}+6/n_{a})p_{\max}$
Same source, different destinations	$I_{s}^{-}(i,j)$	$16p_{\max}/n_{s}$
Forward chain	$J_{a}(j,i)$	$5p_{\max}/n_{a}$
Swap	$J_{s}(i,j)+J_{a}(j,i)$	$5(1/n_{s}+1/n_{a})p_{\max}$
Different sources, same destination	$I_{a}^{+}(i,j)$	$6p_{\max}/n_{a}$
Backward chain	$J_{s}(i,j)$	$5p_{\max}/n_{s}$

Appendix

Table 7: Interactions and their upper bounds after radial relaxation, for the anchor

i:s\to a

. The last column bounds the excess over the two one-point gains.

Dataset	$n$	$d$	$k$	Representation
Yale †	165	1,024	15	Pixels
ORL †	400	1,024	40	Pixels
USPS-1000 †	1,000	256	10	Pixels
COIL20 †	1,440	1,024	20	Pixels
MNIST-2000 †	2,000	784	10	Pixels
COIL100	7,200	1,024	100	Pixels

Appendix

Table 8: Datasets used in the clustering-quality comparison. Here

n

d

, and

k

denote the number of samples, feature dimension, and number of reference classes, respectively. Dataset names link to benchmark source pages; the listed dimensions describe the representations used here. Datasets marked with

\dagger

are used for the descent-trajectory study and heuristic comparisons.

Dataset	Lloyd	lite-KM	Hartigan	D-LO	Min-D-LO	BPS-2PLS
Yale	$38.67\pm 2.35$	$38.67\pm 2.35$	$41.39\pm 3.82$	$41.03\pm 2.46$	$41.03\pm 2.46$	$\mathbf{41.64}\pm 2.16$
ORL	$48.85\pm 4.12$	$48.85\pm 4.12$	$54.73\pm 3.09$	$53.88\pm 4.36$	$53.88\pm 4.36$	$\mathbf{54.88}\pm 3.09$
USPS-1000	$45.69\pm 3.06$	$45.69\pm 3.06$	$\mathbf{47.67}\pm 2.34$	$47.59\pm 2.74$	$47.59\pm 2.74$	$47.48\pm 2.54$
COIL20	$60.76\pm 5.20$	$60.76\pm 5.20$	$61.69\pm 4.33$	$61.88\pm 5.64$	$61.88\pm 5.64$	$\mathbf{62.03}\pm 5.46$
MNIST-2000	$54.78\pm 3.31$	$54.78\pm 3.31$	$54.29\pm 3.50$	$\mathbf{56.12}\pm 4.21$	$\mathbf{56.12}\pm 4.21$	$56.07\pm 4.13$
COIL100	$46.11\pm 1.68$	$46.11\pm 1.68$	—	$46.72\pm 1.64$	$46.72\pm 1.64$	$\mathbf{46.99}\pm 1.80$

Appendix

Table 9: ACC (%) for the sample-reassignment methods over ten shared initializations. Entries are mean

\pm

standard deviation; higher is better. Bold marks the highest available mean on each dataset, including ties at the displayed precision. Dashes mark Hartigan runs omitted because of excessive runtime.

Dataset	Lloyd	lite-KM	Hartigan	D-LO	Min-D-LO	BPS-2PLS
Yale	$47.06\pm 2.14$	$47.06\pm 2.14$	$50.32\pm 2.52$	$49.95\pm 1.84$	$49.95\pm 1.84$	$\mathbf{50.56}\pm 2.74$
ORL	$72.92\pm 2.10$	$72.92\pm 2.10$	$77.08\pm 1.23$	$76.63\pm 1.92$	$76.63\pm 1.92$	$\mathbf{77.29}\pm 1.87$
USPS-1000	$44.10\pm 1.22$	$44.10\pm 1.22$	$\mathbf{46.94}\pm 2.76$	$46.12\pm 2.30$	$46.12\pm 2.30$	$46.06\pm 2.00$
COIL20	$76.51\pm 1.65$	$76.51\pm 1.65$	$76.77\pm 1.62$	$77.28\pm 2.11$	$77.28\pm 2.11$	$\mathbf{77.38}\pm 2.04$
MNIST-2000	$48.75\pm 2.65$	$48.75\pm 2.65$	$49.01\pm 2.76$	$\mathbf{49.49}\pm 2.80$	$\mathbf{49.49}\pm 2.80$	$49.45\pm 2.76$
COIL100	$76.26\pm 0.34$	$76.26\pm 0.34$	—	$76.74\pm 0.36$	$76.74\pm 0.36$	$\mathbf{76.92}\pm 0.38$

Appendix

Table 10: NMI (%) for the sample-reassignment methods over ten shared initializations. Entries are mean

\pm

standard deviation; higher is better. Bold marks the highest available mean on each dataset, including ties at the displayed precision. Dashes mark Hartigan runs omitted because of excessive runtime.

Method	Variant	Yale	ORL	USPS-1000	COIL20	MNIST-2000
Min-D-LO	Baseline	$41.03\pm 2.46$	$53.88\pm 4.36$	$\mathbf{47.59}\pm 2.74$	$61.88\pm 5.64$	$\mathbf{56.12}\pm 4.21$
	Inner	$\mathbf{41.64}\pm 2.16$	$\mathbf{54.88}\pm 3.09$	$47.48\pm 2.54$	$\mathbf{62.03}\pm 5.46$	$56.07\pm 4.13$
LocalSearch	Baseline	$38.61\pm 3.16$	$53.00\pm 2.39$	$37.59\pm 2.64$	$57.72\pm 3.64$	$45.31\pm 3.42$
	Post	$\mathbf{43.15}\pm 4.33$	$\mathbf{57.55}\pm 2.79$	$\mathbf{47.86}\pm 2.07$	$\mathbf{63.01}\pm 4.79$	$\mathbf{53.05}\pm 3.99$
MultiSwap	Baseline	$36.97\pm 3.64$	$52.10\pm 0.78$	$46.14\pm 2.92$	$58.34\pm 4.89$	$\mathbf{54.65}\pm 5.78$
	Post	$40.61\pm 4.13$	$\mathbf{56.33}\pm 1.63$	$\mathbf{48.12}\pm 2.92$	$\mathbf{59.63}\pm 4.25$	$54.62\pm 5.33$

Appendix

Table 11: ACC (%) for heuristic enhancement: mean

\pm

standard deviation over ten initializations. Bold marks the highest mean within each method and dataset, including ties at the displayed precision. Dashes mark J-Means runs omitted because of excessive runtime. For Min-D-LO, Inner reports BPS-2PLS.

Method	Variant	Yale	ORL	USPS-1000	COIL20	MNIST-2000
Min-D-LO	Baseline	$49.95\pm 1.84$	$76.63\pm 1.92$	$\mathbf{46.12}\pm 2.30$	$77.28\pm 2.11$	$\mathbf{49.49}\pm 2.80$
	Inner	$\mathbf{50.56}\pm 2.74$	$\mathbf{77.29}\pm 1.87$	$46.06\pm 2.00$	$\mathbf{77.38}\pm 2.04$	$49.45\pm 2.76$
LocalSearch	Baseline	$45.08\pm 1.87$	$72.18\pm 1.06$	$33.05\pm 1.68$	$70.72\pm 1.04$	$35.68\pm 3.67$
	Post	$\mathbf{50.86}\pm 2.88$	$\mathbf{77.83}\pm 1.26$	$\mathbf{46.60}\pm 1.99$	$\mathbf{77.46}\pm 1.10$	$\mathbf{47.72}\pm 2.19$
MultiSwap	Baseline	$45.47\pm 2.79$	$74.48\pm 0.62$	$44.44\pm 2.65$	$75.72\pm 1.50$	$49.44\pm 2.84$
	Post	$49.44\pm 2.87$	$\mathbf{77.68}\pm 0.90$	$\mathbf{46.65}\pm 2.67$	$\mathbf{76.61}\pm 1.19$	$\mathbf{49.79}\pm 2.52$

Appendix

Table 12: NMI (%) for heuristic enhancement: mean

\pm

Dataset	Total	Descent	Screen	Search	Other
Yale	$0.268\pm 0.013$	$\mathbf{0.149}\pm 0.025$ $(55.5\%)$	$0.077\pm 0.016$ $(28.9\%)$	$0.036\pm 0.015$ $(13.5\%)$	$0.006\pm 0.002$ $(2.1\%)$
ORL	$2.090\pm 0.538$	$0.695\pm 0.091$ $(33.3\%)$	$\mathbf{1.269}\pm 0.442$ $(60.7\%)$	$0.100\pm 0.054$ $(4.8\%)$	$0.026\pm 0.010$ $(1.2\%)$
USPS-1000	$1.032\pm 0.582$	$\mathbf{0.923}\pm 0.582$ $(89.4\%)$	$0.080\pm 0.043$ $(7.7\%)$	$0.021\pm 0.011$ $(2.0\%)$	$0.009\pm 0.005$ $(0.9\%)$
COIL20	$1.256\pm 0.546$	$\mathbf{1.055}\pm 0.584$ $(84.0\%)$	$0.176\pm 0.062$ $(14.0\%)$	$0.005\pm 0.001$ $(0.4\%)$	$0.020\pm 0.008$ $(1.6\%)$
MNIST-2000	$1.681\pm 0.565$	$\mathbf{1.590}\pm 0.544$ $(94.6\%)$	$0.059\pm 0.025$ $(3.5\%)$	$0.015\pm 0.006$ $(0.9\%)$	$0.018\pm 0.010$ $(1.1\%)$
COIL100	$76.852\pm 6.224$	$\mathbf{40.142}\pm 3.222$ $(52.2\%)$	$34.024\pm 4.556$ $(44.3\%)$	$2.027\pm 0.842$ $(2.6\%)$	$0.660\pm 0.087$ $(0.9\%)$

Appendix

Table 13: Runtime by phase with screening enabled. Times are in seconds, reported as mean

\pm

standard deviation across three initializations after averaging five repetitions within each initialization. Parentheses give each phase’s percentage of the mean total time; percentages below

0.1\%

are shown to two decimal places. Bold marks the largest mean phase time for each dataset.

Explore similar work

Jul 28, 2026stat.ML

Lloyd's $K$ -Means Clustering Algorithm Is Frank-Wolfe in Disguise

Lloyd's

K

-means algorithm, also known as naïve

K

-means, is a widely used ad hoc optimization heuristic, designed to minimize the sum of squared errors (SSE) across all

K

-partitions of a dataset via iterative cluster refinement. In this work, we establish a novel connection between Lloyd's algorithm and the Frank-Wolfe (FW) algorithm, a prominent first-order method for projection-free optimization. We demonstrate that Lloyd's algorithm is a special case of FW. Leveraging recent advances in FW methods for concave objectives, we derive a non-asymptotic

\mathcal{O}(1/t)

convergence rate to a local minimum of the SSE objective. To account for empty clusters, an outcome possible under Lloyd's greedy assignment, we develop an FW variant for semismooth objectives while retaining the same convergence rate that is solely controlled by the initial SSE value. We illustrate our findings with a simulation study for spherical Gaussian mixtures and a real-world image segmentation dataset.

Michael Pokojovy, J. Marcus Jobe, Simon Lacoste-Julien

Department of Mathematics & Statistics and School of Data Science Old Dominion University Norfolk, VA 23529, USA · Information Systems & Analytics Department Miami University Oxford, OH 45056, USA · Canada CIFAR AI Chair Department of Computer Science and Operations Research & Mila Université de Montréal Montréal (QC) Canada

K-Means Convex Optimization Clustering Optimal Sample Complexity Gaussian Mixture Models

Ask PDF

Jul 17, 2026cs.LG

Data-Native Global Optimization for Big Data K-means Clustering

Big data clustering remains challenging: the Minimum Sum-of-Squares Clustering (MSSC) problem underlying K-means is NP-hard, and existing methods either reach poor local minima or require prohibitive metaheuristic hybrids. We target arbitrarily tall data: a fixed feature space may contain arbitrarily many, possibly infinitely many, observations, while the algorithm accesses only finite random samples. We propose Big-means++, an algorithm achieving scalability and global-search quality by curating inputs to MSSC optimization on big data. It orchestrates local K-means refinements into a data-native global search for big data clustering. Rather than optimizing the full-data MSSC objective, Big-means++ traverses sample-induced surrogate landscapes. Each sample defines a distinct empirical MSSC approximation with a perturbed local-optimum structure, turning sample-to-sample variation into a global-search mechanism. Unlike Big-means, a flowing-incumbent strategy propagates centroid state across empirical landscapes through K-means refinements on fresh samples without rollback to a best-so-far solution. This increases mobility and favors stable, high-quality configurations across approximations of the full-data structure. A new shaking mechanism varies sample size geometrically, broadening the surrogate landscapes explored across resolution scales, accounting for cluster imbalance, and improving solution quality. A competitive multi-agent system asynchronously explores independent sampled landscapes, transforming diverse stochastic trajectories into collective search intelligence. Automatic convergence detection stops each agent after attaining a high-quality solution but before further search risks degrading it, while providing a universal speed-quality control. Experiments on 22 datasets against 11 competing algorithms demonstrate the effectiveness, efficiency, and robustness of Big-means++.

Ravil Mussabayev, Rustam Mussabayev, Zukhra Yerdaliyeva +1

AI Research Lab, Satbayev University, Almaty, 050000, Kazakhstan · Laboratory for Analysis and Modeling of Information Processes, Institute of Information and Computational Technologies, Almaty, 050010, Kazakhstan

K-Means Clustering Data Lakes Maximum Mean Discrepancy Unsupervised

Ask PDF

Dec 5, 2025stat.ML

BalLOT: Balanced $k$ -means clustering with optimal transport

We consider the fundamental problem of balanced

k

-means clustering. In particular, we introduce an optimal transport approach to alternating minimization called BalLOT, and we show that it delivers a fast and effective solution to this problem. We establish this with several theoretical guarantees and a variety of numerical experiments. On the theory front, we first prove that for generic data, BalLOT produces integral couplings at each step. Next, we perform a landscape analysis to provide theoretical guarantees for both exact and partial recoveries of planted clusters under the stochastic ball model. We also propose initialization schemes that achieve one-step recovery of planted clusters. To conclude, we present numerical experiments that corroborate our theoretical results.

Wenyan Luo, Dustin G. Mixon

Department of Mathematics, The Ohio State University, Columbus, Ohio, USA · Translational Data Analytics Institute, The Ohio State University, Columbus, Ohio, USA

K-Means Optimal Transport Approach Cluster Majority Voting

Ask PDF

Abstract

Explore similar work

Lloyd's KKK-Means Clustering Algorithm Is Frank-Wolfe in Disguise

Data-Native Global Optimization for Big Data K-means Clustering

BalLOT: Balanced kkk-means clustering with optimal transport

Lloyd's $K$ -Means Clustering Algorithm Is Frank-Wolfe in Disguise

BalLOT: Balanced $k$ -means clustering with optimal transport