# Full segmentation of a linear space using a clustering algorithm

Once again, I am not entirely sure how to describe what I am looking for, hence I have a hard time finding answers using Google or any other literal search method.

Let's say I have a time series:

$$(x_1, y_1), (x_2, y_2), (x_3, y_3), (x_4, y_4),\ldots, (x_n, y_n)$$

Let's say I also have 2 key values: $V_1$ and $V_2$

I want to find a continuous decomposition of the time series such that each segment's average of its $y$ values is either closest to $V_1$ or to $V_2$.

Here is a trivial example to illustrate the point:

Time series:

$$(0, 1), (1, 1), (2, 1), (3, 4), (4, 4), (5, 4), (6, 4)\\ V_1 = 1,\quad V_2 = 4$$

The best continuous segmentation of the time series is:

$$s_1=\begin{bmatrix}(0, 1) & (1, 1) & (2, 1)\end{bmatrix}\\ s_2=\begin{bmatrix}(3, 4) & (4, 4) & (5, 4)&(6, 4)\end{bmatrix}$$

Since the average value of $s_1$'s $y$ values is $1$ and the "distance" between $1$ and $V_1$ is zero. Similarly, the average value of $s_2$'s $y$ values is $4$ and the "distance" between $4$ and $V_2$ is zero.

Generalizing, each segment would either be determined to be closer to $V_1$ or closer to $V_2$, the algorithm would minimize the sum of distances (between each segment's average $y$ values and either $V_1$ or $V_2$) and the sequence of segments would equal the time series without gaps or overlaps.

Ouch!

Is that even comprehensible?

Can you think of an algorithm making that possible? (regardless of complexity for now)

Segmentation could include as many segments as necessary to minimize the overall distance. Another trivial example to illustrate this:

$$(0, 1), (1, 1), (2, 1), (3, 4), (4, 4), (5, 4), (6, 4), (7, 1), (8, 1)$$ $$V_1 = 1 \quad V_2 = 4$$

In this case, the following segmentation is perfect:

$$s_1=\begin{bmatrix}(0, 1) & (1, 1) & (2, 1)\end{bmatrix}\\ s_2=\begin{bmatrix}(3, 4) & (4, 4) & (5, 4) & (6, 4)\end{bmatrix}\\ s_3=\begin{bmatrix}(7, 1) & (8, 1)\end{bmatrix}$$

The overall distance between each segment's average value and its proximal key is zero.

• would you specify the number of segments you will chop it into as an input? – spektr Dec 15 '16 at 3:04
• No, the number of segments would not be an input. – L_R_T Dec 15 '16 at 3:19
• It may be an optimal stopping problem. One could start with 1 segment (the entire series) and compute the distance to V1 and V2. Then figure out if the average distance can be reduced with 2 segments. Then, figure out if the average distance can be reduced with 3 segments, etc. – L_R_T Dec 15 '16 at 18:07
• That sounds like a reasonable approach. I was thinking about how it could be formulated as a sort of integer optimization problem, but that approach would better handle the implicit desire to minimize the number of segments used. – spektr Dec 15 '16 at 18:14
• @choward I still need to wrap my head around whether I can assume that if the average distance was not improved with n segments, it won't be improved with n + i segments. Workout through some examples. It does feel like a linear optimization problem, but with a hierarchical twist. – L_R_T Dec 15 '16 at 19:06

It's been a while since this question was posted and looking at it again, I have some new insights into tackling the problem. Let us first define the set $$V = \lbrace V_1, V_2, \cdots, V_k \rbrace$$ as the set of possible key values we care to consider.

Let us also assume the time series data is sorted in time so we can model it as some sequence $$D(n): \mathbb{N}^{+} \rightarrow \mathbb{R}$$ where we have $$M$$ pieces of data. Let us also compute the prefix sum of the time series and denote it $$p(n) = \sum_{k=1}^n D(k) = p(n-1) + D(n)$$, with $$p(1) = D(1)$$, which can be computed in $$O(M)$$ time.

We will now compute the minimum error, using the key set, for some subset of the data ranging from the indices $$i$$ to $$j$$ in the following manner:

\begin{align} e(i,j) &= \min_{v \in V} d\left(p(j) - p(i) + D(i), (j-i+1) v \right) \end{align}

where $$d(u,v)$$ is some metric, perhaps $$d(u,v) = |u-v|$$. This above computation must be done for all $$1 \leq i < j$$, so this computation is ultimately $$O\left(k M^2\right)$$. With this information, we then wish to try and segment the data in a manner that minimizes the error across each segment of the data.

Let us define some positive scalar $$\tau \in \mathbb{R}$$ that we will use to softly penalize each segment we generate, meaning that the larger $$\tau$$ is, the less segments we will strive to have to balance that penalty with the error of the segments. Using this, we can write the following recursive relationship:

\begin{align} f(j) = \min_{1 \leq i \leq j} \left(e(i,j) + \tau + f(i-1)\right) \end{align}

where $$f(0) = 0$$. It is trivial to see we can find $$f(j)$$ in increasing order of $$j=1$$ to $$j=M$$ where each minimization will take $$O(M)$$ computation, making the computation of $$f(j)$$ an $$O\left(M^2\right)$$ process.

To back out the segment indices and keys associated with each segment, we can use the following algorithm:

\begin{align} &S = \emptyset\\ &j^* \leftarrow M\\ &\textbf{while}(j^* \geq 1)\\ &\hspace{0.5cm} \cdot i^* \leftarrow \arg \min_{1 \leq i \leq j^*} \left(e(i,j^*) + \tau + f(i-1)\right)\\ &\hspace{0.5cm} \cdot v^* \leftarrow \arg \min_{v \in V} d\left(p(j^*) - p(i^*) + D(i^*), (j^* - i^* + 1) v \right)\\ &\hspace{0.5cm} \cdot S \leftarrow S \cup \left\lbrace \left(i^*,\; j^*,\; v^*\right) \right\rbrace\\ &\hspace{0.5cm} \cdot j^* \leftarrow \left(i^* - 1\right) \\ \end{align}

which is also at worst case $$O\left(M^2 + k M\right)$$. Looking at the overall algorithm, we can see that the time complexity is $$O \left(k M^2 \right)$$ and the space complexity is $$O \left( M^2 \right)$$. In this algorithm, $$\tau$$ can be tuned to help impact how large of segments you'll allow, though this may take some experimentation.

To validate that this algorithm described above works, below is a code to perform the computation in C++. The primary source code for the algorithm is shown below. I should note that to minimize space requirements down to $$O(M)$$, I slightly modified how some of the above quantities are stored and computed such that the time complexity is still $$O\left(k M^2 \right)$$.

#include <iostream>
#include <stdint.h>
#include <vector>
#include <cmath>
#include <limits.h>

namespace opt {

struct segment {
segment(size_t s, size_t e, double k):start(s),end(e),key(k) {}
size_t start, end;
double key;
void print() const {
std::cout << "(" << start << ", " << end << ", " << key << ")" << std::endl;
}
};

using segment_list = std::vector<struct segment>;

void comp_prefixsum(const std::vector<double>& D,
std::vector<double>& p)
{
p[0] = D[0];
for(size_t i = 1; i < D.size(); ++i){ p[i] = p[i-1] + D[i]; }
}

double err_metric(double u, double v){
return std::abs(u - v);
}

double min_error(size_t i, size_t j,
const std::vector<double>& D,
const std::vector<double>& p,
const std::vector<double>& V )
{
double min_err = std::numeric_limits<double>::max();
for(auto v: V){
double err = err_metric(p[j] - p[i] + D[i],
static_cast<double>(j-i+1)*v );
if( err < min_err ){ min_err = err; }
}// end for v

return min_err;
}

double argmin_error(size_t i, size_t j,
const std::vector<double>& D,
const std::vector<double>& p,
const std::vector<double>& V )
{
double argmin_key, min_err = std::numeric_limits<double>::max();
for(auto v: V){
double err = err_metric(p[j] - p[i] + D[i],
static_cast<double>(j-i+1)*v );
if( err < min_err ){ min_err = err; argmin_key = v; }
}// end for v

return argmin_key;

}

void segmentation(const std::vector<double>& time_series,
const std::vector<double>& key_set,
double tau,
std::vector<segment>& soln_segmentation )
{
// compute size of data
size_t M = time_series.size();

// allocate data arrays
std::vector<double> p(M,0.0), f(M,0.0);
std::vector<size_t> arg_fmin(M, 0);

// compute prefix sum
comp_prefixsum(time_series, p);

// compute sub problem values
for(size_t j = 0; j < M; ++j){
double ft, f_min = std::numeric_limits<double>::max();
for(size_t i = j; i != -1; --i){

// compute minimum key-based error
double e_ij = min_error(i, j, time_series, p, key_set);

// compute recurrence value
if( i != 0 ){   ft = e_ij + tau + f[i-1]; }
else{           ft = e_ij + tau; }

// update minimum recurrence value, if necessary
// and the argmin index of this optimization
if( ft <= f_min ){
f_min = ft;
f[j]  = f_min;
arg_fmin[j] = i;
}

}// end for i
}// end for j

// compute segmentation solution
soln_segmentation.clear();
soln_segmentation.reserve(M);
size_t j = M-1;
while(j != -1){
size_t i_opt = arg_fmin[j];
double v_opt = argmin_error(i_opt, j, time_series, p, key_set);
soln_segmentation.emplace_back(i_opt, j, v_opt);
j = i_opt - 1;
}// end while

}
} // end namespace opt


A sample code to test that the above algorithms work is the following:

int main(int argc, const char** argv) {

// initialize the data
std::vector<double> data { 6, 6, 4, 4, 4, 4, 4, 1, 1, 1, 1 };
std::vector<double> keys { 1, 4, 6 };

// declare the output segment solution list
typename opt::segment_list segments;

// define the tuning tau variable
double tau = 0;

// compute the optimal segmentation
opt::segmentation(data, keys, tau, segments);

// print the segmentation results based on
// 3-tuple defined in mathematical formulation
std::cout << "Print segmentation 3-tuples:" << std::endl;
size_t num_segs = segments.size();
for(size_t i = (num_segs-1); i != -1; --i){
segments[i].print();
}
std::cout << std::endl;

// print the actual segmentations
std::cout << "Print actual segmentations:" << std::endl;
size_t count = 0;
for(size_t i = (num_segs-1); i != -1; --i){
std::cout << "s" << count << " = [ ";
for(size_t j = segments[i].start; j <= segments[i].end; ++j){
std::cout << "(" << j << ", " << data[j] << ") ";
}
std::cout << " ]" << std::endl;
++count;
}

// finish work
return 0;
}


This code prints the $$(i,j,v)$$ tuples and their associated segmentations of the actual data. The above example outputs the following:

Print segmentation 3-tuples:
(0, 1, 6)
(2, 6, 4)
(7, 10, 1)

Print actual segmentations:
s0 = [ (0, 6) (1, 6)  ]
s1 = [ (2, 4) (3, 4) (4, 4) (5, 4) (6, 4)  ]
s2 = [ (7, 1) (8, 1) (9, 1) (10, 1)  ]


Thus, this algorithm is correct and runs in polynomial time!

• Fantastic! Still validating your solution, but I think you nailed it. – L_R_T Mar 2 '19 at 18:09

This sounds hard when you ask it like that, but it actually isn't. Due to the continuity, one could always do a linear search over the whole array. At each iteration a new value from the sequence is picked and the means of both segments are updated. Because the mean could drastically change by inclusion/exclusion of one element, we cannot early stop and have to traverse the whole array. Therefore, we exhaustively iterate over all the elements (all possible continuous segmentations) and compute the mean and the distance from $V_1$ and $V_2$. At the end, we take the one which minimizes the sum of these distances.

Because running mean could be computed in $O(1)$ time, the computational cost of this is $O(N)$. Please correct me if I'm wrong.

• Thanks! I am dead tired at the moment, but I did not mean to imply the segmentation should be in 2 segments only it could in fact in n segments if it yields a lower overall distance between segments and their respective values. I'll add another trivial example to illustrate this. – L_R_T Dec 15 '16 at 2:41
• Doesn't matter. In that case you would also need to loop over all the segments, rendering the complexity to be $O(MN)$. – Tolga Birdal Dec 15 '16 at 15:39
• If you have a large data set and m segments, you also have to find optimal size of each segment. This is thousands of combinations to be checked at every iteration since if you want to check whether m+1 segments is better, you have to do it again, and you don't know when to stop unless you know that for some combination the sum of distances would be zero, but I guess in general case you can't assume that. My point is, this approach should work, but I think it's factorial complexity, since you have (n - m) * (n- m - 1) * ... * n combinations to check in each iteration, where n is data set size. – KjMag Aug 3 '17 at 7:27

some times ago I have tried to find an algorithm on a similar problem: find the optimal partitioning of a time series by performing linear regression on each segments, without prior knowledge of the optimal number of segments and their location.

I had searched for papers with keywords as 'Multiple Structural Change', or 'Time serie segmentation'...

I think the main idea is to consider this as a minimal distance problem on a directed network. Each point in time is a node, which is connected to all its successor nodes. We look for the shortest path from the first to the last time. The measure for the "distance" is, I think, the average error to the fit (constant in your case).

The difficulty is maybe to not over-fit the data: i.e. to end up with as many segments as data points. The cost of creating a new segment had to be taken into account. I am not sure about how to do this.

At the time, I didn't look further so I am curious about the right way to do this,

Best regards,

• You can check my answer to see how you can handle this problem. – spektr Dec 25 '18 at 22:17