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.620
1 61428.9 0.810476 238.200 0.631
2 61429.7 2.05045e-02 5.27634 0.643
3 61426.8 1.74886 2.31499 0.651
4 61415.5 2.65657 2.28781e-02 0.662
5 61414.2 0.511604 0.486061 0.670
6 61403.1 0.234781 3.44865e-02 0.680
7 61401.2 0.761006 7.04428 0.688
8 61399.1 0.749367 13.6101 0.695
9 61389.1 0.224549 19.2999 0.702
10 61361.3 0.816281 13.2791 0.709
11 61294.6 1.34010 19.8685 0.716
12 61221.0 1.04345 33.3070 0.724
13 61040.4 0.890976 22.6221 0.732
14 60753.6 0.736901 46.7827 0.740
15 60435.0 0.646638 52.0580 0.751
16 59991.0 0.447014 163.133 0.765
17 59690.2 1.39782 226.868 0.779
18 59382.5 0.195866 158.893 0.801
19 59299.2 0.890367 242.353 0.822
20 57809.0 44.3912 722.797 0.842
30 55747.8 0.298749 1215.42 0.955
40 53940.9 0.493559 1919.67 1.021
50 46380.1 3.33733 1411.18 1.098
60 43147.7 1.16083 768.696 1.166
70 42322.8 1.00012 79.6175 1.250
80 40446.4 0.293124 8.79460 1.332
90 39382.8 0.123083 104.402 1.428
100 39018.0 1.29656e-02 23.8057 1.524
110 38545.9 3.44844e-02 14.0002 1.615
120 38259.9 4.64016e-02 34.5724 1.688
130 37883.0 5.14092e-02 25.1694 1.773
140 37556.1 0.125260 26.9713 1.863
150 37389.8 1.30155e-03 18.2803 1.935
160 37154.6 0.100262 53.2789 2.011
170 36806.2 3.41522e-02 32.0399 2.098
180 36535.6 1.12986e-02 33.6759 2.188
190 36153.8 5.54266e-05 40.5962 2.293
200 35682.6 1.05951e-02 67.3326 2.394
300 30399.3 1.37069e-03 62.1902 3.272
400 29392.8 2.11691e-02 15.4583 3.970
500 28137.4 8.21879e-03 28.6982 4.743
600 26592.3 4.85259e-03 22.3457 5.498
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 0 26239.0 LEAF 26239.0 0.00% 9.862
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 (9.257s)
# 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 (9.257s)
Total program time (secs) = 9.86244 (8.713 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.264s
===========================================================================
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.167
1 61429.6 0.254031 227.584 0.179
2 61430.4 1.08545e-02 4.81686 0.191
3 61429.4 0.686077 1.12078 0.199
4 61429.0 0.537131 3.87969 0.209
5 61421.2 3.97732 1.56933 0.218
6 61408.1 3.09100 25.7263 0.227
7 61369.6 2.06009 49.1843 0.236
8 61300.3 1.51677 14.8594 0.246
9 61125.6 1.82015 32.8951 0.255
10 61027.1 0.164290 57.1776 0.269
11 60657.7 7.25064 245.033 0.284
12 59708.7 14.7642 534.892 0.295
13 59147.8 33.5905 771.494 0.302
14 59366.5 7.94245 748.651 0.314
15 59384.7 3.93828 967.490 0.325
16 58909.6 11.6697 1050.07 0.332
17 58873.7 11.5071 1015.03 0.340
18 58891.8 11.4921 1027.73 0.350
19 58712.6 10.0673 952.556 0.359
20 58415.2 7.63765 149.433 0.371
30 56017.0 1.09818 1175.92 0.482
40 49673.1 0.754167 606.409 0.583
50 47424.6 0.176186 242.233 0.700
60 47146.8 8.67371e-02 778.502 0.786
70 46817.7 0.400611 1163.51 0.852
80 46321.9 0.231256 1049.62 0.902
90 45069.4 0.136318 817.868 0.994
100 44272.3 7.71488e-02 37.1663 1.092
110 43488.4 1.53221e-02 34.1009 1.192
120 42921.5 1.34353 316.894 1.297
130 41462.3 0.261109 65.7680 1.405
140 40620.2 0.102787 54.5341 1.527
150 39585.8 0.373669 61.5926 1.615
160 39173.4 0.355983 10.1364 1.707
170 38885.8 1.90147e-02 37.0122 1.823
180 38608.7 0.192131 71.3066 1.948
190 38153.7 6.99029e-03 20.0785 2.072
200 37563.3 1.17738e-02 16.8752 2.175
300 33221.5 9.59040e-02 36.7366 3.188
400 31211.1 0.144352 41.5202 4.024
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 2 30829.6 FCRD -inf 4.646
3 0 30829.6 30829.6 -0.00% 8.250
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.525s)
# 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.525s)
Total program time (secs) = 8.25056 (13.286 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.007s
MPEC heuristic = 0 / 0 / 0.000s
Local search heuristic = 6 / 0 / 8.431s
===========================================================================
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.142
1 1489.96 1101.54 3.28216 0.157
2 2554.87 1040.53 3.98401 0.165
3 6337.87 811.462 4.15084 0.173
4 8170.08 690.225 14.3058 0.180
5 12108.0 432.143 22.7216 0.191
6 17872.5 0.127741 9.02348 0.199
7 17665.1 21.3171 20.0914 0.207
8 17541.7 2.45595 0.934700 0.220
9 17459.6 1.77129 0.290171 0.231
10 17288.9 3.85752 25.3903 0.241
11 16757.3 6.41186 17.0476 0.249
12 16564.8 0.180832 30.8332 0.259
13 15899.6 1.84878 110.815 0.269
14 14706.3 2.59603 126.495 0.282
15 14198.1 2.48437 169.424 0.293
16 14075.3 2.52406 140.075 0.307
17 13783.5 2.52502 229.922 0.317
18 13708.6 2.17280 256.247 0.330
19 13134.1 1.86695 301.440 0.340
20 12893.6 1.82505 277.736 0.351
30 12391.7 1.84491e-03 22.1956 0.474
40 11840.0 3.34248e-02 18.7916 0.585
50 11508.3 3.20043e-02 5.08204 0.696
60 11151.1 3.00640e-02 3.05351 0.803
70 10870.9 3.41429e-02 2.13531 0.901
80 10477.9 0.224573 18.8921 1.004
90 10158.6 0.533185 20.3222 1.107
100 9904.70 0.106237 14.8696 1.211
110 9506.61 6.67822e-02 17.6674 1.313
120 8984.51 2.33523e-04 18.6721 1.479
130 8464.28 1.67936e-04 49.7324 1.659
140 8206.58 1.91186e-02 8.60207 1.826
150 7806.42 5.89944e-02 10.4629 1.986
160 7447.03 1.83919e-02 21.3800 2.162
170 7225.31 1.68538e-02 25.8184 2.298
180 6758.03 0.503025 55.1525 2.402
190 6733.69 0.376650 1189.40 2.520
200 7499.78 0.160282 662.631 2.557
300 3394.17 0.152750 9.67668 3.649
400 2446.06 0.673513 47.4696 4.755
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 0 2327.18 LEAF 2327.18 0.00% 5.073
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 (5.063s)
# 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 (5.063s)
Total program time (secs) = 5.07288 (5.531 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.531s
===========================================================================
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.106
1 1366.59 1105.27 1.70860 0.123
2 6125.38 809.372 17.3254 0.134
3 9822.37 562.605 17.5182 0.146
4 12164.7 408.387 36.7584 0.165
5 18127.2 19.7107 20.9620 0.176
6 17792.3 5.14678 6.71419 0.186
7 17644.0 3.34848 4.80460 0.201
8 17536.9 2.46527 8.10628 0.215
9 17335.5 1.52050 15.6926 0.228
10 17165.6 1.14244 13.5905 0.242
11 17132.9 0.117888 20.8888 0.256
12 16989.3 0.551612 37.9685 0.269
13 16892.3 0.674729 56.6391 0.282
14 16550.2 1.50266 66.7281 0.295
15 16489.1 2.89946 5.38732 0.310
16 16211.2 2.34878 27.2421 0.323
17 15845.1 2.22856 30.4512 0.332
18 15537.7 2.48475 29.6480 0.345
19 15375.6 2.40860 33.5132 0.358
20 15098.1 2.67755 34.5911 0.372
30 11743.9 3.28544 43.3686 0.501
40 9399.02 0.620069 19.2005 0.631
50 7443.80 0.723993 23.9865 0.756
60 5896.24 1.30940e-02 21.6349 0.893
70 5510.80 6.48169e-04 17.1411 1.080
80 5151.30 3.37327e-04 16.4899 1.228
90 4882.84 1.46259e-03 21.7494 1.382
100 4627.98 0.170842 68.7438 1.534
110 4460.35 4.15459e-02 21.1606 1.683
120 4369.30 2.02778e-02 19.6876 1.837
130 4134.38 0.227629 16.0811 1.982
140 3817.88 0.388453 46.5484 2.126
150 3642.13 0.167777 37.6463 2.269
160 3538.97 0.131606 49.4239 2.412
170 3424.70 5.32309e-02 18.2683 2.573
180 3068.80 0.106758 13.8749 2.729
190 2944.03 0.337083 7.50018 2.879
200 2876.32 0.158691 7.60475 3.012
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 0 2591.68 LEAF 2591.68 0.00% 3.680
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 (3.671s)
# 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.671s)
Total program time (secs) = 3.68033 (5.036 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.245s
===========================================================================
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.