# JuMP kmeans with constraints

**URL:** <https://discourse.julialang.org/t/jump-kmeans-with-constraints/89906>\
**Category:** Optimization (Mathematical)\
**Tags:** question, jump, optimization\
**Created:** [November 7, 2022, 8:19pm UTC](https://discourse.julialang.org/t/jump-kmeans-with-constraints/89906 "2022-11-07T20:19:45Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![monty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/monty/32/44234_2.png) [@monty](https://discourse.julialang.org/u/monty)\
**Post date:** [November 7, 2022, 8:19pm UTC](https://discourse.julialang.org/t/jump-kmeans-with-constraints/89906/1 "2022-11-07T20:19:45Z")

</div>

I am solving k-means using JuMP.

```julia

X = [
	2.0 2.0;
	2.5 2.1;
	7.0 7.0;
	2.2 2.3;
	6.8 7.0;
	7.2 7.5;]
k=2
m, n = size(X)
W = exp.(-Distances.pairwise(Distances.Euclidean(), X'))
I = Matrix(1.0 * LinearAlgebra.I, m, m)
model = JuMP.Model(SCS.Optimizer)
JuMP.unset_silent(model)
@JuMP.variable(model, Z[1:m, 1:m], PSD)
@JuMP.constraint(model, Z .>= 0)
@JuMP.constraint(model, sum(Z; dims=2) .== ones(m))
@JuMP.constraint(model, LinearAlgebra.tr(Z) == k)
@JuMP.objective(model, Min, LinearAlgebra.tr(W * (I - Z)))
JuMP.optimize!(model)
Z_value = JuMP.value.(Z)

```

Everything works great.

However, I want to impose constraints based on `Z` value.

I wanted to add something like this:

```julia
JuMP.register(model, :proc_matrix, m * m, proc_matrix; autodiff = true)
@JuMP.variable(model, distances[1:m, 1:m] <= 0.1)
@JuMP.constraint(model, distances .== proc_matrix(Z...))

function proc_matrix(a...)
	m = round(Int, floor(sqrt(length(a))))
	M = reshape([a...], m, m)
	Distances.pairwise(Distances.Euclidean(), M)
end

```

However, even this does not work.

Any suggestions?!

I also want to make `proc_matrix` and the constraint more complicated but first I want to figure out how to process a variable matrix and return external output back to JuMP. Thank you!

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [November 7, 2022, 9:35pm UTC](https://discourse.julialang.org/t/jump-kmeans-with-constraints/89906/2 "2022-11-07T21:35:51Z")

</div>

A couple of quick points:

- `JuMP.register` is for registering user-defined nonlinear functions. Currently, you must use them in `@NL...` macros (but we’re working on fixing this)
- `SCS` is a conic optimizer; it does not support arbitrary nonlinear functions.
- User-defined functions must return a scalar; they cannot return a matrix

If you want to keep using a conic optimizer, you’ll have to reformulate the problem into a conic program. If you want to impose arbitrary nonlinear restrictions on the distances, you’ll need a nonlinear solver that also supports PSD constraints, and then I think you’re out of luck?

Assuming I didn’t make any mistakes, I think you’re after something like this:

```julia
# JuMP.register(model, :proc_matrix, m * m, proc_matrix; autodiff = true)
# @JuMP.variable(model, distances[1:m, 1:m] <= 0.1)
# @JuMP.constraint(model, distances .== proc_matrix(Z...))
# function proc_matrix(a...)
# m = round(Int, floor(sqrt(length(a))))
# M = reshape([a...], m, m)
# Distances.pairwise(Distances.Euclidean(), M)
# end

# is equivalent to
# @constraint(
# model, 
# [row in 1:m, col in 1:m], 
# Distances.evaluate(Distances.Euclidean(), Z[row, :], Z[: col]) <= 0.1,
# )

# is equivalent to
# JuMP.@constraint(
# model,
# [row in 1:m, col in 1:m],
# sqrt(sum((Z[row, i] - Z[i, col])^2 for i in 1:m)) <= 0.1
# )

# is equivalent to
JuMP.@constraint(
	model,
	[row in 1:m, col in 1:m],
	[0.1; Z[row, :] .- Z[:, col]] in JuMP.SecondOrderCone(),
)

```

---

<div class="post-metadata">

**Author:** ![monty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/monty/32/44234_2.png) [@monty](https://discourse.julialang.org/u/monty)\
**Post date:** [November 8, 2022, 2:40pm UTC](https://discourse.julialang.org/t/jump-kmeans-with-constraints/89906/3 "2022-11-08T14:40:00Z")

</div>

@odow this is a great example. Now I understand better how it works.

Ultimately I want to get based on the distances an estimate of the probability of each point being associated with a given cluster. There are 2 clusters in this case. I want to create a 6x2 matrix with the probability of each point belonging to each cluster. I expect the sum of the rows of this matrix to add to 1.

I know how to code this outside jump based on the Z values. However, I do not know how to implement this inside jump. i want to add penalties based on the values of this 6x2 matrix.

Are there tricks in SCS or JuMP to do this? I checked the documentation and the source code of SCS, and I could find a clear solution. Thank you!

---

<div class="post-metadata">

**Author:** ![monty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/monty/32/44234_2.png) [@monty](https://discourse.julialang.org/u/monty)\
**Post date:** [November 8, 2022, 3:32pm UTC](https://discourse.julialang.org/t/jump-kmeans-with-constraints/89906/4 "2022-11-08T15:32:04Z")

</div>

I want to be able to do something like this:

```julia
X = [
	2.0 2.0;
	2.5 2.1;
	7.0 7.0;
	2.2 2.3;
	6.8 7.0;
	7.2 7.5;]
C = [10, 5, 5, 10, 1, 1]
k=2
m, n = size(X)
W = exp.(-Distances.pairwise(Distances.Euclidean(), X'))
I = Matrix(1.0 * LinearAlgebra.I, m, m)
model = JuMP.Model(SCS.Optimizer)
JuMP.unset_silent(model)
@JuMP.variable(model, Z[1:m, 1:m], PSD)
@JuMP.constraint(model, Z .>= 0)
@JuMP.constraint(model, sum(Z; dims=2) .== ones(m))
@JuMP.constraint(model, LinearAlgebra.tr(Z) == k)
@JuMP.variable(model, capacity[1:k] == 10)
@JuMP.constraint(model, sum(unique(abs.(trunc.(Z; digits = 1)) .* C; dims=2); dims=1) .== capacity)
@JuMP.objective(model, Min, LinearAlgebra.tr(W * (I - Z)))
JuMP.optimize!(model)

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [November 8, 2022, 8:03pm UTC](https://discourse.julialang.org/t/jump-kmeans-with-constraints/89906/5 "2022-11-08T20:03:47Z")

</div>

I don’t know if you’re on the right track. There _are_ some reformulations that could achieve what you’re looking for, but they get complicated. And they introduce binary and integer variables, so now you’ll need a mixed-integer conic solver, and the only mixed-integer PSD solver in JuMP’s [list of supported solvers](https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers) is Pajarito.

Stepping back, why solve K-means with a conic solver? What are the alternatives?

`unique` and `trunc` give you trouble. The key for `trunc` with `digits = 1` is to see that `Z[i,j]` is between `0` and `1`, so that the trunc will be an integer between `0` and `9` (followed a division by 10).  
So that gives:

```julia
JuMP.@variable(model, Z_trunc[1:m, 1:m] >= 0, Symmetric, Int)
JuMP.@constraint(model. Z_trunc .<= 10 * Z .<= Z_trunc + 1)
JuMP.@constraint(model, sum(unique(abs.(Z_trunc / 10) .* C; dims=2); dims=1) .== capacity)

```

then we can drop the `abs` (why is it needed?) because `Z_trunc` cannot be negative.

```julia
JuMP.@constraint(model, sum(unique(Z_trunc / 10 .* C; dims=2); dims=1) .== capacity)

```

I’m a but confused by the usage of `unique`. I invented some fake data:

```nohighlight
julia> Z_trunc
6×6 Matrix{Float64}:
 0.2 0.1 0.0 0.0 0.1 0.2
 0.2 0.0 0.2 0.2 0.1 0.0
 0.0 0.2 0.1 0.1 0.1 0.2
 0.1 0.1 0.2 0.1 0.0 0.2
 0.2 0.1 0.1 0.0 0.2 0.1
 0.2 0.0 0.2 0.2 0.1 0.1

julia> unique(Z_trunc, dims = 2)
6×6 Matrix{Float64}:
 0.2 0.1 0.0 0.0 0.1 0.2
 0.2 0.0 0.2 0.2 0.1 0.0
 0.0 0.2 0.1 0.1 0.1 0.2
 0.1 0.1 0.2 0.1 0.0 0.2
 0.2 0.1 0.1 0.0 0.2 0.1
 0.2 0.0 0.2 0.2 0.1 0.1

```

Did you want to count the number of unique entries in each row?

If we omit the `unique`, then the last part is:

```julia
Y = Z_trunc .* C / 10
JuMP.@constraint(model, sum(Y; dims=1) .== 10)

```

Puttig it all together

```julia
JuMP.@variable(model, Z_trunc[1:m, 1:m] >= 0, Symmetric, Int)
JuMP.@constraint(model. Z_trunc .<= 10 * Z .<= Z_trunc + 1)
JuMP.@expression(model, Y, Z_trunc .* C / 10)
JuMP.@constraint(model, sum(Y; dims=1) .== 10)

```