A new numerical scheme for discrete constrained total variation flows and its convergence

In this paper, we propose a new numerical scheme for a spatially discrete model of total variation flows whose values are constrained to a Riemannian manifold. The difficulty of this problem is that the underlying function space is not convex; hence it is hard to calculate a minimizer of the functional with the manifold constraint even if it exists. We overcome this difficulty by “localization technique” using the exponential map and prove a finite-time error estimate. Finally, we show a few numerical results for the target manifolds S2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S^2$$\end{document} and SO(3).


Introduction
We are concerned with a numerical scheme for solving a spatially discrete constrained total variation flow proposed by Giga and Kobayashi [17], which is designated as (DTVF GK ; u 0 ) and called the discrete Giga-Kobayashi (GK) model in the present paper (see Definition 1).
A general constrained total variation flow (constrained TV flow for short) for u : Ω× [0, T ) → M is given as where Ω ⊂ R k (k ≥ 1) is a bounded domain with Lipschitz boundary ∂Ω, M a manifold embedded into R ( ≥ 1), u 0 : Ω → M an initial datum, π p the orthogonal projection from the tangent space T p R (= R ) to the tangent space T p M(⊂ R ) at p ∈ M, ν Ω the outer normal vector of ∂Ω and T > 0. If π u is absent, (TVF; u 0 ) is the standard vectorial total variation flow regarded as the L 2 -gradient flow of the isotropic total variation of vector-valued maps: The introduction of π u means that we impose a restriction on the gradient of total variation so that u always takes value in M. The constrained TV flow is also called the "M-valued TV flow" or "1-harmonic map flow". In some literature, Eq. (1.1) is called the constrained TV flow equation or system while the constrained TV flow itself means its solution. However, in this paper we do not distinguish "flow" and "flow equations". The discrete GK model is the constrained gradient flow of the total variation in the space of piecewise constant mappings defined in a given partition of Ω. Originally, [17] proposed the one space dimensional model. Later, the high-dimensional case is studied in [36]. It is formally a system of ordinary differential equations, but the velocity is singular with respect to its variables. Moreover, because of the presence of manifold constraints, it is impossible to interpret the problem as a gradient flow of a convex function. The goal of this paper is to give a time-discrete scheme which is not only practically easy to calculate but also converges to the solution of the discrete GK model.
The readers might be interested whether this (spatially) discrete GK model converges to the constrained TV flow (1.1) if the space grid tends to zero. This problem is widely open and is not simple even for unconstrained case as discussed at the end of Sect. 3.

Applications in science and engineering
Constrained TV flows (TVF; u 0 ) have applications in several fields. The first application of the flow of this type appears in [37], where the authors consider the two-dimensional sphere as the target manifold M in order to denoise color images while preserving brightness. The system with the target manifold M of all threedimensional rotations SO (3), is an important prototype of the continuum model for time-evolution of grain boundaries in a crystal, proposed in [25,26]. Moreover, Eq. (1.1), where the target manifold M is the space of all symmetric positive definite three-dimensional matrices S P D (3), is proposed for denoising MR diffusion tensor image [4,9,32,39].

Mathematical analysis
Despite its applicability, the mathematical analysis of the manifold-constrained TV flows is still developing. Two difficulties lie in mathematical analysis: One is the singularity of the system when ∇u vanishes; the other is the constraint with values of flows in a manifold. Many studies can be found to overcome the first difficulty. In order to explain the second difficulty, we distinguish two types of solutions: "regular solution" and "irregular solution".
We mean by "regular solution" a solution without jumps. In [20], the existence of local-in-time regular solution was proved when Ω is the k-torus T k , the manifold M is an ( −1)-sphere S −1 and the initial datum u 0 is sufficiently smooth and of small total variation. Recently, this work has been improved significantly in [15]. In particular, the assumption has been weakened to convex domain Ω and Lipschitz continuous initial data u 0 . Moreover, in [15], the existence of global-in-time regular solution and its uniqueness have been proved when the target manifold M has non-positive curvature, and the initial datum u 0 is small.
In [18], it has been proved that rotationally symmetric solutions may break down, that is, lose their smoothness in finite time when Ω is the two-dimensional unit disk and M = S 2 . Subsequently, in [5], the optimal blowup criterion for the initial datum given in [18] was found, and it was proved that the so-called reverse bubbling blowup might happen.
For "irregular solution", a solution that may have jumps, two notions of a solution are proposed, depending on the choice of the distance to measure jumps of the function. These choices lead to different notions of a solution of a constrained gradient system, which may not coincide.
Weak solutions derived from "extrinsic distance" or "ambient distance", the distance of the Euclidean space in which the manifold M is embedded, were studied in [17,19]. In [17], the existence and the uniqueness of global-in-time weak solution in the space of piecewise constant functions, have been established when M is compact, the domain Ω is an interval with finite length, and the initial datum u 0 is piecewise constant. Moreover, a finite-time stopping phenomenon of S 1 -valued TV flows was also proved. On the contrary, for S 2 -valued TV flows, an example that does not stop in finite time was constructed in [19], which will be reproduced numerically by our new numerical scheme and used for its numerical verification in this paper.
Weak solutions derived from "intrinsic distance", the geodesic distance of the target manifold M, were studied in [6,13,14]. In [13], the existence and the uniqueness of global-in-time weak solution have been proved when Ω is a bounded domain with Lipschitz boundary, M = S 1 , and the initial datum u 0 has finite total variation and does not have jumps greater than π . These arguments and results were extended in [6] when the target manifold M is a planar curve. As for higher dimensional target manifolds, the existence of global-in-time weak solution was proved in [14] when the target manifold M is a hyperoctant S −1 + of the ( − 1)-sphere. If one takes the intrinsic distance, the uniqueness of solution may fail to hold as pointed out in [36]. This is one reason why we adopt the extrinsic distance.
In [38], discrete models of S 1 -valued TV flows and S 2 -valued TV flows based on the finite difference method were studied. More precisely, a numerical scheme was proposed, and numerical computations were performed.
In [3,11], S −1 -valued regularized TV flows based on the finite element method were studied. In these works, the existence and uniqueness of global-in-time solution of the discrete models were established, and numerical computations were performed. We remark that the convergence of the discrete model to the original model was also studied. However, its argument has some flaws, which were pointed out in [13].
In [17,20,21,36], discrete models of constrained TV flows on the space of piecewise constant functions were studied. Such discrete models can be directly applied to denoising of manifold-valued digital images. For one-dimensional spatial domains, solutions of the discrete models coincide with irregular solutions of the corresponding original model derived by ambient distance but it is not the case for higher dimensional domains. In one-dimensional spatial domain, the existence and the uniqueness of global-in-time solution to the discrete models were established in [17]. Moreover, numerical computations of S 1 -valued discrete models were performed. These discrete models are formulated as ordinary differential inclusions; i.e., differential equation with multi-valued velocity. There are two key ideas in [17] to solve them. The first one is computation of the canonical restriction of the multi-valued velocity. The second one is to use facet-preserving property of flows. We emphasize that these two key ideas do not work in dimensions higher than one. In the higher dimensional case, the existence and the uniqueness of global-in-time solution of the discrete model were established in [20,21,36]. Numerical computations of these discrete models, however, were not performed.

Contribution of this paper
This paper is dedicated to the study of a numerical scheme for simulation of the discrete model studied in [20,21,36], of constrained TV flow. In particular, we propose a new numerical scheme and show its convergence. We also perform numerical simulations based on the proposed scheme. We outline the three main contributions below:

New numerical scheme
Constrained discrete TV flow is formulated as a gradient flow in a suitable manifold. Hence, one can use the minimizing movement scheme (see [2]) to simulate it. It is summarized as follows: Let H := (H , ·, · H ) be a real Hilbert space, E be a submanifold of H , F be a real-valued functional on H allowing the value +∞, I := [0, T ] be a time interval and τ > 0 be a step size. We first consider a sequence {u (n) τ } in E generated by a simple minimizing movement scheme for F , which will be referred as (MM; τ , u 0 ) in Sect. 3. In this scheme, {u (n) τ } is determined successively by taking a minimizer in E of a functional Note that the uniqueness of minimizers is not guaranteed since E is not a convex constraint. In general, if F were geodesically convex and coercive on E, then the piecewise linear interpolation (Rothe interpolation) of {u (n) τ } would converge to the gradient flow of F (see [2]) and this general theory could be applied. However, this is not the case in our situation. In this simple scheme, we need to solve an optimization problem at each step. Each optimization problem is classified as Riemannian optimization problem, i.e., an optimization problem with Riemanninan manifold constraint. Theory of smooth Riemannian optimization, that is, Riemannian optimization problem whose objective function is smooth is well-studied, and we refer to the systematically summarized book [1]. However, our energy F is the total variation energy so it is not smooth. Unfortunately, non-smooth Riemannian optimization is only scarcely explored as a subfield of the theory of Riemannian optimization. We refer to [27,28,40,41] as its references. Moreover, the theory of non-smooth Riemannian optimization is generalized to the theory of non-smooth and non-convex optimization with a separable structure in [30]. Although there are several studies in this direction, it is under development as compared with the linearly constrained problem. For these reasons in this paper, we propose a new minimizing movement scheme which includes only a linearly constrained optimization problem, which is referred as (MM loc ; τ, u 0 ) in Sect. 3. In this scheme, instead of minimizing defined for X which is an element of the tangent space T u (n−1) τ E contained in H . This is a linearly constrained (convex) optimization problem, which will be referred as (VP loc ; u (n−1) τ ). There is a unique minimizer. Let X (n−1) τ be the minimizer. We then determine u

Stability and convergence
In this paper, we study stability and convergence of the scheme from two points of view. Those are "discrete energy dissipation" and "error estimate" for the proposed scheme. Since discrete TV flows have the structure of gradient flows, the proposed scheme should inherit the properties of the gradient flow. Therefore, in this paper, we show that if τ is sufficiently small, the proposed scheme satisfies the energy dissipation inequality, which is one of the properties of the gradient flow.
Moreover, we prove that the sequence {u (n) τ } generated by the modified minimizing movement scheme converges to the original gradient flow as τ → 0. The convergence result of the minimizing movement scheme appears in [2,10,34]. In the classical work [10] on the convergence analysis in Banach space, error estimate of order O( √ τ ) was obtained, and it is improved to the order O(τ ) in [34]. A similar estimate to [10] in general metric space is shown in [2]. Since our scheme contains the "localization" process, we cannot apply these previous works to our scheme directly. We can show, however, the error estimate locally in time between the piecewise linear interpolation, the so-called Rothe interpolation, of {u (n) τ }, and the original flow. This estimate, which will be presented in Sect. 3.2.3, states that the error of the Rothe interpolation is O( √ τ ) as τ → 0, and corresponds to those in [2,10]. We remark that the idea using the exponential map appears in the optimization problem in matrix manifolds and in numerical computation of regularizing flows like constrained heat flows (see [1,7,8]). We are concerned, however, with the approximation of the gradient flow in addition to converging to a minimizer. As far as we know, there is no previous rigorous result on the convergence of the flow itself.

Numerical simulations
The proposed scheme is not enough to simulate constrained discrete TV flow since we need to solve the linearly constrained convex but non-smooth optimization problem (VP loc ; u (n−1) τ ) at each step. We overcome this situation by rewriting (VP loc ; u (n−1) τ ) as an iteration, and adopt alternating split Bregman iteration, proposed by [23], which is adequate to the optimization problem including total variation. We refer to [31,33] for examples of the application of this iteration to calculate the mean curvature flow numerically. One can find a proof of convergence of this iteration in [35].
In this paper, we provide a numerical analysis of three aspects of constrained TV flows. The first one is a property of S 2 -valued TV flow, discovered in [19], of not stopping in finite time. The second one is an error estimate of the proposed scheme; actually, the example in [19] can be rewritten as a simple ordinary differential equation, which can be accurately solved by an explicit scheme. Therefore, we can confirm the theoretical convergence rate by comparing the results obtained by the two numerical methods. The third one is the numerical observation that a facet is preserved in most of evolution. We simulate S 2 -valued TV flows with one spatial dimension constructed in [19], and SO(3)-valued TV flows in two spatial dimensions.

Advantages
The proposed scheme has five advantages: First, this scheme does not restrict the target manifold M. In many previous studies, the target manifold M is fixed in advance, to the sphere, for example. Our method, however, can be applied to any Riemannian manifold as target manifold M.
Second, this scheme is not restricted to constrained TV flows and can be applied to more general constrained gradient flows. Indeed, as we can see from the scheme (MM loc ; τ, u 0 ), it is possible to construct a numerical solution {u Third, the proposed scheme can describe facet-preserving phenomena of constrained TV flows. In the numerical calculation of the TV flow, we should pay attention to whether a numerical scheme can adequately simulate the evolution of facets. Many schemes that have been proposed so far cannot capture this phenomenon since the energies are smoothly regularized. On the other hand, our scheme is capable of preserving facets since the energies are only convexified, not regularized.
Fourth, the proposed scheme is numerically practical. Especially, if the exponential map and the orthogonal projection π of the target manifold M can be calculated easily, the practical advantage of our scheme is clear. If M is in the class of orthogonal Stiefel manifolds, its orthogonal projection π and its exponential map can be written explicitly. See also [1]. Besides, as mentioned above, our method does not use the projection onto the target manifold, which is sometimes hard to calculate.
Finally, the proposed scheme is well-defined, and we shall prove its convergence together with convergence rate.

Organization of this paper
The plan of this paper is as follows: In Sect. 2, we recall notion and notations for describing constrained discrete TV flows we study in this paper. More precisely, we recall the notion of manifolds and define a space of piecewise constant functions, discrete total variation, and a discrete model of constrained TV flows.
In Sect. 3, we propose a new numerical scheme and provide its analysis. Namely, in Sect. 3.1, we derive a new numerical scheme for constrained discrete TV flows, starting from the minimizing movement scheme for constrained discrete TV flows. In Sect. 3.2.1, we explain Rothe interpolation. This interpolation is useful for establishing the energy dissipation inequality and error estimate of the proposed scheme. In Sect. 3.2.2, we prove that the proposed scheme satisfies energy dissipation inequality if the step size is sufficiently small. In Sect. 3.2.3, we establish an error estimate of the proposed scheme. We see that the error estimate implies that the proposed scheme converges to constrained discrete TV flows as step size tends to zero. The key of the proof is to establish the evolutionary variational inequality, used in [2], for Rothe interpolation of the proposed scheme.
In Sect. 4, we propose a practical algorithm to perform numerical computations. More precisely, we rewrite the scheme into a practical version with alternating split Bregman iteration, and we use it to simulate S 2 -valued and SO(3)-valued discrete TV flows in Sect. 5. In "Appendix A", we explain why the constant C M defined by (2.2) in Sect. 2 is bounded by explicit quantities associated to the submanifold M.

Preliminaries
Here and henceforth, we fix the bounded domain Ω in R k with Lipschitz boundary and denote by ν Ω the outer normal vector of ∂Ω.

Notion and notations of submanifolds in Euclidean space
Let M be a C 2 -submanifold in R ( ≥ 2). For p ∈ M, we denote by T p M and T ⊥ p M the tangent and the normal spaces of M at p, respectively, and write π p and π ⊥ p for the the orthogonal projections from R to T p M and T ⊥ p M, respectively. We denote by II p : respectively. If M is compact, then Diam(M) and Curv(M) are finite. Given a point p ∈ M and velocity V ∈ T p M, we consider the ordinary differential equation in R , the so-called geodesic equation in M, for γ : [0, ∞) → R : If M is compact, then the Hopf-Rinow theorem implies that there exists a unique curve γ p,V : [0, ∞) → M satisfying (2.1). The curve γ p,V is called geodesic with initial point p and initial velocity V . Then the exponential map exp p : If M is path-connected and compact, then the constant C M is finite. In fact, as explained in Appendix A, an explicit bound for C M can be obtained.

Finite-dimensional Hilbert spaces
First, we define partitions with rectangles. Let Δ be a finite set of indices and let Ω be a bounded domain in R k . A family Here, by a rectangular domain, we mean that Although our theory works for rather a general partition, we consider only the rectangular partition above for practical use. We We denote by e(Δ) the set of edges associated with Δ defined as where H m denotes the m-dimensional Hausdorff measure. For γ := {α, β} ∈ e(Δ), we take a bijection Sign : γ → {±1} and set E γ := ∂Ω α ∩ ∂Ω β . Moreover, we define the set EΩ Δ of interior edges in Ω associated with Ω Δ and the symbol Sign e(Δ) on e(Δ) by respectively. Subsequently, assume that the partition Ω Δ of Ω is given. Then, set We regard the space H Δ as a closed subspace of the Hilbert space L 2 (Ω; R ) endowed with the inner product defined as For U ∈ H Δ , we denote the facet of U by Further, we set We regard the non-convex space M Δ as a submanifold of H Δ in which the function takes values in M. Subsequently, for u ∈ M Δ , we denote by T u M Δ the tangent space of M Δ at u, i.e., Moreover, we denote by H EΩ Δ the space of piecewise constant R -valued maps on

Discrete constrained TV flows
The problem (TVF; u 0 ) is formally regarded as the gradient system of (isotropic) total variation: The spatially discrete problems we consider in this paper are regarded as the gradient system of discrete (isotropic) total variation. Let us begin with the definition of discrete total variation associated with Ω Δ . Let Ω Δ := {Ω α } α∈Δ be a rectangle partition of Ω. Then the discrete total variation functional T V Δ : H Δ → R associated with Ω Δ is defined as follows: This definition is easily deduced from the original definition of T V (u) when u is a piecewise constant function associated with Ω Δ . We remark that the functional T V Δ is convex on H Δ but not differentiable at a point u whose facet Facet(u) is not empty. The next proposition is an immediate conclusion of the definition of T V Δ . Hence, we state it without proof: The following statements hold: Constrained discrete TV flow is the constrained L 2 -gradient flow of spatially discrete total variation. The gradient ∇ M Δ T V Δ of constrained discrete TV flow is, formally, given by (2.5) Here, ∇ T V Δ denotes a formal gradient in H Δ and P u denotes the orthogonal projection from H Δ to T u M Δ naturally induced from π p : The discrete total variation T V Δ , however, is not differentiable. Therefore, we use the subgradient ∂ T V Δ in H Δ instead of the gradient ∇ T V Δ : Accordingly, the spatially discrete problem we consider is as follows: (2.7) Here, we note that the existence and the uniqueness of global-in-time solution of discrete GK model have been proved in [17,20,36]: [17,20,36] Let u 0 ∈ M Δ , I := [0, T ) and M be a C 2 -compact submanifold in R . Then, there exists a solution u ∈ W 1,2 (0, T ; M) to the discrete GK model of (TVF; u 0 ). Moreover, assuming that M is path-connected, u is the unique solution.

Derivation
We seek a suitable time discretization of (DTVF GK ; u 0 ), so that we fix the step size τ > 0 and denote by N (I , τ ) the maximal number of iterations, that is, the minimal integer greater than T /τ . Moreover, we define the time nodal points We focus on the gradient structure of (DTVF GK ; u 0 ) to use the minimizing movement scheme in [2]. Then we obtain the following numerical scheme of (DTVF GK ; u 0 ): We determine u Here exp x is the exponential map of the Riemannian manifold M. Since Exp u (n−1) Here, we emphasize that Φ τ loc (·; u Subsequently, we consider the optimization problem (VP loc ; u (n−1) τ ): The problem (VP loc ; u Remark 1 (Why the proposed scheme describes facet-preserving phenomena) In the numerical calculation of (constrained) TV flows, we should pay attention to whether the scheme can adequately simulate the evolution of facets. To see why the scheme has the facet-preserving property, note that given u ∈ M Δ and X ∈ T u M Δ , the total variation of u + X is decomposed into Hence, the minimizer X * ∈ T u M Δ of the optimization problem tends to have the same facet as u, that is, Facet(u) = Facet(X * ). Therefore, Exp u (X * ) also tends to have the same facet as u.
This scheme is always well-defined since the minimizer is unique. Here is the statement:  in the above modified minimizing movement scheme are well-defined, and satisfy for all n = 0, . . . , N (I , τ ) − 1.

Rothe interpolation
We consider the Rothe interpolation of sequences generated by (MM loc ; τ, u 0 ). This interpolation is useful to prove energy dissipation in Proposition 5 and error estimate in Theorem 1. Set the time interpolation functions { (n) τ } N (I ,τ )−1 n=0 , τ : I → [0, 1] as follows: n=0 be the sequence generated by the modified minimizing movement scheme (MM loc ; τ, u 0 ). Then, we define two interpolations u τ , u τ : I → M Δ as follows: for t ∈ I . In particular, u τ is called the Rothe interpolation of {u

Remark 2
By definition, the Rothe interpolation u τ is also represented as where X τ := . Especially, u τ is continuous in I . has the following properties:

the acceleration of u τ is bounded by the speed of u τ , that is,
Moreover, Since the speed of the geodesic is constant, we have Here, in the last inequality, we used Proposition 3. 3. Since u (n−1) τ and u (n) τ are joined by the exponential map, u τ | [t (n−1) ,t (n) ) satisfies the geodesic equation (2.1), that is, Hence, we have Since the geodesic has constant speed, Moreover, Proposition 3 implies that

Discrete energy dissipation property
The total variation dissipates along corresponding constrained TV flows. Hence, it is desirable that the proposed scheme also has this property. Indeed, the proposed scheme has this property if the step size is small enough.
The above formula and the triangle inequality imply Since X (n) τ is the minimizer of (VP loc ; u (n) τ ), we have By applying Lipschitz continuity of T V Δ and the Minkowski inequality for integrals

Error estimate
Here, we establish an error estimate between the sequence generated by (MM loc ; τ, u 0 ) and the solution to (DTVF GK ; u 0 ) when u 0 ∈ M Δ is given.
for all t ∈ I , where (3.7) Here, we again note that the constant C M which appears in (2.2) can be bounded by quantities only associated with the manifold M; that is, all the constants C 0 , C 1 and C 2 are independent of time t and step size τ .

Remark 3
In this theorem, we cannot remove the exponentially growing term e C 0 t from the right hand side of (3.5), because the functional T V Δ defined by (2.4) is not in general convex but merely semi-convex.
We immediately see that if u 1 0 = u 2 0 in (3.5), then we have and thus we conclude that u τ converges to u in C(I ; M Δ ) as τ → 0.
The key estimates to prove Theorem 1 are the evolutionary variational inequalities for u τ and u.
for all t ∈ I \I (τ ), where C 0 , C 1 , C 2 are the same constants as in Theorem 1, and I (τ ) is defined in (3.1).
Proof We will compute 1 2 d dt u τ − v 2 H Δ by splitting it into a semi-monotone term and an error term.
Expanding du τ / dt in Taylor series implies that and inserting this into (3.10) Moreover, we split the term I into a monotone term and a non-monotone term. Proposition 3 implies that Self-adjointness of P u τ and the relation P u τ = I − P ⊥ u τ imply that (3.11) We plug (3.11) into the identity (3.10) to obtain We shall estimate I 1 , I 2 and I I, respectively. First, we estimate the term I 1 : Here, we apply the triangle inequality, x, y ∈ H Δ , to the underlined part in the above inequality to obtain The triangle inequality for T V Δ implies that

Proposition 4 implies
(3.13) Next, we estimate the term I 2 . The Cauchy-Schwarz inequality and Proposition 1 imply that (3.14) Here, we claim that Taking the norm · H Δ in the above equation and applying the triangle inequality and the Minkowski inequality for integrals yield As the second term on the right-hand side is estimated in the same way as in (3.13), we focus on the term Moreover, for sequences it is clear that Thus, for f ∈ H Δ , and we obtain that where The Cauchy-Schwarz inequality implies that

Proposition 4 implies
Plugging (3.19) into (3.16) yields the claim (3.15), and taking into account (3.14), we obtain (3.20) Next, we estimate the term I I. The Cauchy-Schwarz inequality and the Minkowski inequality for integrals imply
The solution u of (DTVF G K ; u 0 ) satisfies the following evolutionary variational inequality which is obtained by an argument similar to that in the proof of Proposition 6; for the proof, see [36].
for a.e. t ∈ (0, T ), where C 0 is the same constant as in Proposition 6.
Now, we will finish the proof of Theorem 1. t ∈ (0, T ). By substituting v = u τ (t) into (3.22) and v = u(t) into (3.9) and adding these two inequalities, we obtain
Remark 4 (Convergence of spatially discrete TV flow as the mesh tends to zero) This problem is difficult and is studied only for unconstrained problems. We consider the case for the TV flow of scalar functions. If one uses a rectangular partition of Ω, then the solution of discrete model solves an anisotropic 1 -total variation flow, which is the gradient flow of . This is known for the one-dimensional case for a long time ago [16] and for higher dimensional case by [29]. Thus, a discrete solution u h converges to a solution u of converges to the continuum initial data u 0 in L 2 (Ω). More precisely, holds since the solution semigroup is a contraction semigroup. For constrained cases, it is expected, but so far, there is no literature stating this fact.

Practical algorithms
In order to implement the proposed scheme, we need to solve a minimization problem in each iteration. We simplify this optimization problem by applying alternating split Bregman iterations. We replace the minimization problem (VP loc ; u (n−1) τ ) with alternating split Bregman iteration which is proposed in [23] to solve the L 1 regularization problem efficiently.
First we apply a splitting method to (VP loc ; u Here, we note that (i) both X (k) and Z (k) in the above iterations are solved explicitly when the orthogonal projection and the exponential map in M have explicit formulae. This is the case, for example, for the class of orthogonal Stiefel manifolds, which includes the important manifolds S 2 and SO(3). (ii) X (k) converges to the minimizer of (VP loc ; u (n−1) τ ) in H 0 thanks to Corollary 2.4.10 in [35].
We explain details concerning the first point (i). In step (a), we seek the global minimizer X (k) of the function Differentiating this function with respect to X , we see that the global minimizer can be characterized as the solution to the corresponding Poisson equation, which is a strictly diagonally dominant system and can be efficiently solved by using well-known solvers such as the Gauss-Seidel method. In step (b), we first seek Z (k) 0 which is the global minimizer in H EΩ Δ of the function Since there are no interactions between {Z γ 0 } γ ∈e(Δ) in the above equation, the global minimizer can be obtained by computing it componentwise. Thus, we need to obtain the global minimizer x * ∈ R of the function τ x R + (ρ/2) x − y 2 R l (y ∈ R ), and we can explicitly write down the global minimizer x * ∈ R as where shrink(x, γ ) (x ∈ R , γ ∈ R) is the shrinkage operator defined by As a result, each (Z (k) 0 ) γ is given by We finally compute Z (k) 1 , which is the global minimizer of the function The global minimizer for this equation can be given by using the orthogonal projection as follows: Summarizing the above, in the alternating split Bregman iteration, we do not need to solve minimization problems directly, and each step requires only solving a well- Remark 5 As we have explained in the above, the computational cost of our method based on alternating split Bregman iteration is cheap. Indeed, as pointed out in [23], if we choose parameters properly, the number of iterations in Algorithm (VP loc,split ; u (n−1) τ ) becomes small; however, as far as we know, there are no mathematical guidelines on optimal choice of parameters (see [22,23] for detailed explanation on the choice of the initial values Z (0) and B (0) and the parameter ρ).
For constrained TV flows, several numerical methods have been proposed based on the regularization of the total variation energy or Lagrange multipliers. If we impose a regularization of the total variation energy, the problem on singularity disappears but it becomes impossible to capture the steep edge structure and the facet-preserving phenomenon. When the target manifold M is the sphere S l−1 , then it is not so difficult to obtain the unconstrained problem by introducing the Ginzburg-Landau functionals or Lagrange multipliers. However, if we focus on more complicated manifolds such as SO (3), introducing these kinds of regularizations is not straightforward; also, it is not clear whether or not we can construct a practical numerical scheme to simulate it. Our approach does not regularize but merely convexify the total variation energy, and we adopted the orthogonal projection and the exponential map to obtain the solution at next time step. Namely, our method does not restrict the target manifold a priori and can be applied to a broad class of target manifolds.

Numerical results
In this section, we use the above scheme to simulate S 2 and SO(3) valued TV flows, respectively. Throughout the numerical experiments, we choose initial values Z (0) and B (0) in algorithm (VP loc,split ; u (n−1) τ ) as zeros, and set the parameter ρ to be equal to 0.1 (see [23] for detailed explanation on the choices of the initial values Z (0) and B (0) and the parameter ρ). Moreover, we terminate the iteration for computing X (n−1) τ when the relative error becomes less than 10 −4 : H Δ , where we have defined X (0) as zero.

The tangent spaces, orthogonal projections and exponential maps of S 2
We regard the 2-sphere as Then the tangent spaces, their orthogonal projections and exponential maps in S 2 are given by the following explicit formulae: where I 3 denotes the identity matrix in R 3 and exp denotes the matrix exponential; here x = (x 1 , x 2 , x 3 ) is a column vector and x denotes its transpose.

Euler angles
Vectors in S 2 have three parameters. Euler angle representation is beneficial to reduce parameters of S 2 . Given γ := (x, y, z) ∈ S 2 , its Euler angle representation is given as follows:

Counterexample to finite-time stopping phenomena
In [19], an example of constrained TV flow which does not reach the stationary point in finite time is shown. Here is the statement.
In this theorem, h(t) = (h 1 (t), 0, h 3 (t)) satisfies the following system of differential equations: which can be calculated numerically. Here c = 2 − 1 . Therefore we use it as a benchmark for the validation of our algorithm.
Remark 6 (Dirichlet problem) So far, in this paper, we have considered the Neumann problem of constrained TV flows, while this example solves the Dirichlet problem. Therefore, we can not apply the proposed scheme directly. However, we can derive the Dirichlet problem version of the proposed scheme just by replacing For more on the Dirichlet problem, see [20].  for X ∈ T x (SO (3)) so x (3), where L x denotes the left translation of SO(3) given by g → xg. In other words, the exponential map is given in the following simple form exp x (X ) = x exp(x X ), X ∈ so x (3).
Since the decomposition