Generalized Shrinkage

1 Authors
2 Publication
3 Introduction
4 Total Variation (TV) Minimization
5 The Split Bregman Algorithm
6 Mean Curvature Flow via TV minimization
7 The basic shrinkage problem
8 Generalized Shrinkage for polyhedral norms
9 The Polyhedral Shrinkage Algorithm
10 Numerical Results for Polyhedral Shrinkage
11 Numerical Results for Mean Curvature Flow using Generalized Shrinkage

Authors

This is a joint project between

Adam Oberman Main Page
Ryo Takei Students

This project involves a key step (generalized shrinkage) which is part of a larger project with

Stan Osher people
Richard Tsai people

for anisotropic motion by mean curvature.

Write up by Adam Oberman and Ryo Takei.

Publication

This project has been completed. The resulting publication is

(with Stan Osher, Ryo Takei, and Richard Tsai) Numerical Methods for Smooth and Crystalline Mean Curvature Flow

Cam Report

Introduction

Here we explain the steps in mean curvature motion by bregman iteration and thresholding. Thresholding with a different HJ equation gives anisotropic motion by mean curvature The main ingredient is the shrinkage step, with an anisotropic norm. There are known formulas for shrinkage for the euclidian and L1 norm, we investigate shrinkage using other norms. (Can do L infinity, get rotated L1) We do polyhedral norms, and use convex analysis to work out the mapping.

Total Variation (TV) Minimization

A basic model in image denoising is the TV minimization with quadratic penalty (or simply TV minimization): given a norm $\gamma : \R^2 \rightarrow \R$ , and function $f:\R^2\rightarrow\R$ and a parameter $λ > 0$ , compute

$u^* = \arg\min_u \int \gamma(\nabla u) + \frac{\lambda}{2}\int|u - f|^2.$

For the isotropic case $\gamma(\cdot) = ||\cdot||_2$ , this is the Rudin-Osher-Fatemi image denoising model(cite). Esedoglu and Osher (cite) showed that TV minimization with anisotropic $γ$ is more suited for denoising arbitrary convex shapes. For example, a unit square is denoised better using the L1-norm $γ(x 1, x 2) = | x 1 | + | x 2 |$ rather than standard (isotropic) L2-norm (See example below).

L1 anisotropic ROF denoised image.

Isotropic ROF denoised image

Salt & Pepper noised image.

In general, for TV minimization of a convex shape $W$ containing the origin, $γ$ should be chosen as the Legendre transform, or the dual, of the indicator function of $W$ .

The Split Bregman Algorithm

The Split Bregman Algorithm of Goldstein and Osher (cite) provides a fast way of solving the TV minimization problem. The core of the algorithm is the alternating minimizations:

$u^{k+1} = \arg\min_u \frac{\lambda}{2}||u-f||_2^2 + \frac{\mu}{2}||d_1^k-u_x-b_1^k||_2^2 + \frac{\mu}{2}||d_2^k-u_y-b_2^k||_2^2$

$(d_1^{k+1},d_2^{k+1}) = \arg\min_{(d_1,d_2)} \gamma(d_1,d_2) + \frac{\mu}{2}||d_1-u_x^k-b_1^k||_2^2 + \frac{\mu}{2}||d_2-u_y^k-b_2^k||_2^2$

where $μ$ is a parameter and the $b_i^k$ 's are the Bregman updates,

$(b_1^{k+1},b_2^{k+1})= (b_1^k,b_2^k) + (u_x^k - d_1^{k+1}, u_y^k - d_2^{k+1})$ .

The first of the two minimizations can be solved efficiently by a few iterations of the Gauss-Seidel method. The second minimization is the anisotropic shrinkage problem for $γ$ . The speed of the algorithm depends heavily on how fast this problem can be solved.

Later, we will consider this problem where $γ$ is a polyhedral norm.

Mean Curvature Flow via TV minimization

The following connection between TV minimization and mean curvature flow was observed by Almgren, Wang and Taylor (cite): given a curve $C k = {x : u k (x) = 0}$ , compute

$u^{k+1} = \arg\min_u \int|\nabla u| + \frac{1}{2h}\int |u-SDF(u^k)|$

where $h > 0$ and $S D F (u k)$ is the signed distance function to the zero level set of $u k$ . Then $C k + 1 = {x : u k + 1 (x) = 0}$ is an approximation to the mean curvature motion of $C k$ after time $h$ . In general, replacing the first integrand by $\gamma(\nabla u)$ gives the anisotropic mean curvature motion. For example, if $γ(x 1, x 2) = | x 1 | + | x 2 |$ we get mean curvature motion by square anisotropy.

To conclude, anisotropic mean curvature flow can be generated by solving a sequence of anisotropic TV minimizations (and solving signed distance functions). Furthermore, anisotropic TV minimization reduces to solving anisotropic shrinkage, via the Split Bregman Algorithm. Hence, we focus on the shrinkage problem where $γ$ is anisotropic. In particular, we consider the case where $γ$ is a polyhedral norm.

The basic shrinkage problem

Let $x = (x 1, x 2)$ be a vector in the plane and let $N (x)$ be a norm. The basic problem is to find the mapping $x = x (y)$

$x = \arg \min \left ( N(x) + \frac \mu 2 \| x - y \|^2 \right )$

This involves the minimization of a homogeneous order one function (the norm) and the standard quadratic.

In the case of the L1 norm, we get standard shrinkage. The Euclidean norm is also easily solved. For the L infinity norm, we can obtain the answer by rotation, but it's useful as practice for the general case to work it out, using subdifferentials, i.e. non-smooth convex optimization.

Dual norm formulae

Generalized Shrinkage for polyhedral norms

Let $N(x) = max_i (n_i \cdot x)$ be a polyhedral norm. Consider the shrinkage problem

$x^* = \arg\min_x \{N(x) + \frac{\lambda}{2}||x-y||^2\}$ .

Denote by $\partial N(x)$ the subdifferential of $N (x)$ . Then, by convex analysis, the minimizer $x *$ should satisfy the inclusion $0\in\partial N(x^*) + \lambda (x^*-y)$ , or equivalently $y \in x^* + \frac{1}{\lambda}\partial N(x^*)$ . Note that the inclusion defines a unique $x *$ for a given $y$ , but not vice versa.

If $N$ is smooth at $x *$ , we have $\partial N(x^*) = n_i$ for some $i$ , so we simply have $x^* = y - \frac{1}{\lambda}n_i$ . If $N$ is not smooth at $x *$ , then necessarily $n_i\cdot x^* = n_{j}\cdot x^*$ for some $i, j$ . In this case, $\partial N(x^*) = conv\{n_i, n_j\}$ . Therefore $y\in x^* + \frac{1}{\lambda}conv\{n_i, n_j\}$ . The last inclusion implies that, in a region near the ray $R_{i,j} = \{x : n_i\cdot x = n_j\cdot x \ge 0\}$ , multiple values of $y$ can shrink to the same minimizer $x^*\in R_{i,j}$ .

(Note, the reason for the calculation of the subgradient of the norm function is given by Danskin's Theorem

The above analysis gives rise to a formula to compute $x *$ given $y$ and set of normals ${n 1, n 2,..., n k}$ ordered in a clockwise direction. First consider the following formula: if $y - (n_i\phi_i + n_{i+1}(1-\phi_i))/\lambda \in R_{i,i+1}$ , then,

$\phi_i = (\lambda(n_i - n_{i+1})\cdot y - n_i\cdot n_{i+1} + ||n_{i+1}||^2 )/||n_i - n_{i+1}||^2.$

Thus, $x^*=y - (n_i\phi_i + n_{i+1}(1-\phi_i))/\lambda \in R_{i,i+1}$ if and only if $\phi_i\in [0,1 ]$ . Furthermore, if $φ i < 0$ and $φ i + 1 > 1$ , then $N$ is smooth at $x *$ with subgradient ${n i + 1}$ , so $x^* = y - \frac{1}{\lambda}n_{i+1}$ . In conclusion, the values of $φ i$ characterizes the direction which $y$ should be shrunk, and immediately gives the value of $x *$ .

The Polyhedral Shrinkage Algorithm

Input: $y\in \R^n$ and an ordered set of normals $\{n_i\}_{i=1}^k$ (assume $n k + 1 = n 1$ ).

Output: minimizer $x^*(y)\in \R^n$ of the shrinkage problem.

1. For each $i = 1,2,..., k$ , compute $φ i$ defined as above.

2. If for some $j$ , $\phi_j \in [0,1]$ , then $x * = y - (n j φ j + n j + 1 (1 - φ j)) / λ$ .

3. If for some $j$ , $φ j < 0$ and $φ j + 1 > 1$ , then $x^* = y - \frac{1}{\lambda}n_{j+1}$ .

The algorithm requires a few operations and comparisons for each $i = 1,2,..., k$ .

Numerical Results for Polyhedral Shrinkage

Surface plots of $x^*(y) = (x^*_1(y), x^*_2(y))$ and the unit balls of the corresponding $N (x)$ are shown.

Max Norm

Skewed Hexagon Norm Shrinkage

Octagon Norm Shrinkage

Max Norm Shrinkage

Skewed Hexagon Norm

Octagon Norm

Numerical Results for Mean Curvature Flow using Generalized Shrinkage

Numerical results for anisotropic mean curvature flow using the Split Bregman Algorithm and anisotropic shrinkage are shown below. The anisotropies tested were those used in the previous section, as well as the L1-norm (the square anisotropy) and a triangular anisotropy. Parameters used were $h = 1,μ = 2$ . The distance function was computed using a combination of the fast marching method and a few iterations of the discretization of the PDE $u_t = sgn(u)(1-|\nabla u|)$ .