# Numerical Implementation of Variational Morphological Operators on Graphs

Miguel Amorim, Jean-Christophe Pesquet, Tristan Portugues, Antonio

Silveti-Falls, Hugues Talbot

► To cite this version:

Miguel Amorim, Jean-Christophe Pesquet, Tristan Portugues, Antonio Silveti-Falls, Hugues Talbot. Numerical Implementation of Variational Morphological Operators on Graphs. DGMM 2025 - 4th International Conference on Discrete Geometry and Mathematical Morphology, Nov 2025, Groningen, Netherlands, Netherlands. hal-05421047

HAL Id: hal-05421047

https://inria.hal.science/hal-05421047v1

Submitted on 17 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical Implementation of Variational Morphological Operators on Graphs

Miguel Amorim, Jean-Christophe Pesquet, Tristan Portugues, Antonio Silveti-Falls, Hugues Talbot

CVN, CentraleSupélec, Inria, Université Paris-Saclay
tristan.portugues@centralesupelec.fr

**Abstract.** In this work, we study several implementations of the variational formulation of classical morphological operators such as erosion, dilation, opening, and closing, as well as arbitrary rank filters using convex optimization over graphs. Since the structuring element can naturally vary along nodes in the graph, this formulation allows for an effective implementation of spatially-variant mathematical morphology. It is particularly efficient on massively parallel architectures. We focus on solving the associated optimization problems using various algorithms from convex analysis. The problem being non-smooth, we consider algorithms based on the proximal operator, such as ADMM and PPXA+. We tested these approaches on both images and graph data. Our results show that proximal algorithms can handle the problems efficiently. We believe that this variational framework is a promising direction to combine mathematical morphology with modern computational tools, and is a direct avenue for solving inverse problems with mathematical morphology operators in the formulation.

## 1 Introduction

Mathematical morphology (MM) is a fundamental framework in image and signal processing, providing nonlinear operators such as erosion, dilation, opening and closing [12,7]. A recent variational framework has reformulated these operators, as well as arbitrary rank filters using convex optimization principles, making them well defined on arbitrary graphs [8]. This formulation opens the door to new computational strategies, parameter tuning, and potential integration with machine learning architectures, e.g. to solve inverse problems in a different way than with classical linear operators. In this work, we investigate the practical resolution of this variational formulation using methods suited to the structure of the problem.

In Section 2, we propose zeroth-order methods for solving the variational erosion/dilation problem. Section 3 addresses the more challenging task of implementing variational opening and closing. Two proximal splitting approaches are shown to yield effective solutions. Finally, conclusions and future works are presented in Section 4.
**Notation and Graph Framework.** Let $G = (V, E)$ be a directed, weighted, reflexive graph with $|V| = n$. Each node $v_i \in V$ indexed by $i \in \{1, \dots, n\}$, is associated with a scalar value $y_i \in \mathbb{R}$. The neighborhood of node $i$ is $N_i = \{j | (i, j) \in E\}$ and its reciprocal neighborhood is $\tilde{N}_i = \{j | (j, i) \in E\}$.

## 2 Variational Formulation of Erosion and Dilation

We consider the variational formulation of the morphological erosion and dilation operators, as introduced in [8], which we will now describe. Given an input graph with node values $y \in \mathbb{R}^n$, the variational erosion/dilation are defined by the following convex optimization problem

$$ \min_{x \in \mathbb{R}^n} \sum_{i=1}^{n} \sum_{j \in U_i} \varphi_i (\sigma_{\Omega_i}(x_i - y_j - w_{ij})) \quad (1) $$

where $U_i$, a nonempty subset of $V$, is the set of neighbors of node $i$ and $(w_{ij})_{j \in U_i} \in \mathbb{R}^{|U_i|}$ is a vector of shift parameters. The support function $\sigma_{\Omega_i}$ is defined for an interval $\Omega_i = [\alpha_i^-, \alpha_i^+]$, $\alpha_i^- < 0 < \alpha_i^+$, as

$$ (\forall x \in \mathbb{R}) \quad \sigma_{\Omega_i}(x) = \begin{cases} \alpha_i^- x & \text{if } x < 0 \\ 0 & \text{if } x = 0 \\ \alpha_i^+ x & \text{if } x > 0. \end{cases} $$

We define on $\mathbb{R}_+$ the convex function $\varphi_i: x \mapsto x^p$ for $p \in [1, 2]$. The parameters of these two functions $\sigma_{\Omega_i}$ and $\varphi_i$ can depend on the node $i$ and therefore be different at each node of the graph. This property of the variational formulation allows for adaptive morphology operations. For $p = 1$ and $p = 2$, closed-form solutions can be derived from [8] for the formulation of erosion and dilation.

The choice of neighborhood $U_i$ for each node $i$ is totally flexible. The solution of the optimization problem remains well defined for arbitrary neighborhoods, although varying choices of neighborhoods can lead to significantly different filtering behaviors, as it is well-known in MM. Despite these differences, the computational process for solving the variational problem remains unchanged, making the approach both general and adaptable. Like the parameters of the support and convex cost function, the neighborhoods can differ from node to node, allowing the implementation of structuring functions.

### 2.1 Closed-Form Results from the Variational Formulation

First results found in [8], give closed-form solutions to the specific cases of the variational problem. The morphological erosion and dilation can be recovered by setting $\alpha_i^- / \alpha_i^+$ take infinity values, while a rank filter and a weighted mean are found for $p = 1$ and $p = 2$. For other values of $p \in ]1, 2[$, we have an interpolation between the non-differentiable solution of $p = 1$ and the continuous and differentiable solution of $p = 2$.
Let $U_i = N_i$ be the neighborhood of node $i$ and $(z_{ij})_{1\le j\le |U_i|}$ the coefficients $(y_j + w_{ij})_{j\in U_i}$ indexed in increasing order.

For $\alpha_i^+ = +\infty$, we have

$$x_i = \min_{j \in U_i} y_j + w_{ij} \quad (\text{Erosion}).$$

For $\alpha_i^- = -\infty$, we have

$$x_i = \max_{j \in U_i} y_j + w_{ij} \quad (\text{Dilation}).$$

For $p_i = 1$, let $j_i = \lfloor \frac{|U_i|}{1-\alpha_i^+ / \alpha_i^-} \rfloor + 1$, the optimal $x_i$ is given by

$$\hat{x}_i = z_{ij_i}$$

if $j_i - 1 \notin \mathbb{N}$, otherwise $\hat{x}_i$ can be chosen arbitrarily in $[z_{ij_{i-1}}, z_{ij_i}]$. This corresponds to a weighted rank filter (as defined in [8, Prop. 3]).

For $p_i = 2$, the optimal $x_i$ is given by

$$x_i = \frac{(\alpha_i^+)^2 \sum_{j=1}^{j_i-1} z_j + (\alpha_i^-)^2 \sum_{j=j_i}^{|U_i|} z_j}{(j_i - 1)(\alpha_i^+)^2 + (|U_i| - j_i + 1)(\alpha_i^-)^2}$$

where $j_i$ is the rank index at which the minimizer lies (as defined in [8, Prop. 2]).

## 2.2 General Case for $p > 1$: Root-Finding Approach

For most values $p \in ]1, 2[$, except for a few special cases, the variational formulation of erosion and dilation does not admit a closed-form solution. However, due to the separable structure of the objective, the problem reduces to a set of independent scalar convex minimization problems, one for each pixel or graph node $i \in \{1, \dots, n\}$.

For each pixel value $x_i$, we must solve the following scalar convex minimization problem

$$\min_{x_i \in \mathbb{R}} F_i(x) = \min_{x_i \in \mathbb{R}} \sum_{j \in U_i} \varphi_i (\sigma_{\Omega_i}(x_i - y_j - w_{ij})) \quad (2)$$

For finite values of $\alpha_i^-$ and $\alpha_i^+$, the derivative of the objective function $F_i$ at $x \in \mathbb{R}$ can be computed directly by

$$F'_i(x) = \sum_{j \in U_i} \begin{cases} p\alpha_i^- (\alpha_i^- (x - y_j - w_{ij}))^{p-1} & \text{if } x < y_j + w_{ij}, \\ 0 & \text{if } x = y_j + w_{ij}, \\ p\alpha_i^+ (\alpha_i^+ (x - y_j - w_{ij}))^{p-1} & \text{if } x > y_j + w_{ij}. \end{cases} \quad (3)$$
This derivative $F'_i$ is continuous and strictly increasing, $F_i$ is strictly convex. Hence the minimizer $x_i^*$ is unique and must satisfy

$$F'_i(x_i^*) = 0.$$

We solve this root-finding problem directly for each node $i$. Due to nonsmoothness and lack of global Lipchitz continuity, classical first or second-order solvers are not efficient root finders. Instead, zeroth-order methods can be applied to find the unique zero of $F'$. Due to the convexity and strict monotonicity of the derivative, our root-finding approach is guaranteed to solve the variational erosion and dilation problem for all $p \in [1, 2]$. Since the problem is separable and each output node can be solved in parallel, the proposed solution is efficient and numerically stable, and well-suited for implementation on large graphs or image grids.

We tested different zeroth-order methods to solve the problem such as the Bisection, Golden search, Ternary search and Interpolate-Truncate-Project methods [2]. While all these methods share the same worst-case complexity, the Golden search and bisection methods show better results over larger images (see Figure 1).

<figure><img src="image_2.png" /><figcaption>Fig. 1. Performance against image size for different zeroth-order methods</figcaption></figure>

The results applied to pictures and 3-D meshes show a clear interpolation between the two traditional morphological operators (Erosion and Dilation) when $-\alpha_+/\alpha_-$ goes from zero to large values (see Figure 2 and 3). When using different neighborhood structures (also known as structuring elements), the results show a similar interpolation behavior while producing images with properties that follow the neighborhood shape. We show the resulting images using a diagonal cross neighborhood instead of a square of radius 1 (see Figure 4). One can see the effect of the neighborhood shape by recognizing the diagonal cross or square over the solutions.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Results of variational morphology for different ratios of $\alpha^+/\alpha^-$ and $p = 1.2$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Slice of the first result images for different ratios $\alpha^+/\alpha^-$</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Results for variational erosion and dilation for 3D graph meshes (Top left : original image, Top right : median filter ($-\alpha^+/\alpha^- = 1$), Bottom left : variational erosion, Bottom right : variational dilation</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Results for variational erosion and dilation with a diagonal cross structuring element</figcaption></figure>
### 3 Variational Formulation of Opening and Closing

In this section, we consider the numerical resolution of the variational formulation of the morphological opening and closing operators proposed in [8]:

$$ \min_{(x,t) \in (\mathbb{R}^n)^2} \sum_{i=1}^{n} \left( \sum_{j \in N_i} \varphi_i(\sigma_{\Omega_i}(t_i - y_j - w_{ij})) + \sum_{j \in \tilde{N}_i} \psi_i(\sigma_{\Lambda_i}(x_i - t_j - \nu_{ij})) \right). \quad (4) $$

This formulation is inspired by classical operators opening and closing, which are compositions of the classical erosion and dilation. A naive solution to problem (4) was first used in [8] composing the result of variational erosion and dilation to find an approximation of variational opening and closing. This solution first solves for $t$, then for $x$ ignoring the coupling of the variables. While this approach gives a result close to the optimal solution and is theoretically founded under some asymptotic conditions [8], it cannot be considered a proper variational opening/closing. We will refer to this naive approach to compare the solutions computed in this section (see Figures 5 and 6).

The operators are defined through a more complex convex problem involving two variables per node. Specifically, we solve the following optimization problem for variables $x = (x_i)_{i=1}^n$ and $t = (t_i)_{i=1}^n$. We remark that this formulation's objective function, while convex, couples the variables $x$ and $t$, in contrast to (1). We can rewrite our problem in the form:

$$ \min_{(x,t) \in (\mathbb{R}^n)^2} \bar{F}(t) + \bar{G}(x,t) \quad (5) $$

with

$$ \bar{F}(t) = \sum_{i=1}^{n} \hat{F}_i(t_i) \\ (\forall i \in \{1, \dots, n\}) \quad \hat{F}_i(t_i) = \sum_{j \in N_i} \underbrace{\varphi_i(\sigma_{\Omega_i}(t_i - y_j - w_{ij}))}_{F_{ij}(t_i)} $$

and

$$ \bar{G}(x,t) = \sum_{i=1}^{n} \sum_{j \in \tilde{N}_i} \psi_i(\sigma_{\Lambda_i}(x_i - t_j - \nu_{ij})). $$

While the variational formulation of erosion/dilation (1) leads to an optimal-ity condition which is completely separable, the structure of the objective func-tion in (5) is not as simple. Nonetheless, each function $F$ and $G$ is a sum that is amenable to proximal splitting methods [9]. Indeed, due to the nature of our choice for $\Phi_i$ and $\sigma_{\Omega_i}$, we can compute the proximal operator associated to the sum over each neighborhood or reciprocal neighborhood. This allows us to consider proximal splitting algorithms [4,5] such as Alternating Direction Method
of Multipliers (ADMM) [11] and Parallel Proximal Point Algorithm (PPXA+) [10] to solve (4). The theory of monotone operators which encompasses optimization with proximal operators can be found in [1]. A comprehensive repository of proximal operators can be found in [3].

Since each $\hat{F}_i$ is a function of only $t_i$, the proximal operator associated to $\bar{F}$ can be computed by computing the proximal operator associated to each $\hat{F}_i$ node by node,

$$ (\forall i \in \{1, \dots, n\}) \quad [\text{prox}_{\lambda \bar{F}}(t)]_i = \text{prox}_{\lambda \hat{F}_i}(t_i). $$

The proximal operator associated to each $\hat{F}_i$ can itself be computed by solving a scalar optimization problem, which amounts to solving the following inclusion problem

$$ \begin{aligned} \text{prox}_{\lambda \hat{F}_i}(t_i) = z_i &\iff 0 \in \lambda \partial \hat{F}_i(z_i) + z_i - t_i = 0 \\ &\iff 0 \in \lambda \left( \sum_{j \in \tilde{N}_i} \partial F_{ij}(z_j) \right) + z_i - t_i = 0. \end{aligned} \quad (6) $$

For $\bar{G}$, such a splitting is no longer possible because of the coupling between the variables $t$ and $x$. In light of this, we can define, for all $i \in \{1, \dots, n\}$ and $j \in \tilde{N}_i$, $u_{ij} = x_i - t_j$ and

$$ G_{ij}(u_{ij}) = \psi_i(\sigma_{\Lambda_i}(u_{ij} - \nu_{ij})). $$

Since each $G_{ij}$ only depends on $u_{ij}$, we can rewrite $\bar{G}$ in the following, separable way,

$$ \bar{G}(x, t) = G(u) = \sum_{i=1}^{n} \sum_{j \in \tilde{N}_i} G_{ij}(u_{ij}). $$

From this, we can deduce the following expression of the proximal operator associated to $G$ given componentwise, for all $i \in \{1, \dots, n\}$ and $j \in \tilde{N}_i$, by

$$ [\text{prox}_{\lambda G}(u)]_{ij} = \text{prox}_{\lambda G_{ij}}(u_{ij}). $$

As is the case for $\text{prox}_{\lambda \hat{F}_i}$, each proximal operator $\text{prox}_{\lambda G_{ij}}(u)$ can be computed by solving the following inclusion problem,

$$ \text{prox}_{\lambda G_{ij}}(u_{ij}) = z_{ij} \iff 0 \in \lambda \partial G_{ij}(u_{ij}) + u_{ij} - z_{ij}, \quad (7) $$

which we will describe in more detail later in the paper.

As we have remarked earlier, these different splitting strategies are necessary because of the coupling induced by the function $\bar{G}$. We can rewrite (4) with respect to the lifted variables $(X, U) \in (\mathbb{R}^n)^2 \times \mathbb{R}^{n \times m}$,

$$ \min_{X \in (\mathbb{R}^n)^2} F(X) + G(\underbrace{LX}_{U}), \quad (8) $$
where $X := (x, t) \in (\mathbb{R}^n)^2$, $F(X) = \bar{F}(t)$, $U = (u_{ij})_{i,j} \in \mathbb{R}^{n \times m}$, $L : (\mathbb{R}^n)^2 \to \mathbb{R}^{n \times m}$ is defined by the differences $x_i - t_j$ for every $i \in \{1, \dots, n\}$ and $j \in \tilde{N}_i$, and $m = \max_{i \in \{1, \dots, n\}} |\tilde{N}_i|$. The proximal splitting algorithms [4] that we employ will make use of the adjoint operator associated to the linear operator $L$. This adjoint operator $L^* : \mathbb{R}^{n \times m} \to (\mathbb{R}^n)^2$ maps a vector $u = (u_{ij})_{i,j} \in \mathbb{R}^{n \times m}$ to a pair $(x^*, t^*) \in (\mathbb{R}^n)^2$ as

$$ (\forall i \in \{1, \dots, n\}) \quad \begin{cases} x_i^* = \sum_{j \in \tilde{N}_i} u_{ij} \\ t_j^* = - \sum_{i \mid j \in \tilde{N}_i} u_{ij}. \end{cases} \qquad (9) $$

Once the problem is rewritten in this way, applying proximal splitting methods like ADMM or PPXA+ becomes straightforward as long as we can compute the linear operators $L$ and $L^*$, and the proximal operators $\text{prox}_{\lambda G_{ij}}$ and $\text{prox}_{\lambda \hat{F}_i}$, which we will describe in the following subsection.

## 3.1 Explicit Form of the Proximal Operator Computation

We have seen that using the proximal splitting methods requires two different splitting strategies on $\bar{F}$ and $\bar{G}$. At the end, the computation of the proximal operators $\text{prox}_{\lambda G_{ij}}(u_{ij})$ and $\text{prox}_{\lambda \hat{F}_i}(v_i)$ reduces to solving two independent root-finding equations, in the case when $p > 1$, $\psi_i = \varphi_i$, and $\Lambda_i = \Omega_i$.

Using (7), we have

$$ x = \text{prox}_{\lambda G_{ij}}(u_{ij}) \iff 0 = \frac{1}{\lambda}(x - u_{ij}) + \begin{cases} p_i \alpha_i^{-} (\alpha_i^{-} (x - \nu_{ij}))^{p-1} & \text{if } x < \nu_{ij} \\ p_i (\alpha_i^{+})^{p} (x - \nu_{ij})^{p-1} & \text{if } x > \nu_{ij} \\ 0 & \text{if } x = w_{ij}. \end{cases} $$

Using (6), we have

$$ x = \text{prox}_{\lambda \hat{F}_i}(t_i) \iff 0 = \frac{1}{\lambda}(x - t_i) + \sum_{j \in \tilde{N}_i} \begin{cases} p_i \alpha_i^{-} (\alpha_i^{-} (x - w_{ij}))^{p-1} & \text{if } x < w_{ij} \\ p_i (\alpha_i^{+})^{p} (x - w_{ij})^{p-1} & \text{if } x > w_{ij} \\ 0 & \text{if } x = w_{ij}. \end{cases} $$

Since $G_{ij}$ and $\hat{F}_i$ are differentiable and convex for $p > 1$, these two problems admit a unique minimizer. They are similar in form to the subproblems solved for the variational erosion/dilation in (1). These root-finding problems will be handled via a zeroth-order method, exactly as described in the computation of the variational formulation of erosion and dilation, since the functions encountered are of the same form. The proximal operators being computed term-by-term for the functions $\bar{F}$ and $\bar{G}$, we can now make use of proximal splitting algorithms to solve the variational opening/closing problem (4).
## 3.2 ADMM Algorithm

To solve the variational opening/closing problem (8) using ADMM, we introduce the auxiliary variable $U \in \mathbb{R}^{n \times m}$ and reformulate the problem as

$$ \min_{\substack{X \in (\mathbb{R}^n)^2 \\ U \in \mathbb{R}^{n \times m}}} F(X) + G(U) \quad \text{subject to} \quad LX - U = 0. \qquad (10) $$

This is now a linearly constrained convex optimization problem to which we apply the linearized-ADMM algorithm as detailed in [11]; for simplicity we will refer to the linearized-ADMM as ADMM for the rest of the paper. We present the algorithm in its entirety below.

<caption>**Algorithm 1 ADMM for Opening/Closing**</caption>

1: Initialize: $X_0 = (x_0, t_0) \in (\mathbb{R}^n)^2$, $U_0 \in \mathbb{R}^{n \times m}$, $z_0 \in \mathbb{R}^{n \times m}$
2: **for** $k = 0, 1, \dots$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;$X_{k+1} = \text{prox}_{\mu F}(X_k - \frac{\mu}{\lambda}L^*(z_k + LX_k - U_k))$
4: &nbsp;&nbsp;&nbsp;&nbsp;$U_{k+1} = \text{prox}_{\lambda G}(LX_{k+1} + z_k)$
5: &nbsp;&nbsp;&nbsp;&nbsp;$z_{k+1} = z_k + LX_{k+1} - U_{k+1}$.
6: **end for**

Given step size parameters $\mu > 0$ and $\lambda > 0$ and a dual variable $z_k \in \mathbb{R}^{n \times m}$, the algorithm updates the variables $X$ and $U$ by applying, at each iteration, the proximal operator associated to either $F$ or $G$, respectively. This algorithm is guaranteed to converge under standard assumptions; more specifically, we can ensure convergence under appropriate step sizes $\mu > 0$ and $\lambda > 0$ since the constraint $LX - U = 0$ is linear and nonempty and the functions $F$ and $G$ are proper, lower semicontinuous, and convex.

The choice of step sizes $\lambda$ and $\mu$ influences the convergence speed. One choice for which convergence is guaranteed is given by

$$ 0 < \mu < \frac{\lambda}{\|L\|^2} $$

where $\|L\|$ is the operator norm of $L$. In our case, $L$ is a sparse difference operator coupling the $x_i$ and $t_j$ variables. While computing the exact norm is costly, an upper bound can be estimated from the maximum number of neighbors for a node.

ADMM offers a flexible approach to solving the variational formulation of opening/closing. Its convergence is stable under our hypotheses, and its iterations are simple to implement due to the proximal calculus and separable structure of the functions involved. When the graph domain is fixed, such as in the case of images of a given size, the linear operator $L$ and its adjoint $L^*$ depend only on the chosen neighborhood (or structuring element), so that these operators (and potentially their factorizations) can be precomputed and stored to facilitate their use in ADMM.
**Algorithm 2 PPXA+ (as in [10, Prop. 5.2])**

1: Initialize: $x_0 \in (\mathbb{R}^n)^2$, $z_{0,1} \in (\mathbb{R}^n)^2$, $z_{0,2} \in \mathbb{R}^{n \times m}$

2: Set:

$$x_0 \in \operatorname{argmin}_{y \in (\mathbb{R}^n)^2} \sum_{i=1}^{2} \omega_i \|L_i y - z_{0,i}\|^2$$

3: **for** $n = 0, 1, 2, \dots$ **do**

4: &nbsp;&nbsp;&nbsp;&nbsp; $p_{n,1} = \operatorname{prox}_{F/\omega_1}(z_{n,1})$

5: &nbsp;&nbsp;&nbsp;&nbsp; $p_{n,2} = \operatorname{prox}_{G/\omega_2}(z_{n,2})$

6: Compute:

$$c_n \in \operatorname{argmin}_{y \in (\mathbb{R}^n)^2} \sum_{i=1}^{2} \omega_i \|L_i y - p_{n,i}\|^2$$

7: **for** $i = 1, 2$ **do**

8: &nbsp;&nbsp;&nbsp;&nbsp; $z_{n+1,i} = z_{n,i} + \lambda_n (L_i(2c_n - x_n) - p_{n,i})$

9: **end for**

10: $x_{n+1} = x_n + \lambda_n(c_n - x_n)$

11: **end for**

## 3.3 PPXA+ Algorithm

Another proximal splitting algorithm that is amenable to our problem is the so-called PPXA+ algorithm. This algorithm, introduced in [10], is designed to solve convex optimization problems of the form

$$\min_{X \in (\mathbb{R}^n)^2} F(L_1 X) + G(L_2 X) \qquad (11)$$

where $F, G$ are proper lower-semicontinuous convex functions, and $L_1, L_2$ are bounded linear operators. The method computes proximal operators of each term in parallel and aggregates their results via an update step. To apply PPXA+ to (4), we must first rewrite it by using the lifted variable $X := (x, t) \in (\mathbb{R}^n)^2$ and introducing the bounded linear operators $L_1 := \mathrm{Id}$ and $L_2 := L$. Then, the variational formulation for opening and closing in (8) can be rewritten as (11).

Given a sequence of relaxation parameters $(\lambda_n)_{n \in \mathbb{N}}$, with $\lambda_n \in (0, 2)$, satisfying $\lambda_{n+1} \le \lambda_n$ and strictly positive weights $\omega_1, \omega_2$, with $\omega_1 + \omega_2 = 1$, the algorithm computes the prox for $F$ and $G$ separately at each iteration and then performs a sort of averaging step by solving a least-squares problem.

One of the main differences between ADMM and PPXA+ iterations is the computation of $c_n$ that must be performed in PPXA+. This step requires solving a least squares problem of the form

$$c_n = \operatorname{argmin}_{z \in (\mathbb{R}^n)^2} \omega_1 \|z - p_{n,1}\|^2 + \omega_2 \|L_z - p_{n,2}\|^2.$$

This is a strictly convex quadratic problem, so its solution satisfies the normal equations

$$(\omega_1 \mathrm{Id} + \omega_2 L^* L) c_n = \omega_1 p_{n,1} + \omega_2 L^* p_{n,2}.$$
We solve this system using a linear solver (e.g., conjugate gradient) or, in smaller-
scale cases, a direct method using sparse linear algebra. The operator $\omega_1 \mathrm{Id} + \omega_2 L^* L$
is symmetric positive definite by construction, ensuring that the system is well
conditioned and uniquely solvable. The complexity of solving this linear system
can worsen the performance of PPXA+ relatively to ADMM. However, we have
found that this difference is not substantial for large images since the time re-
quired by both algorithms is dominated by the cost of the proximal operators
involved.

In PPXA+, all proximal subproblems are decoupled and computed indepen-
dently, making the algorithm naturally parallelizable between pixels or graph
elements. Furthermore, PPXA+ offers a way to parallelize both proximal com-
putations at the same time and can potentially be faster than ADMM.

Convergence for PPXA+ is guaranteed under the considered assumptions.
Since $L$ is sparse and $L^*$ has a known structure (see (9)), the matrix vector prod-
ucts involved in the solver are efficient. Moreover, for fixed graph/neighborhood
structure, $L^*L$ can be precomputed (and possible factorized) or reused across
iterations.

Now, we show a comparison in performance for both algorithms and with
the naive solution. We can see that both algorithms converge to an objective
value that is lower than, but still close to, the naive solution. However the pro-
cessed images are visually similar for both solutions of (4), although different
from the naive image (see Figures 5 and 6). This shows that even though the
naive method can come really quickly to an almost optimal objective value, it
produces images that are distinct from the optimal solution. We observe partic-
ularly that the naive method images are less smooth and noisier than the ones
generated by ADMM and PPXA+ algorithms.

4 Conclusion

We have detailed the implementation and numerical resolution of a variational
framework for morphological operators on graphs. For the variational erosion
and dilation, we use the separability of the variational formulation to decom-
pose its solution into the solutions of scalar convex problems, solved efficiently
using zeroth-order root-finding methods. For the variational opening and clos-
ing, we show that proximal operator splitting techniques are an effective ap-
proach, specifically demonstrating the PPXA+ and ADMM algorithms, as they
are adapted to the structure of the problem. The methods presented here pro-
vide flexible and efficient tools for using variational morphological operators on
graphs. They can be easily implemented as layers in deep neural networks, and
can thus open up possibilities for further research, such as integrating adaptive
or learned parameters. Another interesting perspective is the possibility that
this approach be used as part of more sophisticated processing pipeline, where
parameters can be adjusted or learned depending on the context. In future work
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Objective functions in the resolution of a variational closing processing for ADMM, PPXA+, and naive algorithms</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Results of variational opening for ADMM, PPXA+, and naive algorithms</figcaption></figure>
we will also use these implementations of morphological operators for solving
inverse problems using MM operators via variational methods [6].

References

1. Bauschke, H.H., Combettes, P.L.: Convex Analysis and Monotone Operator Theory in Hilbert Spaces. Springer, 2 edn. (2017)

2. Burden, R.L., Faires, J.D.: Numerical Analysis. The Prindle, Weber and Schmidt Series in Mathematics, PWS-Kent Publishing Company, Boston, fourth edn. (1989)

3. Chierchia, G., Chouzenoux, E., Combettes, P.L., Pesquet, J.C.: The proximity operator repository (2025), http://proximity-operator.net/download/guide.pdf

4. Combettes, P.L., Pesquet, J.C.: Proximal splitting methods in signal processing. Fixed-point algorithms for inverse problems in science and engineering pp. 185–212 (2011)

5. Combettes, P.L., Pesquet, J.C.: Fixed point strategies in data science. IEEE Transactions on Signal Processing pp. 3878–3905 (2021)

6. Gharbi, M., Villa, S., Chouzenoux, E., Pesquet, J.C., Duval, L.: Unrolled deep networks for sparse signal restoration in analytical chemistry. In: 2024 IEEE 34th International Workshop on Machine Learning for Signal Processing (MLSP). pp. 1–6. IEEE (2024)

7. Najman, L., Talbot, H. (eds.): Mathematical Morphology: from theory to applications. ISTE-Wiley, London, UK (September 2010), ISBN 978-1848212152

8. Najman, L., Pesquet, J.C., Talbot, H.: When convex analysis meets mathematical morphology on graphs. In: Mathematical Morphology and Its Applications to Signal and Image Processing. vol. 9082, pp. 473–484. Springer LNCS (2015)

9. Parikh, N., Boyd, S.: Proximal algorithms. Found. Trends Optim. 1(3), 127–239 (Jan 2014)

10. Pesquet, J.C., Pustelnik, N.: A parallel inertial proximal optimization method. Pacific Journal of Optimization 7(2), 273–305 (2011)

11. Ryu, E.K., Yin, W.: Large-Scale Convex Optimization: Algorithms and Analyses via Monotone Operators. Cambridge University Press (2022)

12. Serra, J.: Image analysis and mathematical morphology. Academic Press, Inc. (1983)