A Partially Derivative-Free Proximal Method for Composite Multiobjective Optimization in the H\"older Setting

TL;DR

This paper introduces a derivative-free proximal method for multiobjective composite optimization problems under Hölder conditions, providing convergence guarantees and demonstrating practical effectiveness through numerical experiments.

Contribution

It proposes a novel partially derivative-free algorithm with convergence analysis for composite multiobjective optimization in the Hölder setting.

Findings

Converges to a weakly ε-approximate Pareto point within O(ε^{-(β+1)/β}) iterations.

Balances computational accuracy and efficiency using gradient approximations and a scaling matrix.

Shows competitive performance on benchmark problems.

Abstract

This paper presents an algorithm for solving multiobjective optimization problems involving composite functions, where we minimize a quadratic model that approximates $F(x) - F(x^k)$ and that can be derivative-free. We establish theoretical assumptions about the component functions of the composition and provide comprehensive convergence and complexity analysis. Specifically, we prove that the proposed method converges to a weakly $\varepsilon$-approximate Pareto point in at most $\mathcal{O}\left(\varepsilon^{-\frac{\beta+1}{\beta}}\right)$ iterations, where $\beta$ denotes the H\"{o}lder exponent of the gradient. The algorithm incorporates gradient approximations and a scaling matrix $B_k$ to achieve an optimal balance between computational accuracy and efficiency. Numerical experiments on a collection of benchmark problems are presented, illustrating the practical behavior of the proposed approach and its competitiveness with existing composite algorithms.

Tables

Table 1. Table 1: Number of solutions found by PDFPM (200 runs per setting).

Problem	$\delta =0.0$	$\delta =0.02$	$\delta =0.05$	$\delta =0.1$
AAS1	196 (98%)	199 (99.5%)	199 (99.5%)	195 (97.5%)
AAS2	199 (99.5%)	200 (100%)	200 (100%)	200 (100%)

Table 2. Table 2: List of test problems used in the comparison.

Problem Source $n$ $m$ AP1 [ansary] 2 3 AP2 [ansary] 1 2 AP3 [ansary] 2 2 AP4 [ansary] 3 3 BK1 [reviewproblems] 2 2 DD1 [dd1] 5 2 DGO1 [reviewproblems] 1 2 DGO2 [reviewproblems] 1 2 DTLZ1 [Deb2005] 7 3 DTLZ2 [Deb2005] 7 3 DTLZ3 [Deb2005] 7 3 DTLZ4 [Deb2005] 7 3 DTLZ5 [Deb2005] 9 5 FA1 [reviewproblems] 3 3 Far1 [reviewproblems] 2 2 FDS [fliegeanddrummond2009] 5 3 FF1 [reviewproblems] 2 2 Hil1 [hill] 2 2 IKK1 [reviewproblems] 2 3 IM1 [reviewproblems] 2 2 JOS1 [10.5555/2955239.2955427] 100 2 JOS4 [10.5555/2955239.2955427] 100 2 KW2 [kw2] 2 2 LE1 [reviewproblems] 2 2 Lov1 [doi:10.1137/100784746] 2 2 Lov2 [doi:10.1137/100784746] 2 2 Lov3 [doi:10.1137/100784746] 2 2 Lov4 [doi:10.1137/100784746] 2 2 Lov5 [doi:10.1137/100784746] 3 2 Lov6 [doi:10.1137/100784746] 6 2 LTDZ [laumanns2002combining] 3 3 MGH9 [moretest] 3 15 MGH16 [moretest] 4 5 MGH26 [moretest] 4 4 MGH33 [moretest] 10 10

Problem Source $n$ $m$ MHHM1 [reviewproblems] 1 3 MHHM2 [reviewproblems] 2 3 MLF1 [reviewproblems] 1 2 MLF2 [reviewproblems] 2 2 MMR1 [italianos] 2 2 MMR2 [italianos] 2 2 MMR3 [italianos] 2 2 MMR4 [italianos] 3 2 MOP2 [reviewproblems] 2 2 MOP3 [reviewproblems] 2 2 MOP5 [reviewproblems] 2 3 MOP6 [reviewproblems] 2 2 MOP7 [reviewproblems] 2 3 PNR [pnr] 2 2 QV1 [reviewproblems] 10 2 SD [doi:10.2514/5.9781600866234.0209.0249] 4 2 SK1 [reviewproblems] 1 2 SK2 [reviewproblems] 4 2 SLCDT1 [slcdt] 2 2 SLCDT2 [slcdt] 10 3 SP1 [reviewproblems] 2 2 SSFYY2 [reviewproblems] 1 2 TKLY1 [reviewproblems] 4 2 Toi4 [tointtest] 4 2 Toi8 [tointtest] 4 3 Toi9 [tointtest] 4 4 Toi10 [tointtest] 4 3 VU1 [reviewproblems] 2 2 VU2 [reviewproblems] 2 2 ZDT1 [ZDT] 30 2 ZDT2 [ZDT] 30 2 ZDT3 [ZDT] 30 2 ZDT4 [ZDT] 30 2 ZDT6 [ZDT] 10 2 ZLT1 [reviewproblems] 10 5

Equations

$$\underset{x\in {\mathbb{R}}^n}{\min }F(x)$$

$$F_j:=f_j+h_j,j\in \mathcal{J}:=\{1,\dots ,m\},$$

$$\underset{x\in \mathcal{C}}{\min }(f_1(x),\dots ,f_m(x)).$$

$${\Phi }_{x^k}(x)+\frac{{\sigma }_k}{2}\Vert x-x^k{\Vert }^2$$

$${\Phi }_{x^k}(x):=\underset{j\in \mathcal{J}}{\max }[⟨g_{f_j}(x^k,{\lambda }_k)+\frac{1}{2}{B}_j^k(x-x^k),x-x^k⟩+h_j(x)-h_j(x^k)],$$

$${\mathbb{R}}_{+}^m:=\{u\in {\mathbb{R}}^m\mid u_j\ge 0,\forall j\in \mathcal{J}\}\text{and}{\mathbb{R}}_{++}^m:=\{u\in {\mathbb{R}}^m\mid u_j>0,\forall j\in \mathcal{J}\}.$$

$$F(x)\le F(x^{\ast })\text{mbox}andF(x)\ne F(x^{\ast }).$$

$$0\in \sum _{j=1}^m{\gamma }_j(\nabla f_j(x^{\ast })+\partial h_j(x^{\ast })),$$

$$\partial h_j(\bar{x})=\{y\in {\mathbb{R}}^n:h_j(z)\ge h_j(\bar{x})+⟨y,z-\bar{x}⟩,\forall z\in {\mathbb{R}}^n\}.$$

$$\delta_{j}^{\mathcal{C}}(x)\!=\!\left\{\hskip-5.0pt\begin{tabular}[]{ll}$0$,&\hskip-12.0pt if $x\in\mathcal{C}$,\\ $+\infty$,&\hskip-12.0pt otherwise,\end{tabular}\right.\,\,j\in\mathcal{J}.$$

$$N_{\mathcal{C}}(x):=\{z\in {\mathbb{R}}^n:⟨z,y-x⟩\le 0,\forall y\in \mathcal{C}\}.$$

$$\Vert \nabla \phi (y)-\nabla \phi (x)\Vert \le L\Vert y-x{\Vert }^{\theta },\forall x,y\in {\mathbb{R}}^n.$$

$$g_{f_j}:{\mathbb{R}}^n\times [0,1]\to {\mathbb{R}}^n\text{where}\underset{\lambda \to 0}{\lim }{g}_{f_j}(x,\lambda )=\nabla f_j(x).$$

$$g_{f_j}(x,\lambda )=[\frac{f_j(x+\lambda e_1)-f_j(x)}{\lambda },\cdots ,\frac{f_j(x+\lambda e_n)-f_j(x)}{\lambda }];$$

$$g_{f_j}(x,\lambda )=[\frac{f_j(x)-f_j(x-\lambda e_1)}{\lambda },\cdots ,\frac{f_j(x)-f_j(x-\lambda e_n)}{\lambda }];$$

$$g_{f_j}(x,\lambda )=[\frac{f_j(x+\lambda e_1)-f_j(x-\lambda e_1)}{2\lambda },\cdots ,\frac{f_j(x+\lambda e_n)-f_j(x-\lambda e_n)}{2\lambda }].$$

$$0<{\lambda }_k\le \frac{ϵ}{{\sigma }_k\sqrt{n}}$$

$$\underset{x\in {\mathbb{R}}^n}{\min }{\Phi }_{x^k}(x)+\frac{{\sigma }_k}{2}\Vert x-x^k{\Vert }^2,$$

$$F_j({\bar{x}}^k)\le F_j(x^k)-\frac{\alpha {ϵ}^2}{2{\sigma }_k},\forall j\in \mathcal{J},$$

$${\Psi }_{{\sigma }_k}^{x^k}(x)={\Phi }_{x^k}(x)+\frac{{\sigma }_k}{2}\Vert x-x^k{\Vert }^2$$

$$\underset{\Vert x\Vert \to +\infty }{\lim }{\Psi }_{{\sigma }_k}^{x^k}(x)=+\infty ,$$

$$f_j(y)\le f_j(x)+⟨\nabla f_j(x),y-x⟩+\frac{L_j}{2}\Vert y-x{\Vert }^2+\frac{M_j}{{\beta }_j+1}\Vert y-x{\Vert }^{{\beta }_j+1},\forall x,y\in {\mathbb{R}}^n.$$

$$s_j(y)\le s_j(x)+⟨\nabla s_j(x),y-x⟩+\frac{L_j}{2}\Vert y-x{\Vert }^2,\forall x,y\in {\mathbb{R}}^n,j\in \mathcal{J},$$

$$r_j(y)\le r_j(x)+⟨\nabla r_j(x),y-x⟩+\frac{M_j}{{\beta }_j+1}\Vert y-x{\Vert }^{{\beta }_j+1},\forall x,y\in {\mathbb{R}}^n,j\in \mathcal{J}.$$

$$\Vert \nabla f_j(x)-g_{f_j}(x,\lambda )\Vert \le \lambda \frac{\sqrt{n}{L}_j}{2}+\sqrt{n}\frac{M_j}{{\beta }_j+1}{\lambda }^{{\beta }_j}.$$

$$\Vert s_j(y)-s_j(x)-\nabla s_j(x{)}^{\top }(y-x)\Vert \le \frac{L_j}{2}\Vert y-x{\Vert }^2$$

$$\Vert r_j(y)-r_j(x)-\nabla r_j(x{)}^{\top }(y-x)\Vert \le \frac{M_j}{1+{\beta }_j}\Vert y-x{\Vert }^{{\beta }_j+1}.$$

$$|\frac{s_j(x+\lambda e_i)-s_j(x)}{\lambda }-\frac{\partial s_j(x)}{\partial x_i}|$$

$$|\frac{r_j(x+\lambda e_i)-r_j(x)}{\lambda }-\frac{\partial r_j(x)}{\partial x_i}|$$

$$\Vert \nabla s_j(x)-g_{s_j}(x,\lambda )\Vert \le \frac{\sqrt{n}{L}_j}{2}\lambda \text{and}\Vert \nabla r_j(x)-g_{r_j}(x,\lambda )\Vert \le \frac{\sqrt{n}{M}_j}{{\beta }_j+1}{\lambda }^{{\beta }_j}.$$

Full text

MSC Classification]90C29, 65K05, 90C30

1 Introduction

Multiobjective optimization is characterized by the simultaneous minimization (or maximization) of several objective functions. In general, these objectives are conflicting, which prevents the existence of a single point capable of optimizing all objective functions at once. Therefore, the concept of Pareto optimality is employed to characterize solutions of multiobjective problems. Problems that deal with multiple objectives simultaneously are common in various fields, including engineering, industrial processes, telecommunications, and game theory, as shown in the works [BhaskarGuptaRay, BeltonStewart2002, Miettinen1998, Stewart2008, Martinez-MunozMartiYepes2025, Rgamemontc].

There are several iterative approaches for solving multiobjective optimization (MO) problems. To the best of our knowledge, there are three principal approaches: scalarization [Eichfelder2008], heuristic [laumanns2002combining], and direct methods, the latter also being considered an extension of scalar-valued methods to multiobjective (vector) optimization. This approach was initiated in [Fliege2000], where the authors extended the steepest descent method; see also [fliegeanddrummond2009]. Over the past two decades, various direct methods have sought to extend methods such as Newton [chuong2013newton, fliegeanddrummond2009, mauricionewtonvetorial, wang2019extended], quasi-Newton [ansary, math8040616, Mahdavi-AmiriSadaghiani2020, Morovati2017, POVALEJ2014765, QuGohChan2011, QU2014503, Prudente2022, Prudente2024], conjugate gradient [Goncalves2019, cgvo], projected gradient [luis-jef-yun, fazzio2019convergence, fukuda2011convergence, fukuda2013inexact, mauricio&iusem], proximal methods [bonnel2005proximal, ceng2010hybrid, ceng2007approximate, chuong2011generalized, chuong2011hybrid, fukuda2019, fukuda2023, Nonmonotone-Peng2025, CHEN2025116422, Bello-Cruz2025], (generalized) conditional gradient [Assuncao2021, Assuncao2023, Goncalves2025], and Newton-type proximal gradient [Ansary04052023].

In many practical optimization problems, direct access to the derivatives of the functions involved is often unavailable, which makes their computation challenging. Additionally, evaluating the objective function and its derivatives can be computationally expensive. As a result, there is growing interest in methods that offer low computational cost while maintaining good worst-case complexity in terms of function evaluations and iterations. In this context, research on derivative-free methods has advanced significantly for scalar problems. For further details, see [Amaral2025, Grapiglia2023, Grapiglia2024] and the references therein.

The objective of this work is to develop a partially derivative-free method for solving the nonsmooth multiobjective optimization problem defined as follows:

[TABLE]

where $F:\mathbb{R}^{n}\to\mathbb{R}^{m}$ is a vector-valued function with $F(x)=(F_{1}(x),\ldots,F_{m}(x))$, and each component is given by

[TABLE]

with $f_{j}:\mathbb{R}^{n}\to\mathbb{R}$ differentiable and $h_{j}:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$ convex, but not necessarily differentiable. The problem of minimizing a function expressed as the sum of two functions has several applications; see [fukuda2019, Section 5.2]. Let $\mathcal{C}\subset\mathbb{R}^{n}$ be a compact and convex set such that, for all $j\in\mathcal{J}$, $h_{j}(x)=0$ for all $x\in\mathcal{C}$ and $h_{j}(x)=+\infty$ otherwise; then $h_{j}$ is the indicator function of $\mathcal{C}$. Thus, the problem in (1) is equivalent to solving the following problem:

[TABLE]

In addition, the problem in (1) can be applied to robust optimization problems, which involve uncertainties in their parameters and essentially consist of solving the problem in the worst-case scenario; see [Chen2019, fukuda2019] and their references.

Recently, composite problems such as in (1) have been addressed using various methods, including proximal gradient [fukuda2019, Nonmonotone-Peng2025, Bello-Cruz2025, CHEN2025116422], generalized conditional gradient [Assuncao2023], and Newton-type proximal gradient [Ansary04052023]. In all references, the component functions $f_{j}:\mathbb{R}^{n}\to\mathbb{R}$ for $j\in\mathcal{J}$ are continuously differentiable. In [fukuda2023, Assuncao2023], the authors further assumed that $f_{j}$ have Lipschitz continuous gradients, while the functions $h_{j}:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$ are proper and convex. Notably, under the assumption that $h_{j}$ are lower semicontinuous, Assunção et al. [Assuncao2023] proposed a generalized version of the conditional gradient method (Frank-Wolfe) for multiobjective composite optimization. In [Ansary04052023, Bello-Cruz2025], the functions $f_{j}$ are assumed to be convex, and $h_{j}$ are considered continuous and convex; Bello-Cruz [Bello-Cruz2025] additionally assumed that $h_{j}$ are proper. In contrast, [CHEN2025116422, Nonmonotone-Peng2025] addressed cases where $f_{j}$ are continuously differentiable but not necessarily convex, while $h_{j}$ are convex (and proper, lower semicontinuous in [CHEN2025116422]).

Given an iteration $x^{k}$, our approach finds $x^{k+1}$ such that $F_{j}(x^{k+1})$ is sufficiently smaller than $F_{j}(x^{k})$ for all $j\in\mathcal{J}$. As discussed by Calderón et al. [Lizet], this can be done by defining an approximate model for each $F_{j}(x)-F_{j}(x^{k})$ around $x^{k}$ and approximately minimizing the maximum of these models. In our case, to approximate

[TABLE]

where

[TABLE]

$g_{f_{j}}(x^{k},\lambda_{k})$ is an approximation of $\nabla f_{j}(x^{k})$ and $B_{j}^{k}$ is a symmetric positive semidefinite matrix. To our knowledge, Fliege and Svaiter [Fliege2000] were the first to use this approach to solve multiobjective optimization problems.

One of the novelties of our approach lies in the extension of the derivative-free method described by Grapiglia in [Grapiglia2024] to multiobjective optimization. To the best of our knowledge, our work is one of the pioneers in applying this type of method to multiobjective optimization. In [Grapiglia2024], Grapiglia proposed a quadratic regularization method to minimize a function $f:\mathbb{R}^{n}\to\mathbb{R}$, where he used finite-difference approximations for the gradients and proved that, when the objective function is bounded below and has Lipschitz gradients, the method requires at most $\mathcal{O}(\epsilon^{-2})$ iterations to generate an approximate stationary point with $\epsilon$ accuracy.

In our work, we assume a Hölder-Lipschitz composition for the gradient of $f_{j}$ and convexity only in $h_{j}$. We show that the proposed method performs at most $\mathcal{O}\left(\epsilon^{-\frac{\beta+1}{\beta}}\right)$ iterations to find an $\epsilon$-approximate Pareto critical point for the problem in (1), where $\beta=\min\{\beta_{j}\mid j=1,\dots,m\}$ and $\beta_{j}$ are the exponents of the assumed Hölder conditions. Our complexity bound is order-wise consistent with the bound established by Calderón et al. [Lizet] for the first-order version of their method of order $p$, as well as with the complexity order obtained by Grapiglia et al. [PinheiroGrapiglia2024]. Furthermore, when $\beta_{1}=\dots=\beta_{m}=1$, i.e., when only Lipschitz conditions are assumed for the gradients of $f_{j}$, our bound becomes $O\left(\epsilon^{-2}\right),$ which is consistent with previously established bounds; see, for example, Grapiglia et al. [Grapiglia2015]. In the case $m=1$, our complexity order coincides with that obtained by Martinez [Martinez2017] for the first-order version of their method of order $p$. Finally, numerical experiments on robust biobjective problems with Lipschitz and Hölder-gradient components are reported. These experiments illustrate the method’s ability to approximate the Pareto front under different levels of uncertainty and to consistently recover diverse solutions, even in challenging settings where only Hölder continuity is assumed for the gradients.

The remainder of the article is structured as follows. First, in Section 2, we present the notation and important results for understanding the work. Next, in Section 3, we provide a detailed description of the proposed proximal point method and demonstrate its effectiveness. In Section 4, we analyze the iteration complexity of the method in the worst-case scenario. To further demonstrate the applicability of our approach, we provide numerical experiments in Section 5. Finally, in Section 6, we present the concluding remarks of the paper.

2 Notations and Auxiliary Results

We now introduce some notations, definitions, and results that will be used throughout the paper. We denote by $\mathbb{N}\coloneqq\{0,1,2,\ldots\}$ the set of nonnegative integers, and by $\mathbb{N}^{*}\coloneqq\{1,2,3,\ldots\}$ the set of positive integers. The sets $\mathbb{R}$, $\mathbb{R}_{+}$, and $\mathbb{R}_{++}$ represent the set of real numbers, the set of nonnegative real numbers, and the set of positive real numbers, respectively. Let $\mathcal{J}=\{1,\ldots,m\}$, we define:

[TABLE]

If $u,v\in{\mathbb{R}}^{m}$, then $v\leq u$ means that $u-v\in{\mathbb{R}}^{m}_{+}$, and if $v<u$ then $u-v\in{\mathbb{R}}^{m}_{++}$. The symbol $\langle\cdot,\cdot\rangle$ is the usual inner product in $\mathbb{R}^{n}$ and $\|\cdot\|$ denotes the Euclidean norm.

Definition 1**.**

We say that a point $x^{*}\in\mathbb{R}^{n}$ is a Pareto optimal point for the Problem (1) if there does not exist $x\in\mathbb{R}^{n}$ such that

[TABLE]

The Fritz John conditions for multiobjective optimization were extended in [Miettinen1998, Section 3.2.1 (First-Order Conditions)]. If $x^{*}$ is Pareto optimal for the Problem (1), then there exists $\gamma\in\mathbb{R}^{m}_{+}$ with $\sum_{j=1}^{m}\gamma_{j}=1$ such that

[TABLE]

where $\partial h_{j}$ is the subdifferential of the convex function $h_{j}:\mathbb{R}^{n}\to\mathbb{R}$ defined by

[TABLE]

Definition 2**.**

For a convex, closed and, non-empty set $\mathcal{C},$ we use $\delta^{\mathcal{C}}:\mathbb{R}^{n}\to\mathbb{R}^{m}$ as its indicator(vector) function, with $\delta^{\mathcal{C}}(x)=(\delta_{1}^{\mathcal{C}}(x),\ldots,\delta_{m}^{\mathcal{C}}(x))$ where,

[TABLE]

It is important to note that $\delta_{j}^{\mathcal{C}}(x)$ is a convex, proper, and lower semicontinuous function. Another important fact is that $\partial\delta_{j}^{\mathcal{C}}(x)(\cdot)=N_{\mathcal{C}}(\cdot),$ where

[TABLE]

Let $\varphi:\mathbb{R}^{n}\to\mathbb{R}$ be a differentiable function. We say that $\nabla\varphi$ is Hölder continuous with constant $L>0$ if there exists $\theta\in(0,1]$ such that

[TABLE]

If $\theta=1$, we say that $\nabla\varphi$ is Lipschitz continuous with constant $L>0$.

3 Proposed Method

For the development of the method proposed in this work, for each $j\in\mathcal{J}$, we consider the function

[TABLE]

We observe that the function in (4) can be defined as $g_{f_{j}}=\nabla f_{j}$. However, when the calculation of the gradient of $f_{j}$ is prohibitive, this definition is not recommended. In this case, it is preferable to define $g_{f_{j}}$ as an approximation of the gradient $\nabla f_{j}$ that is less expensive to obtain. Three possible definitions of $g_{f_{j}}$ are presented in the following remark.

Remark 3.1**.**

If $f_{j}$ is differentiable and $e_{i},\,\,i=1,\ldots,n$ denotes the $i$-th canonical vector of $\mathbb{R}^{n}$, then $g_{f_{j}}$ can be defined in one of the following ways:

• (i)

Forward difference:

[TABLE]

• (ii)

Backward difference:

[TABLE]

• (iii)

Central difference:

[TABLE]

It is important to note that, since $f_{j}$ is differentiable at $x$, each of the approximations defined above converges to the gradient $\nabla f_{j}(x)$ as $\lambda\to 0$. More precisely, we have $\lim_{\lambda\to 0}g_{f_{j}}(x,\lambda)=\nabla f_{j}(x),$ in all three cases (forward, backward, or central differences).

3.1 The Method

In this subsection, we present a partially derivative-free method, which can be made fully derivative-free, for solving the Problem (1).

The parameter $\alpha$ controls the reduction level of $F_{j}$ in (6), where values close to $1$ result in a more aggressive reduction, which may cause the penalty parameter $\sigma_{k}$ to grow unnecessarily. Therefore, it is crucial to tune $\alpha$ carefully. The descent test in (6) shows that small values of $\sigma_{k}$ tend to produce larger steps, which can speed up the convergence of the method. As will be shown in Remark 3.5, the criterion established in Step 3 of Algorithm 1 is appropriate.

The remainder of this subsection is devoted to establishing the well-definedness of the steps of Algorithm 1. To ensure the correct formulation of Step 2, we consider the following function:

[TABLE]

with $B^{k}_{j}$ symmetric positive semidefinite matrices. Since $h_{j}$ is convex, then there exists a linear function $\overline{h}_{j}(x)=\langle a,x\rangle+b,\,a,x\in\mathbb{R}^{n},\,b\in\mathbb{R}$ such that $h_{j}(x)\geq\overline{h}_{j}(x)$ for all $x\in\mathbb{R}^{n}$, hence we have that

[TABLE]

in other words, $\Psi_{\sigma}^{x^{k}}(x)$ is coercive. Furthermore, it is easy to see that the function $\Psi_{\sigma_{k}}^{x^{k}}$ is strongly convex, and, as a consequence, the Problem (5) has a unique solution, ensuring the well-definedness of Step 2 in Algorithm 1.

It is important to note that the matrix $B_{j}^{k}$ used in Step 2 of Algorithm 1 does not need to be an approximation of the Hessian $\nabla^{2}f_{j}(x^{k})$, as the expression in (2) might suggest. It can even be defined as the null matrix, relying only on first-order information, which is useful when the Hessian of $f_{j}$ is expensive to compute. However, second-order information, when available, can improve the efficiency of the step. The flexibility in the choice of $B_{j}^{k}$ allows us to use the Hessian, a cheap approximation (e.g., quasi-Newton), or the null matrix, as long as it is uniformly bounded, as presented in Assumption 1.

Next, we analyze the well-definedness of Step 4. For this purpose, we consider the following assumptions.

Assumption 1**.**

There exists $\overline{B}\geq 0$ such that $\|B_{j}^{k}\|\leq\overline{B}$ for all iterations $k$ of Algorithm 1 and $j\in\mathcal{J}$.

Assumption 2**.**

For each $j\in\mathcal{J}$, there exist $L_{j},M_{j}\in\mathbb{R}_{++}$ and $\beta_{j}\in(0,1]$ such that

[TABLE]

In Remark 3.2, we present a specific situation in which Assumption 2 is satisfied, thereby illustrating a concrete scenario in which this assumption is applicable.

Remark 3.2**.**

If $f_{j}=s_{j}+r_{j}$, where $s_{j}$ has a continuously differentiable gradient that is $L_{j}$-Lipschitz and $r_{j}$ has a continuously differentiable gradient that is $M_{j}$-Hölder continuous with exponent $\beta_{j}$, we have that

[TABLE]

and

[TABLE]

See, e.g., [Yashtini2016, Lemma 1], which directly implies (7).

Assumption 3**.**

For $L_{j}$ and $M_{j}$ given in Assumption 2, the following inequality holds for each $j\in\mathcal{J}$:

[TABLE]

In the following, we present situations that illustrate how Assumptions 2 and 3 are not restrictive and can be satisfied by a broad class of functions.

Remark 3.3**.**

It is important to highlight that, in Step 1, any function $g_{f_{j}}$ satisfying (4) and Assumption 3 may be used; it does not need to be one of the specific constructions presented in Remark 3.1.

In the following Lemma 3.1, we will demonstrate that, if $f_{j}$ satisfies the condition presented in Remark 3.2 and if we consider $g_{f_{j}}$ as one of the three definitions presented in Remark 3.1, then Assumption 3 holds.

Lemma 3.1**.**

Suppose $f_{j}=s_{j}+r_{j}$, where $s_{j}$ has a gradient that is $L_{j}$-Lipschitz continuous and $r_{j}$ has a gradient that is $M_{j}$-Hölder continuous with exponent $\beta_{j}$. Then the condition (8) is satisfied when considering $g_{f_{j}}(x,\lambda)$ in any of the definitions given in Remark 3.1.

Proof.

Considering $s_{j}$ with a gradient that is $L_{j}$-Lipschitz continuous and $r_{j}$ with a gradient that is $M_{j}$-Hölder continuous with exponent $\beta_{j}$, by [Yashtini2016, Lemma 1] we can conclude that

[TABLE]

and

[TABLE]

Using $y=x+\lambda e_{i}$ in (9) and (10), we have:

[TABLE]

and

[TABLE]

this implies that

[TABLE]

Using $y=x-\lambda e_{i}$ in (9) and (10), we have:

[TABLE]

and

[TABLE]

this implies that

[TABLE]

From the definitions given in Remark 3.1, we obtain the following:

[TABLE]

In an analogous way, we obtain

[TABLE]

this implies that

[TABLE]

Therefore, using the triangle inequality, we obtain

[TABLE]

concluding the proof of the lemma. ∎

Remark 3.4**.**

If there exists $\mathcal{L}\subset\mathcal{J}$ such that $\nabla f_{j}$ is Lipschitz for all $j\in\mathcal{L}$ and $\nabla f_{j}$ is Hölder for all $j\in\mathcal{J}\setminus\mathcal{L}$, then $f_{j}$ satisfies the conditions of Lemma 3.1. In this case, it suffices to define

[TABLE]

A scenario that satisfies Remark 3.4 is given below. For each $j\in\mathcal{J}$, let us consider $f_{j}:\mathbb{R}^{n}\to\mathbb{R}$ and a $\mathcal{L}\subset\mathcal{J}$ where

[TABLE]

where $A_{j}\in\mathbb{R}^{m_{j}\times n}$ is a positive semidefinite matrix, $b_{j}\in\mathbb{R}^{m_{j}}$, $D_{j}\in\mathbb{R}^{k_{j}\times n}$ is a linear operator, $\mu_{j}>0$, and $1<p<2$. It is straightforward to verify that the quadratic functions $f_{j}(x)=\tfrac{1}{2}\|A_{j}x-b_{j}\|_{2}^{2}$ possess Lipschitz continuous gradients. Furthermore, it is established in [Bernigaud2024] that the functions $f_{j}(x)=\tfrac{\mu_{j}}{p}\|D_{j}x\|_{p}^{p}$ have $(p-1)$-Hölder continuous gradients. The following Example 3.1 presents a particular instance of the function defined in (11).

Example 3.1**.**

Consider the function $f:\mathbb{R}\to\mathbb{R}^{2}$ defined as in (11), with $j=1,2$, where we take $A_{1}=I$, $b_{1}=0$, $D_{2}=I$, $\mu_{2}=1$, and $p=\tfrac{3}{2}$. In this setting, $f_{1}(x)=|x|^{2}$ and $f_{2}(x)=|x|^{\tfrac{3}{2}}$. We claim that $\nabla f_{2}$ is Hölder continuous but not Lipschitz. Indeed, for $x\in\mathbb{R}$,

[TABLE]

Taking $x=0$ and $y=\tfrac{1}{n}$, where $n\in\mathbb{N}$, we obtain

[TABLE]

If $\nabla f_{2}$ were Lipschitz, there would exist a constant $M>0$ such that

[TABLE]

so we would have

[TABLE]

This would imply that the set of natural numbers is bounded above, which is absurd. Therefore, we conclude that the gradient $\nabla f_{2}$ is not Lipschitz.

The Assumptions 2 and 3 provide greater modeling flexibility than traditional approaches based solely on Lipschitz or Hölder continuity. Specifically, they allow each function $f_{j}$ to be decomposed into a sum of two components: one with a gradient satisfying the Lipschitz condition and the other with a gradient satisfying the Hölder condition. This decomposition encompasses a broader class of problems in which the smoothness of each $f_{j}$ is non-uniform. As a result, the method remains valid and effective regardless of whether each $f_{j}$ has a purely Lipschitz gradient, a purely Hölder gradient, or any structured combination of both. It is important to highlight that for each $f_{j}$ that has a gradient satisfying only Lipschitz conditions or only Hölder conditions, Assumption 2 is also satisfied.

The next lemmas are useful in demonstrating that the sufficient descent condition at Step 4 is satisfied for sufficiently large $\sigma_{k}$, as stated in Theorem 3.4.

Lemma 3.2**.**

If $\bar{x}^{k}$ is a solution of Problem (5), then there exist $\gamma\in~\mathbb{R}_{+}^{m},\,\displaystyle\sum_{j=1}^{m}\gamma_{j}=1$ and $w_{j}^{k}\in\partial h_{j}(\bar{x}^{k}),$ such that

[TABLE]

for all $j\in\mathcal{J}$, and

[TABLE]

Proof.

Let $\bar{x}^{k}$ be the solution of Problem (5), then $\Psi_{\sigma}^{x^{k}}(\bar{x}^{k})\leq\Psi_{\sigma}^{x^{k}}(x)$ for all $x\in\mathbb{R}^{n}$, whence it follows that $\Psi_{\sigma}^{x^{k}}(\bar{x}^{k})\leq\Psi_{\sigma}^{x^{k}}(x^{k})=0$, that is,

[TABLE]

which implies

[TABLE]

This proves condition (12).

To prove condition (13), we follow the same steps taken by Calderón [teselizet] in the proof of Theorem 2.1. We can see that solving Problem (5) is equivalent to solving the following problem:

[TABLE]

If $\bar{x}^{k}$ is a solution of (14), then by the Karush–Kuhn–Tucker (KKT) conditions [martinez1], there exist $\gamma\in\mathbb{R}_{+}^{m}$ and $w_{j}^{k}\in\partial h_{j}(\bar{x}^{k})$ such that

[TABLE]

From which it follows that

[TABLE]

The condition (13) follows directly from (15). This concludes the proof of Lemma 3.2. ∎

Lemma 3.3**.**

Let $\bar{x}^{k}$ be the point computed in Step 2. Suppose that Assumptions 1, 2 and 3 hold. If

[TABLE]

then

[TABLE]

where $s^{k}\coloneqq\bar{x}^{k}-x^{k}$.

Proof.

By Assumption 2, we have

[TABLE]

Using (12), Assumption 3, and (16), we obtain

[TABLE]

By $0<\lambda_{k}\leq\dfrac{\epsilon}{\sigma_{k}\sqrt{n}}$ and the last inequality we have that

[TABLE]

Using $\sigma_{k}\|s^{k}\|\geq\epsilon$ in the last inequality we obtain

[TABLE]

This concludes the proof of Lemma 3.3. ∎

In Theorem 3.4 below we prove that there exists a $\sigma_{k}$ such that the decrement given in (6) is achieved, thus ensuring the well-definedness of Step 4.

Theorem 3.4**.**

Let $\bar{x}^{k}$ be the point computed in Step 2 of Algorithm 1. Suppose that Assumptions 1, 2, and 3 hold. If

[TABLE]

and

[TABLE]

then

[TABLE]

and

[TABLE]

Proof.

We define $s^{k}\coloneqq\bar{x}^{k}-x^{k}$. In Lemma 3.3, we obtain

[TABLE]

Thus, taking this inequality into consideration, we can conclude that to prove inequality (18) it is necessary to show that

[TABLE]

which is equivalent to proving the inequality

[TABLE]

By (17), we have that $\dfrac{\epsilon}{\sigma_{k}}\leq\|s^{k}\|$, which implies

[TABLE]

This completes the proof of (18). To obtain (19), we combine (3.3) with (18). ∎

Theorem 3.4 guarantees that, as long as the test point $\bar{x}^{k}$ does not satisfy the stopping criterion described in Step $3$, there will be a decrease in the value of the objective function, with a reduction established in Step $4$. This reduction is fundamental for obtaining the complexity results, which determine the maximum number of iterations and evaluations necessary to reach an approximate solution with the given tolerance.

The following remark provides an estimate for $\lambda_{k}$ such that the stopping criterion established in Step $3$ provides a reasonable approximation to a stationary point as defined in (3).

Remark 3.5**.**

In Lemma 3.2, we showed that if $\bar{x}^{k}$ is a solution to (5), then there exist $\gamma\in\mathbb{R}_{+}^{m}$ with $\sum_{j=1}^{m}\gamma_{j}=1$ and $w_{j}^{k}\in\partial h_{j}(\bar{x}^{k})$, for $j\in\mathcal{J}$, such that

[TABLE]

In particular, if $\bar{x}^{k}=x^{k}$, then $\|\sum_{j=1}^{m}\gamma_{j}\big[g_{f_{j}}(x^{k},\lambda_{k})+w_{j}^{k}\big]\|=0.$ Hence,

[TABLE]

If $\lambda_{k}\in\left[0,\,\displaystyle\min_{1\leq j\leq m}\left\{\dfrac{\epsilon}{\sqrt{n}L_{j}},\;\left(\dfrac{(1+\beta_{j})\epsilon}{2\sqrt{n}M_{j}}\right)^{\tfrac{1}{\beta_{j}}},\;\dfrac{\epsilon}{\sqrt{n}\sigma_{k}}\right\}\right],$ then

[TABLE]

This implies that the distance from [math] to the set $\sum_{j=1}^{m}\gamma_{j}\big(\nabla f_{j}(x^{k})+\partial h_{j}(x^{k})\big)$ is at most $\epsilon$. Thus, it is reasonable to stop the method when $\bar{x}^{k}$ is close to $x^{k}$. Therefore, the stopping criterion adopted in Step 3 is appropriate.

The Section 4 is dedicated to establishing upper bounds on the number of iterations and the number of function evaluations required to reach an approximate solution within a predefined tolerance.

4 Complexity Analysis

In Theorem 3.4 we proved that the reduction in $F_{j}$, for each $j\in\mathcal{J}$, required in Step 4 of Algorithm 1 is achieved for all sufficiently large $\sigma_{k}$. In Theorem 4.1, we then prove that the expected reduction in each $F_{j}$ is a constant factor, depending only on the constants given in the algorithm and on Assumptions 1, 2, and 3.

Theorem 4.1**.**

Suppose that Assumptions 1, 2, and 3 are satisfied, and let $x^{k+1}$ be the output of Algorithm 1. Then the minimum reduction in the objective function is given by

[TABLE]

where

[TABLE]

Proof.

It is worth remembering that if

[TABLE]

and the method does not stop at Step 3, then the decrease in the function value at each $F_{j}$ is always achieved.

Therefore, we can conclude that $\sigma_{k}\leq\sigma_{\max}$, where

[TABLE]

This implies that $\rho_{\max}=c\epsilon^{\tfrac{\beta-1}{\beta}}$. By substituting $\rho_{\max}$ in (19), we conclude that

[TABLE]

This completes the proof of (20). ∎

In the inequality (20) of Theorem 4.1, we can observe that, as $\beta$ approaches zero, the expected reduction in each $F_{j}$ decreases. Therefore, it is expected that the smaller $\beta$ is, the greater will be the number of iterations required to reach an approximate solution with a given tolerance.

Next, we demonstrate one of the main results of this section, showing that the number of iterations required to reach an approximate solution with a given tolerance is bounded by a multiple of $\epsilon^{-\tfrac{\beta+1}{\beta}}$. Theorems 4.2 and 4.3, presented below, are adaptations of Theorems 4.2 and 4.3 presented in [amaral2025complexity].

Theorem 4.2**.**

Assume that Assumptions 1, 2, and 3 hold. Given $\overline{F}\in\mathbb{R}^{m}$, the maximum number of iterations for which $F(x^{k})>\overline{F}$ and $\sigma_{k}\|\bar{x}^{k}-x^{k}\|\geq\epsilon$ is

[TABLE]

where $c$ and $\beta$ are the same as those given in Theorem 4.1.

Proof.

The case in which $\overline{F}$ does not satisfy $F(x^{0})>\overline{F}$ is trivial. It remains to analyze the case in which $F(x^{0})>\overline{F}$. In Theorem 4.1, we proved that if Algorithm 1 does not meet the established stopping criterion at iteration $k$, then the condition

[TABLE]

holds for every $i=1,2,\ldots,k$. Hence, if $F(x^{k})>\overline{F}$, we have that

[TABLE]

[TABLE]

[TABLE]

[TABLE]

This inequality implies that

[TABLE]

This concludes the proof of the theorem. ∎

It is important to note that, in Theorem 4.2, we do not consider the number of test points rejected in Step 4, that is, the number of evaluations of the function $F$. This aspect will be properly addressed in Theorem 4.3 below, which provides an upper bound for the number of evaluations of $F$ required for the method to reach an approximate solution with a given tolerance.

Theorem 4.3**.**

Suppose that the conditions of Theorem 4.2 hold. Algorithm 1 uses at most the following number of evaluations of $F_{j}$ and its subdifferential:

[TABLE]

where

[TABLE]

with $c$ and $\beta$ the same as in Theorem 4.1.

Proof.

Since $\sigma\leq\sigma_{\max}$, the parameter $\sigma$ is increased a finite number of times. Let $r$ be the maximum number of times that $\sigma$ is increased. Then, we have:

[TABLE]

This implies that:

[TABLE]

From Theorem 4.2, we know that the maximum number of iterations is upper bounded by:

[TABLE]

Therefore, the number of evaluations of $F_{j}$ and its subdifferentials employed by Algorithm 1 is upper bounded by:

[TABLE]

This completes the proof. ∎

It is relevant to compare Theorems 4.2 and 4.3 with other complexity results in the literature, in which Hölder continuity is imposed on the gradients. In this case, the complexity obtained in Theorem 4.2 is $\mathcal{O}(\epsilon^{-\frac{\beta+1}{\beta}})$ for the number of iterations. This is the same complexity as the method presented in [Lizet] for $p=1$ (although this method was also designed for constrained problems), as shown in [Lizet, Theorem 2.3]. When $m=1$, we recover the complexity established in the theory for models with quadratic regularizations and Hölder conditions imposed on the gradient of the objective function, see for example [Martinez2017, Yashtini2016]. In summary, in this work we obtain complexity results compatible with those established in both multiobjective optimization theory and scalar optimization.

5 Applications and Numerical Results

In this section, we present numerical experiments to illustrate the applicability of Algorithm 1 (PDFPM) to multiobjective optimization problems. We are particularly interested in assessing the effectiveness of the algorithm when applied to composite problems (1), where $f_{j}:\mathbb{R}^{n}\to\mathbb{R}$ exhibits Hölder regularity for some objective $j\in\mathcal{J}$, and $h_{j}:\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$ is convex but not necessarily differentiable. A direct application of composite problems is robust multiobjective optimization, which will be discussed in Subsection 5.1. For further details on the wide range of applications involving the minimization of composite problems, see [fukuda2019, Chen2019] and the references therein.

The numerical experiments, which will be presented in Subsection 5.2, were conducted using the Julia programming language (version 1.11.3) on a machine equipped with an AMD Ryzen 5 processor and 32 GB of RAM. We used standard packages for nonlinear and linear optimization. The source code is available at https://github.com/Souza-DR/AAS2025-PDFreeMO.

5.1 Application to Robust Multiobjective Optimization

Let $f_{j}:\mathbb{R}^{n}\to\mathbb{R}$ be differentiable for all $j\in\mathcal{J}$. The multiobjective optimization problem with uncertainty is defined as follows:

[TABLE]

where, for each $j\in\mathcal{J}$, $z\in\mathcal{Z}_{j}$ is a deterministic uncertain parameter, $\mathcal{Z}_{j}$ is an uncertainty set that can be estimated, and $\hat{h}_{j}:\mathbb{R}^{n}\times\mathbb{R}^{n}\to\mathbb{R}\cup\{+\infty\}$ is convex with respect to the first argument. The robust version of problem (21) is given by:

[TABLE]

If we define $h_{j}(x)\coloneqq\max_{z\in\mathcal{Z}_{j}}\hat{h}_{j}(x,z)$ for all $j\in\mathcal{J}$, then Problem (1) becomes a robust optimization problem in which the objectives are affected by uncertainties. For further details on the robust formulation of uncertain problems, see [Chen2019, fukuda2019] and the references therein.

In the proof of Lemma 3.2, we saw that the solution to the subproblem (5) is equivalent to solving the following constrained optimization problem:

[TABLE]

Note that, for $j\in\mathcal{J}$, even if $\hat{h}_{j}(\cdot,\cdot)$ is differentiable, $h_{j}(\cdot)$ is not necessarily differentiable. In general, the term $h_{j}(\cdot)$ is not easily computable, which makes solving subproblem (23) particularly challenging. However, when $\hat{h}_{j}(\cdot,\cdot)$ and $\mathcal{Z}_{j}$ have a special structure, it is possible to simplify the resolution of subproblem (23) by applying duality theory.

According to the framework presented in [fukuda2019], for each $j\in\mathcal{J}$, suppose that $\hat{h}_{j}(x,z)\coloneqq\langle x,z\rangle$ (linear with respect to the first argument) and $\mathcal{Z}_{j}=\{z\in\mathbb{R}^{n}\mid A_{j}z\leq b_{j}\}$ (a polyhedron), where $A_{j}\in\mathbb{R}^{d\times n}$ and $b_{j}\in\mathbb{R}^{d}$. In this case, we can rewrite $h_{j}(x)=\max_{z\in\mathcal{Z}_{j}}\langle x,z\rangle$ as the following linear programming problem:

[TABLE]

Its dual problem is given by

[TABLE]

By applying strong duality theory, we conclude that solving Problem (5) is equivalent to solving Problem (23), which in turn can be reformulated as the following quadratic programming problem:

[TABLE]

Thus, for specific structures of $h_{j}(\cdot)$ and $\mathcal{Z}_{j}$, the subproblem arising in Algorithm 1 becomes tractable. Moreover, this reformulation highlights how duality theory provides not only a theoretical characterization but also a practical computational advantage, since quadratic programming is a well-studied class of problems for which efficient solvers are available. This observation reinforces the applicability of our method in robust multiobjective optimization settings, closing the gap between abstract formulations and implementable algorithms.

5.2 Numerical Results

In the tests conducted with the PDFPM method, we used the following parameters: $\alpha=0.1$, $\epsilon=10^{-4}$, $\sigma_{0}=1.0$, and $B_{j}^{0}=I_{n}$ for all $j\in\mathcal{J}$. We also employed quasi-Newton matrices $B_{j}^{k}$ updated by BFGS-type schemes, one matrix per objective function. In the scalar case, it is well known (see [nocedal]) that, given a symmetric matrix $B_{j}^{k}>0$, the classical BFGS update preserves positive definiteness, i.e., $B_{j}^{k+1}>0$, provided that the curvature condition $(s^{k})^{\top}y_{j}^{k}>0$ holds, where $s^{k}\coloneqq x^{k+1}-x^{k}$ and $y_{j}^{k}\coloneqq\nabla f_{j}(x^{k+1})-\nabla f_{j}(x^{k})$ for all $j\in\mathcal{J}$. Such a condition can be enforced, for instance, if $f_{j}$ is strictly convex or if the step length satisfies Wolfe conditions. In the multiobjective setting there is no exact counterpart: if all components $f_{j}$ are strictly convex, then $(s^{k})^{\top}y_{j}^{k}>0$ holds for each $j$, but without convexity assumptions, even a multiobjective Wolfe line search [PerezPrudente2019] does not guarantee that the curvature condition is satisfied for every objective. Motivated by this difficulty, [Prudente2022] proposed a modified BFGS-type update that preserves positive definiteness when the step length is obtained by a Wolfe line search, without assuming convexity.

Since Algorithm 1 does not employ a line search, we combine this modified BFGS update with a cautious strategy in order to maintain $B_{j}^{k}$ positive definite. For each objective $j\in\mathcal{J}$, if the curvature condition $(s^{k})^{\top}y_{j}^{k}>0$ holds, we perform the classical BFGS update

[TABLE]

If $(s^{k})^{\top}y_{j}^{k}\leq 0$ but

[TABLE]

then we update $B_{j}^{k}$ by the modified BFGS rule

[TABLE]

If neither $(s^{k})^{\top}y_{j}^{k}>0$ nor $\rho_{j}^{k}>0$ holds, we simply set $B_{j}^{k+1}\coloneqq B_{j}^{k}$, i.e., no update is performed. It is worth noting that the vectors $y_{j}^{k}$ are computed from the gradients of $f_{j}$ in the Composite Problem (1), and whenever $\nabla f_{j}$ is not available in closed form these gradients are approximated by finite differences, for all $j\in\mathcal{J}$, see Remark 3.1. The maximum number of iterations allowed was set to $100$, and a run was considered unsuccessful whenever no point satisfying the stopping criterion was obtained within this limit; see Step 3 of Algorithm 1.

Considering the discussion in Subsection 5.1, we now apply the PDFPM method to a robust optimization problem. For each test problem and each objective $j$, we generate a nonsingular matrix $\tilde{A}_{j}\in\mathbb{R}^{n\times n}$ with entries sampled uniformly in $[0,1]$ and keep it fixed across all uncertainty levels. Then, for a given $\delta\geq 0$, we define

[TABLE]

where $e=(1,\ldots,1)^{\top}\in\mathbb{R}^{n}$. By taking $A_{j}\coloneqq[\tilde{A}_{j};-\tilde{A}_{j}]\in\mathbb{R}^{2n\times n}$ and $b_{j}\coloneqq(\delta,\ldots,\delta)^{\top}\in\mathbb{R}^{2n}$, we obtain a robust optimization problem as defined in (22), where $\delta$ controls the level of uncertainty in the problem. In particular, $\delta=0$ corresponds to the nominal (non-robust) case without uncertainty.

The presence of the pointwise maximization in (25) makes the function $h_{j}$, for each $j\in\mathcal{J}$, convex but not necessarily differentiable. Consequently, at each iteration of the algorithm we solve, for every objective $j\in\mathcal{J}$, one linear program of the form (25) in order to evaluate $h_{j}(x)$, together with a single quadratic program of the form (24). To address these computational challenges, we employ HiGHS [HiGHS] to solve the linear programs (25) and Ipopt [Ipopt] to solve the quadratic subproblem (24), an interior-point optimizer that efficiently handles nonlinear constraints through a primal–dual method with a filter line-search approach.

5.2.1 Hölder Setting

To illustrate the applicability of the proposed algorithm, we define two biobjective functions with the following properties: the first vector function has one objective with a Lipschitz continuous gradient, while the other has a $(p-1)$-Hölder continuous gradient; the second vector function has both objectives with $(p-1)$-Hölder continuous gradients. See Examples 5.1 and 5.2 below. Note that, for $1<p<2$, the mapping $x\mapsto\|x\|_{p}^{p}$ has a gradient that is $(p-1)$-Hölder continuous. Consequently, in Example 5.1 the gradient of $f_{2}$ is $(p-1)$-Hölder continuous, whereas in Example 5.2 the gradients of $f_{1}$ and $f_{2}$ are $(p_{1}-1)$- and $(p_{2}-1)$-Hölder continuous, respectively. Hence both AAS1 and AAS2 satisfy the Hölder-type assumptions adopted in this work.

Example 5.1** (AAS1).**

Let $F:\mathbb{R}^{2}\to\mathbb{R}^{2}$ be defined by $F(x)=(f_{1}(x),f_{2}(x))$ with

[TABLE]

where $1<p<2$. The first objective exhibits a quadratic structure with strong convexity, while the second objective represents a $p$-norm regularization term with Hölder continuous gradients. The problem parameters are specified as follows:

[TABLE]

Example 5.2** (AAS2).**

Consider the following biobjective function:

[TABLE]

where $1<p_{1},p_{2}<2$. Both objective functions are convex and possess Hölder continuous gradients. The problem parameters are specified as follows:

[TABLE]

[TABLE]

In the experiments, we considered different uncertainty scenarios in the minimization of the biobjective functions AAS1 and AAS2, namely $\delta\in\{0.0,~0.02,~0.05,~0.1\}$. For each fixed value of $\delta$, we generated $m$ matrices $A_{j}\in\mathbb{R}^{2n\times n}$, with $j\in\mathcal{J}$, thereby defining our robust Problem (22). For each problem instance, the algorithm was executed $200$ times with randomly selected starting points $x^{0}\in\mathcal{X}$, where $\mathcal{X}\coloneqq[-2,2]\times[-2,2]$ for AAS1 and $\mathcal{X}\coloneqq[-5,5]\times[-5,5]$ for AAS2.

Figures 1 and 2 correspond to the case $\delta=0$, that is, a scenario without uncertainty. The right panel illustrates the convergence behavior of the algorithm, highlighting its ability to approximate the true Pareto front. We also plot the initial point (in gray) and the final point (in red) of the optimization, connected by a dashed line. In addition, the left panel in both figures shows the complete set of objective function values obtained by evaluating the objectives on a uniform grid of points, thus providing a comprehensive view of the objective space landscape. Figure 3 presents the comparison of the reconstructions obtained for the other values of $\delta$. On the left, the approximations of the Pareto fronts obtained for AAS1 are shown, and on the right, those for AAS2. One observes that as the level of uncertainty increases, the approximation of the Pareto front produced by PDFPM deviates further from the expected front.

To quantify the reliability of the PDFPM solver in the Hölder setting, we recorded, for each problem and each value of $\delta$, how many of the $200$ independent runs converged to a point satisfying the stopping criterion. Table 1 reports the results for the AAS1 and AAS2 problems. For AAS1, between $195$ and $199$ runs were successful, which corresponds to success rates in the interval $97.5\%\text{ to }99.5\%$ across all uncertainty levels. For AAS2, the performance is even better: between $199$ and $200$ runs converged successfully for every tested value of $\delta$, yielding success rates between $99.5\%\text{ and }100\%$. Although problems with Hölder continuous gradients are typically harder to solve than their Lipschitz counterparts, these high success rates, with averages of approximately $98.6\%$ for AAS1 and $99.9\%$ for AAS2, suggest that PDFPM handles different uncertainty levels in a stable way and performs reliably on composite multiobjective problems with Hölder continuous objectives.

5.2.2 Comparison with ProxGrad and CondG

To assess the performance of Algorithm 1 on standard multiobjective test problems, we compare PDFPM with two first-order methods designed for composite multiobjective optimization. The first one is the proximal gradient method proposed in [fukuda2019, Algorithm 3.1], which we refer to as ProxGrad. The second is the conditional gradient method introduced in [Assuncao2023, Algorithm 1], denoted by CondG. Although ProxGrad, CondG, and Algorithm 1 are all tailored to composite multiobjective problems, there are important differences in the underlying subproblems, in the corresponding search directions, and in the stopping criteria.

As discussed in Section 5.1, Step 2 of Algorithm 1 requires the solution of the subproblem (5), which is equivalent to (23). In an analogous way, ProxGrad requires the solution of the following subproblem

[TABLE]

where $\ell>0$ is an algorithmic parameter. In our implementation the solution of (26), denoted by $\bar{x}^{k}_{pg}$, is computed with Ipopt [Ipopt]. The corresponding search direction is then given by $d(x^{k})=\bar{x}^{k}_{pg}-x^{k},$ and we denote the optimal value of (26) at $\bar{x}^{k}_{pg}$ by $\theta_{pg}(x^{k})$.

For CondG the subproblem is

[TABLE]

and the search direction is defined as $d(x^{k})=\bar{x}^{k}_{cg}-x^{k},$ where $\bar{x}^{k}_{cg}$ is a solution of (27). The optimal value of (27) at $\bar{x}^{k}_{cg}$ is denoted by $\theta_{cg}(x^{k})$. In our experiments, (27) is solved with HiGHS [HiGHS]. For both ProxGrad and CondG, a line search is performed along the direction $d(x^{k})$. If $x^{k}$ is noncritical and $d(x^{k})$ is a descent direction, we look for $\eta_{k}\in(0,1]$ such that

[TABLE]

where $\zeta=10^{-4}$ and $\theta(x^{k})$ denotes the optimal value of the corresponding subproblem (that is, $\theta_{pg}$ for ProxGrad and $\theta_{cg}$ for CondG). After the line search, the new iterate is given by $x^{k+1}=x^{k}+\eta_{k}\,d(x^{k})$.

In our tests we use, for each method, the stopping criterion proposed in the corresponding reference. For ProxGrad, the algorithm is designed to stop when $\bar{x}^{k}_{pg}=x^{k}$, which is equivalent to $\theta_{pg}(x^{k})=0$. For CondG, the authors require $\theta_{cg}(x^{k})=0$. Accordingly, in the numerical comparison we declare success for ProxGrad and CondG whenever the associated optimal value satisfies $|\theta(x^{k})|<\varepsilon$. For PDFPM, see Step 3 of Algorithm 1, the stopping criterion is $\sigma_{k}\|\bar{x}^{k}-x^{k}\|<\varepsilon$, where $\bar{x}^{k}$ is the solution of the subproblem (5). We set $\varepsilon=10^{-4}$ for all methods and use $\ell=1$ in (26). The details of the parameter choices for PDFPM were given at the beginning of this section; see Subsection 5.2. In all experiments reported in this subsection, the maximum number of iterations for each algorithm is set to $100$, and a run is declared unsuccessful whenever no point satisfying the corresponding stopping criterion is obtained within this limit.

The test problems are listed in Table 2 and are frequently used to validate multiobjective optimization methods, see for instance [Assuncao2021, Prudente2022] and the references therein. As in Subsection 5.2.1, we consider different uncertainty levels in the minimization of the functions in Table 2, namely $\delta\in\{0.0,0.02,0.05,0.1\}.$ For each fixed value of $\delta$, we generate $m$ matrices $A_{j}\in\mathbb{R}^{2n\times n}$, with $j\in\mathcal{J}$, thus defining the robust problem (22), as discussed previously. For each problem instance, the algorithm under consideration is executed $200$ times with randomly selected starting points. The same collection of $200$ starting points is used for PDFPM, ProxGrad, and CondG to ensure a fair comparison.

The numerical results are presented by means of performance profile plots [Dolan2002], which are a standard tool for comparing several methods on a large collection of test problems. To compare the behavior of each algorithm, we consider all $56{,}000$ instances arising from $70$ problems, $4$ uncertainty levels, and $200$ starting points. Figures 4 and 5 display performance profiles for PDFPM, ProxGrad, and CondG with respect to the number of function evaluations, the number of gradient evaluations, the iteration count, and the total execution time. The plots show that PDFPM is more efficient and more robust than ProxGrad and CondG when function and gradient evaluations are used as performance measures. When iteration count and execution time are considered, ProxGrad is typically the most efficient method, but PDFPM is the most robust algorithm overall, in the sense that it attains the highest fraction of successfully solved instances across all performance measures.

6 Conclusions

In this work, we addressed a multiobjective optimization problem involving the composition of a differentiable function with a convex one, under the assumption that the gradient of each differentiable component can be decomposed into parts with Lipschitz and Hölder regularity. This flexible modeling covers both classical cases and mixed scenarios with different regularities across objectives. One of the main contributions of the method is its flexibility in handling approximations of the differentiable part, allowing exact derivatives, first-order approximations (such as finite differences), or second-order information through a matrix $B_{j}^{k}$, which may be either the exact Hessian or an approximation obtained via methods such as BFGS, secant Newton, or finite differences, provided that $\{B_{k}\}$ is a bounded sequence of positive semidefinite matrices. This makes the method applicable to a wide range of problems, including situations where derivative information is only partially or not directly available. Moreover, large values of the regularization parameter $\sigma$ may lead to small steps and slow convergence, so it is advisable to start with a small $\sigma$. The method does not require prior knowledge of the Lipschitz–Hölder constants, and its complexity analysis ensures convergence to an $\epsilon$-approximate solution in a finite number of iterations, with order $\mathcal{O}(\epsilon^{-\frac{\beta+1}{\beta}})$, consistent with the best-known complexity results for this type of regularity.

The numerical experiments focused on robust multiobjective optimization problems in the Hölder setting. In the biobjective examples AAS1 and AAS2, we observed that, as the level of uncertainty increases, the Pareto fronts approximated by PDFPM move away from the nominal front corresponding to the case without uncertainty, as expected from the robust formulation. At the same time, PDFPM maintained very high success rates across all uncertainty levels, including the more challenging instances with Hölder continuous gradients, which suggests that the method handles uncertainty in a stable way on these problems. In a larger benchmark study involving standard multiobjective test problems and several uncertainty levels, a performance-profile comparison with two first-order composite methods indicated that PDFPM is particularly robust in terms of the fraction of successfully solved instances, while remaining competitive with respect to the usual performance measures.

We developed a partially derivative-free method, and for future work, it would be interesting to consider a fully derivative-free algorithm. Additionally, other situations could be explored, such as vector environments and closed pointed convex cones.

\bmhead

Acknowledgements This work has been partially supported by CNPq grant 407147/2023-3 (V. S. Amaral).

\bmhead

Data availability We do not analyze or generate any datasets, as our work proceeds within a theoretical approach.

Declarations

The authors have no conflicts of interest to declare that are relevant to the content of this article.