Python sequence cluster exercise

Question

I am working through an exercise in my textbook and implementing the code in Python to practice dynamic programming. I feel like I am right on the edge of figuring it out, but after many hours, I come here for help.

Basically, my code is going through a list of values x, and given a k, breaks that list into k clusters based on calculating the minimum sum of squared errors (SSE) for a particular cluster.

The code creates a table, calculating the SSE for 1 cluster, 2 clusters, ..., k clusters, if we were to put the cluster parentheses around all variations of values within list[0:1], list[0:2], list[0:3], ..., list[0:n], and choosing the minimum SSE for that particular step in the table.

For example:

Given x= [7,6,9,15,18,17,30,28,29] and k=3

we would return clusters (7,6,9)(15,18,17)(30,28,29), which would translate to sum of squared error equal to (4.666)(4.666)(2) for each cluster. So our max SSE would be 4.666 for that clustering on that list.

Now when I try it on my second list x = [52, 101, 103, 101, 6, 5, 7], I should get clustering (52)(101, 103, 101)(6, 5, 7), which should give (0)(2.666)(2) or a max of 2.666, but am not getting that. I believe the error lives in the def f(s, j_down, t) for the 2nd return statement, and how I increment s and t. Hopefully, I have not made a silly mistake!

def mean(numbers):
    return float(sum(numbers)) / max(len(numbers), 1)

def sum_square(x):
    if isinstance(x, (int,)):
        return 0
    w = 0
    for i in x:
        w += (i - mean(x))**2
    return w

def f(s, j_down, t):
    if not r[s][j_down] and r[s][j_down] != 0:
        return sum_square(x[:t - s])

    return max(r[s][j_down], sum_square(x[:t-s]))

def get_min_f_and_s(j_down, t):
    """ range s from 1 to t-1 and set s to minimize f(s)
    """
    items = [(s, f(s, j_down, t)) for s in range(t)]
    s, min_f = min(items, key=lambda x:x[1])
    return s, min_f

def seq_out(n,k):
    for j in range(k):
        if j == 0:
            for t in range(n):
                r[t][j] = sum_square(x[:t+1])

                c[t][j] = x[:t+1]
        else:
            for t in range(1, n):
                s, min_f = get_min_f_and_s(j - 1, t)
                r[t][j] = min_f
                c[t][j] = [c[s][j - 1]] + x[s+1:t+1]

    print('the max SSE is: {}'.format(r[-1][-1]))
    print('the cluster centers are: {}'.format(c[-1][-1]))

#x = [7,6,9,15,18,17,30,28,29]    
x = [52, 101, 103, 101, 6, 5, 7]
k = 3
n = len(x)

r = [[[] for _ in range(k)] for _ in range(n)]
c = [[[] for _ in range(k)] for _ in range(n)]

print(seq_out(n,k))
print(r)
print(c)

Edit: Question layout

Given a sequence X = [x_1, x_2, ... x_n] and integer k > 1, partition X into clusters C_1,..., C_k of sizes n_1, ..., n_k, so that the sum of squared errors is minimized.

Answer 1

Until you define what recursion you used, it's hard to really be sure what you're trying to implement. I took a look at the fundamental problem you are solving, though, and defined a recurrence that will work. I assume here we have some sequence $d(i)$ with $n$ components that starts with indexing at $1$. Running with this, the minimum cost clustering such that you get exactly $k$ clusters can be expressed with the recurrence

\begin{align} f(i,j) &= \begin{cases} \infty & (i=0 \land j \neq k+1) \lor (i \neq 0 \land j = k+1) \\ 0 & (i = 0 \land j = k+1) \\ \text{cost}(1,i) & j = k\\ \min_{1 \leq l \leq i} f(l-1, j+1) + \text{cost}(l,i) & \text{otherwise} \end{cases} \end{align}

where the $\text{cost}(i,j)$ cost function returns the cluster cost of viewing elements $i$ to $j$ in the sequence as a cluster. We can define the cluster cost as

\begin{align} \text{cost}(i,j) &= \sum_{s = i}^j \left(d(s) - \mu(i,j)\right)^2 \\ \text{where }\mu(i,j) &= \frac{1}{j - i + 1}\sum_{t = i}^j d(t) \end{align}

Using these definitions, one can implement the recurrence however they want and to find $f$. The simplest is a recursion with memoization built in but I tend to look for where the know elements are in $f(i,j)$ off the bat and using those to fill in the other quantities more efficiently. Regardless of how you do it, you can fill in $f$ and then use the table to find the optimal $k$ clusters by using $\arg \min_{1 \leq l \leq i} f(l-1, j+1) + \text{cost}(l,i)$ with the starting indices $i = n$ and $j = 1$.

I wrote a code in Python that implements this solution, which you can compare against if desired

#import useful libs
import numpy as np

# problem is found here: 
# https://scicomp.stackexchange.com/questions/31244/python-sequence-cluster-exercise

# define efficient representation of cluster cost
class cluster_cost:
    def __init__(self, data):
        self.n = len(data)
        self.ss_x = np.cumsum(data)    # subset sum of x
        self.ss_x2= np.cumsum(data**2) # subset sum of x^2

    def cost(self, i, j):
        if i == 0:
            dx = self.ss_x[j]
            dx2= self.ss_x2[j]
        else:
            dx = self.ss_x[j] - self.ss_x[i-1]
            dx2= self.ss_x2[j] - self.ss_x2[i-1]
        mean = dx/(j-i+1)
        return dx2  + (j-i+1)*mean**2 - 2.0*mean*dx



def dp_soln(data, k):

    # define cluster cost function
    cc = cluster_cost(data)

    # get number of elements in data
    n = cc.n

    # init the DP cost table
    f = np.full((n+1,k+2), np.inf)

    # construct f[i,j] for all i and j
    # init for j = k+1
    f[0,k+1] = 0

    # init for j <= k
    for j in range(k, 0, -1):
        for i in range(1,n+1):

            # get the min cost option
            best = np.inf
            for l in range(1, i+1):
                best = np.min([best, f[l-1,j+1] + cc.cost(l-1,i-1)])

            # set best cost to DP table
            f[i,j] = best

    # solve for the optimal clustering
    # using cost table from DP algo
    partition_indices = []
    partition = []

    # set initial state indices to (n,1)
    i = n
    j = 1

    # loop until k clusters found
    while j <= k:

        # find argmin that matches DP value
        best = np.inf
        best_idx = 1
        for l in range(1, i+1):
            bestn = np.abs( f[i,j] - f[l-1,j+1] - cc.cost(l-1,i-1) )
            if( bestn < best ):
                best_idx = l
                best = bestn

        # get the next partition based on argmin
        partition_indices.append((i-1, best_idx-1))
        partition.append(data[(best_idx-1):i])

        # update indices
        i = best_idx - 1
        j = j + 1

    # return the partition indices and partition
    return (partition_indices, partition)

if __name__ == "__main__":
    k = 3
    data = np.array([7, 6, 9, 15, 18, 17, 30, 28, 29])
    (indices, part) = dp_soln(data, k)
    print(indices)
    print(part)

The above sample $k$ and $\text{data}$ return the expected output when run:

[(8, 6), (5, 3), (2, 0)]
[array([30, 28, 29]), array([15, 18, 17]), array([7, 6, 9])]