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.258
    1        61426.5         0.630342           281.628       0.265
    2        61426.8      6.50194e-02           18.0260       0.273
    3        61417.5         0.794125           1.50530       0.280
    4        61410.5          2.36487           13.0541       0.288
    5        61384.1          1.21386       5.09325e-02       0.298
    6        61357.5          1.98839           32.8442       0.306
    7        61280.9          2.84366           51.3784       0.312
    8        61050.1          15.4803           142.905       0.319
    9        61028.9          7.53393           215.335       0.328
   10        60839.7          5.56279          0.230524       0.336
   11        60802.7          5.55543           3.08086       0.344
   12        59885.4          52.4410           130.276       0.351
   13        59575.5          47.6464           21.1059       0.360
   14        58946.9          45.8579           151.715       0.366
   15        58604.3          13.0982           257.044       0.377
   16        58064.4         0.427745          0.417459       0.385
   17        57775.7         0.375513          0.456829       0.393
   18        57144.1          1.47808           2261.32       0.400
   19        56592.6          16.3205          0.523538       0.409
   20        56065.6          4.18539           1439.51       0.416
   30        44616.1          4.98346           271.086       0.489
   40        43262.7         0.491788           17.1441       0.565
   50        41972.6         0.397145           143.737       0.632
   60        40981.0         0.127645           367.458       0.699
   70        40551.5         0.165141           100.189       0.755
   80        40252.1         0.167335           201.469       0.798
   90        39926.0         0.119874           49.1061       0.854
  100        39590.7      1.83578e-02           19.2105       0.931
  110        39166.0      1.68466e-02           15.4386       1.017
  120        38735.6      1.12765e-02           34.5056       1.089
  130        38083.7      7.26469e-02           30.2949       1.169
  140        37359.6      3.06462e-02           31.7484       1.255
  150        36265.2         0.408907           31.7567       1.334
  160        35448.6      5.86283e-04           47.8987       1.427
  170        34746.9      1.16792e-02           12.0904       1.508
  180        33672.6      3.35414e-02           37.0736       1.596
  190        32751.3      1.75575e-02           27.1773       1.677
  200        31941.0      2.92960e-02           44.3317       1.771
  300        26377.0         0.157874           24.7651       2.529

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

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

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 (2.620s)
# 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 (2.620s)
Total program time (secs)           =  2.72549 (4.010 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.373s

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

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.026
    1        61426.7         0.591051           197.582       0.224
    2        61427.1      6.63340e-02           16.3770       0.232
    3        61417.7          1.71448           3.56013       0.238
    4        61413.4          2.12257           9.83269       0.246
    5        61399.2         0.937200       4.25888e-02       0.259
    6        61383.0          2.10165           10.6803       0.267
    7        61366.7          2.28822           20.2986       0.273
    8        61318.3         0.215271           39.3928       0.283
    9        61123.0          1.23451          0.151873       0.293
   10        60950.0          4.07360           114.586       0.302
   11        60683.5          2.45776           327.595       0.312
   12        60213.7          12.3542           448.880       0.318
   13        60213.8          12.3542           448.862       0.323
   14        60261.7          5.47322           444.808       0.334
   15        60255.8          4.19062           428.692       0.342
   16        58752.7          12.3165           733.436       0.348
   17        58370.2          2.76969           175.789       0.359
   18        58030.4          1.09956           301.288       0.370
   19        57512.0          1.07736           434.060       0.382
   20        57230.1         0.220811           384.733       0.395
   30        52926.0         0.531627           6.67588       0.513
   40        45938.2         0.739153           173.419       0.604
   50        43967.7         0.283920          0.484802       0.712
   60        42412.1      1.90211e-02           3.95295       0.809
   70        41201.1      2.92444e-03           6.03693       0.901
   80        40263.1         0.113854           11.8735       0.986
   90        39733.5      2.30983e-02           13.8659       1.073
  100        39210.1      4.35561e-02           45.6406       1.157
  110        38628.0      7.21905e-03           27.6717       1.264
  120        38318.1      2.10533e-02           271.726       1.351
  130        37809.2         0.159179           29.5659       1.441
  140        37347.9      1.72047e-04           47.3631       1.556
  150        37116.4      1.04859e-02           54.2490       1.629
  160        36882.3      2.91801e-02           58.1421       1.696
  170        36458.6      8.47973e-02           44.1092       1.791
  180        36162.4      2.23908e-02           7.91962       1.899
  190        35916.8         0.113543           177.697       1.966
  200        35806.1      3.21231e-02           1.38222       2.034
  300        35370.9         0.701185           798.616       2.917
  400        31427.0      9.85838e-02           26.1264       3.617

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

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

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.144s)
# 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.144s)
Total program time (secs)           =  9.65895 (17.851 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.005s
MPEC heuristic                      =  0 / 0 / 0.000s
Local search heuristic              =  5 / 0 / 7.663s

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

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.029
    1        1487.26          1102.39           4.45510       0.040
    2        2058.04          1070.27           3.91727       0.047
    3        6849.19          782.921           4.36437       0.056
    4        7889.53          715.945           75.7832       0.064
    5        17362.3          45.2973           56.0640       0.082
    6        17890.7         0.662585           8.37312       0.091
    7        17853.7          2.14573           7.52516       0.100
    8        17762.1          6.62611           14.2644       0.110
    9        17568.2          10.1849           17.9898       0.119
   10        16967.5          15.0375           49.8952       0.129
   11        16014.5          11.4053           168.231       0.139
   12        15675.1          10.7008           175.034       0.148
   13        15783.5          10.4709           151.589       0.159
   14        15031.9          5.44026           175.122       0.168
   15        14506.5          3.00338           124.122       0.177
   16        13368.8          3.87339           166.787       0.187
   17        13135.4          4.23828           148.172       0.196
   18        12937.6          6.08570          0.699952       0.207
   19        12812.7          5.84701           19.4389       0.217
   20        12780.2          5.76639           44.1726       0.226
   30        12750.5          1.63942           53.9187       0.321
   40        12524.0         0.225108           20.7614       0.414
   50        12259.5      9.13712e-02           9.04623       0.508
   60        11941.3         0.268268           10.9666       0.617
   70        11645.8         0.136398           6.35913       0.720
   80        11389.7      8.05797e-02           19.4013       0.820
   90        11218.8         0.531764           18.0356       0.926
  100        10994.7         0.638335           18.4313       1.022
  110        10783.5         0.444057           31.9235       1.122
  120        10564.1         0.227965           45.0711       1.226
  130        10279.7         0.195787           58.4761       1.327
  140        9947.10      7.86821e-02           16.0948       1.425
  150        9753.03      2.00208e-02           22.6803       1.525
  160        9293.74         0.263786           5.37519       1.622
  170        9142.45      1.81938e-03           36.8759       1.730
  180        8860.88      1.27845e-03           36.3541       1.831
  190        8495.04      6.10454e-02           52.5508       1.935
  200        7857.72      1.23733e-02           14.6896       2.043
  300        5158.68         0.137093           28.0698       3.118
  400        4237.96         0.557513           21.6649       4.159
  500        3248.17         0.260605           16.6118       5.306
  600        2671.96      7.20648e-03          0.115830       6.335

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

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

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.663s)
# 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.663s)
Total program time (secs)           =  6.66978 (7.265 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.591s

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

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.021
    1        1366.12          1106.47           1.31740       0.028
    2        1950.78          1072.77           1.38279       0.034
    3        8407.03          674.042           4.34476       0.041
    4        9956.41          569.081           9.38570       0.047
    5        12942.4          367.850           9.26963       0.060
    6        14810.0          227.544           14.5533       0.066
    7        16497.8          103.282           14.0713       0.079
    8        17816.1      3.20291e-02           3.29643       0.089
    9        17721.4          1.09937           6.65690       0.099
   10        17660.3         0.761767           11.6332       0.112
   11        17535.3         0.798871           23.7737       0.122
   12        17501.3         0.803505           16.8992       0.136
   13        17458.0          1.35236           18.2236       0.153
   14        17455.6          1.24801          0.592996       0.169
   15        17416.5          1.16457           7.61737       0.183
   16        17300.6          1.41027           13.0000       0.193
   17        17050.6          1.63572           8.92913       0.206
   18        16363.7          6.71881           24.3250       0.219
   19        16327.8          11.5606          0.333615       0.235
   20        16321.4          9.15179          0.258204       0.249
   30        15088.8      5.63175e-02           43.0317       0.383
   40        12022.5          2.31670           70.0070       0.505
   50        10476.6         0.969247           11.9341       0.641
   60        9866.46         0.360787           36.1910       0.773
   70        9482.59         0.309431           36.2022       0.902
   80        9147.15         0.448380           35.8501       1.033
   90        8828.91         0.281554           19.5151       1.167
  100        8301.66      2.60765e-03           64.8484       1.304
  110        7856.34      1.10141e-02           12.7694       1.456
  120        7302.88      3.04567e-04           23.1308       1.601
  130        6812.44      8.85710e-05           12.1087       1.746
  140        6379.82      4.53837e-02           9.84716       1.897
  150        5777.31      1.50649e-02           6.60727       2.039
  160        5484.62      2.97521e-02           9.54750       2.181
  170        5160.35      5.70674e-02           25.8640       2.325
  180        4833.28      7.69890e-04           12.7415       2.476
  190        4623.19      6.94168e-02           10.0239       2.619
  200        4238.98         0.179234           9.02433       2.759
  300        2626.49      1.91158e-08       7.36914e-07       4.113

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

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

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.109s)
# 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.109s)
Total program time (secs)           =  4.73344 (6.776 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.003s
MPEC heuristic                      =  0 / 0 / 0.000s
Local search heuristic              =  6 / 0 / 2.272s

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

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.