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 knitrominimize_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 add_constraint takes them as written.
def minimize_distance(
positions,
num_clusters,
*,
must_link=(),
cannot_link=(),
min_size=None,
max_size=None,
):
n = len(positions)
prob = knitro.Problem()
observations, clusters = range(n), range(num_clusters)
pairs = [(i, k) for i in observations for k in clusters]
a = {pair: prob.add_variable(vtype=knitro.KN_VARTYPE_BINARY) for pair in pairs}
d = {pair: prob.add_variable(lb=0) for pair in pairs}
cx = {k: prob.add_variable() for k in clusters}
cy = {k: prob.add_variable() for k in clusters}
prob.add_objective(prob.nsum(a[pair] * d[pair] for pair in pairs))
for i in observations:
prob.add_constraint(prob.nsum(a[i, k] for k in clusters) == 1)
for k in clusters:
size = prob.nsum(a[i, k] for i in observations)
total_x = prob.nsum(a[i, k] * positions[i][0] for i in observations)
total_y = prob.nsum(a[i, k] * positions[i][1] for i in observations)
prob.add_constraint(cx[k] * size == total_x)
prob.add_constraint(cy[k] * size == total_y)
for i in observations:
x, y = positions[i]
for k in clusters:
prob.add_constraint(d[i, k] == (x - cx[k]) ** 2 + (y - cy[k]) ** 2)
for group in must_link:
for i in group[1:]:
for k in clusters:
prob.add_constraint(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:
prob.add_constraint(a[cannot_link[x], k] + a[cannot_link[y], k] <= 1)
for k in clusters:
if max_size is not None:
prob.add_constraint(prob.nsum(a[i, k] for i in observations) <= max_size)
if min_size is not None:
prob.add_constraint(prob.nsum(a[i, k] for i in observations) >= min_size)
prob.solve()
labels = [max(clusters, key=lambda k: a[i, k].value) for i in observations]
centers = [(cx[k].value, cy[k].value) for k in clusters]
total_distance = prob.get_attr(knitro.KN_ATTR_OBJ_VALUE)
return Clusters(total_distance, labels, centers)Comparison with scikit-learn
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.
concurrent_evals 0
feastol 1e-06
feastol_abs 1e-06
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 expressions: 10829 | 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, 4e+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, 1e+00] | [1e+00, 2e+03]
Root node relaxation
--------------------
Iter Objective Feasibility Optimality Time
error error (secs)
---- --------- ----------- ---------- ------
0 3695.50 1109.15 0.280016 0.150
1 61427.1 1.46486 110.560 0.157
2 61427.8 4.43006e-02 5.14429 0.164
3 61420.7 2.21936 4.06169 0.170
4 61418.2 1.98961 7.23102 0.177
5 61410.1 0.502792 3.09268e-02 0.187
6 61402.1 0.788797 4.11273 0.194
7 61392.4 1.17870 4.30275 0.201
8 61109.0 2.80582 65.0264 0.208
9 61035.3 1.92984 57.0354 0.215
10 60749.6 0.281936 34.1836 0.222
11 60560.3 0.204399 86.0664 0.232
12 60243.9 5.12656 250.824 0.242
13 59967.6 2.79784 0.286995 0.250
14 59319.6 16.8188 554.049 0.258
15 58880.8 10.9504 553.930 0.265
16 59068.7 3.97989 626.526 0.273
17 59159.4 0.971826 811.995 0.282
18 59176.5 0.590024 852.490 0.288
19 59201.4 0.178497 914.084 0.297
20 59173.5 0.709994 991.480 0.305
30 57938.6 2.52692 1464.59 0.372
40 55428.6 5.78292 1132.77 0.454
50 48255.2 0.235207 756.542 0.541
60 40811.6 1.76421 2048.20 0.593
70 40591.5 2.91997e-03 41.2780 0.692
80 39836.4 8.32420e-03 18.9218 0.783
90 39241.5 2.42676e-03 30.5086 0.871
100 38607.6 2.60814e-03 38.9491 0.955
110 38071.8 1.27730e-03 8.86080 1.050
120 37526.3 7.76244e-03 19.9145 1.140
130 36959.6 0.145227 33.9095 1.230
140 36408.4 3.97384e-02 33.7348 1.305
150 35763.1 3.03368e-02 45.1469 1.412
160 35390.9 4.28555e-04 51.7920 1.512
170 34894.7 0.132845 47.3950 1.604
180 34565.5 9.12914e-03 56.4190 1.713
190 34043.2 2.20938e-02 57.4717 1.810
200 33528.7 6.11846e-03 23.6403 1.890
300 30904.2 7.29125e-02 28.4540 2.707
400 29248.2 1.21879e-03 4.05594 3.506
500 27987.3 1.74104e-02 98.3979 4.277
600 26273.3 1.76612e-05 4.71021 5.049
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 0 26240.7 LEAF 26240.7 0.00% 5.219
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.62406575903795310e+04
Final bound value = 2.62406575903795310e+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 (5.207s)
# 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.207s)
Total program time (secs) = 5.22294 (6.490 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.650s
===========================================================================
Show the full outputHide the full output
report_clusters(positions, clusters)cluster observations centre x centre y distance share
--------------------------------------------------------------
1 82 20.10 4.89 6981.86 26.6%
2 118 1.48 4.56 9569.17 36.5%
3 100 11.15 21.20 9689.63 36.9%
Total squared distance 26240.66, cluster sizes 82 to 118draw_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
)=======================================
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.01s.
concurrent_evals 0
feastol 1e-06
feastol_abs 1e-06
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 expressions: 10937 | 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, 4e+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, 1e+00] | [1e+00, 2e+03]
Root node relaxation
--------------------
Iter Objective Feasibility Optimality Time
error error (secs)
---- --------- ----------- ---------- ------
0 3718.77 1109.15 0.279452 0.083
1 61428.2 0.752658 115.521 0.094
2 61428.9 3.59998e-02 4.62970 0.107
3 61411.8 1.23668 3.77991 0.117
4 61384.1 3.32263 3.54977 0.130
5 61292.6 2.49705 1.27337 0.143
6 60996.2 13.1232 25.6338 0.156
7 60924.2 14.0042 64.1590 0.170
8 60636.9 6.12117 0.247255 0.191
9 60506.5 5.47249 58.1188 0.204
10 60208.3 4.69138 68.7643 0.215
11 59817.4 3.78799 95.7534 0.228
12 59028.5 5.06467 170.077 0.243
13 58865.1 2.04944 42.3302 0.262
14 58545.8 1.98730 73.8716 0.277
15 58436.5 0.894407 78.7006 0.297
16 57432.0 11.1264 289.960 0.308
17 56391.6 16.8890 398.427 0.318
18 55984.5 5.86100 865.482 0.335
19 55755.5 5.78590 933.245 0.348
20 55703.7 4.04436 1030.66 0.369
30 51193.8 0.532337 478.584 0.512
40 50376.7 2.66935e-02 41.2278 0.625
50 50047.0 1.06420e-02 54.3047 0.733
60 48469.0 7.75799e-02 26.9710 0.836
70 47332.9 2.28698e-02 147.696 0.945
80 46012.6 4.61655e-02 45.8507 1.058
90 43933.9 1.12477e-02 228.413 1.170
100 42816.7 0.129803 124.479 1.262
110 42391.2 0.356528 113.863 1.344
120 41919.2 0.135571 72.6286 1.437
130 41264.7 0.100665 34.7707 1.524
140 40861.6 7.49320e-03 32.2713 1.640
150 40448.2 2.48372e-02 31.5220 1.757
160 40058.1 1.53845e-03 78.3904 1.934
170 39725.1 2.11126e-02 84.4898 2.045
180 39308.0 1.93312e-02 21.7517 2.156
190 37136.5 1.23251 117.711 2.255
200 36454.6 6.98424e-03 35.2640 2.391
300 31603.1 1.83374e-04 50.0684 3.284
400 30829.6 1.57699e-07 11.6604 4.041
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 2 30829.6 FCRD -inf 4.100
3 0 30829.6 30829.6 0.00% 8.086
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.08296499547825333e+04
Final optimality gap (abs / rel) = 1.75700e-05 / 5.69906e-10 (0.00%)
# of root cutting plane rounds = 1
# of restarts = 0
# of nodes processed = 3 (4.043s)
# 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.043s)
Total program time (secs) = 8.08739 (14.420 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 / 10.068s
===========================================================================
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.
concurrent_evals 0
feastol 1e-06
feastol_abs 1e-06
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 expressions: 10262 | 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, 4e+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, 1e+00] | [1e+00, 2e+03]
Root node relaxation
--------------------
Iter Objective Feasibility Optimality Time
error error (secs)
---- --------- ----------- ---------- ------
0 1353.86 1109.12 0.255056 0.020
1 1408.83 1105.95 2.77395 0.026
2 2652.19 1034.35 7.47763 0.031
3 6845.10 776.437 11.0477 0.037
4 7481.17 735.078 11.9023 0.042
5 16440.7 113.768 26.4435 0.053
6 17832.6 3.02890 19.4946 0.061
7 17720.8 1.83506 117.490 0.070
8 17380.8 8.73919 109.398 0.078
9 17260.1 7.84965 97.2386 0.088
10 17117.0 0.216981 0.881562 0.100
11 17008.1 2.85782 6.93578 0.110
12 16790.6 1.85313 0.693047 0.122
13 16184.9 13.7071 23.8737 0.132
14 15650.1 12.1496 25.2712 0.139
15 15181.2 0.494600 28.1361 0.150
16 14728.0 1.05001 0.783843 0.161
17 14411.0 0.362034 0.924467 0.173
18 14181.4 0.121771 1.06941 0.183
19 13390.7 0.556190 66.0772 0.191
20 13180.7 0.556865 24.4949 0.199
30 12653.3 0.766207 16.0930 0.296
40 12338.4 0.201530 4.75078 0.395
50 11919.5 0.677553 36.7197 0.495
60 11253.2 1.54308 70.0946 0.592
70 10735.9 0.103027 6.92797 0.693
80 10003.9 1.92900e-02 21.3788 0.796
90 9136.95 0.372813 31.5789 0.903
100 8425.64 0.267365 23.1098 1.007
110 7825.28 3.01980e-03 30.4255 1.116
120 7337.86 6.46777e-02 33.8307 1.223
130 6880.39 7.13774e-02 62.5398 1.325
140 6714.58 2.61975e-02 72.6677 1.427
150 6389.23 4.47745e-02 18.5179 1.528
160 6229.16 7.31402e-03 17.7292 1.628
170 5962.68 1.14265e-02 24.1955 1.743
180 5775.81 3.12063e-04 26.7308 1.842
190 5495.03 2.34993e-02 5.66067 1.952
200 5307.96 4.99778e-02 5.22593 2.055
300 4252.44 1.59170e-04 3.38690 3.198
400 3224.25 0.532993 32.7134 4.281
500 3097.29 0.393481 7.80063 5.208
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 0 2216.00 LEAF 2216.00 0.00% 5.468
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.21599548133198414e+03
Final bound value = 2.21599548133198414e+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.462s)
# 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.462s)
Total program time (secs) = 5.46850 (6.107 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.613s
===========================================================================
Show the full outputHide the full output
report_clusters(positions_2, unbalanced)cluster observations centre x centre y distance share
--------------------------------------------------------------
1 12 4.28 -3.10 363.83 16.4%
2 9 21.96 -0.36 346.84 15.7%
3 12 16.32 14.02 233.73 10.5%
4 10 3.47 14.18 107.94 4.9%
5 4 21.25 30.16 218.62 9.9%
6 5 -9.11 3.37 206.05 9.3%
7 11 2.29 6.02 184.16 8.3%
8 8 11.95 7.85 84.15 3.8%
9 11 7.57 23.11 150.44 6.8%
10 8 28.42 12.19 320.23 14.5%
Total squared distance 2216.00, cluster sizes 4 to 12draw_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)=======================================
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
feastol 1e-06
feastol_abs 1e-06
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 expressions: 10262 | 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, 4e+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, 1e+01] | [1e+00, 2e+03]
Root node relaxation
--------------------
Iter Objective Feasibility Optimality Time
error error (secs)
---- --------- ----------- ---------- ------
0 1290.90 1109.12 0.259250 0.032
1 1317.08 1107.51 1.34656 0.047
2 4567.51 905.105 12.4905 0.059
3 8146.36 669.513 26.3210 0.069
4 9577.93 570.390 32.6265 0.078
5 12411.4 384.825 15.7704 0.092
6 17787.7 0.186804 7.79838 0.103
7 17669.2 0.808774 3.70947 0.113
8 17517.9 0.724213 5.80388 0.126
9 17392.2 0.335118 0.305357 0.143
10 17290.1 0.230633 3.28107 0.156
11 17074.1 0.935114 4.89089 0.166
12 16928.1 0.710185 5.50899 0.179
13 16781.0 0.761296 10.8371 0.193
14 16247.8 2.47114 5.63185 0.202
15 15609.7 3.57355 8.46401 0.216
16 14968.3 4.25448 19.1055 0.229
17 13993.0 5.08790 33.5284 0.243
18 13117.5 5.26773 27.2724 0.256
19 12375.7 5.00643 38.8465 0.270
20 12165.2 4.83416 38.8589 0.283
30 9968.46 4.43857 28.0790 0.420
40 8717.25 3.12368 17.4706 0.555
50 7777.59 1.79943 14.5063 0.705
60 6765.33 0.869958 7.20423 0.851
70 6059.06 3.61800e-03 8.77838 0.993
80 5515.16 2.08992e-04 14.2629 1.138
90 5123.48 1.00147e-02 15.8843 1.275
100 4537.94 6.40318e-03 18.6595 1.416
110 4198.10 0.182432 9.07008 1.558
120 3838.68 1.55369e-02 7.01654 1.695
130 3639.15 6.25700e-02 12.3614 1.835
140 3367.10 0.375746 18.5609 1.975
150 3144.98 0.390591 9.83157 2.111
160 3040.63 0.464382 14.9895 2.236
170 2952.00 0.243887 26.3686 2.358
180 2856.68 0.185019 25.5800 2.498
190 2679.96 7.19675e-04 19.9252 2.631
200 2502.63 0.249926 9.08155 2.765
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 2 2440.50 FCRD -inf 2.945
3 0 2440.50 2440.50 -0.00% 3.769
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.44050077279231937e+03
Final bound value = 2.44050086853016546e+03
Final optimality gap (abs / rel) = -9.57378e-05 / -3.92288e-08 (-0.00%)
# of root cutting plane rounds = 1
# of restarts = 0
# of nodes processed = 3 (2.924s)
# 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 (2.924s)
Total program time (secs) = 3.76966 (5.112 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.004s
MPEC heuristic = 0 / 0 / 0.000s
Local search heuristic = 6 / 0 / 1.726s
===========================================================================
Show the full outputHide the full output
report_clusters(positions_2, balanced)cluster observations centre x centre y distance share
--------------------------------------------------------------
1 8 27.57 10.22 240.86 9.9%
2 8 20.86 -2.12 332.24 13.6%
3 8 18.04 10.06 108.93 4.5%
4 8 11.20 14.07 62.54 2.6%
5 10 7.08 23.17 123.30 5.1%
6 9 2.12 13.74 91.66 3.8%
7 11 1.94 -2.79 253.39 10.4%
8 8 -5.82 5.25 403.75 16.5%
9 8 20.62 25.25 608.19 24.9%
10 12 8.29 5.11 215.64 8.8%
Total squared distance 2440.50, 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.