% TEX 7(ascii) bits
\documentclass{article}
\usepackage[fleqn]{amsmath}
\usepackage{amssymb}
\usepackage{epsfig}
\usepackage{epsf}
\usepackage{hero}
%\usepackage{siamptex}


%% change sizes and margins
\topmargin      -16.0mm
\oddsidemargin  -11.0mm
\evensidemargin -11.0mm
\textheight     243.5mm
\textwidth      190.0mm
\columnsep        6.0mm
\parindent        1.0em
\headsep          6.3mm
\headheight        12pt
\lineskip           1pt
\normallineskip     1pt
\def\baselinestretch{1}

%\newtheorem{proof}{Proof}
\newtheorem{definition}{Definition}
\newtheorem{corollary}{Corollary}
\newtheorem{theorem}{Theorem}
\newtheorem{proposition}{Proposition}
\newtheorem{lemma}{Lemma}
\newtheorem{algorithm}{Algorithm}
\newtheorem{procedure}{Procedure}
\newtheorem{example}{Example}
\newtheorem{assumptions}{Assumptions}
\newtheorem{assumption}{Assumption}
\title{Generalized Proximal Point Algorithms and Bundle Implementations}
%\author{\thanks{Dpt. of Electrical Engineering and 
%Computer Science, The University of Michigan, 1301 Beal Avenue, Ann Arbor, 
%Michigan 48109-2122, USA. (chretien@eecs.umich.edu and hero@eecs.umich.edu).
%This work has been partially supported by AFOSR F49620-96-0028.} 
%\and Alfred O. Hero$^*$}

\author{St\'ephane Chr\'etien and Alfred O.  Hero\\
\mbox{}
\\
COMMUNICATIONS \& SIGNAL PROCESSING LABORATORY
\\
Department of Electrical Engineering and Computer Science
\\
The University of Michigan 
\\
Ann Arbor, MI 48109-2122
\\
\mbox{}
\\
{\bf Technical Report No. 316}, Aug. 1998}



\begin{document}

\maketitle

\vskip 1in

\bc
{\bf ABSTRACT}
\ec

  In this paper, we present a study of the {\em proximal point} algorithm
  using very general regularizations for minimizing possibly
  nondifferentiable and nonconvex locally Lipschitz functions. We deduce from
  the {\em proximal point} scheme simple and implementable bundle methods for
  the convex and nonconvex cases. The originality of our bundle method is
  that the bundle information incorporates the subgradients of both the
  objective and the regularization function. The resulting method opens up a
  broad class of regularizations which are not restricted to quadratic,
  convex or even differentiable functions.


\bt
{\bf Keywords}:
mathematical programming, proximal point, bundle methods, nonsmooth regularization
\et


This work was partially supported by AFOSR F49620-96-0028.

\newpage

\newpage

%\begin{AMS}
%26B05, 90C26, 90C30, 90C90, 65U05
%\end{AMS}

\pagestyle{myheadings}
\markboth{CSPL-315  Aug. 1998}{Chretien and Hero, ``Generalized proximal
  methods \ldots''}


\section{Introduction}
In this paper, we address the problem of minimizing a locally Lipschitz
possibly nondifferentiable and nonconvex function $f(x)$ on $\mathbb R^n$, i.e.
\begin{equation}\label{minimini}
\min_{x \in \mathbb R^n} f(x).
\end{equation}
One of the most widely studied methods for solving nondifferentiable optimization 
problems is the bundle method first proposed by Lemarechal \cite{Lemarechal:1975} and
Wolfe \cite{Wolfe:1975} for convex minimization and further developed by Mifflin 
\cite{Mifflin:1977,
Mifflin:1982} and Kiwiel \cite{Kiwiel:1985,Kiwiel:1986,Kiwiel:1996} for the nonconvex case; 
see also \cite{Bihain:1984}, \cite{Lemarechal:1981}, \cite{Schramm:1992} and the references 
therein. The bundle method can be interpreted as a {\em cutting plane} algorithm stabilized by a
quadratic penalty or regularization. In its simplest form, for $f$ convex, 
the bundle method generates a sequence of iterates, starting from $x^1$ and defined by 
\begin{equation}\label{bund}
x^{k+1}={\rm argmin}_{x\in \mathbb R^n}\{\hat{f}(x)+\frac{\lambda}{2} |x^k-x|^2\}.
\end{equation}
where 
\begin{equation}
\hat{f}(x)=\max_{j \in J^k} \{f(y_j)+\langle g(y_j),x-y_j\rangle \}
\nonumber\end{equation}
is a piecewise linear approximation called the 
{\em cutting plane} model, $y_j$, $j\in J^k$ are some points in a neighborhood of the 
current iterate $x^k$ and $g(y)$ is a subgradient of $f$ at the point $y$. In the case 
where $f$ is nonconvex, the following polyhedral 
approximation is usually chosen, as in \cite{Kiwiel:1985,Kiwiel:1986,Kiwiel:1996},
\begin{equation}\label{polyhf}
\hat{f}(x)=f(x^k)+
\max_{j\in J^k}\{-\alpha_j^k+\langle g(y_j),x-x^k\rangle,\;\;\;j\in J^k \},
\end{equation}    
where
\begin{align*}
& \alpha^k_j =\alpha(x^k,y^j) \\
& \alpha(x,y) =|f(x)-f(y)-\langle g(y),x^k-y\rangle|.
\end{align*}
One fruitful interpretation of the bundle method is to consider iteration (\ref{bund})
as an implementable approximation of the well known {\em proximal point} algorithm using a
{\em cutting plane} model $\hat{f}$ of the objective function. In the original form 
\cite{Martinet:1970,Rockafellar:1976}, the {\em proximal point} algorithm is 
defined by the reccurence
\begin{equation}
x^{k+1}={\rm argmin}_{y\in \mathbb R^n}\{f(y)+\frac{\lambda_k}{2} |x-y|^2\}.
\nonumber\end{equation}
where $(\lambda_k)_{k \in \mathbb N}$ is a sequence of positive relaxation parameters.
The proximal point algorithm and the bundle method share the same property of solving
a sequence of minimization subproblems incorporating a quadratic penalty, also denoted 
Moreau-Yosida regularization \cite{Lemarechal:1997}. 

In this paper, we address the study of the proximal point algorithm and its
bundle implementations using very general nonquadratic penalty functions. 
In particular, we establish convergence for a class of locally Lipschitz 
regularizations without any convexity nor differentiability assumptions. The utility
of such nonquadratic regularization is motivated by the following examples.

\begin{example}[EM-type algorithms for maximum likelihood estimation]
{\rm
We show here that the case of Kullback regularization results in a 
proximal point method which is a generalization of the well known Expectation 
Maximization (EM) algorithm for maximum likelihood estimation \cite{Dempster:1977}. 
Consider as in \cite{Dempster:1977} the sample spaces 
$\Omega_1$ and $\Omega_2$ on which one defined the random variables $V_1$ and $V_2$ with respective 
probability densities $p_1(v_1;x)$ 
and $p_2(v_2;x)$, both indexed by an unknown parameter $x\in \mathbb R^n$ to be estimated. 
Assume that $V_2$ is obtained from $V_1$ 
through a many-to-one mapping $V_1\rightarrow V_2=h(V_1)$ and define 
$p(v_1\mid v_2;x)=\frac{p_1(v_1;x)}{p_2(v_2;x)}$ the  
density of $V_1$ conditioned on $V_2=v_2$. Then, the Kullback information measure between
$p(v_1\mid v_2;x)$ and $p(v_1 \mid v_2;y)$ for two parameter values $x$ and $y$ 
takes the following form
\begin{equation}\label{Kull}
I(x,y \mid v_2)=\int_{\{h^{-1}(v_2)\}}
\log \Big(\frac{p(v_1\mid v_2;y)}{p(v_1\mid v_2 ;x)}\Big) \;\;p(v_1\mid v_2;y)dv_1
\end{equation}
Now, consider the following proximal point
algorithm for maximizing the {\em likelihood} function $f(x)=\log p_2(v_2;x)$. 
\begin{equation}\label{truc}
x^{k+1}={\rm argmin}_{y\in \mathbb R^n}\{-\log p_2(v_2;y)+\lambda_k I(x^k,y\mid v_2)\}.
\end{equation}
Using the facts that 
\begin{equation}
p(v_1\mid v_2;y)=\frac{p(v_1;y)}{p(v_2;y)}
\nonumber
\end{equation}
and
\begin{equation}
\int p(v_1\mid v_2 ; y)dv_1=1
\nonumber
\end{equation}
it can be shown that (\ref{truc}) takes the equivalent form
\begin{align}
\label{ulmit}
x^{k+1}= & {\rm argmax}_{y \in \mathbb R^n} \Big\{ -(1-\lambda_k) \log p_2(v_2;y)\\
& \hspace{1cm} +\lambda_k {\sf E}[\log p_1(v_1;y) \mid V_2=v_2;x^k]\Big\}\nonumber
\end{align}
where ${\sf E}[\log p_1(v_1;y)\mid V_2=v_2;x^k]=\int \log p_1(v_1;y) p(v_1\mid v_2;x^k)dv_1$
denotes conditional expectation. When $\lambda_k=1$, recursion (\ref{ulmit}) is identical to
the EM algorithm introduced in \cite{Dempster:1977} where $V_2$ is
the incomplete data and $V_1$ is the complete data. 
Furthermore, as implied by Lemma 4.1 below the recursion (\ref{ulmit}) monotonically increases the
log-likelihood $\log p_2(v_2;y)$ as does the standard EM algorithm of \cite{Dempster:1977}. 
A special case of (\ref{truc}) is
the case of Laplacian data,
\begin{equation}
p_2(v_2;x)=\frac{x}{2}\exp(-x|v_2|),\:\: x\geq 0.
\nonumber
\end{equation}
When the complete data $V_1$ is also chosen as Laplacian, it is easy to show that the
Kullback regularization given by (\ref{Kull}) is nonsmooth and nonconvex. 
}
\end{example}

\begin{example}[Methods of multipliers]
{\rm
In \cite{Rockafellar:1976b} Rockafellar shows that the proximal point approach 
can be applied to the dual of a constrained optimization problem to yield interesting
classes of multiplier methods. Subsequent studies \cite{Teboulle:1992,Teboulle:1997,
Iusem:1994} have demonstrated the benefit of using nonquadratic regularization functions.
Among the possible choices for regularization functions proposed in \cite{Teboulle:1992} is 
the $\phi-$divergence
\begin{equation}
d_{\phi}(x,x^{\prime})=\sum_{j=1}^nx_j 
\phi(\frac{x^{\prime}_j}{x_j}).
\nonumber\end{equation} 
where, in \cite{Teboulle:1992}, $\phi$ was assumed strictly convex. 
In particular, consider the convex program
\begin{equation}\label{const}
\min_{x \in \mathbb R^n} f(x)
\hspace{.5cm}\text{ subject to }\hspace{.5cm}
g_i(x)\leq 0, \;\;\; i\in \{1,\ldots,m\}
\end{equation}
where $f$ and $g_1,\ldots,g_m$ are convex functions.
The proximal point algorithm applied to the dual takes the form (see \cite[section 6]{Iusem:1994})
\begin{align*}
& p^0 \in \mathbb R^m_+ \\
& p^{k+1}={\rm argmax}_{p\geq 0} \{c(p)-\lambda_k d_{\phi}(p^k,p)\}
\end{align*}
with $c(p)$ being the dual functional defined by $\inf_{x\in \mathbb R^n} L(x,p)$ where
$L(x,p)$ is the Lagrangian 
\begin{equation*}
L(x,p)=\begin{cases}
\begin{array}{ll}
f(x)+\sum_{i=1}^m p_i g_i(x) & \text{if } p_i \geq 0,\:\: \forall i \in \{1,\ldots,m\}\\
-\infty & \text{otherwise}
\end{array}
\end{cases}
\end{equation*}
Thus, one obtains convergence of $(x^k)_{k\in \mathbb N}$ and $(p^k)_{k \in \mathbb N}$ to the
solution of problem (\ref{const}) (for instance, see \cite{Iusem:1994}).
Using different choices for the function $\phi$, some well known multiplier methods
can be recovered. Our generalization of $\phi$ to nonconvex functions opens up 
many new possibilities. 
}
\end{example}

Further examples of nonsmooth and nonconvex regularizations have also recently been
studied in the context of inverse problems in \cite{Nikolova:1997} and 
\cite{Nikolova:1997b}.

The outline of the paper is the following. In Section \ref{prox} the generalized
proximal point algorithm is introduced for a wide class of possible
regularizations. In Section \ref{NDOpt}, the fixed points of the method are studied. In 
particular an analysis of nondifferentiable nonconvex regularization is provided 
which seems to have no precedent in the literature. In Section \ref{Conv}, global convergence
of the method is established. We then demonstrate
local convergence when $f$ is strictly convex in
an open neighborhood of an accumulation point of the minimizing sequence.
In sections \ref{bundle1} and \ref{bundle2}, we turn to implementable bundle methods for 
approximation of the generalized proximal point iterations.
The case of convex function is discussed first in Section \ref{bundle1} for the sake of clarity in the
exposition. Then, algorithmic refinements, including a linesearch, are introduced in section \ref{bundle2}
in order to accomodate the nonconvex case. 
 
We recall the following notations and definitions \cite{Clarke:1983}. 
The inner product on $\mathbb R^n$ is denoted by $\langle.,.\rangle$ and the associated
norm is $|.|$. The convex hull of a set $S$ is denoted ${\rm conv}(S)$. The function
$f$: $\mathbb R^n\rightarrow R$ is locally Lipschitz if for any bounded set 
$B\subset \mathbb R^n$
there exists a constant $L<\infty$ such that $|f(x)-f(y)|\leq L|x-y|$ 
$\forall x,y \in B$. The effective domain of $f$ is
$C_{f}=\{x\;\;\mid\;\; f(x)<\infty\}$.
The generalized derivative of $f$ at $x$ in the direction of $d$ is
\begin{equation}
f^{\prime}(x,d)=
\underset{\underset{t\downarrow 0}{h\rightarrow 0}}{\limsup} \;\frac{f(x+h+td)-f(x)}{t}
\nonumber\end{equation}
The subdifferential of $f$ at $x$ is the set
\begin{equation}
\partial f(x)={\rm conv}\{\lim_{x^i\rightarrow x} \nabla f(x^i)
\;\;\mid\;\;\nabla f(x^i) \text{ exists}\},
\nonumber\end{equation}
Where $\nabla f$ denotes the gradient of $f$.
The subdifferential has the property that it is a closed and convex set.
An equivalent definition of $\partial f(x)$ which will be useful is
\begin{equation}\label{secdef}
\partial f(x)=\{g \in \mathbb R^n \;\;\mid\;\; \langle g,d\rangle\leq f^{\prime}(x,d)\;\;
\forall d \in \mathbb R^n\}.
\end{equation}
The multivalued function $x \mapsto \partial f(x)$ is upper semicontinuous and locally
bounded. Notice that if a point $x$ is a local minimum of the function $f$ 
we have $f^{\prime}(x,d)\geq 0$ for all $d$. In this case, following the second 
definition (\ref{secdef}) of the subdifferential, $0\in \partial f(x)$. 
More generally, we say that $x$ is a stationary point for $f$ if $0\in \partial f(x)$. 
For convex functions $f(x)$ the subdifferential is equivalently defined by
\begin{equation}\label{troizedef}
\partial f(x)=\{g \in \mathbb R^n \;\;\mid\;\; \langle g,y-x\rangle\leq f(y)-f(x)\;\;
\forall y \in \mathbb R^n\}.
\end{equation}
An extention of the subdifferential, called the $\epsilon$-subdifferential,
can be defined for convex function as follows
\begin{equation}\label{quarte}
\partial_\epsilon f(x)=\{g \in \mathbb R^n \;\;\mid\;\; 
f(y)\geq f(x)+\langle g,y-x\rangle-\epsilon\;\;
\forall y \in \mathbb R^n\}.
\end{equation}
In the nonconvex case, we will need to introduce another type of approximate subdifferential
called the Goldstein $\epsilon$-subdifferential.
For any $x\in \mathbb R^n$ and any $\epsilon \geq 0$ the Goldstein 
$\epsilon$-subdifferential \cite{Kiwiel:1985} of $f$ at $x$ is the set 
\begin{equation}
\partial^{\epsilon} f(x)={\rm conv}\{\partial f(y) \;\;\mid\;\; |y-x|\leq \epsilon \}.
\nonumber\end{equation}
The multivalued function $(x,\epsilon)\mapsto \partial^{\epsilon} f(x)$ is locally bounded and
upper semicontinuous, i.e., $x^k\rightarrow x$, $\epsilon^k\rightarrow \epsilon$, 
$p^k\in \partial_{\epsilon^k} f(x^k)$
and $p^k\rightarrow p$ imply $p \in \partial_\epsilon f(x)$.
We will also need the notion of {\em weak upper semismoothness} introduced by Mifflin
\cite{Mifflin:1977}. A function $f$ : $\mathbb R^n \rightarrow \mathbb R$ is said to
be weakly upper semismooth at $x \in \mathbb R^n$ if

a. f is Lipschitz on a ball about $x$ and

b. for each $d\in \mathbb R^n$ and for any sequences $\{t^k\}\subset \mathbb R_+$ and
$\{g^k\}\subset \mathbb R^n$ such that $\{t^k\}\downarrow 0$ and $g^k \in
\partial f(x+t^kd)$ it follows that 
\begin{equation}
\liminf_{k\rightarrow\infty}\langle g^k,d\rangle\geq \limsup_{t\downarrow 0}
\frac{f(x+td)-f(x)}{t}.
\nonumber\end{equation}
Finally, for any function $F(x,x^{\prime})$ satisfying the local 
Lipschitz property in $x^{\prime}$, we will use the notation 
$\partial F(x,x^{\prime})$ for the subdifferential in the
second variable at $x^{\prime}$. 

\section{Generalized Proximal Point Algorithms}\label{prox} 
In what follows, the regularization function 
is denoted $\Psi(.,.)$. We first state some assumptions on the
objective function $f$ and the regularization function $\Psi$.

\begin{assumptions}[Objective function]
\label{ass1}
(i) {\rm $f$ is inf-compact, i.e. the $\alpha-$level sets 
$\mathcal L_f(\alpha)=\{x \in \mathbb R^n \;\;\mid\;\; f(x)\leq \alpha\}$ are bounded for any 
$\alpha \in \mathbb R$}  \\
(ii) {\rm $f$ is locally lipschitz over $\mathbb R^n$}.
\end{assumptions}

Notice that Assumptions \ref{ass1} (i) and (ii) together imply that $f$ is
bounded from below. 
\begin{assumptions}[Regularization function]
\label{ass2}
(i) {\rm (Positivity)} $\Psi(x,y)\geq 0$ {\rm for all} $x$, $y\in \mathbb R^n$. \\
(ii) {\rm (Identifiability)} $\Psi(x,y)=0 \Leftrightarrow x=y$. \\ 
(iii) $\Psi(x,y)$ {\rm is locally lipschitz over} $\mathbb R^n\times \mathbb R^n$.\\
(iv) {\rm The effective domain} $C_{\Phi(x,.)}=\mathbb R^n$ {\rm for any} $x\in \mathbb R^n$.
\end{assumptions}
A usefull property of the regularizing function $\Psi$ is now given
\begin{lemma}\label{heub}
Assume that $\Psi(x,y)$ satisfies Assumptions 2. Then,
\begin{equation}
\partial \Psi(x,x)=\{0\},
\nonumber\end{equation}
for all $x$ in $\mathbb R^n$. 
\end{lemma}
\noindent{\it Proof}:
Fix $x$ in $\mathbb R^n$. Due to Assumptions 2, $\Psi(x,y)>\Psi(x,x)$ for all
$y$ in $\mathbb R^n$. Therefore, the first order optimality condition gives
\begin{equation}
0 \in \partial \Psi(x,x),
\nonumber\end{equation}
as desired. 
%\end{proof}

We now define the {\em generalized proximal point} algorithm 

\begin{definition}
Assume that $f$ and $\Psi(.,.)$ satisfy assumptions 1 and 2 respectively. 
Then, the {\em generalized proximal point} algorithm is defined by 
the following recursion starting at $x^1$ 
\begin{equation}\label{genprox}
x^{k+1}\in {\rm argmin}_{y\in \mathbb R^n}\{f(y)+\lambda_k \Psi(x^k,y)\}.
\end{equation}
\end{definition}

The following result proves that recursion (\ref{genprox})
definition is well defined. 
\begin{proposition}
The set of minimizers of (\ref{genprox}) is nonempty for any $x^k$ in $\mathbb R^n$. 
\end{proposition}
\noindent{\it Proof}:
This result is a straightforward consequence of Assumptions 1 and 2. Indeed,
the function $f(x)+\lambda_k \Psi(x^{\prime},x)$ inherits the inf-compactness
and local Liptschitz properties of $f$ and $\Psi$ and thus possesses at least one bounded minimizer. 
%\end{proof}

Notice that, due to nonconvexity of the functions involved, the minimum
in (\ref{genprox}) may not be unique. In such cases, 
$x^{k+1}$ in (\ref{genprox}) can be arbitrarily chosen among the set of minimizers of (\ref{genprox}).

For the sake of notational convenience, in the remainder of this paper, the regularized 
objective function with relaxation parameter $t$ 
will be denoted 
\begin{equation}\label{regfunc}
F_t(x,y)=f(y)+t\Psi(x,y)
\end{equation}
 and the {\em generalized proximal} operator will be 
defined by 
\begin{equation}
P_t(x)={\rm argmin}_{y \in \mathbb R^n} F_t(x,y).
\nonumber\end{equation}
\section{Optimality for nonsmooth regularization}\label{NDOpt} 
We first investigate theoretical difficulties concerning
optimality conditions for nonsmooth regularization. Indeed the problem which
one encounters in nonsmooth situations is that a stationary
point of the regularized objective function is not in general a stationary point
of the unregularized objective function itself. Thus, it is important to establish
sufficient conditions which guarantee that stationary points of the 
objective function coincide with those of the regularized problem. For that
purpose, we introduce the condition of {\em subdifferential domination}.
Under this condition, we show that the fixed
points of the generalized proximal point mapping are
locally optimal. 

\subsection{Fixed points and nonsmooth regularization}
Consider the regularized function (\ref{regfunc}) where, with no loss of generality, 
$\lambda$ is set to $1$, i.e.
\begin{equation}
F(x,y)\triangleq F_1(x,y)=f(y)+\Psi(x,y),
\nonumber\end{equation}
and the associated proximal operator
\begin{equation}
P(x)={\rm argmin}_{y\in \mathbb R^n} F(x,y).
\nonumber\end{equation}
The success of the proximal point approach relies on the condition that fixed points
of $P$ also be stationary points of $f$.   
In the case where $\Psi$ is quadratic, say $\Psi(x,y)=\frac12 |x-y|^2$, this property
is well known and straightforward. In the case 
of the general possibly nonconvex regularization considered in this paper, we only have
\begin{equation}
x^*={\rm argmin}_{y\in \mathbb R^n} F(x^*,y),
\nonumber\end{equation}
which is equivalent to 
\begin{equation}\label{etalors}
0 \in \partial F(x^*,x^*).
\end{equation}
Using the calculus of subdifferentials, we have
\begin{equation}
\partial F(x,y)\subset \partial f(y)+\partial \Psi(x,y),
\nonumber\end{equation}
and thus, (\ref{etalors}) implies
\begin{equation}\label{optim}
0 \in \partial F(x^*,x^*)\subset \partial f(x^*)+\partial \Psi(x^*,x^*)
\end{equation}
Hence, consistence of the generalized proximal point approach reduces to the question
of knowing in which circumstances a point $x^*$ satisfying (\ref{optim}) is a 
stationary point of $f$. 
 
In the case where $\Psi(x^*,y)$ is smooth in the second variable, 
e.g. $\Psi(x^*,y)=|x^*-y|^2$, then $\partial \Psi(x^*,x^*)$ reduces  to $\{0\}$,
and we deduce from (\ref{optim}) that $0\in \partial f(x^*)$ which
proves that $x^*$ is a stationary point of the objective function. 
For $\Psi(x^*,y)$ nonsmooth in the
second variable, $x^*$ is no longer guaranteed to be a stationary
point of $f$. To illustrate, consider the following one dimensional example,
\begin{align*}
& f(y) =
\begin{cases}
x^*-y, \text{ if } y \leq x^*\\ 
\frac12(x^*-y), \text{ if } y > x^*
\end{cases}\\
& \Psi(x^*,y) =|x^*-y|.
\end{align*}
We then have $\partial f(x^*)=[-1,-\frac12]$ and $\partial
\Psi(x^*,x^*)=[-1,1]$. On the other hand, $\partial F(x^*,x^*)=[-2,\frac12]$.
Thus, in this case, we have $0\in \partial F(x^*,x^*)$ whereas
$x^*$ is not a stationary point of $f$ (indeed, $f$ has no stationary
point).

\subsection{Subdifferential domination}
The reason that condition (\ref{optim}) fails in the latter example is that the
subdifferential of the regularization $\Psi$ is ``bigger'' than the
subdifferential of the objective $f$. To overcome this difficulty,
we will impose that the objective
have ``greater'' variation than the regularization. The following
definition makes precise the idea of ``greater'' variation.

\begin{definition}Let $f_1$ and $f_2$ be locally Lipschitzian functions.
The function $f_1$ subdifferentially dominates $f_2$ at the point $x$ if
\begin{equation}\label{subdom}
|f^{\prime}_1(x,d)| > |f^{\prime}_2(x,d)|\;\;\; \forall d \in \mathbb R^n.
\end{equation} 
\end{definition}

An important consequence is that if $f_1$ 
subdifferentially dominates $f_2$ at $x$, then $f_1$ subdifferentially
dominates $tf_2$ for any $t$ satisfying $0\leq t \leq 1$ at $x$. The following 
lemma establishes that the conditions given by (\ref{optim}) specify a
stationary point of the objective function under the subdifferential domination
hypothesis.

\begin{lemma}
\label{lemsum}
Let $F(x,y)=f(y)+\Psi(x,y)$ and assume the existence of 
a point $x^*$ satisfying the following stationarity condition
\begin{equation}\label{optim2}
0 \in \partial f(x^*)+\partial \Psi(x^*,x^*), 
\end{equation}
with $\Psi$ satisfying Assumptions 2. If $f$ subdifferentially dominates $\Psi(x^*,.)$ at
$x^*$, then $x^*$ is a stationary point of $f$.
\end{lemma}

\noindent{\it Proof}:
For simplicity denote $\Psi(x^*,y)$ by $\Psi_{x^*}(y)$.
Since $0 \in \partial \Psi_{x^*}(x^*)$, due to Lemma \ref{heub}, 
we have $\Psi^{\prime}_{x^*}(x^*,d)\geq 0$ for
all $d\in \mathbb R^n$. Due to the first equation in (\ref{optim2}), we 
deduce the existence of two vectors $g_1$ and $g_2$ such that
\begin{align*}
& g_1+g_2=0, \\
& g_1 \in \partial f(x^*),\\
& g_2 \in \partial \Psi_{x^*}(x^*),
\end{align*} 
and thus
\begin{equation*}
\langle g_1, d\rangle+\langle g_2, d\rangle=0,
\end{equation*}
for all $d\in \mathbb R^n$. 
Since, by definition, $f^{\prime}(x,d)\geq \langle g_1, d\rangle$ and 
$\Psi_{x^*}(x^*,d)\geq \langle g_2, d\rangle$, we have
$f^{\prime}(x^*,d)+\Psi^{\prime}_{x^*}(x^*,d)\geq 0$. Thus, we obtain 
\begin{equation}\label{renoir}
f^{\prime}(x^*,d)\geq -\Psi^{\prime}_{x^*}(x^*,d).
\end{equation} 
We now proceed by contradiction. Suppose that $f^{\prime}(x^*,d)<0$
for  some $d \in \mathbb R^n$. In this case, as $\Psi^{\prime}_{x^*}(x^*,d)\geq 0$, 
equation (\ref{renoir}) implies
\begin{equation}
|f^{\prime}(x^*,d)|\leq |\Psi^{\prime}_{x^*}(x^*,d)|.
\nonumber\end{equation} 
Since this last equation contradicts the domination assumption, we
deduce that $f^{\prime}(x^*,d)\geq 0$ for all $d \in \mathbb R^n$. Using the
(second) definition (\ref{secdef}) of the subdifferential, we conclude that $0 \in
\partial f(x^*)$.
%\end{proof}

With this result in hand, we are now ready to discuss the case of the generalized
proximal point operator. Using Lemma \ref{lemsum} we first deduce the 
optimality of the fixed points of the generalized proximal point mapping 
$P_t$ with relaxation parameter $t\geq 0$ in the following straightforward lemma. 

\begin{lemma}\label{rums}
\label{ptfix}
Assume that $f(x)$ subdifferentially dominates $t\Psi(x, x)$ at each
fixed point $x^*$ of the {\em generalized proximal} mapping $P_t$ with
relaxation  parameter $t \geq 0$.  Then any fixed point $x^*$ of $P_t$
is a stationary point of $f(x)$.
\end{lemma} 

A last question of computational importance remains to be discussed. In real
life situations, one may not know whether a given regularization is
subdifferentially dominated by the function to be minimized. This problem
is easily overcome by forcing the relaxation parameter $t$ towards zero in the generalized
proximal point operator $P_t$. Indeed, the definition of subdifferential
domination and Lemma \ref{rums} prescribe that 
\begin{equation}
\label{urbule}
|f^{\prime}(x^*,d)| > t |\Psi^{\prime}_{x^*}(x^*,d)|,\;\;\; \forall d \in \mathbb R^n,
\nonumber\end{equation}  
where, as above, $x^*$ is a fixed point of $P_t$ and 
$\Psi_{x^*}^{\prime}(x^*,d)$ is the directional derivative of  
$\Psi(x^*,.)$ in the direction of $d$ at $x^*$. Therefore, one easily checks
that, given $f$ and $\Psi$, it is sufficient to take a small enough 
relaxation parameter $t$ to guarantee (\ref{urbule}).  
As a consequence, we may conclude that a safe strategy, when performing the 
generalized proximal point algorithm with a nonsmooth regularization, is to 
take a sequence of relaxation parameters $t=\lambda_k$ indexed by iteration $k$ converging
towards a sufficiently small value\footnote{We will require nevertheless that 
$(\lambda_k)_{k\in \mathbb N}$ does not vanish to push through our convergence
analysis}. 
\section{Convergence analysis}\label{Conv}
In this section, we give an asymptotic analysis of the 
{\em generalized proximal point} algorithm which does not require differentiability nor
convexity assumptions. A
Lyapunov method is the guideline of the proof where the Lyapunov function
is simply the objective $f(x)$. We show that the accumulation points of the sequence
defined by (\ref{genprox}) are locally optimal. 
Under a strict local convexity assumption convergence is established.
\subsection{Main results}
We start with the following monotonicity result.

\begin{lemma}
\label{truit} Let $f$ and $\Psi$ satisfy Assumptions 1 and 2. 
For any iteration $k >1$, the sequence $(x^k)_{k\in \mathbb N}$ satisfies
\begin{equation}\label{majgrad}
f(x^{k+1})-f(x^k)\leq -\lambda_k \Psi(x^k,x^{k+1})\leq 0.
\end{equation} 
\end{lemma}

\noindent{\it Proof}:
Iteration (\ref{genprox}) implies that $f(x^{k+1})+\lambda_k \Psi(x^k,x^{k+1})
\leq f(x^k)+\lambda_k\Psi(x^k,x^k)$. Recall that $\Psi(x^k,x^k)=0$ due
to identifiability assumption \ref{ass2} $(ii)$ and $\Psi(x^k,x^{k+1})\geq 0$
by positivity assumption \ref{ass2} $(i)$. Thus (\ref{majgrad}) follows.
%\end{proof}

We next deduce the following important property which is sometimes referred
to as ``asymptotic regularity''\cite{Bauschke:1996}.

\begin{lemma}
\label{asymreg}
Assume that there exists a real number $\lambda$ such that $\lambda_k\geq \lambda >0$.
Then, the sequence of iterates $(x^k)_{k \in \mathbb N}$ satisfies
$\lim_{k\rightarrow \infty}|x^k-x^{k+1}|=0$.
\end{lemma}

\noindent{\it Proof}:
By lemma \ref{truit}, $(f(x^k))_{k \in \mathbb N}$ is decreasing and thus
the fact that $f$ is bounded from below implies that $(f(x^k))_{k \in \mathbb N}$ converges.
Thus, the left hand side of (\ref{majgrad}) tends towards zero and
since  $\lambda_k \geq \lambda >0$, for all $k$, we have $\lim_{k\rightarrow
\infty}\Psi(x^k,x^{k+1})=0$.

We now prove that $\lim_{k\rightarrow \infty}|x^k-x^{k+1}|=0$ by
contradiction.  Assume that there exits a subsequence
$(x^{\sigma(k)})_{k\in \mathbb N}$ such that $|x^{\sigma(k)}-x^{\sigma(k)+1}|\geq 3\epsilon$  
for some $\epsilon>0$ and for all large $k$. 
The fact that $(f(x)^{\sigma(k)})_{k\in \mathbb N}\leq f(x^1)$ implies that 
$(x^{\sigma(k)})_{k\in \mathbb N}$ is bounded due to the 
inf-compactness assumption \ref{ass1} $(ii)$. Since $(x^{\sigma(k)})_{k\in \mathbb N}$ is bounded,
one can extract a convergent subsequence, and thus, we may assume
without any loss of generality that $(x^{\sigma(k)})_{k\in \mathbb N}$ is
convergent with limit $x^{\prime}$.
Using the triangle inequality, we have
$|x^{\sigma(k)}-x^{\prime}|+|x^{\prime}-x^{\sigma(k)+1}|\geq 3\epsilon$. Since
$(x^{\sigma(k)})_{k\in \mathbb N}$ converges to $x^{\prime}$, there exists a integer
$K$ such that $k\geq K$ implies $|x^{\sigma(k)}-x^{\prime}|\leq \epsilon$.
Thus for $k\geq K$ we have $|x^{\prime}-x^{\sigma(k)+1}|\geq 2\epsilon$. Now extract
from $(x^{\sigma(k)+1})_{k\geq K}$ a convergent subsequence
$(x^{\sigma(\gamma(k))+1})_{k\geq K}$ with limit $x^{\prime\prime}$. Then, using the 
same arguments as above, we
obtain $|x^{\prime}-x^{\prime\prime}|\geq \epsilon$. Finally, recall that 
$\lim_{k\rightarrow\infty}\Psi(x^k,x^{k+1})=0$. We thus have
$\lim_{k\rightarrow\infty}\Psi(x^{\sigma(\gamma(k))},x^{\sigma(\gamma(k))+1})=0$,
and, due to the fact that the sequences are bounded and $\Psi(.,.)$ is
locally Lipschitz (and therefore continuous in both variables), we have
$\Psi(x^{\prime},x^{\prime\prime})=0$. Thus assumption \ref{ass2} (ii) implies
that $|x^{\prime}-x^{\prime\prime}|=0$ and we obtain a contradiction. Hence, 
$\lim_{k\rightarrow \infty}|x^k-x^{k+1}|=0$ as claimed. 
%\end{proof}

We are now ready to establish our global convergence theorem.

\begin{theorem}
\label{theo}
Assume that the sequence $(\lambda_k)_{k\in \mathbb N}$ is bounded and satisfies
$\lambda^k \geq \lambda >0$, $\forall k \in \mathbb N$ for a given $\lambda$. 
Define $\lambda^+=\limsup_{k\rightarrow \infty}\lambda_k$. If
$f$ subdifferentially dominates $\lambda^+ \Psi(x^*,.)$ at any fixed
point $x^*$ of the operator $P_{\lambda^+}$, then every accumulation point
of  the sequence $(x^k)_{k \in \mathbb N}$ is a stationary point of $f(x)$.
\end{theorem}

\noindent{\it Proof}:
Take a convergent subsequence $(x^{\sigma(k)})_{k\in \mathbb N}$ of
$(x^k)_{k\in \mathbb N}$ with limit point $x^*$. Lemma \ref{asymreg} implies
that  $(x^{\sigma(k)+1})_{k\in \mathbb N}$ also converges to $x^*$. In
accordance with Lemma \ref{ptfix}, we need to prove the existence
of a real  $\lambda^+\geq t\geq 0$, such that $F_t(x^*,x^*)\leq F_t(x^*,x)$ for all $x \in
\mathbb R^n$. In the following, we prove that this result holds with
$t=\lambda^+$. By definition of $(x^k)_{k\in \mathbb N}$,
$F_{\lambda_k}(x^{\sigma(k)},x^{\sigma(k)+1}) \leq
F_{\lambda_k}(x^{\sigma(k)},x)$ for all $x \in  \mathbb R^n$. Therefore,
for all $x$
\begin{gather}\label{eq0}
F_{\lambda^+}(x^{\sigma(k)},x^{\sigma(k)+1})
+(\lambda_k-\lambda^+) \Psi(x^{\sigma(k)},x^{\sigma(k)+1}) \notag\\
\hspace{3cm} \leq F_{\lambda^+}(x^{\sigma(k)},x)+(\lambda_k-\lambda^+)
\Psi(x^{\sigma(k)},x).
\end{gather}
Since $\lim_{k\rightarrow
\infty}\Psi(x^{\sigma(k)},x^{\sigma(k)+1})=0$, for any $\epsilon>0$,
there exits an integer $K_1$ such that
$\Psi(x^{\sigma(k)},x^{\sigma(k)+1})\leq \epsilon$ for all $k\geq
K_1$.  On the other hand, $F_{\lambda^+}$ is continuous in both
variables, due to the locally Lipschitz property.  Fix $x \in
\mathbb R^n$. By continuity in the first and second arguments of $F_{\lambda^+}(.,.)$, 
respectively, we have for any $\epsilon > 0$ there exists $K_2\in \mathbb N$  such that 
for all $k\geq K_2$ 
\begin{equation}\label{eq1}
F_{\lambda^+}(x^*,x)\geq F_{\lambda^+}(x^{\sigma(k)},x)-\epsilon,
\end{equation}
and
\begin{equation}\label{eq2}
F_{\lambda^+}(x^*,x^*)\leq
F_{\lambda^+}(x^{\sigma(k)},x^{\sigma(k)+1})+2\epsilon.
\end{equation}
Combining equations (\ref{eq1}) and (\ref{eq2}) with (\ref{eq0}), we obtain
\begin{equation}
F_{\lambda^+}(x^*,x^*)\leq F_{\lambda^+}(x^*,x)-(\lambda^+-\lambda_k)
(\Psi(x^{\sigma(k)},x)-\Psi(x^{\sigma(k)},x^{\sigma(k)+1})) +3\epsilon.
\nonumber\end{equation}
Now, since $\lambda^+=\limsup_{k\rightarrow \infty}\lambda_k$, there
exists an integer $K_3$ such that $\lambda^+-\lambda_k\geq -\epsilon$ for
all $k\geq K_3$.
Then for all $k\geq \max \{K_1,K_2,K_3\}$, we easily obtain
\begin{equation*}
F_{\lambda^+}(x^*,x^*)\leq
F_{\lambda^+}(x^*,x)+\epsilon\Psi(x^{\sigma(k)},x)-\epsilon^2+3\epsilon.
\end{equation*}
Since $(x^{\sigma(k)})_{k \in \mathbb N}$ is bounded, and due to
Assumption \ref{ass2}(iv) and Assumption \ref{ass2}(iii),
there exists an upper bound $C$ such that $\Psi(x^{\sigma(k)},x)\leq C$
for all $k\geq \max \{K_1,K_2,K_3\}$. Thus, we have
\begin{equation*}
F_{\lambda^+}(x^*,x^*)\leq
F_{\lambda^+}(x^*,x)+\epsilon C-\epsilon^2+3\epsilon.
\end{equation*} 
Since no assumption was made on $x$, this holds for any
$x \in  \mathbb R^n$. Thus, letting $\epsilon$ tend to zero, we see that 
$x^*$ is a fixed point of ${\rm argmin}_{x \in \mathbb R^n}F_{\lambda^+}(x^k,x)$. 
Furthermore, recall that $f(x)$ subdifferentially dominates $\lambda^+ \Psi(x,x)$ at the point
$x^*$. Therefore, Lemma \ref{ptfix} implies that $x^*$ is a
stationary point of $f(x)$. 
%\end{proof}

\subsection{Convergence to local minima under additional convexity assumptions}
Let $S^*$ be the set of accumulation
points of the sequence $(x^k)_{k\in \mathbb N}$. We first establish the
following lemma

\begin{lemma}
Let $f$ and $\Psi$ satisfy Assumptions 1 and 2. Then, for a given starting point $x^1$,
the set $x^*$ of accumulation points of the sequence $(x^k)_{k \in \mathbb N}$ is compact 
and connected.
\end{lemma}

\noindent{\it Proof}:
This result follows directly from the fact that
$\lim_{k\rightarrow\infty}|x^{k+1}- x^k|=0$ and from \cite[Theorem
28.1]{Ostrowski:1966}.
%\end{proof}

\begin{corollary}
Suppose, in addition to the assumptions of Theorem \ref{theo}, that
$f(x)$ is strictly convex in an open neighborhood
$\mathcal N$ of an  accumulation point $x^*$ of
$(x^k)_{k\in \mathbb N}$. Then the sequence $(x^k)_{k\in \mathbb N}$ converges to a local
minimizer of $f(x)$.
\end{corollary}

\noindent{\it Proof}:
We obtained in Theorem \ref{theo} that every accumulation point of
$(x^k)_{k\in \mathbb N}$  is a stationary point of $f(x)$. Since $f(x)$ is
strictly convex over $\mathcal N$, the set of stationary points of
$f(x)$ belonging to $\mathcal N$ reduces to singleton. Thus $x^*$ is
the unique stationary point in $\mathbb N$ of $f(x)$, and {\em a fortiori}, the unique
accumulation point of $(x^k)_{k\in \mathbb N}$ belonging to $\mathcal N$. To
complete the proof, it remains to show that there is no accumulation
point in the exterior of  $\mathcal N$. For that purpose, consider an
open ball $\mathcal B$ of center $x^*$ and radius $\epsilon$  included
in $\mathcal N$. Then, $x^*$ is the unique accumulation point in
$\mathcal B$.  Moreover, any accumulation point $x^{\prime}$, lying in
the exterior of  $\mathcal N$ must satisfy $|x^*-x^{\prime}|\geq
\epsilon$, and we obtain a contradiction with the fact that $S^*$ is
connected.  Thus every accumulation point lies in $\mathcal N$, from
which we conclude that $x^*$ is the only accumulation point of
$(x^k)_{k\in \mathbb N}$  or, in other words, that $(x^k)_{k\in \mathbb N}$
converges towards $x^*$.  Finally, notice that the strict convexity
of $f(x)$ over $\mathcal N$ implies  that $x^*$ is a local minimizer
and the proof is completed.
%\end{proof}

\section{Bundle implementations: the convex case} \label{bundle1}
The study of the generalized proximal point algorithm gives an elegant framework
for the exploration of a large class of regularizations. 
In the two next sections, we use the bundle framework introduced by Lemar\'echal
\cite{Hirriart:1993} to make
this approach tractable in applications. The present section introduces the main
ideas governing incorporation of nonquadratic regularization functions into
the bundle mechanism in the case where the functions are convex. More precisely
we will require that 
$f(x)$ be convex and $\Psi(x^{\prime},x)$ be convex with respect to $x$ 
for every $x^{\prime}$ in $\mathbb R^n$. 
The case where the functions $f$ and $\Psi$ are only required to be 
locally Lipschitz will be discussed in the next section. 

\subsection{Background}
Bundle methods have been widely recognized as a very efficient technique
for minimization of nondifferentiable functions. They can be interpreted
as implementations of stepwise approximations of the proximal point algorithm 
\cite{Lemarechal:1997,Kiwiel:1990}.
However, to our knowledge, bundle implementations of general nonquadratic 
regularizations have not been considered in the litterature. This may be 
due to be fact that bundle methods are essentially
sequential quadratic programs, which seems to exclude the possibility 
of more general regularization functions. 

The main idea behind application of bundle methods to nonquadratic 
proximal algorithms is the following. The generalized proximal iteration
is approximated using subgradient information with respect to $f$ and $\Psi$
about the current iterate $x^k$. Let $\{y^j\}_{J^k}$ be a set of points
in a neighborhood of $x^k$, indexed by $j\in J^k$. For any point $y
\in \mathbb R^n$, let $g(y)$ (resp. $h_k(y)$) denote an arbitrary
subgradient in $\partial f(y)$ (resp. $\partial \Psi(x^k,y)$). Hence,
any point $y^j$, $j \in J^k$ allows to define approximate models
\begin{equation}
\hat{f}_j(x)=f(y^j)+\langle g(y^j),x-y^j\rangle
\nonumber\end{equation}
and 
\begin{equation}
\hat{\Psi}_j(x^k,x)=\Psi(x^k,y^j)+\langle h_k(y^j),x-y^j\rangle
\nonumber\end{equation}
of $f$ and $\Psi$ respectively. To stabilize the model
$\hat{\Psi}_j(x^k,x)$ a quadratic term
$\frac12|x-x^k|^2$ is added and the following approximation is 
obtained for the function
$F_{\lambda_k}(x^k,x)$ in the proximal point algorithm
\begin{align*}
\hat{F}_{\lambda_k,j}(x^k,x)= & f(y^j)+
\langle g(y^j), x-y^j\rangle+\lambda_k\big(\Psi(x^k,y^j)+
\langle h_k(y^j),x-y^j\rangle \\
& +\frac12|x-x^k|^2\big).
\end{align*}
Note that convex functions are bounded from below by their
local first order approximations and that the quadratic
stabilization term is independent from $y^j$. 
Hence, the best approximation of $F_{\lambda_k}(x^k,x)$ can be
obtained by gathering the information at each $y^j$ in the neighborhood
of $x^k$ via the following max-function,
\begin{align*}
\hat{F}_{\lambda_k}(x^k,x)= & \max_{j\in J^k} \{f(y^j)+
\langle g(y^j),x-y^j\rangle + \lambda_k
\big(\Psi(x^k,y^j)+\langle h_k(y^j),x-y^j\rangle\\
& +\frac12 |x-x^k|^2 \},
\end{align*}
Equivalently, 
$\hat{F}_{\lambda_k}(x^k,x)$ can be written
\begin{equation}
\hat{F}_{\lambda_k}(x^k,x)=\max_{j\in J^k} \{
f(x^k)-\alpha_j^k+\langle s^k(y^j),x-x^k \rangle\} +\frac12
\lambda_k |x-x^k|^2,
\nonumber\end{equation}
where $\alpha_j^k$ is the accumulated linearization error due to the
fact that the subgradients are computed at the points $y^j$ instead of 
$x^k$, i.e.
\begin{equation}
\label{urlub}
\alpha_j^k=f(x^k)-f(y^j)-\langle g(y^j),x^k-y^j \rangle
-\Psi(x^k,y^j)-\langle h_k(y^j),x^k-y^j\rangle, 
\end{equation}
and $s(y^j)$ is the accumulated subgradient 
\begin{equation}
s^k(y^j)=g(y^j)+\lambda_k h_k(y^j).
\nonumber\end{equation}

With this model in hand, we are now ready to introduce the bundle mechanism. 
First, an approximate proximal step is taken which defines a candidate for
the next iterate, 
\begin{equation}\label{cand}
y^{k+1}={\rm argmin}_{x \in \mathbb R^n} \hat{F}_{\lambda_k}(x^k,x),
\end{equation}
or, equivalently,
\begin{equation}
y^{k+1}={\rm argmin}_{x \in \mathbb R^n}\max_{j\in J^k} \{
f(x^k)-\alpha_j^k+\langle s(y^j),x-x^k \rangle\} +\frac{\lambda_k}2 |x-x^k|^2\}.
\nonumber\end{equation}
The bundle approach is a strategy
in order to decide wether the accumulated subgradient information at every point
of the set $\{y^j\}$, $j\in J^k$, is sufficiently accurate so that a reliable proximal 
step can be achieved, i.e. $x^{k+1}=y^{k+1}$. It is well known that iteration
(\ref{cand}) may even not provide a descent step with regard to $f$ if the 
subgradient information is innaccurate Hence, a reasonable selection of candidates
$\{y^{k+1}\}$
for a proximal step should be based on a test of descent in the objective function
$f$. In order to implement such a test, a parameter $\delta_k$
is computed, representing the expected decrease given the model $\hat{F}_{\lambda_k}$,
\begin{equation}
\delta_k=f(x^k)-\hat{F}_{\lambda_k}(y^{k+1},x^k),
\nonumber\end{equation}
or equivalently
\begin{equation}
\delta_k=f(x^k)-\max_{j\in J^k} \{f(x^k)-\alpha_k^j
+\langle s^k(y^j),y^{k+1}-x^k\rangle\}-\frac{\lambda_k}2 |y^{k+1}-x^k|^2.
\nonumber\end{equation}
The decrease obtained at $y^{k+1}$ is then compared to a fraction $m \in (0,1)$ of
the expected decrease $\delta_k$, following the rule
\begin{itemize}
\item if $f(x^k)-f(y^{k+1})\geq m \delta_k$, then a descent step is taken,
i.e. $x^{k+1}=y^{k+1}$,
\item otherwise $x^{k+1}=x^k$.
\end{itemize}
In either case, the subgradient information at $y^{k+1}$ is 
collected and is incorporated to the polyhedral approximation of $f$ and 
$\Psi$ at iteration $k+1$. In this manner the accuracy of the approximation 
$\hat{F}(x^k,x)$ to $f(x)+\lambda_k\Psi(x^k,x)$ improves at each iteration.

\subsection{A generalized proximal bundle method}
In this section, we present the details of our bundle implementation. 

\begin{algorithm}[Generalized Proximal Bundle Method for Convex Functions] 
\label{bund1}

{\bf Step 0} {\rm (Initialization)} Choose a final accuracy parameter $\delta_s >0$ and
a parameter $m \in (0,1)$. Choose the starting point $x^1 \in \mathbb R^n$. Set
$k=1$, $y^1=x^1$, $J^1=\{1\}$, $g^1=g(y^1)$, $h^1=h_1(y^1)$, $s^1=
g^1+h_1^1$ and $\alpha_1^1=0$.

{\bf Step 1} {\rm (Proximal step)} Compute
\begin{equation}\label{proxstep}
y^{k+1}={\rm argmin}_{x \in \mathbb R^n}\max_{j\in J^k} \big\{
f(x^k)-\alpha_k^j+\langle s^j,x-x^k \rangle\} +\frac{\lambda_k}2 |x-x^k|^2 \big\}.
\end{equation}

{\bf Step 2} {\rm (Descent test)} Set 
\begin{equation}
\delta_k=f(x^k)-\max_{j\in J^k} \{f(x^k)-\alpha_k^j
+\langle s^j,y^{k+1}-x^k\rangle\}-\frac{\lambda_k}2 |y^{k+1}-x^k|^2.
\nonumber\end{equation}
\begin{itemize}
\item if $f(x^k)-f(y^{k+1})\geq m \delta_k$, then set $x^{k+1}=y^{k+1}$ \hfill 
{\rm ({\bf descent step})},
\item otherwise set $x^{k+1}=x^k$ \hfill {\rm ({\bf null step})}.
\end{itemize}

{\bf Step 3} {\rm (Stopping criterion)} If $\delta_k < \delta_s$, then stop.

{\bf Step 4} {\rm (Variable updating)} Set $J^{k+1}=J^k\cup \{k+1\}$, choose
$\lambda_{k+1}$, set $g^{k+1}=g(y^{k+1})$, $h^{k+1}=h_{k+1}(y^{k+1})$, 
$s^{k+1}=g^{k+1}+\lambda_{k+1}h^{k+1}$, set $\alpha_{k+1}^j=
\alpha_k^j$, for $j \in \{1,\ldots,k\}$,
$\alpha_{k+1}^{k+1}=f(x^{k+1})-\langle s^{k+1},x^{k+1}-y^{k+1}\rangle$.
Finally, increase $k$ by 1 and go to Step 1.
\end{algorithm} 

We now establish some preliminary results concerning this algorithm. The
following lemmas are standard in the analysis of bundle methods. 
In particular, the following lemma expresses proximal step (\ref{proxstep}) in the form of a dual 
quadratic program:
\begin{lemma}{\rm \cite{Bonnans:1997}}
\label{benlem}
Consider the constrained minimization problem
\begin{align}\label{dualquad}
\bar{u} = & {\rm argmin}_{u\in \mathbb R^{k}}
\Big\{ \frac12 |\sum_{j \in J^k} u_j s^j|^2+\lambda_k \sum_{j \in J^k} u_j \alpha_k^j\Big\}\\
\text{ subject to } & \\
        & u \in \Delta_k=\big\{ z\in (0,1)^k \mid \sum_{j \in J^k} z_j=1 \big\}\nonumber.
\end{align}
Then, 
\begin{equation}\label{pstepdual}
{\rm (i)} \hspace{.5cm} y^{k+1}=x^k-\frac1{\lambda_k} \sum_{j \in J^k} \bar{u}_j s^j,
\end{equation}
where $y^{k+1}$ solves the proximal step (\ref{proxstep}) in Step 1 of Algorithm 1, and the polyhedral
component $\Phi_k(x)=\max_{j\in J^k} \{f(x^k)-\alpha_k^j+\langle s^j,x-y^j \rangle \}$ of 
$\hat{F}_{\lambda_k}(x^k,x)$ satisfies
\begin{equation} 
{\rm (ii)} \hspace{.5cm} \Phi_k(y^{k+1})=f(x^k)- \frac1{\lambda_k} 
\big|\sum_{j \in J^k} u_j s^j\big|^2-
\sum_{j \in J^k} \bar{u}_j \alpha_k^j,
\nonumber\end{equation}.
\end{lemma}

In the sequel, we adopt the notation
\begin{equation}\label{zob}
p^k=\sum_{j\in J^k}\bar{u}_j s^j.
\end{equation}
We must also introduce another dual
variable which plays an important role in the bundle method. For this purpose,
notice that the quadratic program given by (\ref{dualquad}) is equivalent to
\begin{align}\label{duald}
\bar{u}= & {\rm argmin}_{u \in \mathbb R^k} \frac12 \big|\sum_{j \in J^k} u_j s^j\big|^2\\
\text{subject to } & \\
         & u \in \Delta_k \nonumber \\
         & \sum_{j\in J_k} u_j\alpha_j \leq \epsilon_k, \label{urdate}
\end{align}
where $\lambda_k$ in (\ref{dualquad}) is identified with the Lagrange multiplier associated with the 
contraint (\ref{urdate}). Using lemma \ref{benlem} (ii), and the definition of $\delta_k$ in step 2
of Algorithm 1 
\begin{equation}\label{tord}
f(x^k)-\Phi(y^{k+1})-\frac{\lambda_k}2|y^{k+1}-x^k|^2=
\frac1{2\lambda_k}\big|\sum_{j\in J^k}\bar{u}_js^j\big|^2+\sum_{j\in J^k}\bar{u}\alpha_k^j.
\end{equation}
Recalling that $\bar{u}$ is solution to the quadratic program (\ref{duald}), we have
\begin{equation}
\label{osher}
\epsilon_k=\sum_{j \in J^k}\bar{u}_j\alpha_k^j. 
\end{equation}
Therefore (\ref{tord}) is equivalent to (see also \cite{Bonnans:1997})
\begin{equation}\label{turl}
\delta_k=\epsilon_k+\frac1{2\lambda_k}|p^k|^2.
\end{equation}
With this result in hand, we obtain the following lemma.
\begin{lemma}\label{proson}
Let $p^k$ be the direction defined by (\ref{zob}). Then, 
\begin{equation}
p^k \in \partial_{\epsilon_k} \big(f(x^k)+\lambda_k \Psi(x^k,x^k)\big).
\nonumber\end{equation}
\end{lemma}
\noindent{\it Proof}:
The polyhedral function $\Phi_k(x)$ may be written
\begin{equation}
\Phi_k(x)=\max_{j \in J^k}\{f(y^j)+\langle g^j,x-y^j\rangle +
\lambda_k\big( \Psi(x^k,y^j)+\langle h^j,x-y^j\rangle \big) \}.
\nonumber\end{equation}
Then, due to convexity, we have
\begin{equation}
f(x)\geq f(y^j)+\langle g^j,x-y^j \rangle
\nonumber\end{equation}
and
\begin{equation}
\Psi(x^k,x)\geq \Psi(x^k,x^k)+\langle h^j, x-y^j \rangle.
\nonumber\end{equation}
Therefore,
\begin{equation}\label{tutute}
f(x)+\lambda_k \Psi(x^k,x) \geq \Phi_k(x).
\end{equation}
On the other hand, optimality in the proximal step (\ref{proxstep}) (Step 1 of Algorithm 1), gives
\begin{equation}
0 \in \partial \Phi_k(y^{k+1})+\lambda_k (y^{k+1}-x^k).
\nonumber\end{equation}
Noticing that $\lambda_k(y^{k+1}-x^k)=p^k$, we thus obtain that
\begin{equation}
p^k=\partial \Phi_k(y^{k+1}).
\nonumber\end{equation}
Hence, using the subgradient inequality, we obtain
\begin{equation}\label{tote}
\Phi_k(x)\geq \Phi_k(y^{k+1})+\langle p^k,x-y^{k+1}\rangle.
\end{equation}
Recalling that, due to Lemma \ref{benlem} (ii)
\begin{equation}
\Phi_k(y^{k+1})=f(x^k)-\frac1{\lambda_k}|p^k|^2-\sum_{j\in J_k}\bar{u}_j \alpha_k^j,
\nonumber\end{equation}
and combining this result with (\ref{tote}), (\ref{tutute}) becomes
\begin{equation}
f(x)+\lambda_k \Psi(x^k,x)\geq f(x^k)-\frac1{\lambda_k}|p^k|^2
-\sum_{j\in J^k} \bar{u}_j \alpha_k^j + \langle p^k,x-y^{k+1}\rangle.
\nonumber\end{equation} 
Since $-\frac1{\lambda_k}|p^k|^2=\langle p^k,y^{k+1}-x^k\rangle$, we obtain
\begin{equation}
f(x)+\lambda_k \Psi(x^k,x)\geq f(x^k)+\langle p^k,x-x^k\rangle -
\sum_{j \in J^k} \bar{u}_j \alpha_k^j.
\nonumber\end{equation}
Furthermore, recalling that 
\begin{equation}
\epsilon_k=\sum_{j\in J^k} \bar{u}_j \alpha_k^j,
\nonumber\end{equation}
and that $\Psi(x^k,x^k)=0$, we easily obtain
\begin{equation}
f(x)+\lambda_k\Psi(x^k,x) \geq f(x^k)+\lambda_k \Psi(x^k,x^k)+
\langle p^k,x-x^k\rangle-\epsilon_k,
\nonumber\end{equation}
which is equivalent to 
\begin{equation}
p^k \in \partial_{\epsilon_k} \big(f(x^k)+\lambda_k \Psi(x^k,x^k) \big). 
\nonumber\end{equation}

%\end{proof}
This last result shows in particular the well known fact that bundle methods may be 
interpreted as $\epsilon-$subgradient methods (see \cite{Hirriart:1993}). The main 
feature of bundle methods is therefore
the control of the parameter $\epsilon_k$ via (\ref{turl}) and the descent 
step/null step strategy using the expected decrease $\delta_k$. We now discuss 
convergence of this method to a minimizer of $f$. 

\subsection{Convergence}
For the convex case, the proof of convergence of the generalized proximal bundle method is similar
to the proof of convergence of the standard form of bundle methods. In the
sequel, $K$ denotes the set of indices $k$ where a descent step is taken.
We start with the following lemma.
\begin{lemma}\label{propsur}
Consider the following convergent sequences, $x^k\rightarrow \bar{x}$,
$y^k\rightarrow \bar{y}$, $\epsilon_k\rightarrow \bar{\epsilon}$,
$\lambda^k\rightarrow \bar{\lambda}$ and 
$p^k\rightarrow \bar{d}$. Assume that 
\begin{equation}
p^k\in \partial_{\epsilon_k}\big(f(y^k)+\lambda_k\Psi(x^k,y^k)\big)
\nonumber\end{equation}
for all $k\in \mathbb N$. Then,
$\bar{d} \in \partial_{\bar{\epsilon}}\big(f(\bar{y})+\bar{\lambda}\Psi(\bar{x},\bar{y})\big)$.
\end{lemma}
\noindent{\it Proof}:
Fix $y$ in $\mathbb R^n$. The fact that $p^k\in \partial_{\epsilon_k}\big(f(y^k)+
\lambda_k\Psi(x^k,y^k)\big)$
implies that for fixed $y$
\begin{equation}
f(y)+\lambda_k\Psi(x^k,y)\geq f(y^k)+\lambda_k\Psi(x^k,y^k)+
\langle p^k,y-y^k\rangle-\epsilon_k.
\nonumber\end{equation}
Therefore,
\begin{align}\label{peute}
f(y)+\lambda_k\Psi(x^k,y)\geq & f(y^k)
+\lambda_k\Psi(x^k,\bar{y})+\lambda_k\big(\Psi(x^k,y^k)-\Psi(x^k,\bar{y})\big)\\
 & +\langle p^k,y-y^k\rangle-\epsilon_k. \nonumber
\end{align}
Now, since $\Psi^k$ is assumed locally Lipschitz over $\mathbb R^n\times \mathbb R^n$, 
and since $(x^k)_{k\in \mathbb N}$ and $(y^k)_{k\in \mathbb N}$ are bounded
due to convergence, there exists a constant $C$ such that
\begin{equation}
|\Psi(x^k,y^k)-\Psi(x^k,\bar{y})|\leq C \sqrt{|x^k-x^k|^2+|y^k-\bar{y}|^2}
= C |y^k-\bar{y}|.
\nonumber\end{equation}
Therefore, as $y^k$ converges to $\bar{y}$
\begin{equation}
\lim_{k\rightarrow\infty} \Psi(x^k,y^k)-\Psi(x^k,\bar{y})=0.
\nonumber\end{equation}
Using a similar argument one easily deduces that $\Psi(x^k,y)$ converges
to $\Psi(\bar{x},y)$ and that $\Psi(x^k,\bar{y})$ converges to $\Psi(\bar{x},\bar{y})$.
Furthermore, Assumption \ref{ass1} (ii) implies (lower semi-) continuity  of $f$ and
thus, passing to the limit in (\ref{peute}) gives
\begin{equation}
f(y)+\bar{\lambda}\Psi(\bar{x},y)\geq f(\bar{y})+\bar{\lambda}\Psi(\bar{x},\bar{y})+
\langle \bar{d},y-\bar{y}\rangle-\bar{\epsilon},
\nonumber\end{equation}
which proves the desired result.
%\end{proof}

We will also need a more technical result the proof of which can be found in \cite{Bonnans:1997}.
\begin{lemma}\label{ealor}
Let $f$ and $\Psi$ satisfy assumptions \ref{ass1} and \ref{ass2}. Let
$K$ be the set of indices where a descent step is taken. 
Let $f_*=\lim_{k\rightarrow\infty, k\in K}f(x^k)$. Then
\begin{equation}
\sum_{k\in K} \delta_k \leq \frac{f(x^1)-f_*}{m}.
\nonumber\end{equation}
\end{lemma}

The following convergence analysis is divided into two parts. In the first part, we
consider the case where an infinite number of descent steps are taken. Then, we will turn to
the finite case. 
\begin{lemma}\label{lemprem}
Assume that $f$ and $\Psi$ satisfy Assumptions \ref{ass1} and \ref{ass2}.
Assume that $f$ subdifferentially dominates $\lambda^+ \Psi$ at every point $x^*$
such that $0\in \partial f(x^*)+\lambda^+ \partial \Psi(x^*,x^*)$ for some 
real number $\lambda^+>0$. Assume in addition that 
$\limsup_{k\rightarrow \infty} \lambda_k \leq \lambda^+$
and that $\liminf_{k\rightarrow \infty}\lambda_k \geq \lambda^->0$ for some real number $\lambda^-$ and for
all $k$ in $\mathbb N$. Finally, assume that an infinite number of descent
steps is taken. Then, every accumulation point of $(x^k)_{k\in \mathbb N}$
is a minimizer of $f$. 
\end{lemma}
\noindent{\it Proof}:
Using (\ref{turl}), we obtain
\begin{equation}
\frac1{2\lambda_k} |p^k|^2=\delta_k-\epsilon_k \geq \delta_k,
\nonumber\end{equation}
for all $k \in K$. Hence, we have
\begin{equation}
\label{turbule}
\sum_{k\in K}|p^k|^2 \leq 2\lambda_k \sum_{k\in K} \delta_k.
\end{equation}
Recall that Assumption \ref{ass1} implies that $f$ is bounded from below. Thus,
Lemma \ref{ealor} implies that $\lim_{k\rightarrow\infty,
k\in K} \delta_k=0$, and therefore, (\ref{turl}) implies that
$ \lim_{k\rightarrow\infty, k\in K} \epsilon_k=0$. 
The assumptions on the sequence $(\lambda_k)_{k\in \mathbb N}$ imply boundedness 
of $\lambda_k$ and thus, equation (\ref{turbule}) gives
\begin{equation}
\sum_{k\in K} |p^k|^2<+\infty,
\nonumber\end{equation}
which proves that $\lim_{k\rightarrow\infty, k\in K} |p^k|=0$. 
On the other hand, inf-compactness of $f$ and the fact that 
$(f(x^k))_{k\in K}$ is strictly decreasing imply that the sequence 
$(x^k)_{k\in K}$ is bounded. Now, take a convergent subsequence 
$(x^{\sigma(k)})_{k\in K}$ and let $\bar{x}$ be its limit. Since $(\lambda_k)_{k\in K}$
is bounded, a convergent subsequence can be taken, say 
$(\lambda_{\sigma(\gamma(k))})_{k\in K}$, tending towards $\bar{\lambda}$. 
Then, Lemma \ref{proson}
and Lemma \ref{propsur} together imply
\begin{equation}
0\in \partial \big(f(\bar{x})+\bar{\lambda} \Psi(\bar{x},\bar{x})\big),
\text{ and }
0\in \partial f(\bar{x})+\bar{\lambda}\partial \Psi(\bar{x},\bar{x}).
\nonumber\end{equation}
Finally, since every accumulation point $\bar{\lambda}$ satisfies 
$\bar{\lambda}\leq \lambda^+$, subdifferential domination and lemma \ref{lemsum} imply that 
$0\in \partial f(\bar{x})$. Finally, convexity implies that $\bar{x}$ is a minimizer
of $f$ is convex, which finishes the proof.
%\end{proof}

We now discuss the case of a finite number of descent steps.

\begin{lemma}
Let the functions $f$ and $\Psi$ satisfy the assumptions in Lemma \ref{lemprem}
Assume that only a finite number of descent
steps is taken and let $k_0$ denote the index of the last descent step. 
Then, $x^{k_0}$ is a minimizer of $f$. 
\end{lemma}
\noindent{\it Proof}:
The proof is easily adapted from \cite{Lemarechal:1993}. The definition
of $\delta_k$ gives 
\begin{equation}
f(x^k)-\delta^k=\max_{j\in J^k}\{f(y^j)+\lambda_k\Psi(x^k,y^j)
+\langle s^j,y^{k+1}-y^j \rangle \}+\frac{\lambda_k}{2}|y^{k+1}-x^k|^2.
\nonumber\end{equation}
Take $k\geq k_0$. Since $x^k=x^{k_0}$ for all $k\geq k_0$, we obtain
\begin{equation}\label{pulpot}
f(y^j)+\lambda_k\Psi(x^{k_0},y^j)+\langle s^j,y^{k+1}-y^j\rangle
+\frac{\lambda_k}2|y^{k+1}-x^{k_0}|^2\leq f(x^{k_0})-\delta^k,
\end{equation}
for all $j$ in $J^k$. On the other hand, nondescent steps imply
\begin{equation}
f(x^{k})-m\delta^j \leq f(y^j),
\nonumber\end{equation}
for $k\geq j\geq k_0$. Thus, as $f(x^k)=f(x^{k_0})$ for $k\geq k_0$, (\ref{pulpot}) gives
\begin{equation}
\lambda_k\Psi(x^{k_0},y^j)+\langle s^j,y^{k+1}-y^j\rangle
+\frac{\lambda_k}{2}|y^{k+1}-x^{k_0}|^2\leq m\delta^j-\delta^k,
\nonumber\end{equation}
for $k\geq j\geq k_0$,
which implies, due to positivity of both $\Psi$ and the quadratic term, that
\begin{equation}\label{polpot}
\langle s^j,y^{k+1}-y^j\rangle \leq m\delta^j-\delta^k,
\end{equation}
for $k\geq j\geq k_0$.
The case $j=k_0$ gives $y^j=y^{k_0}=x^{k_0}$, the last equality coming
from the fact that $k_0$ is the index of the last descent step. Using this fact along
with (\ref{pulpot})
\begin{equation}
\langle s^{k_0},y^{k+1}-x^{k_0}\rangle +\frac{\lambda_k}2|y^{k+1}-x^{k_0}|^2
\leq -\delta_k\leq 0.
\nonumber\end{equation}
The left term of this last equation is a quadratic form in the variable $y^{k+1}$,
which cannot take negative values except on a bounded neighborhood of $x^{k_0}$. Hence, 
the sequence $(y^k)_{k\in \mathbb N}$ is bounded. As a consequence, 
$s^j$ is also bounded, due to the local boundedness property of the
subgradients of $f(.)$ and $\Psi(x^{k_0},.)$ and the boundedness of $\lambda_k$. 
On the other hand, due to refinement of the polyhedral approximation, and since
$\lambda_k=\lambda_{k_0}$ and $x^{k}=x^{k_0}$ for $k\geq k_0$, 
\begin{equation}
\hat{F}_{\lambda_k}(x^{k},x)\leq \hat{F}_{\lambda_{k+1}}(x^{k+1},x),
\nonumber\end{equation}
for all $x$ in $\mathbb R^n$, 
which proves that $(\delta_k)_{k \geq k_0}$ is decreasing and thus,
has a limit. We now prove that this limit is zero. Indeed, take
a convergent subsequence $(y^{\sigma(k)})_{k\geq k_0}$. Hence, for
any $j$ satisfying $\sigma(j)\geq k_0$, we have, using (\ref{polpot}),
\begin{equation}
\langle s^{\sigma(j)},y^{\sigma(j+1)}-y^{\sigma(j)}\rangle\leq m \delta^{\sigma(j)}
-\delta^{\sigma(j+1)-1}.
\nonumber\end{equation}
Since $\lim_{j\rightarrow\infty}|y^{\sigma(j+1)}-y^{\sigma(j)}|=0$ and due to
boundedness of $s^j$, we have
\begin{equation}\label{pulbot}
\lim_{j\rightarrow\infty} m\delta^{\sigma(j)}-\delta^{\sigma(j+1)-1}\geq 0.
\end{equation}
Recall that $(\delta^k)_{k \in \mathbb N}$ is convergent. Let $\bar{\delta}$ be its
limit. By continuity, (\ref{pulbot}) gives
\begin{equation}
(m-1)\bar{\delta} \geq 0.
\nonumber\end{equation}
Since $0< m < 1$, and $\bar{\delta}\geq 0$ by positivity of $\delta^k$ for all $k
\in \mathbb N$, we have $(m-1)\bar{\delta}=0$ which implies $\bar{\delta}=0$, as desired.
Now recall that
\begin{equation}
\delta^k=\epsilon_k+\frac1{2\lambda_k}|p^k|^2,
\nonumber\end{equation}
to conclude, using lower and upper boundedness of $\lambda_k$, that 
$\lim_{k\rightarrow\infty} \epsilon_k=0$ and $\lim_{k\rightarrow\infty}|p^k|=0$.
As in lemma \ref{lemprem}, we therefore obtain that
\begin{equation}
0\in \partial \big(f(x^{k_0})+\bar{\lambda}\Psi(x^{k_0},x^{k_0})\big)
\nonumber\end{equation}
for every accumulation point $\bar{\lambda}$ of $(\lambda_k)_{k\in \mathbb N}$.
Since every accumulation point $\bar{\lambda}$ satisfies $\bar{\lambda}\leq \lambda^+$, 
subdifferential domination gives
\begin{equation}
0\in \partial f(x^{k_0}),
\nonumber\end{equation}
which proves optimality of $x^{k_0}$ since $f$ is convex. 
%\end{proof}

Combining these two lemmas yields the following convergence
theorem for our bundle implementation of the generalized proximal point
algorithm for convex optimization.
\begin{theorem}
Let the functions $f$ and $\Psi$ satisfy the assumptions in Lemma \ref{lemprem}
Then, every accumulation point of $(x^k)_{k\in \mathbb N}$ defined by
Algorithm \ref{bund1} is a minimizer of $f$. 
\end{theorem}



\section{Bundle implementations: the nonconvex case}
\label{bundle2}
\subsection{Preliminary comments}
In this section, a bundle approach for implementing generalized proximal steps
is developed for the case of nonconvex functions. The main ideas remain the same
as in the convex situation. Nevertheless, several modifications need to be introduced
in order to overcome the difficulties associated with nonconvexity. The 
algorithmic structure used in the sequel is similar to the one proposed by
Kiwiel \cite{Kiwiel:1985} and therefore inherits useful convergence properties. 
This allows us to concentrate on the particular problems
induced by the use of our generalized regularization. We first introduce
the main characteristics of the method.

One important property satisfied  in the convex
case by the polyhedral approximations to $f$ and $\Psi$ is that they are lower approximations,
i.e. they lie below the original functions $f$ and $\Psi$. 
Therefore $\alpha_k^j$ in
(\ref{urlub}) is always positive. This property no longer holds in the nonconvex case. 
Nevertheless, a similar property can be obtained when $\alpha_k^j$ is defined by
\begin{equation}
\alpha_k^j=|f(x^k)-f(y^j)-\langle g^j, x^k-y^j \rangle|+
\lambda_k |\Psi(x^k,y^j)+\langle h^j,x^k-y^j\rangle |
\nonumber\end{equation}
with $g^j \in \partial f(y^j)$ and $h^j \in \partial \Psi(x^k,y^j)$ and
recalling that $\Psi(x^k,x^k)=0$. The approximate model is therefore written in the same
manner as in the convex case, i.e.
\begin{equation}
\hat{F}_{\lambda_k}(x^k,x)=\max_{j\in J^k} \{
f(x)-\alpha_j^k+\langle s^j ,x-x^k \rangle\} +\frac12
\lambda_k |x-x^k|^2,
\nonumber\end{equation}
using the new definition of $\alpha_k^j$ and where as before $s^j$ is defined
by $s^j \in g^j+\lambda_k h^j$ and
$\{y^j\}_{j \in J^k}$ is a collection of points in a neighborhood of $x^k$. 

With this lower approximation in hand, an approximate proximal step is taken, yielding a
precandidate for $x^{k+1}$
\begin{equation}
\label{prurle}
z^{k+1}={\rm argmin}_{z\in \mathbb R^n} \max_{j\in J^k} \{f(x^k)-\alpha_k^j
+\langle s^j,x-x^k\rangle\}+\frac{\lambda_k}2 |x-x^k|^2.  
\nonumber\end{equation}
Let $d^k$ denote the direction given by $d^k=z^{k+1}-x^k$ or, in other words,
\begin{equation}
d^k=\frac1{\lambda_k} \sum_{j\in J^k} \bar{u}_js^j,
\nonumber\end{equation}
and let $p^k$ denote the convex combination of collected subgradients
\begin{equation}
p^k=\sum_{j\in J^k} \bar{u}_js^j.
\nonumber\end{equation}
In the following implementation, $z^{k+1}$ does not directly provide a candidate for a 
descent step. Rather, a linesearch is performed along the direction $d^k$. The goal
of this linesearch is twofold. Firstly, the linesearch provides a improvement of the
result obtained using only first order approximation of $f$ and $\Psi$, and
therefore refines the implementation given in section 5 for the convex case. 
Secondly, if the linesearch only leads to a null step, the linesearch provides a 
systematic refinement of the local subgradient information. Indeed, as discussed
in \cite{Hirriart:1993} and\cite{Mifflin:1977}, under some assumptions of 
weak upper semi-smothness, a new candidate is obtained which satisfies
\begin{equation}
\langle s^{k+1}, d^k\rangle >0
\nonumber\end{equation} 
when a null step is taken. 

The question of optimality is solved in the following manner. The main difference
between the convex case and the nonconvex case is the manner in which an approximate
subdifferential is generated. In the convex case, the same
threshold $\delta_k$ was used to control both descent tests and the size
of the $\epsilon$-subdifferential (recall relation (\ref{osher}) and Lemma \ref{proson}). 
In the nonconvex case, the use of the Goldstein $\epsilon$-subdifferential,
based on the distance from past subgradients to the current iterate, leads to
a different strategy. A locality measure $a^k$ is chosen for bounding from 
above the ``size'' of the Golstein $\epsilon$-subdifferentials. 
In Lemma \ref{loctrut} below, we show that $p^k$ belongs to the sum of Goldstein
$\epsilon$-subdifferentials at $x^k$, whose size is controled by $a^k$. 
For this purpose, we introduce a property called {\em quazi-symmetry}. 
{\em Quazi-symmetry} is then added to the assumptions on $\Psi$.
Along the iterations, past subgradients are deleted from the bundle information 
when the distance from the corresponding neighbor $y^j$ to the current 
iterate becomes excessive. More precisely, deletion occurs when 
\begin{equation}
|p^k|>m_aa^k,
\nonumber\end{equation} 
where $m_a$ is a 
fixed parameter supplied to the algorithm. In other words, $|p^k|$ should not vanish faster
than a fraction of the locaty measure. Convergence of $p^k$ towards zero is obtained
using the same type of proof as in \cite{Kiwiel:1985} which implies convergence of
$a^k$ to zero. 
\subsection{A generalized proximal bundle method for nonconvex functions}

In the sequel, at iteration $k$ we will use the notation
$s(y)$ for any element of $\partial f(y)+\lambda_k\partial \Psi(x^k,y)$ and 
\begin{equation}
\alpha(x,y)=|f(x)-f(y)-\lambda_k\Psi(x,y)-\langle s(y),x-y\rangle |
\nonumber\end{equation}
\begin{algorithm}
\label{bund2}

{\bf Step 0} {\rm (Initialization)} Select the starting point $x^1$ and a final
accuracy tolerance $\epsilon_s>0$. Choose the line search parameters $m_L$ and 
$m_R$ $\in (0,1)$, $m_\alpha \in (0,m_R)$, $\kappa\in (0,1)$ and
a positive reset tolerance $m_a$. Set the threshold stepsize $t^1=1$ and the reset
indicator $r_a^1=1$. Set $J^1=\{1\}$, $y^1=x^1$, $s^1=s(y^1)$, 
and $s_1^1=0$. Set the counters $k=1$, $l=0$ and $k(0)=1$.

{\bf Step 1} {\rm (Proximal step)} Solve the following quadratic subproblem 
\begin{equation}\label{quad}
\begin{array}{l}
\bar{u}={\rm argmin}_{u\in \mathbb R^{\rm card(J^k)}}
\{\frac12 \Big| \sum_{j\in J^k}u_j s^j\Big|^2+\sum_{j\in J^k}u_j \alpha^k_j\}\\
\text{ subject to }\\
u_j\geq 0, \hspace{.5cm} \sum_{j\in J^k}u_j=1.
\end{array}
\end{equation}  
Set 
\begin{align}\label{chtecite}
p^k & =-\sum_{j\in J^k}\bar{u}_j s^j\\
\alpha^k & =\sum_{j\in J^k}\bar{u}_j\alpha_k^j \nonumber\\
v^k & =|p^k|^2+\alpha^k \nonumber
\end{align}
If $r_a^k=1$, set $a^k=\max \{s_j^k \;\;\mid\;\; j\in J^k\}$.

{\bf Step 2} {\rm (Stopping criterion)} If $\max\{|d^k|,m_a a^k \}\leq \epsilon_s$, 
terminate. Otherwise go to Step 3.

{\bf Step 3} {\rm (Resetting test)} If $|d^k|>m_a a^k$ then go to Step 5. Otherwise go to
Step 4.

{\bf Step 4} {\rm (Resetting)} 
(i) If $r_a^k=0$ then set $r_a^k=1$ and go to step 1.

(ii) If $J^k\neq\{k(l)\}$ then delete from $J^k$ one among the elements $j$ with the 
largest $s_j^k>0$ and go to Step 1.

{\bf Step 5} {\rm (Linesearch)} Compute two stepsizes $t_L^k$ and $t_R^k$, 
$t_L^k\leq t_R^k$ using procedure \ref{linsch} below and such that the two corresponding 
points $x^{k+1}=x^k+t^k_Ld^k$ and $y^{k+1}=x^k+t^k_Rd^k$ satisfy either
\begin{itemize}

\item $f(x^k)-f(x^{k+1})-\lambda_k\Psi(x^k,x^{k+1})>m_Lt|v^k|$ 
and 
\begin{equation}
\begin{cases}
t_L=t_R>t^k \text{ or } \\
\alpha(x^k,x^{k+1})>m_{\alpha}|v^k|, 
\end{cases}
\nonumber
\end{equation}
\hfill {\bf (descent step)}

or

\item $-\alpha(x^{k+1},y{k+1})+\langle s^{k+1},d^k \rangle \geq m_R v^k$, 
$t_R^k \leq t^k$ and $t_L^k=0$, \hfill {\bf (null step)}

with 
\begin{equation}
\label{subloc}
s^{k+1}=g(y^{k+1})+\lambda_k h_k(y^{k+1}).
\end{equation}

\end{itemize}

{\bf Step 6} {\rm (Threshold stepsize updating)} If $t_L^k=0$, set $t^{k+1}=\kappa t^k$.
Otherwise, set $t^{k+1}=1$, $k(l+1)=k+1$ and increase $l$ by $1$.

{\bf Step 7} {\rm (Updating)} Update $\lambda_k$. Set
\begin{align}
n_j^{k+1} & =n_j^k+|x^{k+1}-x^k|, \\
n_{k+1}^{k+1} & =|x^{k+1}-y^{k+1}| 
\end{align}
$J^{k+1}=J^k\cup \{k+1\}$, calculate $a^{k+1}=\max\{a^k + |x^{k+1}-x^k|,n_{k+1}^{k+1}\}$
and set $r_a^k=0$.

{\bf Step 8} Increase $k$ by $1$ and go to Step 1.
\end{algorithm}

We now describe the line search procedure.

\begin{procedure}[Linesearch] 
\label{linsch}

$(i)$ Set $t_L=0$, $t=t_U=1$ and $m=(m_R-m_{\alpha}+m_L)$.

$(ii)$ If $f(x^k)-f(x^k+td^k)-\lambda_k\Psi(x^k,x^t+td^k)>mt|v^k|$ set $t_L=t$. 
Otherwise set $t_U=t$.

$(iii)$ If $f(x^k)-f(x^k+td^k)-\lambda_k\Psi(x^k,x^k+td^k)>m_Lt|v^k|$ and 
either $t_L\geq t^k$ or $\alpha(x^k,x^k+td^k>m_\alpha |v^k|$, set $t_L^k=t_R^k=t$, 
$s^{k+1}=s(x^k+t_Ld^k)$ and return.

$(iv)$ If $t<t^k$ and 
\begin{equation}
-\alpha(x^{k},x^k+td^k)+\langle s(x^k+td^k),d^k \rangle \geq m_R v^k,
\nonumber\end{equation}
set $t_R^k=t$, $t_L^k=0$, $s^{k+1}=s(x^k+t_Rd^k)$ and return.

$(v)$ Choose $t=\frac{t_L+t_U}{2}$ and go to $(ii)$.
\end{procedure}

\subsection{Intermediate results}
The convergence analysis for Algorithm 2 will require the additional assumptions below. 
\begin{assumptions}
\label{ass3}
{\rm (i)} The functions $f(y)$ and $\Psi(x,y)$ are semi-smooth in the variable $y$ for
any $x$ in $\mathbb R^n$. 

{\rm (ii)} There exists a function $q(.,.)$, $\mathbb R_+^2\mapsto \mathbb R_+$ satisfying\\
a. $q(.,b)$ is increasing for all $b\in \mathbb R_+$ and $q(a,.)$ is increasing for all 
$a\in \mathbb R_+$,\\
b. $q(.,.)$ is continuous at the point $(0,0)$,\\
c. $\lim_{(a,b)\rightarrow (0,0)} q(a,b)=0$, \\
and such that if $g\in \partial^a \Psi(x,y)$, $a\geq 0$ then 
$g \in -\partial^{q(a,|x-y|)} \Psi(y,x)$. 
\end{assumptions}
We assume that the functions $f$ and $\Psi$ satisfy assumptions 
\ref{ass1}, \ref{ass2} and \ref{ass3}. 
We now prove two essential lemmas. The first one establishes that $p^k$ belongs
to certain Goldstein $\epsilon$-subdifferentials at $x^k$. The second lemma proves
an important continuity result at stationary points. 
\begin{lemma}\label{loctrut}
At any iteration $k$, define $\lambda_k^+=\sup_{j\in J^k}\lambda_j$. Fix a point
$x^*$ in the level set $\{x\;\mid\;f(x)\leq f(x^1)\}$ and define $b^k=|x^k-x^*|$.
Then, under Assumptions \ref{ass3}, for any $k \in \mathbb N$, the direction $p^k$ 
defined in Step 1 of Algorithm \ref{bund2} satisfies
\begin{equation}
p^k \in \partial^{a^k} f(x^k)-\lambda_k^+\partial^{q(a^k+b^k,a^k+b^k)+a^k}
\Psi(x^*,x^k),
\nonumber\end{equation} 
where $a^k$ is the locality measure defined in steps 1 and 6.
\end{lemma}
\noindent{\it Proof}:
Using (\ref{subloc}), we have for any $j\in J^k$,
\begin{equation}
s^j \in \partial f(y^j)+\lambda_j \partial \Psi(x^j,y^j).
\nonumber\end{equation}
Therefore, using the fact that $|x^k-y^j|\leq a^k$ and the definition of the 
Goldstein $\epsilon-$subdifferential
\begin{equation}
s^j \in \partial^{a^k} f(x^k)+\lambda_j \partial^{a^k} \Psi(x^j,x^k),
\nonumber\end{equation}
and 
\begin{equation}
s^j \in \partial^{a^k} f(x^k)+\lambda_j \partial^{a^k+b^k} \Psi(x^j,x^*).
\nonumber\end{equation}
Now, by assumption \ref{ass3} $(ii)$ of quazisymmetry,
\begin{equation}
s^j \in \partial^{a^k} f(x^k)-\lambda_j \partial^{q(a^k+b^k,|x^*-x^j|)} \Psi(x^*,x^j),
\nonumber\end{equation} 
Thus, using the fact that $|x^k-x^j|\leq a^k$, $|x^*-x^j|\leq |x^*-x^k|+a^k$,
and that $q(a,.)$ is increasing for all fixed $a\in \mathbb R_+$ as required by
Assumptions \ref{ass3}(ii) a, we obtain
\begin{equation}
s^j \in \partial^{a^k} f(x^k)-\lambda_j \partial^{q(a^k+b^k,|x^*-x^k|+a^k)+a^k} 
\Psi(x^*,x^k),
\nonumber\end{equation}
and, recalling that $b^k=|x^*-x^k|$,
\begin{equation}\label{flut}
s^j \in \partial^{a^k} f(x^k)-\lambda_j \partial^{q(a^k+b^k,a^k+b^k)+a^k} \Psi(x^*,x^k),
\end{equation}
Since $\lambda_k^+\geq \lambda_j$, we have
\begin{equation}
-\lambda_j \partial^{q(a^k+b^k,a^k+b^k)+a^k} \Psi(x^*,x^k)\subset -\lambda_k^+ 
\partial^{q(a^k+b^k,a^k+b^k)+a^k} \Psi(x^*,x^k)
\nonumber\end{equation}
which, when combined with (\ref{flut}), gives 
\begin{equation}
s^j \in \partial^{a^k} f(x^k)-\lambda_k^+\partial^{q(a^k+b^k,a^k+b^k)+a^k}\Psi(x^*,x^k),
\nonumber\end{equation} 
for any $j \in J^k$. 
Now, due to (\ref{chtecite}), $p^k$ is a convex combination of the subgradients
$s^j$, $j \in J^k$. Recalling that $\epsilon-$Goldstein subdifferentials are closed and convex 
sets, we have
\begin{equation}
p^k \in \partial^{a^k} f(x^k)-\lambda_k^+\partial^{q(a^k+b^k,a^k+b^k)+a^k}
\Psi(x^*,x^k)
\nonumber\end{equation}
as desired.
%\end{proof}

The following lemma will be needed in the sequel.
\begin{lemma}
\label{rizhom}
All the sequences $(x^k)_{k \in \mathbb N}$, $(y^k)_{k \in \mathbb N}$, 
$(s^k)_{k \in \mathbb N}$ and
$(p^k)_{k \in \mathbb N}$ defined in the regularized bundle algorithm are bounded.
\end{lemma} 
\noindent{\it Proof}:
Due to assumption \ref{ass1} (i) of inf-compactness, and the fact that algorithm 
\ref{bund2} 
is a descent method, the sequence $(x^k)_{k \in \mathbb N}$ is bounded. The remainder of
the proof can be easily adapted from \cite[Lemma 3.3]{Kiwiel:1985}.
%\end{proof}

\begin{lemma}
\label{turlop}
Suppose that there exist sequences $(x^k)_{k\in \mathbb N}$, $(p^k)_{k\in \mathbb N}$
and $(a^k)_{k\in \mathbb N}$ and a point $\bar{x}\in \mathbb R^n$ such that 
$\lim_{k\rightarrow \infty}x^k=\bar{x}$, $\lim_{k\rightarrow \infty}p^k=0$, 
$\lim_{k\rightarrow \infty}a^k=0$. Define  
$\lambda_k^+=\max_{j\in J^k}\lambda_j$ for all $k$ in $\mathbb N$.
Assume that $f$ subdifferentially dominates $\lambda^+ \Psi$ at every point $x^*$
such that $0\in \partial f(x^*)+\lambda^+ \partial \Psi(x^*,x^*)$ for some 
real number $\lambda^+>0$. Assume in addition that 
$\limsup_{k\rightarrow \infty} \lambda^+_k \leq \lambda^+$
and that $\liminf_{k\rightarrow\infty}\lambda_k \geq \lambda^->0$ for some 
real number $\lambda^-$ and for every $k\in \mathbb N$. 
Then $\bar{x}$ is a stationary point of $f(x)$.
\end{lemma}

\noindent{\it Proof}:
Lemma \ref{loctrut}
above implies the existence of two sequences $(p_1^k)_{k\in \mathbb N}$ and 
$(p_2^k)_{k\in \mathbb N}$ such that $p^k=p_1^k+p_2^k$
and 
\begin{align*}
p_1^k & \in \partial^{a^k} f(x^k),\\
p_2^k & \in -\lambda_k^+\partial^{q(a^k+b^k,a^k+b^k)+a^k}\Psi(x^*,x^k).
\end{align*}
Due to lemma \ref{rizhom}, $(x^k)_{k \in \mathbb N}$ is bounded and thus the local 
boundedness of the subgradients implies that the sequences $(p_1^k)_{k\in \mathbb N}$ and 
$(p_2^k)_{k\in \mathbb N}$ are bounded. Thus, we can extract convergent subsequences
$(p_1^{\sigma(k)})_{k\in \mathbb N}$, and $(p_2^{\sigma(\gamma(k))})_{k\in \mathbb N}$ with
respective limits $p_1^*$ and $p_2^*$ such that 
\begin{equation}\label{sumzero}
p_1^*+p_2^*=0.
\end{equation} 
On the other hand, using \cite[lemma 3.1]{Kiwiel:1996}, we obtain 
\begin{equation}\label{prems}
p_1^*\in \partial f(\bar{x}).
\end{equation} 
Furthermore, since $\lim_{k\rightarrow \infty}x^k=\bar{x}$, we have 
$\lim_{k\rightarrow \infty}b^k=0$ and thus 
$\lim_{k\rightarrow \infty} q(a^k+b^k,b^k+a^k)+a^k=0$, due to Assumption \ref{ass3}(ii) c. 
Using \cite[lemma 3.1]{Kiwiel:1996} and the fact that 
$\lambda^+=\limsup_{k\rightarrow \infty}\lambda_k^+$ 
one easily obtains
\begin{equation}
p_2^*\in -\lambda^+\partial \Psi(\bar{x},\bar{x}).
\nonumber\end{equation}
Now, due to quazisymmetry Assumption \ref{ass3} (ii) in the case $a=0$ and 
$x=y=\bar{x}$, we have
\begin{equation}\label{deuze}
p_2^*\in \lambda^+\partial \Psi(\bar{x},\bar{x}).
\end{equation}
Finally, (\ref{sumzero}), (\ref{prems}) and (\ref{deuze}) together give 
\begin{equation}\label{plurp}
0 \in \partial f(\bar{x})+\lambda^+\partial \Psi(\bar{x},\bar{x}).
\end{equation}
The proof is then easily completed using Lemma \ref{lemsum}.
%\end{proof}
\subsection{Convergence}
Convergence of Algorithm \ref{bund2} is established using the properties of
Algorithm 2.1 in \cite{Kiwiel:1985}, to which it is structurally very close, and the 
preliminary properties established in the previous section. In order to adapt the study
of \cite{Kiwiel:1985} to our generalized proximal bundle method, we make the following
important observation. 

The behavior of the method principally relies on the output of the linesearch (Procedure
\ref{linsch}). In particular, in the case of descent steps, the linesearch procedure gives
\begin{equation}
f(x^k)-f(y^{k+1})-\lambda_k \geq m_Lt_Lv^k.
\nonumber\end{equation}
Due to positivity of $\Psi$, we obtain
\begin{equation}
\label{pourque}
f(x^k)-f(y^{k+1})\geq m_Lt_Lv^k.
\end{equation}
This last equation corresponds exactly to the output given by Algorithm 2.1 of
\cite{Kiwiel:1985} in the descent case. On the other hand, the case of a null step gives
\begin{equation}
f(x^k)-f(y^{k+1})-\lambda_k \Psi(x^k,y^{k+1}) \leq m_Rt_Rv^k,
\nonumber\end{equation}
which is equivalent to
\begin{equation}
F_{\lambda_k}(x^k,x^k)-F_{\lambda_k}(x^k,y^{k+1})\leq m_Rt_Rv^k. 
\nonumber\end{equation}
This corresponds exactly to the output of the linesearch procedure in Algorithm 2.1
of \cite{Kiwiel:1985} when $f(.)$ is replaced by $F_{\lambda_k}(x^k,.)$. Using these
strong similarities between the two methods, we now establish convergence of our procedure.
The first step is study the linesearch procedure. In particular, we will need
the following simple lemma. 
\begin{lemma}
\label{preced}
Let $f_1$ and $f_2$ be two weakly upper semismooth functions. Then $f_1+f_2$ is also
weakly upper semismooth. 
\end{lemma}
\noindent{\it Proof}:
The proof is an immediate consequence of the definition. Indeed, take 
$p^k \in \partial \big(f_1+f_2\big)(x+t^kd)$. Then, we have
$p^k=p_1^k+p_2^k$ with 
\begin{equation}
p_1^k \in \partial f_1(x+t^kd),
\nonumber\end{equation}
and
\begin{equation}
p_2^k \in \partial f_2(x+t^kd). 
\nonumber\end{equation}
Since $f_1$ and $f_2$ are weakly upper semismooth, we have
\begin{equation}
\liminf_{k\rightarrow\infty} \langle p_1^k,d\rangle\geq \limsup_{t\downarrow 0}
\frac{f_1(x+td)-f_1(x)}{t}
\nonumber\end{equation}
and
\begin{equation}
\liminf_{k\rightarrow\infty} \langle p_2^k,d\rangle\geq \limsup_{t\downarrow 0}
\frac{f_2(x+td)-f_2(x)}{t}.
\nonumber\end{equation}
Therefore, we easily obtain
\begin{equation}
\liminf_{k\rightarrow\infty} \langle p_1^k+p_2^k,d\rangle\geq \limsup_{t\downarrow 0}
\frac{f_1(x+td)+f_2(x+td)-f_1(x)-f_2(x)}{t}.
\nonumber\end{equation}
Recalling that the Lipschitz property is preserved by summation, this last equation
proves the desired result. 
%\end{proof}

\begin{lemma}
The linesearch (Procedure \ref{linsch}) terminates in a 
finite number of steps. 
\end{lemma}
\noindent{\it Proof}:
Due to Lemma \ref{preced} and Assumption \ref{ass3}(i), the function 
$f(.)+\lambda_k \Psi(x^k,.)$ is weakly upper semismooth. Therefore, the proof of Theorem
4.1 (a) in \cite{Mifflin:1977} can be easily adapted to our linesearch procedure. 
%\end{proof}

The following lemma due to Kiwiel \cite{Kiwiel:1983} will be used in the proof
of Theorem \ref{ombre}. 
In particular, given $\epsilon>0$, this lemma implies existence of a finite
number of successive null steps without reset such that $|d^k|<\epsilon$. 
\begin{lemma}\label{trulput}{\rm \cite{Kiwiel:1983,Kiwiel:1985}}
Let $w^k=\frac12 |d^k|^2+\alpha^k$ and assume $t^k_L=0$ and $r_a^k=0$. Then, there exists 
a constant $C$ independent of $k$ such that 
\begin{equation}\label{trucak}
0\leq w^{k+1} \leq w^k -(1-m_R)^2 (w^k)^2/8C^2.
\end{equation}
\end{lemma}
We are now in a position to establish convergence of the generalized proximal bundle method. 
We first study the case of a finite number of descent steps. 
\begin{lemma}
\label{oplote}
Assume that $f$ and $\Psi$ satisfy Assumptions \ref{ass1}, \ref{ass2} and \ref{ass3}.
Assume that $f$ subdifferentially dominates $\lambda^+ \Psi$ at every point $x^*$
such that $0\in \partial f(x^*)+\lambda^+ \partial \Psi(x^*,x^*)$ for some 
real number $\lambda^+>0$. Assume in addition that 
$\limsup_{k\rightarrow \infty} \lambda_k \leq \lambda^+$
and that $\liminf_{k\rightarrow \infty}\lambda_k \geq \lambda^->0$ for some 
real number $\lambda^-$ and for every $k\in \mathbb N$. Assume in addition that
only a finite number of descent step is taken and let $k_0$ denote the index 
of the last descent step. Then $x^{k_0}$ is a stationary point of $f$.
\end{lemma}
\noindent{\it Proof}:
Using the same type of proof as in Lemma 3.4 of \cite{Kiwiel:1985}, we obtain that
\begin{equation}
\lim_{k\rightarrow\infty}p^k=0
\text{ and }
\lim_{k\rightarrow\infty}a^k=0.
\nonumber\end{equation}
Therefore, Lemma \ref{turlop} implies the desired result. 
%\end{proof}

The case of an infinite number of descent steps follows. 
\begin{lemma}
Assume that $f$ and $\Psi$ satisfy the assumptions in lemma \ref{oplote}. 
Assume that a infinite number of descent
steps is taken. Then, every accumulation point of $(x^k)_{k\in \mathbb N}$
is a minimizer of $f$. 
\end{lemma}
\noindent{\it Proof}:
Since the linesearch procedure \ref{linsch} gives equation (\ref{pourque}) in the case of a descent
step, the sequence $(f(x^k))_{k\in \mathbb N}$ is nondecreasing. Furthermore, 
one can easily check that the proof of lemma 3.5 in \cite{Kiwiel:1985} can also be adapted to our
method. Therefore, we obtain 
\begin{equation}
\lim_{k\rightarrow\infty}p^k=0
\text{ and }
\lim_{k\rightarrow\infty}a^k=0.
\nonumber\end{equation}
Now, take any convergent subsequence from $(x^k)_{k\in \mathbb N}$ and let 
$\bar{x}$ denote its limit. Then, Lemma \ref{turlop} implies the desired result. 
%\end{proof}

Combining these two last lemmas, we obtain our convergence theorem for the 
generalized proximal bundle method \ref{bund2}.
\begin{theorem}
\label{ombre}
Let $f$ and $\Psi$ satisfy the assumptions in lemma \ref{oplote}
Then, every accumulation point of $(x^k)_{k\in \mathbb N}$
defined by Algorithm \ref{bund2} is a minimizer of $f$. 
\end{theorem}


\section{Numerical experiments}
\label{simul}
In this section, we illustrate the bundle implementation of Section \ref{bundle2} for
a particular nonsmooth optimization problem arising in adaptive digital signal processing 
for channel equalization. This problem is formulated as follows
\begin{equation}\label{cma}
\min_{x \in \mathbb R^n} f(x),\hspace{2cm}\text{ where } f(x)=\sum_{i=1}^n |1-|a_i^Tx|^p|^{q}. 
\end{equation}
for a given set of $n-$dimensional vectors $(a_i)_{1\leq i \leq n}$. For odd $p$ or $q$, the
function is nondifferentiable.
A popular algorithm for implementing (\ref{cma}) is the well known constant modulus algorithm
(CMA) \cite{Godard:1980,Treichler:1985,Tanrikulu:1997}. The objective function has the
shape represented in Figure \ref{obj}, for the case $n=1$, $p=2$, $q=1$. 

Numerical experiments were performed in order to compare
the computational properties of the three following regularization functions.  
\begin{align*}
\Psi_1(x,y) & =\frac12 |x-y|^2 \text{ (Moreau-Yosida)},\\
\Psi_2(x,y) & =\sum_{i=1}^n|x_i-y_i| \text{ ($l_1$ norm)}\\
\Psi_3(x,y) & =\sum_{i=1}^n \log (1+|x_i-y_i|),
\end{align*}
These three functions offer smooth convex, nonsmooth convex and nonsmooth nonconvex 
alternatives; see Figures \ref{mor}, \ref{l1} and \ref{log}.

We will compare the three regularizations on the basis of the number of iterations required, 
the number of function/subgradient calls and the computed minimum value of the objective 
function. The algorithms were run over a range of the relaxation parameters $\lambda_k=\lambda$
which we kept constant along iterations. The parameters of the regularized bundle method are 
set as in Table \ref{param}.
Our results are reported in Table \ref{tableau}. The symbols $it$, $f^*$
and $nb$ denote, respectively, the minimum number of iterations $N$ necessary to achieve 
$f(x^N)\leq \epsilon_s$, and the total number of function/subgradient calls. 

The table clearly demonstrates that the use of nonsmooth regularization functions
is competitive with the standard Moreau-Yosida regularization. Moreover, the best
candidate often appears to be one of the nonsmooth regularizations $\Psi_2$ and
$\Psi_3$ with regard to the number of function/subgradient calls and the number of 
iterations required for the algorithm to terminate.
 
Algorithmic refinements of the basic regularized bundle method proposed in this paper
will be investigated in future work. Notice that the theoretical results demonstrated
in Section \ref{bundle2} may apply to more sophisticated approaches recently developed
for the standard Moreau-Yosida regularization, including trust region techniques 
\cite{Kiwiel:1996} and second order methods \cite{Lemarechal:1997} in the 
case of convex objective and regularization functions.
 
\section{Conclusion}
This paper developed a convergence analysis for the proximal point algorithm
with general regularization functions. Strong convergence results were obtained
without any assumption on differentiability nor convexity of the regularization. 
A sufficient condition, denoted ``subdifferential domination'', on the regularization 
was defined in order for the method to converge toward a stationary point. Similar 
convergence results were given for regularized bundle implementations of the proximal
method. The resulting bundle algorithms enjoys additional degrees of freedom when compared to the 
standard procedures, since many types of regularization may be
easily implemented in the procedure. Numerical experiments
demonstrate that nonconvex and even nonsmooth regularizations are
competitive with the standard Moreau-Yosida regularization.


\bibliographystyle{plain}
\bibliography{nonsmooth}

\newpage

\begin{figure}[ht] 
\begin{center}
\caption{The objective $f$ and its derivative $f^{\prime}$.}
\epsfig{figure=fcma.eps,width=6cm}
\label{obj}
\end{center} 
\end{figure} 
\begin{figure}
\begin{center}
\caption{ Moreau-Yosida regularization $\Psi_1$}
\epsfig{figure=quadreg.eps,width=5cm}
\label{mor}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\caption{ $l_1$ norm regularization $\Psi_2$}
\epsfig{figure=l1reg.eps,width=5cm}
\label{l1}
\end{center}
\end{figure}
\begin{figure}
\begin{center}
\caption{ $log(1+|.|)$  regularization $\Psi_3$}
\epsfig{figure=logreg.eps,width=5cm}
\label{log}
\end{center}
\end{figure}
\begin{table}
\caption{Parameters of the regularized bundle algorithm}
\begin{center}\footnotesize
\begin{tabular}{|r|l|}
\hline
$n$ & 2 \\
$p$ & 2 \\
$q$ & 1 \\
$\epsilon_s$ & $10^{-5}$ \\
$m_L$ & .1 \\
$m_R$ & .3 \\
$m_{\alpha}$ & .15 \\
$\kappa$ & .8 \\
\hline
\end{tabular}
\label{param}
\end{center}
\end{table}
\begin{table}
\caption{Numerical results for the three regularizations}
\begin{center}\footnotesize
\begin{tabular}{|c|lcr|lcr|lcr|}
\hline
\multicolumn{1}{|l|}{} & \multicolumn{3}{|c|}{$\Psi_1$}
& \multicolumn{3}{|c|}{$\Psi_2$} & \multicolumn{3}{|c|}{$\Psi_3$}\\\hline
$\lambda$ & $it$ & $f^*$ & $nb$ & $it$ & $f^*$  & $nb$ & $it$ & $f^*$  & $nb$ \\\hline
.05       & 155  & 0.0001 & 278  & 183  & 0.0000 & 353   & 203  & 0.0000 & 347  \\\hline
.15       & 198  & 0.0000 & 314  &  96  & 0.0002 & 180   &  90  & 0.0000 & 164  \\\hline
.20       & 220  & 0.0005 & 360  &  73  & 0.0000 & 137   & 121  & 0.0000 & 227  \\\hline
.25       & 182  & 0.0002 & 336  & 102  & 0.0000 & 194   & 145  & 0.0000 & 280  \\\hline
.30       & 138  & 0.0002 & 257  &  90  & 0.0000 & 172   &  79  & 0.0001 & 149  \\\hline
.35       &  79  & 0.0003 & 141  & 124  & 0.0000 & 245   &  88  & 0.0001 & 166  \\\hline
.40       & 136  & 0.0002 & 208  & 136  & 0.0000 & 267   & 118  & 0.0000 & 229  \\\hline
.45       & 148  & 0.0014 & 287  & 112  & 0.0000 & 219   &  78  & 0.0000 & 148  \\\hline
.50       &  74  & 0.0004 & 139  & 129  & 0.0000 & 254   & 147  & 0.0000 & 288  \\\hline
.55       &  83  & 0.0000 & 154  & 115  & 0.0000 & 222   &  84  & 0.0001 & 161  \\\hline
.60       &  99  & 0.0008 & 186  &  65  & 0.0000 & 127   &  75  & 0.0000 & 147  \\\hline
.65       & 123  & 0.0003 & 267  &  96  & 0.0000 & 186   &  96  & 0.0000 & 186  \\\hline
.70       & 165  & 0.0005 & 269  &  84  & 0.0000 & 163   &  66  & 0.0000 & 127  \\\hline
.75       & 106  & 0.0001 & 203  & 102  & 0.0000 & 199   &  58  & 0.0000 & 108  \\\hline
.80       &  82  & 0.0000 & 158  &  74  & 0.0000 & 145   & 101  & 0.0000 & 210  \\\hline
.85       &  95  & 0.0003 & 179  &  66  & 0.0001 & 129   &  67  & 0.0001 & 129  \\\hline
\end{tabular}
\end{center}
\label{tableau}
\end{table}


\end{document}
