k-Means Clustering

Partition observations into k clusters by minimizing the within-cluster squared distance, a nonconvex model that, unlike the classical algorithm, also accepts must-link, cannot-link and cluster-size constraints.

Notebook

Python / Knitro API Python / Pyomo Julia / JuMP

Introduction

\(k\)-means clustering is a classical problem in data science. It consists in finding a partition of a given set of observations into a given number of clusters.

\(k\)-means clustering being a classical problem in data science, algorithms to solve it are already available, and they are fast. The advantage of nonlinear programming shows when additional constraints need to be taken into account. These extensions are usually not supported by the algorithms available in the data science libraries. Here we illustrate three of them and how to model them with nonlinear programming:

• Force some pairs of observations to be in the same cluster
• Force some pairs of observations to be in different clusters
• Impose a minimum and maximum size for clusters

Problem description

Inputs

• \(N\) observations with their coordinates \(X_i\), \(i = 1, \dots, N\) in the considered space
• A number of clusters to build \(K\)

Problem: partition the \(N\) observations into \(K\) clusters such that each observation belongs to exactly one cluster.

Objective: minimize the total sum of Euclidean distances between each observation and the center of its cluster

\[ \min \sum_{k=1}^K \sum_{i \in \mathrm{cluster}_k} \| X_i - c_k \|^2_2 \]

with \(c_k\) the coordinates of the center of cluster \(k\).

Mathematical model

Now, we show how to model the \(k\)-means clustering problem using nonlinear programming.

Variables

• \(a^k_i \in \{0, 1\}\), \(i = 1, \dots, N\), \(k = 1, \dots, K\): a boolean variable which is 1 if the observation \(i\) is in cluster \(k\), else 0
• \(d^k_i \in \mathbb{R}^+\), \(i = 1, \dots, N\), \(k = 1, \dots, K\): a positive variable which represents the Euclidean squared distance between observation \(i\) and cluster \(k\)
• \(c^k \in \mathbb{R}^2\), \(k = 1, \dots, K\): position of the barycenter of the cluster \(k\)

Objective

\[ \min \sum_{i=1}^N \sum_{k=1}^K a^k_i d^k_i \]

Constraints

• Each observation has to be in one cluster

\[ \forall i = 1, \dots, N \qquad \sum_{k=1}^K a^k_i = 1 \]

• Coordinates of the barycenter of the clusters

\[ \forall k = 1, \dots, K \qquad c_k \sum_{i=1}^N a^k_i = \sum_{i=1}^N a^k_i X_i \]

• Square distance formula

\[ \forall i = 1, \dots, N, \quad \forall k = 1, \dots, K \qquad d^k_i = \| X_i - c_k \|^2_2 \]

The products \(a^k_i d^k_i\) and \(c_k a^k_i\) make this a nonconvex mixed-integer nonlinear program.

Model extensions

In some cases, additional constraints must be taken into account. We add the following sets of constraints to the model above.

Cannot-link and must-link observations. In some cases, some pairs of observations might be known to belong to the same cluster or to different clusters. Let \(P^\text{must}\) be the set of sets of observations that must belong to the same cluster, and \(P^\text{cannot}\) the set of pairs of observations that must belong to different clusters:

\[ \forall k = 1, \dots, K, \quad \forall S \in P^\text{must}, \quad \forall i \in S \qquad a^k_i = a^k_{S[0]} \]

\[ \forall k = 1, \dots, K, \quad \forall i_1, i_2 \in P^\text{cannot} \qquad a^k_{i_1} + a^k_{i_2} \leq 1 \]

Constrained clustering. Another common case in clustering is to have lower and/or upper bounds on the cluster sizes. Let \(u_\text{limit}\) be the maximum number of observations in a cluster and \(l_\text{limit}\) the minimum:

\[ \forall k = 1, \dots, K \qquad l_\text{limit} \leq \sum_{i=1}^N a^k_i \leq u_\text{limit} \]

Input data

An instance is the coordinates of the observations and the number of clusters to build. Clusters carries total_distance, the total squared distance the solve reached, along with the cluster label of every observation and the cluster centers.

struct Clusters
    total_distance::Float64
    labels::Vector{Int}
    centers::Vector{Tuple{Float64,Float64}}
end

cluster(clusters::Clusters, k) = findall(==(k), clusters.labels)

The observations are a CSV of X and Y columns under shared/data. link_groups picks the must-link and cannot-link sets out of the positions themselves rather than naming observations by index.

const DATA_DIR = joinpath("shared", "data")

function load_observations(name)
    rows = Tuple{Float64,Float64}[]
    for (index, line) in enumerate(eachline(joinpath(DATA_DIR, name)))
        index == 1 && continue
        isempty(strip(line)) && continue
        x, y = split(strip(line), ',')
        push!(rows, (parse(Float64, x), parse(Float64, y)))
    end
    return [row[axis] for row in rows, axis in 1:2]
end

function link_groups(positions; band_y=20, band_x=(0, 20), split=10)
    left, right = Int[], Int[]
    for i in axes(positions, 1)
        x, y = positions[i, 1], positions[i, 2]
        y >= band_y && band_x[1] < x < band_x[2] && push!(x < split ? left : right, i)
    end
    gap((i, j)) = sum((positions[i, :] - positions[j, :]) .^ 2)
    closest = argmin(gap, [(i, j) for i in left, j in right])
    return [left, right], collect(closest)
end

positions = load_observations("kmeans.csv")

num_clusters = 3

Both axes take one scale, never one each. The model measures a Euclidean distance to a center, and an axis stretched to fill the canvas would draw clusters that do not match the arithmetic underneath them.

draw_observations(positions, "Observations")

A baseline with Clustering.jl

Clustering.jl ships a \(k\)-means implementation, which gives us something to compare against.

using Clustering
using Printf
using Random

function baseline_clusters(positions, k; restarts=10)
    Random.seed!(0)
    data = permutedims(positions)
    best = kmeans(data, k)
    for _ in 2:restarts
        candidate = kmeans(data, k)
        candidate.totalcost < best.totalcost && (best = candidate)
    end
    centers = [(best.centers[1, j], best.centers[2, j]) for j in 1:k]
    total_distance = best.totalcost
    return Clusters(total_distance, collect(best.assignments), centers)
end

baseline = baseline_clusters(positions, num_clusters)
draw_clusters(positions, baseline, "Clustering.jl solution")

@printf("Total squared distance (Clustering.jl): %.2f\n", baseline.total_distance)

Total squared distance (Clustering.jl): 26239.03

Model implementation

using JuMP
using KNITRO

minimize_distance(positions, num_clusters; ...) builds the model, solves it with Knitro, and returns Clusters. The keyword arguments are the side constraints, all optional and all left empty for the first solve.

Both the objective \(\sum_i \sum_k a^k_i d^k_i\) and the barycenter constraint multiply two variables, and the same @constraint macro takes them as written.

function minimize_distance(
    positions,
    num_clusters;
    must_link=Vector{Int}[],
    cannot_link=Int[],
    min_size=nothing,
    max_size=nothing,
)
    n = size(positions, 1)
    clusters = 1:num_clusters

    model = Model(KNITRO.Optimizer)

    @variable(model, a[1:n, clusters], Bin)
    @variable(model, d[1:n, clusters] >= 0)
    @variable(model, c[clusters, 1:2])

    @objective(model, Min, sum(a[i, k] * d[i, k] for i in 1:n, k in clusters))

    @constraint(model, [i in 1:n], sum(a[i, k] for k in clusters) == 1)

    @constraint(
        model,
        [k in clusters, axis in 1:2],
        c[k, axis] * sum(a[i, k] for i in 1:n) ==
            sum(positions[i, axis] * a[i, k] for i in 1:n)
    )

    @constraint(
        model,
        [i in 1:n, k in clusters],
        d[i, k] == (positions[i, 1] - c[k, 1])^2 + (positions[i, 2] - c[k, 2])^2
    )

    for group in must_link, i in group[2:end]
        @constraint(model, [k in clusters], a[i, k] == a[group[1], k])
    end
    for x in eachindex(cannot_link), y in (x + 1):lastindex(cannot_link)
        i, j = cannot_link[x], cannot_link[y]
        @constraint(model, [k in clusters], a[i, k] + a[j, k] <= 1)
    end
    for k in clusters
        cluster_size = sum(a[i, k] for i in 1:n)
        max_size === nothing || @constraint(model, cluster_size <= max_size)
        min_size === nothing || @constraint(model, cluster_size >= min_size)
    end

    optimize!(model)

    labels = [argmax([value(a[i, k]) for k in clusters]) for i in 1:n]
    centers = [(value(c[k, 1]), value(c[k, 2])) for k in clusters]
    total_distance = objective_value(model)
    return Clusters(total_distance, labels, centers)
end

Comparison with Clustering.jl

Let’s solve the same instance with Knitro.

clusters = minimize_distance(positions, num_clusters);

=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro changing mip_method from AUTO to 1.
No start point provided -- Knitro computing one.

Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              1e-06
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001
Knitro changing mip_root_nlpalg from AUTO to 1.
Knitro changing mip_node_nlpalg from AUTO to 1.
Knitro changing mip_branchrule from AUTO to 2.
Knitro changing mip_selectrule from AUTO to 2.
Knitro changing mip_mir from AUTO to 2.
Knitro changing mip_clique from AUTO to 0.
Knitro changing mip_zerohalf from AUTO to 0.
Knitro changing mip_liftproject from AUTO to 0.
Knitro changing mip_knapsack from AUTO to 2.
Knitro changing mip_gomory from AUTO to 0.
Knitro changing mip_cut_flowcover from AUTO to 2.
Knitro changing mip_cut_probing from AUTO to 1.
Knitro changing mip_rounding from AUTO to 3.
Knitro changing mip_heuristic_strategy from AUTO to 1.
Knitro changing mip_heuristic_feaspump from AUTO to 1.
Knitro changing mip_heuristic_misqp from AUTO to 0.
Knitro changing mip_heuristic_mpec from AUTO to 1.
Knitro changing mip_heuristic_diving from AUTO to 1926.
Knitro changing mip_heuristic_fixpropagate from AUTO to 62.
Knitro changing mip_heuristic_lns from AUTO to 0.
Knitro changing mip_heuristic_localsearch from AUTO to 1.
Knitro changing mip_pseudoinit from AUTO to 1.
WARNING: Problem appears to have nonlinear equalities and be non-convex.
         The Knitro mixed integer solver is designed for convex problems.
         For non-convex problems it is only a heuristic, and the reported
         bounds and optimality claims cannot be verified.


Problem Characteristics                     |           Presolved
-----------------------
Problem type: MIQCQP
Objective: minimize / quadratic
Number of variables:                   1806 |                          1806
  bounds:         lower     upper     range |     lower     upper     range
                    900         0       900 |       900         0       900
                             free     fixed |                free     fixed
                                6         0 |                   6         0
                  cont.    binary   integer |     cont.    binary   integer
                    906       900         0 |       906       900         0
Number of constraints:                 1206 |                          1206
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:           300         0         0 |       300         0         0
  quadratic:        906         0         0 |       906         0         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0      5400           |         0      5400          
  quadratic:       1800      3606      2706 |      1800      3606      2706
  total:           1800      5406      2706 |      1800      5406      2706

Knitro using Branch and Bound method with 8 threads.

Initial points
--------------
No initial point provided for the root node relaxation.
No primal point provided for the MIP.

Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [9e-02, 8e+01] |                [9e-02, 8e+01]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 2e+03] |                [1e+00, 2e+03]

Root node relaxation
--------------------

 Iter      Objective      Feasibility        Optimality       Time 
                             error              error        (secs)
 ----      ---------      -----------        ----------      ------
    0        3675.73          1109.80          0.360456       0.240
    1        61426.5         0.630342           281.628       0.246
    2        61426.8      6.50194e-02           18.0260       0.253
    3        61417.5         0.794125           1.50530       0.260
    4        61410.5          2.36487           13.0541       0.267
    5        61384.1          1.21386       5.09325e-02       0.276
    6        61357.5          1.98839           32.8442       0.283
    7        61280.9          2.84366           51.3784       0.289
    8        61050.1          15.4803           142.905       0.301
    9        61028.9          7.53393           215.335       0.316
   10        60839.7          5.56279          0.230524       0.331
   11        60802.7          5.55543           3.08086       0.343
   12        59885.4          52.4410           130.276       0.355
   13        59575.5          47.6464           21.1059       0.370
   14        58946.9          45.8579           151.715       0.380
   15        58604.3          13.0982           257.044       0.398
   16        58064.4         0.427745          0.417459       0.413
   17        57775.7         0.375513          0.456829       0.426
   18        57144.1          1.47808           2261.32       0.437
   19        56592.6          16.3205          0.523538       0.452
   20        56065.6          4.18539           1439.51       0.464
   30        44616.1          4.98346           271.086       0.582
   40        43262.7         0.491788           17.1441       0.705
   50        41972.6         0.397145           143.737       0.819
   60        40981.0         0.127645           367.458       0.930
   70        40551.5         0.165141           100.189       1.015
   80        40252.1         0.167335           201.469       1.090
   90        39926.0         0.119874           49.1061       1.183
  100        39590.7      1.83578e-02           19.2105       1.309
  110        39166.0      1.68466e-02           15.4386       1.450
  120        38735.6      1.12765e-02           34.5056       1.569
  130        38083.7      7.26469e-02           30.2949       1.703
  140        37359.6      3.06462e-02           31.7484       1.846
  150        36265.2         0.408907           31.7567       1.978
  160        35448.6      5.86283e-04           47.8987       2.121
  170        34746.9      1.16792e-02           12.0904       2.263
  180        33672.6      3.35414e-02           37.0736       2.409
  190        32751.3      1.75575e-02           27.1773       2.543
  200        31941.0      2.92960e-02           44.3317       2.698
  300        26377.0         0.157874           24.7651       3.608

Tree search
-----------

       Nodes        Best solution   Best bound      Gap       Time 
   Expl  |  Unexpl      value         value                  (secs)
   ---------------  -------------   ----------      ---      ------
      1       0      26241.0 LEAF      26241.0      0.00%     3.807

EXIT: Optimal solution found (assuming convexity).

HINT: The problem may be a non-convex mixed-integer problem.  Set
      mip_multistart=1 to enable a mixed-integer multistart heuristic,
      which may improve the chances of finding the global solution.

Final Statistics for MIP
------------------------
Final objective value               =  2.62409546350688906e+04
Final bound value                   =  2.62409546350688906e+04
Final optimality gap (abs / rel)    =  0.00000e+00 / 0.00000e+00 (0.00%)
# of root cutting plane rounds      =  0
# of restarts                       =  0
# of nodes processed                =  1 (3.790s)
# of strong branching evaluations   =  0 (0.000s)
# of function evaluations           =  0 (0.000s)
# of gradient evaluations           =  0 (0.000s)
# of hessian evaluations            =  0 (0.000s)
# of hessian-vector evaluations     =  0
# of subproblems processed          =  1 (3.790s)
Total program time (secs)           =  3.80759 (4.926 CPU time)
Time spent in evaluations (secs)    =  0.00000

Cuts statistics (gen / add)
---------------------------
Knapsack cuts                       =  0 / 0
Mixed-integer rounding cuts         =  0 / 0
Flow-cover cuts                     =  0 / 0
Probing cuts                        =  0 / 0

Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump                    =  0 / 0 / 0.000s
Rounding heuristic                  =  0 / 0 / 0.000s
MPEC heuristic                      =  0 / 0 / 0.000s
Local search heuristic              =  4 / 0 / 1.331s

===========================================================================

Show the full outputHide the full output

draw_clusters(positions, clusters, "Knitro solution")

report_clusters(positions, clusters)

cluster  observations   centre x   centre y    distance  share
--------------------------------------------------------------
      1            89      20.86       7.30     8360.74  31.9%
      2            92       9.16      21.32     7742.78  29.5%
      3           119       2.01       3.77    10137.44  38.6%

Total squared distance 26240.95, cluster sizes 89 to 119

The two totals agree to four significant figures, and Knitro takes longer to get there. Neither is proved optimal. The model is nonconvex, so Knitro returns a local optimum, and baseline_clusters keeps the best of ten random restarts for the same reason. On the plain problem there is nothing to gain by switching. The rest of this page is about what the dedicated algorithm cannot do.

Must-link and cannot-link

We single out the observations in the upper band and split them either side of its midline, giving two groups that must each stay together. The pair that straddles the split most closely then has to end up in different clusters, dragging the two groups apart with it.

must_link, cannot_link = link_groups(positions)
draw_link_groups(positions, must_link, cannot_link, "Affected observations")

linked = minimize_distance(
    positions, num_clusters; must_link=must_link, cannot_link=cannot_link
)

=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro changing mip_method from AUTO to 1.
No start point provided -- Knitro computing one.

Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              1e-06
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001
Knitro changing mip_root_nlpalg from AUTO to 1.
Knitro changing mip_node_nlpalg from AUTO to 1.
Knitro changing mip_branchrule from AUTO to 2.
Knitro changing mip_selectrule from AUTO to 2.
Knitro changing mip_mir from AUTO to 2.
Knitro changing mip_clique from AUTO to 0.
Knitro changing mip_zerohalf from AUTO to 0.
Knitro changing mip_liftproject from AUTO to 0.
Knitro changing mip_knapsack from AUTO to 2.
Knitro changing mip_gomory from AUTO to 0.
Knitro changing mip_cut_flowcover from AUTO to 2.
Knitro changing mip_cut_probing from AUTO to 1.
Knitro changing mip_rounding from AUTO to 3.
Knitro changing mip_heuristic_strategy from AUTO to 1.
Knitro changing mip_heuristic_feaspump from AUTO to 1.
Knitro changing mip_heuristic_misqp from AUTO to 0.
Knitro changing mip_heuristic_mpec from AUTO to 1.
Knitro changing mip_heuristic_diving from AUTO to 1926.
Knitro changing mip_heuristic_fixpropagate from AUTO to 62.
Knitro changing mip_heuristic_lns from AUTO to 0.
Knitro changing mip_heuristic_localsearch from AUTO to 1.
Knitro changing mip_pseudoinit from AUTO to 1.
WARNING: Problem appears to have nonlinear equalities and be non-convex.
         The Knitro mixed integer solver is designed for convex problems.
         For non-convex problems it is only a heuristic, and the reported
         bounds and optimality claims cannot be verified.


Problem Characteristics                     |           Presolved
-----------------------
Problem type: MIQCQP
Objective: minimize / quadratic
Number of variables:                   1806 |                          1806
  bounds:         lower     upper     range |     lower     upper     range
                    900         0       900 |       900         0       900
                             free     fixed |                free     fixed
                                6         0 |                   6         0
                  cont.    binary   integer |     cont.    binary   integer
                    906       900         0 |       906       900         0
Number of constraints:                 1314 |                          1314
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:           405         3         0 |       405         3         0
  quadratic:        906         0         0 |       906         0         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0      5616           |         0      5616          
  quadratic:       1800      3606      2706 |      1800      3606      2706
  total:           1800      5622      2706 |      1800      5622      2706

Knitro using Branch and Bound method with 8 threads.

Initial points
--------------
No initial point provided for the root node relaxation.
No primal point provided for the MIP.

Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [9e-02, 8e+01] |                [9e-02, 8e+01]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 2e+03] |                [1e+00, 2e+03]

Root node relaxation
--------------------

 Iter      Objective      Feasibility        Optimality       Time 
                             error              error        (secs)
 ----      ---------      -----------        ----------      ------
    0        3681.31          1109.81          0.367365       0.029
    1        61426.7         0.591051           197.582       0.036
    2        61427.1      6.63340e-02           16.3770       0.044
    3        61417.7          1.71448           3.56013       0.050
    4        61413.4          2.12257           9.83269       0.058
    5        61399.2         0.937200       4.25888e-02       0.072
    6        61383.0          2.10165           10.6803       0.080
    7        61366.7          2.28822           20.2986       0.086
    8        61318.3         0.215271           39.3928       0.095
    9        61123.0          1.23451          0.151873       0.106
   10        60950.0          4.07360           114.586       0.114
   11        60683.5          2.45776           327.595       0.124
   12        60213.7          12.3542           448.880       0.130
   13        60213.8          12.3542           448.862       0.136
   14        60261.7          5.47322           444.808       0.146
   15        60255.8          4.19062           428.692       0.154
   16        58752.7          12.3165           733.436       0.161
   17        58370.2          2.76969           175.789       0.172
   18        58030.4          1.09956           301.288       0.182
   19        57512.0          1.07736           434.060       0.195
   20        57230.1         0.220811           384.733       0.207
   30        52926.0         0.531627           6.67588       0.327
   40        45938.2         0.739153           173.419       0.418
   50        43967.7         0.283920          0.484802       0.539
   60        42412.1      1.90211e-02           3.95295       0.635
   70        41201.1      2.92444e-03           6.03693       0.729
   80        40263.1         0.113854           11.8735       0.818
   90        39733.5      2.30983e-02           13.8659       0.905
  100        39210.1      4.35561e-02           45.6406       0.988
  110        38628.0      7.21905e-03           27.6717       1.092
  120        38318.1      2.10533e-02           271.726       1.179
  130        37809.2         0.159179           29.5659       1.278
  140        37347.9      1.72047e-04           47.3631       1.403
  150        37116.4      1.04859e-02           54.2490       1.476
  160        36882.3      2.91801e-02           58.1421       1.545
  170        36458.6      8.47973e-02           44.1092       1.640
  180        36162.4      2.23908e-02           7.91962       1.756
  190        35916.8         0.113543           177.697       1.824
  200        35806.1      3.21231e-02           1.38222       1.889
  300        35370.9         0.701185           798.616       2.779
  400        31427.0      9.85838e-02           26.1264       3.477

Tree search
-----------

       Nodes        Best solution   Best bound      Gap       Time 
   Expl  |  Unexpl      value         value                  (secs)
   ---------------  -------------   ----------      ---      ------
      1       2      30829.6 FCRD         -inf                4.073
      3       0      30829.6           30829.6     -0.00%     9.504

EXIT: Optimal solution found (assuming convexity).

HINT: The problem may be a non-convex mixed-integer problem.  Set
      mip_multistart=1 to enable a mixed-integer multistart heuristic,
      which may improve the chances of finding the global solution.

Final Statistics for MIP
------------------------
Final objective value               =  3.08296499723525121e+04
Final bound value                   =  3.08296499903364311e+04
Final optimality gap (abs / rel)    =  -1.79839e-05 / -5.83332e-10 (-0.00%)
# of root cutting plane rounds      =  1
# of restarts                       =  0
# of nodes processed                =  3 (4.040s)
# of strong branching evaluations   =  0 (0.000s)
# of function evaluations           =  0 (0.000s)
# of gradient evaluations           =  0 (0.000s)
# of hessian evaluations            =  0 (0.000s)
# of hessian-vector evaluations     =  0
# of subproblems processed          =  2 (4.040s)
Total program time (secs)           =  9.50472 (17.757 CPU time)
Time spent in evaluations (secs)    =  0.00000

Cuts statistics (gen / add)
---------------------------
Knapsack cuts                       =  0 / 0
Mixed-integer rounding cuts         =  0 / 0
Flow-cover cuts                     =  0 / 0
Probing cuts                        =  1944 / 0

Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump                    =  0 / 0 / 0.000s
Rounding heuristic                  =  1 / 1 / 0.004s
MPEC heuristic                      =  0 / 0 / 0.000s
Local search heuristic              =  5 / 0 / 7.579s

===========================================================================

Clusters(30829.649972352512, [2, 3, 2, 3, 1, 3, 1, 2, 1, 2  …  2, 1, 2, 2, 3, 1, 2, 1, 2, 2], [(17.710206281690144, 1.5157433211267606), (2.195431164873417, 9.995688184810126), (18.78723630985915, 19.335682112676054)])

Show the full outputHide the full output

draw_clusters(positions, linked, "Knitro solution with link constraints")

Cluster sizes

A second dataset, split into 10 clusters.

positions_2 = load_observations("kmeans2.csv")
draw_observations(positions_2, "Observations")

Left alone, the clusters come out very uneven:

unbalanced = minimize_distance(positions_2, 10);

=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro changing mip_method from AUTO to 1.
No start point provided -- Knitro computing one.

Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              1e-06
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001
Knitro changing mip_root_nlpalg from AUTO to 1.
Knitro changing mip_node_nlpalg from AUTO to 1.
Knitro changing mip_branchrule from AUTO to 2.
Knitro changing mip_selectrule from AUTO to 2.
Knitro changing mip_mir from AUTO to 2.
Knitro changing mip_clique from AUTO to 0.
Knitro changing mip_zerohalf from AUTO to 0.
Knitro changing mip_liftproject from AUTO to 0.
Knitro changing mip_knapsack from AUTO to 2.
Knitro changing mip_gomory from AUTO to 0.
Knitro changing mip_cut_flowcover from AUTO to 2.
Knitro changing mip_cut_probing from AUTO to 1.
Knitro changing mip_rounding from AUTO to 3.
Knitro changing mip_heuristic_strategy from AUTO to 1.
Knitro changing mip_heuristic_feaspump from AUTO to 1.
Knitro changing mip_heuristic_misqp from AUTO to 0.
Knitro changing mip_heuristic_mpec from AUTO to 1.
Knitro changing mip_heuristic_diving from AUTO to 1926.
Knitro changing mip_heuristic_fixpropagate from AUTO to 62.
Knitro changing mip_heuristic_lns from AUTO to 0.
Knitro changing mip_heuristic_localsearch from AUTO to 1.
Knitro changing mip_pseudoinit from AUTO to 1.
WARNING: Problem appears to have nonlinear equalities and be non-convex.
         The Knitro mixed integer solver is designed for convex problems.
         For non-convex problems it is only a heuristic, and the reported
         bounds and optimality claims cannot be verified.


Problem Characteristics                     |           Presolved
-----------------------
Problem type: MIQCQP
Objective: minimize / quadratic
Number of variables:                   1820 |                          1820
  bounds:         lower     upper     range |     lower     upper     range
                    900         0       900 |       900         0       900
                             free     fixed |                free     fixed
                               20         0 |                  20         0
                  cont.    binary   integer |     cont.    binary   integer
                    920       900         0 |       920       900         0
Number of constraints:                 1010 |                          1010
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            90         0         0 |        90         0         0
  quadratic:        920         0         0 |       920         0         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0      5400           |         0      5400          
  quadratic:       1800      3620      2720 |      1800      3620      2720
  total:           1800      5420      2720 |      1800      5420      2720

Knitro using Branch and Bound method with 8 threads.

Initial points
--------------
No initial point provided for the root node relaxation.
No primal point provided for the MIP.

Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [9e-02, 8e+01] |                [9e-02, 8e+01]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 2e+03] |                [1e+00, 2e+03]

Root node relaxation
--------------------

 Iter      Objective      Feasibility        Optimality       Time 
                             error              error        (secs)
 ----      ---------      -----------        ----------      ------
    0        1364.64          1109.32          0.251704       0.017
    1        1487.26          1102.39           4.45510       0.023
    2        2058.04          1070.27           3.91727       0.028
    3        6849.19          782.921           4.36437       0.033
    4        7889.53          715.945           75.7832       0.038
    5        17362.3          45.2973           56.0640       0.049
    6        17890.7         0.662585           8.37312       0.056
    7        17853.7          2.14573           7.52516       0.065
    8        17762.1          6.62611           14.2644       0.074
    9        17568.2          10.1849           17.9898       0.084
   10        16967.5          15.0375           49.8952       0.093
   11        16014.5          11.4053           168.231       0.103
   12        15675.1          10.7008           175.034       0.112
   13        15783.5          10.4709           151.589       0.123
   14        15031.9          5.44026           175.122       0.132
   15        14506.5          3.00338           124.122       0.141
   16        13368.8          3.87339           166.787       0.151
   17        13135.4          4.23828           148.172       0.160
   18        12937.6          6.08570          0.699952       0.171
   19        12812.7          5.84701           19.4389       0.180
   20        12780.2          5.76639           44.1726       0.190
   30        12750.5          1.63942           53.9187       0.285
   40        12524.0         0.225108           20.7614       0.379
   50        12259.5      9.13712e-02           9.04623       0.472
   60        11941.3         0.268268           10.9666       0.580
   70        11645.8         0.136398           6.35913       0.683
   80        11389.7      8.05797e-02           19.4013       0.783
   90        11218.8         0.531764           18.0356       0.888
  100        10994.7         0.638335           18.4313       0.987
  110        10783.5         0.444057           31.9235       1.086
  120        10564.1         0.227965           45.0711       1.188
  130        10279.7         0.195787           58.4761       1.289
  140        9947.10      7.86821e-02           16.0948       1.387
  150        9753.03      2.00208e-02           22.6803       1.488
  160        9293.74         0.263786           5.37519       1.583
  170        9142.45      1.81938e-03           36.8759       1.692
  180        8860.88      1.27845e-03           36.3541       1.792
  190        8495.04      6.10454e-02           52.5508       1.896
  200        7857.72      1.23733e-02           14.6896       2.002
  300        5158.68         0.137093           28.0698       3.065
  400        4237.96         0.557513           21.6649       4.098
  500        3248.17         0.260605           16.6118       5.235
  600        2671.96      7.20648e-03          0.115830       6.250

Tree search
-----------

       Nodes        Best solution   Best bound      Gap       Time 
   Expl  |  Unexpl      value         value                  (secs)
   ---------------  -------------   ----------      ---      ------
      1       0      2491.35 LEAF      2491.35      0.00%     6.568

EXIT: Optimal solution found (assuming convexity).

HINT: The problem may be a non-convex mixed-integer problem.  Set
      mip_multistart=1 to enable a mixed-integer multistart heuristic,
      which may improve the chances of finding the global solution.

Final Statistics for MIP
------------------------
Final objective value               =  2.49134987285722309e+03
Final bound value                   =  2.49134987285722309e+03
Final optimality gap (abs / rel)    =  0.00000e+00 / 0.00000e+00 (0.00%)
# of root cutting plane rounds      =  0
# of restarts                       =  0
# of nodes processed                =  1 (6.564s)
# of strong branching evaluations   =  0 (0.000s)
# of function evaluations           =  0 (0.000s)
# of gradient evaluations           =  0 (0.000s)
# of hessian evaluations            =  0 (0.000s)
# of hessian-vector evaluations     =  0
# of subproblems processed          =  1 (6.564s)
Total program time (secs)           =  6.56855 (7.176 CPU time)
Time spent in evaluations (secs)    =  0.00000

Cuts statistics (gen / add)
---------------------------
Knapsack cuts                       =  0 / 0
Mixed-integer rounding cuts         =  0 / 0
Flow-cover cuts                     =  0 / 0
Probing cuts                        =  0 / 0

Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump                    =  0 / 0 / 0.000s
Rounding heuristic                  =  0 / 0 / 0.000s
MPEC heuristic                      =  0 / 0 / 0.000s
Local search heuristic              =  4 / 0 / 0.612s

===========================================================================

Show the full outputHide the full output

draw_clusters(positions_2, unbalanced, "Clusters of various sizes")

report_clusters(positions_2, unbalanced)

cluster  observations   centre x   centre y    distance  share
--------------------------------------------------------------
      1             9      14.35      15.27      125.50   5.0%
      2            11       7.14       8.41      163.88   6.6%
      3             7      16.61      28.53      476.18  19.1%
      4             6       9.71      -0.04       64.58   2.6%
      5             7      23.53      -2.19      298.08  12.0%
      6             7      -2.39       8.86      164.22   6.6%
      7            14       4.80      19.08      298.38  12.0%
      8             5      30.05      15.42       88.10   3.5%
      9            11      19.80       8.37      154.93   6.2%
     10            13      -1.15      -1.89      657.50  26.4%

Total squared distance 2491.35, cluster sizes 5 to 14

Bounding every cluster between 8 and 12 observations evens them out:

balanced = minimize_distance(positions_2, 10; min_size=8, max_size=12);

=======================================
          Commercial License
         Artelys Knitro 16.0.0
=======================================

Knitro changing mip_method from AUTO to 1.
No start point provided -- Knitro computing one.

Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.

datacheck                0
feastol                  1e-06
feastol_abs              1e-06
hessian_no_f             1
opttol                   1e-06
opttol_abs               0.001
Knitro changing mip_root_nlpalg from AUTO to 1.
Knitro changing mip_node_nlpalg from AUTO to 1.
Knitro changing mip_branchrule from AUTO to 2.
Knitro changing mip_selectrule from AUTO to 2.
Knitro changing mip_mir from AUTO to 2.
Knitro changing mip_clique from AUTO to 0.
Knitro changing mip_zerohalf from AUTO to 0.
Knitro changing mip_liftproject from AUTO to 0.
Knitro changing mip_knapsack from AUTO to 2.
Knitro changing mip_gomory from AUTO to 0.
Knitro changing mip_cut_flowcover from AUTO to 2.
Knitro changing mip_cut_probing from AUTO to 1.
Knitro changing mip_rounding from AUTO to 3.
Knitro changing mip_heuristic_strategy from AUTO to 1.
Knitro changing mip_heuristic_feaspump from AUTO to 1.
Knitro changing mip_heuristic_misqp from AUTO to 0.
Knitro changing mip_heuristic_mpec from AUTO to 1.
Knitro changing mip_heuristic_diving from AUTO to 1926.
Knitro changing mip_heuristic_fixpropagate from AUTO to 62.
Knitro changing mip_heuristic_lns from AUTO to 0.
Knitro changing mip_heuristic_localsearch from AUTO to 1.
Knitro changing mip_pseudoinit from AUTO to 1.
WARNING: Problem appears to have nonlinear equalities and be non-convex.
         The Knitro mixed integer solver is designed for convex problems.
         For non-convex problems it is only a heuristic, and the reported
         bounds and optimality claims cannot be verified.


Problem Characteristics                     |           Presolved
-----------------------
Problem type: MIQCQP
Objective: minimize / quadratic
Number of variables:                   1820 |                          1820
  bounds:         lower     upper     range |     lower     upper     range
                    900         0       900 |       900         0       900
                             free     fixed |                free     fixed
                               20         0 |                  20         0
                  cont.    binary   integer |     cont.    binary   integer
                    920       900         0 |       920       900         0
Number of constraints:                 1030 |                          1030
                    eq.     ineq.     range |       eq.     ineq.     range
  linear:            90        20         0 |        90        20         0
  quadratic:        920         0         0 |       920         0         0
Number of nonzeros:
              objective  Jacobian   Hessian | objective  Jacobian   Hessian
  linear:             0      7200           |         0      7200          
  quadratic:       1800      3620      2720 |      1800      3620      2720
  total:           1800      7220      2720 |      1800      7220      2720

Knitro using Branch and Bound method with 8 threads.

Initial points
--------------
No initial point provided for the root node relaxation.
No primal point provided for the MIP.

Coefficient range:
  linear objective:          [0e+00, 0e+00] |                [0e+00, 0e+00]
  linear constraints:        [9e-02, 8e+01] |                [9e-02, 8e+01]
  quadratic objective:       [1e+00, 1e+00] |                [1e+00, 1e+00]
  quadratic constraints:     [1e+00, 1e+00] |                [1e+00, 1e+00]
  variable bounds:           [1e+00, 1e+00] |                [1e+00, 1e+00]
  constraint bounds:         [1e+00, 2e+03] |                [1e+00, 2e+03]

Root node relaxation
--------------------

 Iter      Objective      Feasibility        Optimality       Time 
                             error              error        (secs)
 ----      ---------      -----------        ----------      ------
    0        1317.45          1109.32          0.242816       0.020
    1        1366.12          1106.47           1.31740       0.028
    2        1950.78          1072.77           1.38279       0.033
    3        8407.03          674.042           4.34476       0.040
    4        9956.41          569.081           9.38570       0.046
    5        12942.4          367.850           9.26963       0.059
    6        14810.0          227.544           14.5533       0.065
    7        16497.8          103.282           14.0713       0.078
    8        17816.1      3.20291e-02           3.29643       0.087
    9        17721.4          1.09937           6.65690       0.098
   10        17660.3         0.761767           11.6332       0.111
   11        17535.3         0.798871           23.7737       0.121
   12        17501.3         0.803505           16.8992       0.134
   13        17458.0          1.35236           18.2236       0.150
   14        17455.6          1.24801          0.592996       0.166
   15        17416.5          1.16457           7.61737       0.178
   16        17300.6          1.41027           13.0000       0.188
   17        17050.6          1.63572           8.92913       0.201
   18        16363.7          6.71881           24.3250       0.213
   19        16327.8          11.5606          0.333615       0.229
   20        16321.4          9.15179          0.258204       0.243
   30        15088.8      5.63175e-02           43.0317       0.373
   40        12022.5          2.31670           70.0070       0.493
   50        10476.6         0.969247           11.9341       0.626
   60        9866.46         0.360787           36.1910       0.756
   70        9482.59         0.309431           36.2022       0.883
   80        9147.15         0.448380           35.8501       1.013
   90        8828.91         0.281554           19.5151       1.145
  100        8301.66      2.60765e-03           64.8484       1.276
  110        7856.34      1.10141e-02           12.7694       1.426
  120        7302.88      3.04567e-04           23.1308       1.569
  130        6812.44      8.85710e-05           12.1087       1.710
  140        6379.82      4.53837e-02           9.84716       1.857
  150        5777.31      1.50649e-02           6.60727       1.999
  160        5484.62      2.97521e-02           9.54750       2.139
  170        5160.35      5.70674e-02           25.8640       2.279
  180        4833.28      7.69890e-04           12.7415       2.422
  190        4623.19      6.94168e-02           10.0239       2.560
  200        4238.98         0.179234           9.02433       2.698
  300        2626.49      1.91158e-08       7.36914e-07       4.031

Tree search
-----------

       Nodes        Best solution   Best bound      Gap       Time 
   Expl  |  Unexpl      value         value                  (secs)
   ---------------  -------------   ----------      ---      ------
      1       2      2626.49 FCRD         -inf                4.040
      3       0      2626.49           2626.49     -0.00%     4.649

EXIT: Optimal solution found (assuming convexity).

HINT: The problem may be a non-convex mixed-integer problem.  Set
      mip_multistart=1 to enable a mixed-integer multistart heuristic,
      which may improve the chances of finding the global solution.

Final Statistics for MIP
------------------------
Final objective value               =  2.62649190314850193e+03
Final bound value                   =  2.62649190398071187e+03
Final optimality gap (abs / rel)    =  -8.32210e-07 / -3.16852e-10 (-0.00%)
# of root cutting plane rounds      =  1
# of restarts                       =  0
# of nodes processed                =  3 (4.026s)
# of strong branching evaluations   =  0 (0.000s)
# of function evaluations           =  0 (0.000s)
# of gradient evaluations           =  0 (0.000s)
# of hessian evaluations            =  0 (0.000s)
# of hessian-vector evaluations     =  0
# of subproblems processed          =  2 (4.026s)
Total program time (secs)           =  4.64900 (6.632 CPU time)
Time spent in evaluations (secs)    =  0.00000

Cuts statistics (gen / add)
---------------------------
Knapsack cuts                       =  0 / 0
Mixed-integer rounding cuts         =  0 / 0
Flow-cover cuts                     =  0 / 0
Probing cuts                        =  0 / 0

Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump                    =  0 / 0 / 0.000s
Rounding heuristic                  =  1 / 1 / 0.002s
MPEC heuristic                      =  0 / 0 / 0.000s
Local search heuristic              =  6 / 0 / 2.264s

===========================================================================

Show the full outputHide the full output

draw_clusters(positions_2, balanced, "Clusters with restricted sizes")

report_clusters(positions_2, balanced)

cluster  observations   centre x   centre y    distance  share
--------------------------------------------------------------
      1            10       3.47      14.18      107.94   4.1%
      2             9      -0.41       7.84      254.38   9.7%
      3             8       6.50      21.90       41.44   1.6%
      4             8       7.17      -0.42      126.16   4.8%
      5            10      14.03       7.38      151.38   5.8%
      6             8      13.89      14.71       87.36   3.3%
      7             8      16.80      27.43      545.61  20.8%
      8            10      26.80      12.24      429.09  16.3%
      9            11      -2.05      -1.94      582.71  22.2%
     10             8      22.37      -1.05      300.43  11.4%

Total squared distance 2626.49, cluster sizes 8 to 11

Balancing is not free. The bounds pull observations away from their nearest center, so the best balanced clustering costs more than the best unbalanced one. The two solves above need not show that, and may even come out the other way round, since each stops at a local optimum and neither is the best of its problem. What matters is that the requirement can be stated at all, and that is the reason to write the problem out rather than call a library.