# Variational Multiscale Nonparametric Regression: Algorithms and Implementation
## 3 Computational Methods  
In this section, we discuss different algorithmic approaches to compute the MIND \(\hat{f}\) in (3). First, we stress that the optimisation problem (3) is a non-smooth, convex, high-dimensional optimisation problem, which allows in principle to exploit any optimisation method designed for such situations. In the following, we restrict to three different approaches that we found particularly suited for our scenario. To set the notation in this section, we rewrite (3) as  
\[\min_{v\in\mathbb{R}^{n}}J\left(v\right),\qquad J(v):=F(Kv)+G(v), \tag{6}\]  
where \(F\) and \(G\) are lower semi-continuous, proper convex functionals given by  
\[F(w) :=1_{\leq 0}(w-KY-q_{n})+1_{\leq 0}(-w+KY-q_{n})\ \ \text{for}\ w\in \mathbb{R}^{\#\Lambda_{n}},\] \[G(v) :=R(v)\ \ \text{for}\ v\in\mathbb{R}^{n}.\]  
and \(K\) is the linear (bounded) operator that maps from \(\mathbb{R}^{n}\) to \(\mathbb{R}^{\#\Lambda_{n}}\) and is defined by  
\[[Kg_{n}]_{\lambda}:=\left\langle\phi_{\lambda},g_{n}\right\rangle\ \ \text{for any}\ \lambda\in\Lambda_{n}.\]  
Here, and in what follows, \(1_{\leq 0}\) denotes the indicator function of the negative half-space, which is  
\[1_{\leq 0}\left(v\right)= \begin{cases}0&v_{i}\leq 0\ \text{for all possible}\ i,\\ \infty&\text{else}.\end{cases} \]  
Due to the convexity of \(F\) and \(G\), the problem (6) can equivalently be solved by finding a root of the subdifferential \(\partial J(v)\) (see, e.g., [54]), which is a generalisation of the classical derivative. Such roots correspond to stationary points of the evolution equation  
\[\partial v\left(t\right)\in-\partial J\left(v(t)\right),\qquad t>0.\]  
Two fundamental algorithms for the solution of this equation arise from applying the explicit or implicit Euler method, which leads to the steepest decent procedure  
\[v_{k+1}\in\left(I-\lambda_{k}\partial J\right)\left(v_{k}\right),\qquad k\in \mathbb{N}\] (7a) and the so-called proximal point method \[v_{k+1}\in\left(I+\lambda_{k}\partial J\right)^{-1}\left(v_{k}\right),\qquad k \in\mathbb{N} \tag{7b}\]  
respectively, where \(I\) denotes the identity operator and \(\lambda_{k}>0\) is a step size parameter. For (7a), \(\lambda_{k}\) should be chosen such that \(\left(I-\lambda_{k}\partial J\right)\) is a contraction to ensure convergence. For continuously differentiable \(J\), it therefore suffices to choose \(\lambda_{k}<L^{-1}\) with the Lipschitz constant \(L\) of \(\
abla J\); in practice, one often uses \(\lambda_{k}=\left(2L\right)^{-1}\) with an estimator \(L\) of the Lipschitz constant (generated, e.g., by a power method). In the case of (7b), convergence can be obtained for any choice of \(\lambda_{k}>0\), but the computation of the so-called resolvent operator \(\left(I+\lambda_{k}\partial J\right)^{-1}\) is clearly as difficult as the original problem. Motivated by these two methods, different algorithms have been suggested. Therefore, the subdifferential \(\partial J(v)\) of \(J\) is split into the subdifferentials of \(F\) and \(G\), which allow for a much simpler computation of the corresponding resolvent operators. This makes use of the formula  
\[\partial J(v)=K^{*}\partial F\left(Kv\right)+\partial G(v), \tag{8}\]  
which holds true whenever there exists some vector \(v\in\mathbb{R}^{n}\) such that both \(F\) and \(G\) are finite and continuous at \(Kv\) and \(v\) respectively (see e.g., Prop. 5.6 in [55]). In our situation, this is the case whenever \(R\) is continuous at a point in the interior of the feasible set.