2 8 4 Segmentation With Dynamic Programming

We consider the linear model $Y=f^{*}+\epsilon=X\beta^{*}+ \epsilon$ where the $f^{*}_{i}$ are piecewise constants. It corresponds to the variation sparse setting with $p=n$ and $X_{1},...,X_{n}$ the canonical basis of $R^{n}$.

We want to estimate $f^{*}$ by model selection. We will use the family of models $\mathscr{M}=\mathscr{P}([\![1,p]\!] )$, $\{ S_m, m \in \mathscr{M} \}$ with $S_m$ define by $S_m = span(\sum_{j\geq k} X_j, k \in m)$ and $S_m$ gather all the vectors $X\beta$ with $\beta$ fulfilling $\beta_j = \beta_{j-1}$ if $j \neq m$.

As this family is indexed by the $2^n$ subsets of $\![1,...,n]\!$ we can't use a naive minimization of the Criterion $\hat{m} \in \underset{m \in \mathscr{M}}{argmin} \{ \|Y - \hat{f_m}\| + \sigma^{2}pen(m) \}$ . Yet, we can minimize this Criterion much more efficiently thanks to dynamic programing.

We start by introducing some terms we will use in this demonstration.

For $2 \leq k \leq n$

(1)

\begin{align} N^{2}_{k} = \sum_{i=1}^{k-1} Y_{i}^{2}, \quad N^{2}_{1}=0 \quad for \quad 2 \leq k \leq n \end{align}

(2)

\begin{align} \bar{y}_{(k,j)} = \frac{1}{j-k} \sum_{i=k}^{j-1} Y_i, \quad for \quad 1\leq k\leq n+1 \end{align}

and

(3)

\begin{align} R^{2}(k,j)=\sum_{i=k}^{j-1} (Y_i - \bar{y}_{(k,j)})^{2}, \quad for \quad 1\leq k\leq n+1. \end{align}

1. Given $m = \!\{i_1,...,i_D\}\! \subset \!\{1,...,n\}\!$, we want to prove that

(4)

\begin{align} \hat{f_m} = \sum_{q=1}^{D} \bar{y}_{(i_q,i_{q+1})}1_{i_q}^{i_{q+1}} \end{align}

where $i_{D+1}=n+1$ and $1_{k}^{j} = X_k + ... + X_{j-1}$.

By definition $\hat{f_m} = Proj_{S_m}Y$ with $S_m$ as we defined earlier. We want to make an orthogonal projection of $Y$ on $S_m$ which is a subspace in which the vectors have their coordinates equal between $i_k$ and $i_{k+1}-1$ for $1\leq k\leq D$.
The $1_{i_q}^{i_{q+1}}$ for $1\leq q \leq D$ is a basis of $S_m$. It means we can write

(5)

\begin{align} \boxed{Proj_{S_m}Y = \sum_{q=1}^{D} \frac{<Y,1_{i_q}^{i_{q+1}}>}{\|1_{i_q}^{i_{q+1}}\|^{2}}1_{i_q}^{i_{q+1}}} \end{align}

We can easily calculate

(6)

\begin{align} \|1_{i_q}^{i_{q+1}}\|^{2} = i_{q+1} - i_q \end{align}

and

(7)

\begin{align} <Y,1_{i_q}^{i_{q+1}}> = \sum_{k=i_q}^{i_{q+1}-1} Y_k. \end{align}

And finally

(8)

\begin{align} \boxed{Proj_{S_m}Y = \sum_{q=1}^{D} \frac{\sum_{k=i_q}^{i_{q+1}-1} Y_k}{i_{q+1} - i_q} 1_{i_q}^{i_{q+1}}} \end{align}

Which gives us the result we wanted as $\bar{y}_{(i_q,i_{q+1})} = \frac{\sum_{k=i_q}^{i_{q+1}} Y_k}{i_{q+1} - i_q}$.

2. In the following, we'll assume that the probability $\Pi_m$ on $\mathscr{M}$ depends on $m$ only through its cardinality $|m|$. As a result we will write $pen(|m|)$ instead of $pen(m)$ to emphasize the dependence.
We then define for $0 \leq D\leq n$ :

(9)

\begin{align} \hat{m}_D = \underset{m:|m|=D}{argmin} {\|Y-\hat{f_m}\|^{2}}, \end{align}

(10)

\begin{align} C_n(D) = \|Y-\hat{f_{\hat{m}_D}}\|^{2}. \end{align}

We want to prove that the minimizer $\hat{m}$ of the Criterion mentioned earlier is given by $\hat{m} = \hat{m}_{\hat{D}}$ where

(11)

\begin{align} \hat{D} = \underset{D=0,...n}{argmin} \{C_n(D) + \sigma^{2}pen(D)\}. \end{align}

We want to find $\hat{m}$ minimizing $\underset{m}{argmin} {\|Y-\hat{f_m}\| + \sigma^{2}pen(|m|)}$.
The minimization problem is equivalent to :

(12)

\begin{align} \underset{D=0,...,n}{min} \{\underset{m : |m|=D}{min} \{\|Y-\hat{f_m}\| + \sigma^{2}pen(|m|)\}\} \end{align}

(13)

\begin{align} = \underset{D=0,...,n}{min} \{\underset{m : |m|=D}{min} \{\|Y-\hat{f_m}\| + \sigma^{2}pen(D)\}\} \end{align}

And we know

(14)

\begin{align} \underset{m : |m|=D}{argmin} \{\|Y-\hat{f_m}\| + \sigma^{2}pen(D)\} = \underset{m : |m|=D}{argmin} \{\|Y-\hat{f_m}\|\} = \hat{m}_D. \end{align}

Which means, we have to find $D'$ minimizing $\|Y-\hat{f}_{\hat{m}_D}\| + \sigma^{2}pen(D)$, we would then have $\|Y-\hat{f}_{\hat{m}_{D'}}\| + \sigma^{2}pen(D')$ the minimum of $\|Y-\hat{f}_{\hat{m}_D}\| + \sigma^{2}pen(D)$. and by definition it's $\hat{D}$

And at the end

(15)

\begin{align} \boxed{\hat{m}=\hat{m}_{\hat{D}}} \end{align}

3.We now have to prove that $C_n(D)$ and $\hat{m}_D$ are solutions of

(16)

\begin{align} C_n(D) = \underset{1\leq i_1\leq...\leq i_D\leq n}{min} {N^{2}(i_1) + \sum_{q=1}^{D} R^{2}(i_q,i_{q+1})} \end{align}

(17)

\begin{align} \hat{m}_D = \underset{1\leq i_1\leq...\leq i_D\leq n}{argmin} {N^{2}(i_1) + \sum_{q=1}^{D} R^{2}(i_q,i_{q+1})} \end{align}

By definition of $\hat{m}_D$,

(18)

\begin{align} C_n(D) = \|Y - \hat{f}_{\hat{m}_D}\|^{2} = \underset{m : |m|=D}{min} \|Y - \hat{f_m}\|^{2} \end{align}

By using the expression found for $\hat{f_m}$, and since $\underset{m : |m|=D}{min} = \underset{1\leq i_1\leq...\leq i_D\leq n}{min}$ we finally have

(19)

\begin{align} \boxed{C_n(D)} =\underset{1\leq i_1\leq...\leq i_D\leq n}{min} \{\sum_{k=1}^{i_1 – 1} Y_i^{2} + \sum_{q=1}^{D} \sum_{k=i_q}^{i_{q+1}-1} (Y_i - \bar{y}_{(i_q,i_{q+1})})^{2}\} \end{align}

(20)

\begin{align} = \boxed{\underset{1\leq i_1\leq...\leq i_D\leq n}{min} \{N^{2}(i_1) + \sum_{q=1}^{D} R^{2}(i_q,i_{q+1})\}} \end{align}

The same operations give us

(21)

\begin{align} \boxed{\hat{m}_D = \underset{1\leq i_1\leq...\leq i_D\leq n}{argmin} \{\|Y - \hat{f_m}\|\} = \underset{1\leq i_1\leq...\leq i_D\leq n}{argmin} \{N^{2}(i_1) + \sum_{q=1}^{D} R^{2}(i_q,i_{q+1})\}} \end{align}

4. Let us check that for $2\leq D\leq n$ we have

(22)

\begin{align} C_n(D) = \underset{i=D,...,n}{min} \{C_{i-1}(D-1) + R^{2}(i,n+1)\} \end{align}

For $2\leq D\leq n$ :

(23)

\begin{align} \underset{i=D,...,n}{min} \{C_{i-1}(D-1) + R^{2}(i,n+1)\} \end{align}

(24)

\begin{align} = \underset{i=D,...,n}{min} \{\underset{1\leq i_1< … <i_{D-1}\leq i-1}{min} {N^{2}(i_1) + \sum_{q=1}^{D-1}R^{2}(i_q,i_{q+1})} + R^{2}(i,n+1)\} \end{align}

(25)

\begin{align} =\underset{1\leq i_1 < ...< i_{D-1}< i\leq n}{min} \{N^{2}(i_1) + \sum_{q=1}^{D-1}R^{2}(i_q,i_{q+1}) + R^{2}(i,n+1)\} \end{align}

(26)

\begin{align} =\underset{1\leq i_1 < ...< i_{D-1}< i\leq n}{min} \{N^{2}(i_1) + \sum_{q=1}^{D}R^{2}(i_q,i_{q+1})\} =C_n(D) \end{align}

And so

(27)

\begin{align} \boxed{\underset{i=D,...,n}{min} \{C_{i-1}(D-1) + R^{2}(i,n+1)\} = C_n(D), \quad for \quad 2\leq D\leq n} \end{align}

5. We finally want to prove we can compute $\hat{m}$ with $O(n^{3})$ operations.

First, let's check that we can compute all the $N^{2}(k)$ and $R^{2}(k,j)$ for $1\leq k<j\leq n+1$ in $O(n^{3})$ operations.

To compute $N^{2}(k)$ we'll need $O(n)$ operations to have $Y^{2}$ and $O(n)$ operations to sum the coordinates of $Y^{2}$. So to compute the $N^{2}(k)$ for every $k$ we need $O(n^{2})$ operations.

With the same reasoning :
To compute $R^{2}(k,j)$ we first need to compute $\bar{y}_{(k,j)}$, which take $O(n)$ operations. To compute all the $\bar{y}_{(k,j)}$ for $1\leq k<j\leq n+1$ we need $O(n^{3})$ operations. Once we have all the $\bar{y}_{(k,j)}$ for $1\leq k<j\leq n+1$, we can compute $R^{2}(k,j)$ in $O(n)$ operations, and then compute all the $R^{2}(k,j)$ for $1\leq k<j\leq n+1$ in $O(n^{3})$.

We will now propose an algorithm building on the formula

(28)

\begin{align} C_n(D) = \underset{i=D,...,n}{min} \{C_{i-1}(D-1) + R^{2}(i,n+1)\} \end{align}

that computes $C_n(D)$ and $\hat{m}_D$ for $1\leq D \leq n$ in $O(n^{3})$ operations thanks to dynamic programing.

Algorithm :

Compute all the $N^{2}(k)$ and $R^{2}(k,j)$ for $1\leq k<j\leq n+1$ ($O(n^{3})$ operations)
Compute all the $C_k(1)$ for $1\leq k\leq n$ with :
$C = zeros(n,n)$
$m = zeros(n,n)$
for $k$ in $c(1,n)$
$P = zeros(n)$
for $i$ in $c(1,k)$
$P[i] = N^{2}[i] + R^{2}(i,k+1)$
end
$C[k,1] = min(P)$
$m[k,1] = whichmin(P)$
end
($O(n^{2})$ operations)

Compute $C[j,k]$ and $m(j,k)$ for $1\leq k < j\leq n$
for $j$ in $c(1,n)$
for $k$ in $c(1,j)$
$P = zeros(j-k+1)$
for $i$ in $c(k,j)$
$P[i] = C[i-1,k-1] + R^{2}(i,j+1)$
end
$C[j,k] = min(P)$
$m[j,k] = whichmin(P)$
***
end
end
($O(n^{3})$ operations)

Then we have the $C_n(D)$ for $1\leq D\leq n$ with $\boxed{C[n,D]}$ and the $\hat{m}_D$ for $1\leq D\leq n$ with $\boxed{m[c(1:n),D]}$.

Finally to have $\hat{m} = \hat{m}_{\hat{D}}$ we just have to add one line to the previous algorithm (instead of ***)

(29)

\begin{align} F[j,k] = C[j,k] + \sigma^{2}pen(D) \end{align}

where, for a given $K>1$

(30)

\begin{align} pen(D) = K(\sqrt{d_D}+\sqrt{2log(\frac{1}{\pi_D}})^{2}) \end{align}

with $\pi_D = (1+\frac{1}{p})^{-p}p^{-D}$ and $d_D = dim(S_D)$

Then take $whichmin(F[n,c(1:n)])$ to obtain $\hat{D}$ and return $\boxed{m[c(1:n),\hat{D}] = \hat{m}_{\hat{D}} = \hat{m}}$.