from dataclasses import dataclass
@dataclass
class Clusters:
total_distance: float
labels: list
centers: list
@property
def num_clusters(self):
return len(self.centers)
def cluster(self, k):
return [i for i, label in enumerate(self.labels) if label == k]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.
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.
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.
import csv
from pathlib import Path
from plot import draw_clusters, draw_link_groups, draw_observations
from report import report_clusters
DATA_DIR = Path("shared") / "data"
def load_observations(name):
with (DATA_DIR / name).open(encoding="utf-8") as handle:
return [(float(row["X"]), float(row["Y"])) for row in csv.DictReader(handle)]
def link_groups(positions, band_y=20, band_x=(0, 20), split=10):
left, right = [], []
for i, (x, y) in enumerate(positions):
if y >= band_y and band_x[0] < x < band_x[1]:
(left if x < split else right).append(i)
def gap(pair):
i, j = pair
return sum((a - b) ** 2 for a, b in zip(positions[i], positions[j]))
closest = min(((i, j) for i in left for j in right), key=gap)
return [left, right], list(closest)
positions = load_observations("kmeans.csv")
num_clusters = 3Both 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 scikit-learn
Scikit-learn ships a \(k\)-means implementation, so we have something to compare against.
def baseline_clusters(positions, num_clusters):
from sklearn.cluster import KMeans
kmeans = KMeans(n_clusters=num_clusters, random_state=0).fit(positions)
centers = [tuple(center) for center in kmeans.cluster_centers_]
total_distance = kmeans.inertia_
return Clusters(total_distance, list(kmeans.labels_), centers)
baseline = baseline_clusters(positions, num_clusters)
print(f"Total squared distance (scikit-learn): {baseline.total_distance:.2f}")Total squared distance (scikit-learn): 26243.63draw_clusters(positions, baseline, "Scikit-learn solution")Model implementation
import pyomo.environ as pyo
SOLVER_NAME = "knitroampl"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 they go into the model as written.
def minimize_distance(
positions,
num_clusters,
*,
must_link=(),
cannot_link=(),
min_size=None,
max_size=None,
):
n = len(positions)
model = pyo.ConcreteModel()
model.observations = pyo.RangeSet(0, n - 1)
model.clusters = pyo.RangeSet(0, num_clusters - 1)
observations, clusters = model.observations, model.clusters
model.a = pyo.Var(observations, clusters, within=pyo.Binary)
model.d = pyo.Var(observations, clusters, within=pyo.NonNegativeReals)
model.cx = pyo.Var(clusters, within=pyo.Reals)
model.cy = pyo.Var(clusters, within=pyo.Reals)
a, d, cx, cy = model.a, model.d, model.cx, model.cy
cost = pyo.quicksum(a[i, k] * d[i, k] for i in observations for k in clusters)
model.obj = pyo.Objective(expr=cost, sense=pyo.minimize)
def one_cluster_rule(model, i):
return pyo.quicksum(a[i, k] for k in clusters) == 1
model.one_cluster = pyo.Constraint(observations, rule=one_cluster_rule)
def center_x_rule(model, k):
size = pyo.quicksum(a[i, k] for i in observations)
total_x = pyo.quicksum(a[i, k] * positions[i][0] for i in observations)
return cx[k] * size == total_x
model.center_x = pyo.Constraint(clusters, rule=center_x_rule)
def center_y_rule(model, k):
size = pyo.quicksum(a[i, k] for i in observations)
total_y = pyo.quicksum(a[i, k] * positions[i][1] for i in observations)
return cy[k] * size == total_y
model.center_y = pyo.Constraint(clusters, rule=center_y_rule)
def distance_rule(model, i, k):
x, y = positions[i]
return d[i, k] == (x - cx[k]) ** 2 + (y - cy[k]) ** 2
model.distance = pyo.Constraint(observations, clusters, rule=distance_rule)
model.extra = pyo.ConstraintList()
for group in must_link:
for i in group[1:]:
for k in clusters:
model.extra.add(a[i, k] == a[group[0], k])
for x in range(len(cannot_link)):
for y in range(x + 1, len(cannot_link)):
for k in clusters:
model.extra.add(a[cannot_link[x], k] + a[cannot_link[y], k] <= 1)
for k in clusters:
if max_size is not None:
model.extra.add(pyo.quicksum(a[i, k] for i in observations) <= max_size)
if min_size is not None:
model.extra.add(pyo.quicksum(a[i, k] for i in observations) >= min_size)
solver = pyo.SolverFactory(SOLVER_NAME)
solver.solve(model, tee=True)
labels = [max(range(num_clusters), key=lambda k: a[i, k]()) for i in range(n)]
centers = [(cx[k](), cy[k]()) for k in range(num_clusters)]
total_distance = model.obj()
return Clusters(total_distance, labels, centers)Comparison with scikit-learn
Let’s solve the same instance with Knitro.
clusters = minimize_distance(positions, num_clusters)Artelys Knitro 16.0.0:
=======================================
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.
concurrent_evals 0
datacheck 0
feastol 1e-06
feastol_abs 1e-06
findiff_numthreads 1
hessian_no_f 1
hessopt 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 3678.52 1109.81 0.240118 0.210
1 61428.9 0.810476 238.200 0.218
2 61429.7 2.05045e-02 5.27634 0.226
3 61426.8 1.74886 2.31499 0.233
4 61415.5 2.65657 2.28781e-02 0.242
5 61414.2 0.511604 0.486061 0.250
6 61403.1 0.234781 3.44865e-02 0.260
7 61401.2 0.761006 7.04428 0.269
8 61399.1 0.749367 13.6101 0.277
9 61389.1 0.224549 19.2999 0.284
10 61361.3 0.816281 13.2791 0.291
11 61294.6 1.34010 19.8685 0.299
12 61221.0 1.04345 33.3070 0.307
13 61040.4 0.890976 22.6221 0.315
14 60753.6 0.736901 46.7827 0.324
15 60435.0 0.646638 52.0580 0.335
16 59991.0 0.447014 163.133 0.349
17 59690.2 1.39782 226.868 0.362
18 59382.5 0.195866 158.893 0.375
19 59299.2 0.890367 242.353 0.388
20 57809.0 44.3912 722.797 0.400
30 55747.8 0.298749 1215.42 0.504
40 53940.9 0.493559 1919.67 0.573
50 46380.1 3.33733 1411.18 0.654
60 43147.7 1.16083 768.696 0.724
70 42322.8 1.00012 79.6175 0.809
80 40446.4 0.293124 8.79460 0.891
90 39382.8 0.123083 104.402 0.990
100 39018.0 1.29656e-02 23.8057 1.090
110 38545.9 3.44844e-02 14.0002 1.176
120 38259.9 4.64016e-02 34.5724 1.255
130 37883.0 5.14092e-02 25.1694 1.341
140 37556.1 0.125260 26.9713 1.432
150 37389.8 1.30155e-03 18.2803 1.507
160 37154.6 0.100262 53.2789 1.583
170 36806.2 3.41522e-02 32.0399 1.672
180 36535.6 1.12986e-02 33.6759 1.764
190 36153.8 5.54266e-05 40.5962 1.872
200 35682.6 1.05951e-02 67.3326 2.005
300 30399.3 1.37069e-03 62.1902 2.919
400 29392.8 2.11691e-02 15.4583 3.721
500 28137.4 8.21879e-03 28.6982 4.579
600 26592.3 4.85259e-03 22.3457 5.528
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 0 26239.0 LEAF 26239.0 0.00% 6.049
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.62390265181325776e+04
Final bound value = 2.62390265181325776e+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 (6.037s)
# 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.037s)
Total program time (secs) = 6.04903 (7.816 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 = 3 / 0 / 1.202s
===========================================================================
Show the full outputHide the full output
report_clusters(positions, clusters)cluster observations centre x centre y distance share
--------------------------------------------------------------
1 90 20.75 7.27 8452.86 32.2%
2 93 9.14 21.22 7833.92 29.9%
3 117 1.88 3.69 9952.24 37.9%
Total squared distance 26239.03, cluster sizes 90 to 117draw_clusters(positions, clusters, "Knitro solution")The two totals agree to four significant figures, and Knitro takes a few times longer to get there. Neither is proved optimal. The model is nonconvex, so Knitro returns a local optimum, and scikit-learn restarts from several random initializations 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
)Artelys Knitro 16.0.0:
=======================================
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.
concurrent_evals 0
datacheck 0
feastol 1e-06
feastol_abs 1e-06
findiff_numthreads 1
hessian_no_f 1
hessopt 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 3710.84 1109.81 0.245562 0.194
1 61429.6 0.254031 227.584 0.202
2 61430.4 1.08545e-02 4.81686 0.212
3 61429.4 0.686077 1.12078 0.220
4 61429.0 0.537131 3.87969 0.229
5 61421.2 3.97732 1.56933 0.239
6 61408.1 3.09100 25.7263 0.248
7 61369.6 2.06009 49.1843 0.257
8 61300.3 1.51677 14.8594 0.266
9 61125.6 1.82015 32.8951 0.275
10 61027.1 0.164290 57.1776 0.289
11 60657.7 7.25064 245.033 0.303
12 59708.7 14.7642 534.892 0.314
13 59147.8 33.5905 771.494 0.321
14 59366.5 7.94245 748.651 0.334
15 59384.7 3.93828 967.490 0.348
16 58909.6 11.6697 1050.07 0.356
17 58873.7 11.5071 1015.03 0.363
18 58891.8 11.4921 1027.73 0.372
19 58712.6 10.0673 952.556 0.380
20 58415.2 7.63765 149.433 0.392
30 56017.0 1.09818 1175.92 0.501
40 49673.1 0.754167 606.409 0.612
50 47424.6 0.176186 242.233 0.726
60 47146.8 8.67371e-02 778.502 0.813
70 46817.7 0.400611 1163.51 0.883
80 46321.9 0.231256 1049.62 0.937
90 45069.4 0.136318 817.868 1.033
100 44272.3 7.71488e-02 37.1663 1.130
110 43488.4 1.53221e-02 34.1009 1.234
120 42921.5 1.34353 316.894 1.342
130 41462.3 0.261109 65.7680 1.451
140 40620.2 0.102787 54.5341 1.570
150 39585.8 0.373669 61.5926 1.659
160 39173.4 0.355983 10.1364 1.753
170 38885.8 1.90147e-02 37.0122 1.875
180 38608.7 0.192131 71.3066 2.014
190 38153.7 6.99029e-03 20.0785 2.146
200 37563.3 1.17738e-02 16.8752 2.256
300 33221.5 9.59040e-02 36.7366 3.309
400 31211.1 0.144352 41.5202 4.185
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 2 30829.6 FCRD -inf 4.717
3 0 30829.6 30829.6 -0.00% 8.358
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.08296499723525303e+04
Final bound value = 3.08296499820281315e+04
Final optimality gap (abs / rel) = -9.67560e-06 / -3.13841e-10 (-0.00%)
# of root cutting plane rounds = 1
# of restarts = 0
# of nodes processed = 3 (4.680s)
# 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.680s)
Total program time (secs) = 8.35786 (13.476 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 = 6 / 0 / 8.549s
===========================================================================
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)Artelys Knitro 16.0.0:
=======================================
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.
concurrent_evals 0
datacheck 0
feastol 1e-06
feastol_abs 1e-06
findiff_numthreads 1
hessian_no_f 1
hessopt 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 1355.40 1109.29 0.225227 0.139
1 1489.96 1101.54 3.28216 0.150
2 2554.87 1040.53 3.98401 0.156
3 6337.87 811.462 4.15084 0.164
4 8170.08 690.225 14.3058 0.171
5 12108.0 432.143 22.7216 0.184
6 17872.5 0.127741 9.02348 0.194
7 17665.1 21.3171 20.0914 0.204
8 17541.7 2.45595 0.934700 0.221
9 17459.6 1.77129 0.290171 0.237
10 17288.9 3.85752 25.3903 0.250
11 16757.3 6.41186 17.0476 0.260
12 16564.8 0.180832 30.8332 0.274
13 15899.6 1.84878 110.815 0.287
14 14706.3 2.59603 126.495 0.304
15 14198.1 2.48437 169.424 0.318
16 14075.3 2.52406 140.075 0.334
17 13783.5 2.52502 229.922 0.347
18 13708.6 2.17280 256.247 0.364
19 13134.1 1.86695 301.440 0.375
20 12893.6 1.82505 277.736 0.388
30 12391.7 1.84491e-03 22.1956 0.535
40 11840.0 3.34248e-02 18.7916 0.659
50 11508.3 3.20043e-02 5.08204 0.804
60 11151.1 3.00640e-02 3.05351 0.965
70 10870.9 3.41429e-02 2.13531 1.098
80 10477.9 0.224573 18.8921 1.230
90 10158.6 0.533185 20.3222 1.345
100 9904.70 0.106237 14.8696 1.462
110 9506.61 6.67822e-02 17.6674 1.573
120 8984.51 2.33523e-04 18.6721 1.685
130 8464.28 1.67936e-04 49.7324 1.802
140 8206.58 1.91186e-02 8.60207 1.897
150 7806.42 5.89944e-02 10.4629 1.999
160 7447.03 1.83919e-02 21.3800 2.112
170 7225.31 1.68538e-02 25.8184 2.235
180 6758.03 0.503025 55.1525 2.369
190 6733.69 0.376650 1189.40 2.559
200 7499.78 0.160282 662.631 2.615
300 3394.17 0.152750 9.67668 3.624
400 2446.06 0.673513 47.4696 4.654
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 0 2327.18 LEAF 2327.18 0.00% 4.975
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.32718470165902045e+03
Final bound value = 2.32718470165902045e+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 (4.961s)
# 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 (4.961s)
Total program time (secs) = 4.97582 (5.405 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.541s
===========================================================================
Show the full outputHide the full output
report_clusters(positions_2, unbalanced)cluster observations centre x centre y distance share
--------------------------------------------------------------
1 5 21.94 -4.40 138.43 5.9%
2 10 20.65 7.22 125.81 5.4%
3 6 30.26 13.54 195.50 8.4%
4 11 7.57 23.11 150.44 6.5%
5 11 11.05 11.69 132.91 5.7%
6 12 1.51 12.39 176.16 7.6%
7 12 7.56 2.70 274.12 11.8%
8 5 17.88 16.69 39.96 1.7%
9 4 21.25 30.16 218.62 9.4%
10 14 -2.53 -0.83 875.23 37.6%
Total squared distance 2327.18, cluster sizes 4 to 14draw_clusters(positions_2, unbalanced, "Clusters of various sizes")Bounding every cluster between 8 and 12 observations evens them out:
balanced = minimize_distance(positions_2, 10, min_size=8, max_size=12)Artelys Knitro 16.0.0:
=======================================
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.
concurrent_evals 0
datacheck 0
feastol 1e-06
feastol_abs 1e-06
findiff_numthreads 1
hessian_no_f 1
hessopt 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 1300.53 1109.29 0.224282 0.143
1 1366.59 1105.27 1.70860 0.154
2 6125.38 809.372 17.3254 0.162
3 9822.37 562.605 17.5182 0.169
4 12164.7 408.387 36.7584 0.184
5 18127.2 19.7107 20.9620 0.194
6 17792.3 5.14678 6.71419 0.204
7 17644.0 3.34848 4.80460 0.218
8 17536.9 2.46527 8.10628 0.231
9 17335.5 1.52050 15.6926 0.244
10 17165.6 1.14244 13.5905 0.257
11 17132.9 0.117888 20.8888 0.270
12 16989.3 0.551612 37.9685 0.283
13 16892.3 0.674729 56.6391 0.296
14 16550.2 1.50266 66.7281 0.310
15 16489.1 2.89946 5.38732 0.327
16 16211.2 2.34878 27.2421 0.340
17 15845.1 2.22856 30.4512 0.350
18 15537.7 2.48475 29.6480 0.363
19 15375.6 2.40860 33.5132 0.377
20 15098.1 2.67755 34.5911 0.391
30 11743.9 3.28544 43.3686 0.522
40 9399.02 0.620069 19.2005 0.658
50 7443.80 0.723993 23.9865 0.791
60 5896.24 1.30940e-02 21.6349 0.977
70 5510.80 6.48169e-04 17.1411 1.180
80 5151.30 3.37327e-04 16.4899 1.407
90 4882.84 1.46259e-03 21.7494 1.649
100 4627.98 0.170842 68.7438 1.872
110 4460.35 4.15459e-02 21.1606 2.062
120 4369.30 2.02778e-02 19.6876 2.314
130 4134.38 0.227629 16.0811 2.534
140 3817.88 0.388453 46.5484 2.756
150 3642.13 0.167777 37.6463 2.929
160 3538.97 0.131606 49.4239 3.105
170 3424.70 5.32309e-02 18.2683 3.310
180 3068.80 0.106758 13.8749 3.494
190 2944.03 0.337083 7.50018 3.674
200 2876.32 0.158691 7.60475 3.834
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 0 2591.68 LEAF 2591.68 0.00% 4.652
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.59168345860141699e+03
Final bound value = 2.59168345860141699e+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 (4.642s)
# 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 (4.642s)
Total program time (secs) = 4.65203 (5.955 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 = 3 / 0 / 0.234s
===========================================================================
Show the full outputHide the full output
report_clusters(positions_2, balanced)cluster observations centre x centre y distance share
--------------------------------------------------------------
1 8 24.58 1.79 308.66 11.9%
2 8 11.91 -1.36 233.09 9.0%
3 8 26.66 13.79 293.74 11.3%
4 11 5.54 20.30 186.26 7.2%
5 8 13.89 14.71 87.36 3.4%
6 9 14.18 7.84 129.83 5.0%
7 12 3.31 11.18 172.68 6.7%
8 8 -5.82 4.27 395.80 15.3%
9 8 16.80 27.43 545.61 21.1%
10 10 1.79 -2.54 238.65 9.2%
Total squared distance 2591.68, cluster sizes 8 to 12draw_clusters(positions_2, balanced, "Clusters with restricted sizes")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.