Multi-block Bregman proximal alternating linearized minimization and its application to orthogonal nonnegative matrix factorization

TL;DR

This paper introduces BPALM and A-BPALM algorithms using Bregman distances for solving structured nonconvex problems, with convergence analysis and application to orthogonal nonnegative matrix factorization.

Contribution

It proposes novel multi-block proximal algorithms with convergence guarantees for nonconvex problems, applied specifically to ONMF with closed-form solutions.

Findings

Sequences converge to critical points of the objective

Global convergence under KL inequality

Preliminary numerical results on ONMF

Abstract

We introduce and analyze BPALM and A-BPALM, two multi-block proximal alternating linearized minimization algorithms using Bregman distances for solving structured nonconvex problems. The objective function is the sum of a multi-block relatively smooth function (i.e., relatively smooth by fixing all the blocks except one) and block separable (nonsmooth) nonconvex functions. It turns out that the sequences generated by our algorithms are subsequentially convergent to critical points of the objective function, while they are globally convergent under KL inequality assumption. Further, the rate of convergence is further analyzed for functions satisfying the {\L}ojasiewicz's gradient inequality. We apply this framework to orthogonal nonnegative matrix factorization (ONMF) that satisfies all of our assumptions and the related subproblems are solved in closed forms, where some preliminary numerical results is reported.

Tables

Table 1. Table 1 . A comparison among BPALM, A-BPALM1, and A-BPALM2 for the Hubble image, where the algorithms stopped after 90 seconds of running time. Here, F e ​ r ​ r ​ o ​ r subscript 𝐹 𝑒 𝑟 𝑟 𝑜 𝑟 F_{error} stands for the relative fidelity error ‖ X − U k ​ V k ‖ F / ‖ X ‖ F subscript norm 𝑋 superscript 𝑈 𝑘 superscript 𝑉 𝑘 𝐹 subscript norm 𝑋 𝐹 \|X-U^{k}V^{k}\|_{F}/\|X\|_{F} while O e ​ r ​ r ​ o ​ r subscript 𝑂 𝑒 𝑟 𝑟 𝑜 𝑟 O_{error} denotes the orthogonal error ‖ I − V k ​ ( V k ) T ‖ F subscript norm 𝐼 superscript 𝑉 𝑘 superscript superscript 𝑉 𝑘 𝑇 𝐹 \|I-V^{k}(V^{k})^{T}\|_{F} . In each row, the smallest number of F e ​ r ​ r ​ o ​ r subscript 𝐹 𝑒 𝑟 𝑟 𝑜 𝑟 F_{error} and F e ​ r ​ r ​ o ​ r subscript 𝐹 𝑒 𝑟 𝑟 𝑜 𝑟 F_{error} are displayed in bold.

Penalty par.	BPALM		A-BPALM1		A-BPALM2
$\lambda$	$F_{error}$	$O_{error}$	$F_{error}$	$O_{error}$	$F_{error}$	$O_{error}$
1	$6.75\times {10}^{-2}$	$8.28\times {10}^{-2}$	$5.23\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}$	$2.60\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}$	$5.32\times {10}^{-2}$	$3.25\times {10}^{-2}$
10	$1.54\times {10}^{-1}$	$4.35\times {10}^{-2}$	$6.04\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}$	$8.36\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟑}}$	$9.33\times {10}^{-2}$	$2.16\times {10}^{-2}$
100	$2.05\times {10}^{-1}$	$3.28\times {10}^{-2}$	$9.21\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}$	$5.15\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟑}}$	$1.99\times {10}^{-1}$	$1.01\times {10}^{-2}$
1000	$2.07\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏}}$	$3.51\times {10}^{-2}$	$2.09\times {10}^{-1}$	$2.51\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟑}}$	$2.48\times {10}^{-1}$	$6.83\times {10}^{-3}$
10000	$2.09\times {10}^{-1}$	$3.44\times {10}^{-2}$	$2.62\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟐}}$	$1.45\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟑}}$	$2.55\times {10}^{-1}$	$6.33\times {10}^{-3}$

Equations

$$\underset{\mathbf{x}=(x_1,\dots ,x_N)\in {\mathbb{R}}^{{\sum }_i{n}_i}}{\mathrm{minimize}}\text{leavevmode}\text{leavevmode}\phi (\mathbf{x})\equiv f(\mathbf{x})+\sum _{i=1}^N{g}_i(x_i),$$

$$\partial f(x)=$$

$$\hat{\partial }f(x)=$$

$$h_{\mathbf{x}}^i(z):=h(\mathbf{x}+U_i(z-x_i)).$$

$$h_{\mathbf{x}}^i(x_i)\ge h_{\mathbf{x}}^i(y_i)+⟨{\nabla }_{i}h(\mathbf{y}),x_i-y_i⟩+\frac{{\sigma }_h^i}{2}\Vert x_i-y_i{\Vert }^2\forall \mathbf{x},\mathbf{y}\in B({\mathbf{x}}^{\star };\delta );$$

$${\operatorname{D}_{h}(\bm{y},\bm{x}){}\coloneqq{}{\mathopen{}\left\{\begin{array}[]{l @{~~} >{\text{if~}}l<{\text{,}} }h(\bm{y})-h(\bm{x})-\langle{}{\nabla}h(\bm{x}){},{}\bm{y}-\bm{x}{}\rangle\hfil\leavevmode\nobreak\ \leavevmode\nobreak\ &\leavevmode\nobreak\ }\bm{x}\in\operatorname{int}\operatorname{dom}h\\ \infty\hfil\leavevmode\nobreak\ \leavevmode\nobreak\ &\lx@intercol\text{otherwise.}\hfil\lx@intercol\end{array}\right.\mathclose{}}$$

$$D_h(\mathbf{x}+U_i(y_i-x_i),\mathbf{x})$$

$$=h_{\mathbf{x}}^i(y_i)-h_{\mathbf{x}}^i(x_i)-⟨{\nabla }_{i}h(\mathbf{x}),y_i-x_i⟩,$$

$${\varphi }_{\mathbf{x}}^i(z):=L_{i}h(\mathbf{x}+U_i(z-x_i))-f(\mathbf{x}+U_i(z-x_i))$$

$$f(\mathbf{x}+U_i(y_i-x_i))\le f(\mathbf{x})+⟨{\nabla }_{i}f(\mathbf{x}),y_i-x_i⟩+L_i{D}_h(\mathbf{x}+U_i(y_i-x_i),\mathbf{x});$$

$$⟨{\nabla }_{i}f(\mathbf{x})-{\nabla }_{i}f(\mathbf{y}),x_i-y_i⟩\le L_i⟨{\nabla }_{i}h(\mathbf{x})-{\nabla }_{i}h(\mathbf{y}),x_i-y_i⟩;$$

$$L_i{\nabla }_{x_i{x}_i}^{2}h(\mathbf{x})-{\nabla }_{x_i{x}_i}^{2}f(\mathbf{x})⪰0,$$

$${\mathrm{prox}}_{\gamma g}^h(x):=\underset{z\in {\mathbb{R}}^n}{\mathrm{argmin}}\left\{g(z)+\frac{1}{\gamma }{D}_h(z,x)\mid \right\}.$$

$${\mathrm{prox}}_{\gamma g}^h(x)=\left\{y\in \mathrm{dom}g\cap \mathrm{dom}h\text{leavevmode}|\text{leavevmode}g(y)+\frac{1}{\gamma }{D}_h(y,x)=\underset{z}{\min }\left\{g(z)+\frac{1}{\gamma }{D}_h(z,x)\mid \right\}<+\infty \mid \right\},$$

$$g^{\text{nicefrac}h{\gamma }_i}(\mathbf{x}):=\underset{z\in {\mathbb{R}}^{n_i}}{\min }\left\{g(\mathbf{x}+U_i(z-x_i))+\frac{1}{{\gamma }_i}{D}_h(\mathbf{x}+U_i(z-x_i),\mathbf{x})\mid \right\}>-\infty .$$

$${\gamma }_{i,g}^h:=\sup \left\{{\gamma }_i>0\text{leavevmode}|\text{leavevmode}\exists \mathbf{x}\in {\mathbb{R}}^n\text{leavevmode}s.t.\text{leavevmode}{g}^{\text{nicefrac}h{\gamma }_i}(\mathbf{x})>-\infty \mid \right\}.$$

$$g^{\text{nicefrac}h{\gamma }_i}(\mathbf{x})=\sum _{j\ne i}{g}_i(x_i)+\underset{z\in {\mathbb{R}}^{n_i}}{\min }\left\{g_i(z)+\frac{1}{{\gamma }_i}{D}_h(\mathbf{x}+U_i(z-x_i),\mathbf{x})\mid \right\},$$

$$g_i(z)$$

$$\ge g^{\text{nicefrac}h{\gamma }_i}(\mathbf{x})-\sum _{j\ne i}{g}_i(x_i)+\frac{r_i{\gamma }_i-1}{{\gamma }_i}{h}_{\mathbf{x}}^i(z)+\frac{1}{{\gamma }_i}(h(\mathbf{x})+⟨{\nabla }_{i}h(\mathbf{x}),z-x_i⟩)=:{\tilde{g}}_i(z).$$

$$g^{\text{nicefrac}h{\gamma }_i}(\mathbf{x})$$

$$\ge \sum _{j\ne i}{g}_i(x_i)+{\alpha }_i+\underset{z}{\inf }\left\{-r_i{h}_{\mathbf{x}}^i(z)+\frac{1}{{\gamma }_i}{D}_h(\mathbf{x}+U_i(z-x_i),\mathbf{x})\mid \right\}$$

$$\ge \sum _{j\ne i}{g}_i(x_i)+{\alpha }_i-\frac{1}{{\gamma }_i}h(\mathbf{x})+\frac{1}{{\gamma }_i}⟨{\nabla }_{i}h(\mathbf{x}),x_i⟩+\underset{z}{\inf }\left\{\frac{1-{\gamma }_i{r}_i}{{\gamma }_i}{h}_{\mathbf{x}}^i(z)-\frac{1}{{\gamma }_i}⟨{\nabla }_{i}h(\mathbf{x}),z⟩\mid \right\},$$

$$\underset{\Vert z\Vert \to \infty }{\mathrm{liminf}}\frac{g_i(z)}{h_{\mathbf{x}}^i(z)}\ge -r_i+\underset{\Vert z\Vert \to \infty }{\mathrm{liminf}}\frac{{\alpha }_i}{h_{\mathbf{x}}^i(z)}=-r_i>-\infty .$$

$$\underset{\Vert z\Vert \ge M_i}{\inf }{g}_i(z)+r_i{h}_{\mathbf{x}}^i(z)\ge \underset{\Vert x\Vert \ge M_i}{\inf }{h}_{\mathbf{x}}^i(z)({\ell }_i+r_i)>-\infty ,$$

$$M_{\text{nicefrac}h\gamma }(\mathbf{z},\mathbf{x}):=⟨\nabla f(\mathbf{x}),\mathbf{z}-\mathbf{x}⟩+\frac{1}{\gamma }{D}_h(\mathbf{z},\mathbf{x})+\sum _{i=1}^N{g}_i(z_i)$$

$$T_{\text{nicefrac}h{\gamma }_i}^i(\mathbf{x}):=\underset{z\in {\mathbb{R}}^{n_i}}{\mathrm{argmin}}{M}_{\text{nicefrac}h{\gamma }_i}(\mathbf{x}+U_i(z-x_i),\mathbf{x}),$$

$$\phi (\mathbf{x}+U_i(y_i-x_i))$$

$$\le f(\mathbf{x})+⟨{\nabla }_{i}f(\mathbf{x}),y_i-x_i⟩+\frac{1}{{\gamma }_i}{D}_h(\mathbf{x}+U_i(y_i-x_i),\mathbf{x})+g_i(y_i)+\sum _{j\ne i}{g}_j(x_j),$$

$$\phi (\mathbf{x}+U_i(\overline{z}-x_i))\le \phi (\mathbf{x})-\frac{1-{\gamma }_i{L}_i}{{\gamma }_i}{D}_h(\mathbf{x}+U_i(\overline{z}-x_i),\mathbf{x}),$$

$$\begin{array}{rl}{T}_{\text{nicefrac}h{\gamma }_i}^i(\mathbf{x})& =\underset{z\in {\mathbb{R}}^{n_i}}{\mathrm{argmin}}\left\{⟨\nabla f(\mathbf{x}),U_i(z-x_i)⟩+\frac{1}{{\gamma }_i}{D}_h(\mathbf{x}+U_i(z-x_i),\mathbf{x})+\sum _{i=1}^N{g}_i(z)\mid \right\}\\ & =\underset{z\in {\mathbb{R}}^{n_i}}{\mathrm{argmin}}\left\{⟨{\nabla }_{i}f(\mathbf{x}),z-x_i⟩+\frac{1}{{\gamma }_i}{D}_h(\mathbf{x}+U_i(z-x_i),\mathbf{x})+g_i(z)\mid \right\}.\end{array}$$

$$⟨{\nabla }_{i}f(\mathbf{x}),\overline{z}-x_i⟩+\frac{1}{{\gamma }_i}{D}_h(\mathbf{x}+U_i(\overline{z}-x_i),\mathbf{x})+g_i(\overline{z})\le g_i(x_i).$$

$$f(\mathbf{x}$$

$$\le f(\mathbf{x})+L_i{D}_h(\mathbf{x}+U_i(\overline{z}-x_i),\mathbf{x})+g_i(x_i)-g_i(\overline{z})-\frac{1}{{\gamma }_i}{L}_i{D}_h(\mathbf{x}+U_i(\overline{z}-x_i),\mathbf{x})$$

$$=f(\mathbf{x})+g_i(x_i)-g_i(\overline{z})-\frac{1-{\gamma }_i{L}_i}{{\gamma }_i}{D}_h(\mathbf{x}+U_i(\overline{z}-x_i),\mathbf{x}),$$

Full text

1Department of Electrical Engineering (ESAT-STADIUS) – KU Leuven, Kasteelpark Arenberg 10, 3001 Leuven, Belgium

2Department of Mathematics and Operational Research, University of Mons, Belgium

MA and PP acknowledge the support by the Research Foundation Flanders (FWO) research projects G086518N and G086318N; Research Council KU Leuven C1 project No. C14/18/068; Fonds de la Recherche Scientifique - FNRS and the Fonds Wetenschappelijk Onderzoek - Vlaanderen (FWO) under EOS project no 30468160 (SeLMA). NG also acknowledges the support by the European Research Council (ERC starting grant no 679515).

1. Introduction

Consider the structured nonsmooth nonconvex minimization problem

[TABLE]

where we will systematically assume the following hypotheses (see Section 2 for details):

Assumption I** (requirements for composite minimization (1.1)).**

1. a1

$g_{i}:\mathbb{R}^{n_{i}}\rightarrow\overline{\mathbb{R}}\coloneqq\mathbb{R}\cup{\mathopen{}\left\{\infty{}\mid{}{}\right\}\mathclose{}}$* is proper and lower semicontinuous (lsc);* 2. a2

$f:\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}$* is $\mathcal{C}^{1}(\operatorname{int}\operatorname{dom}h)$ and $(L_{1},\ldots,L_{N})$-smooth relative to $h$; here $n=\sum_{i=1}^{N}n_{i}$;* 3. a3

$h:\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}$* is multi-block strictly convex, $1$-coercive and essentially smooth;* 4. a4

$\varphi$* has a nonempty set of minimizers, i.e., $\operatorname*{arg\,min}\varphi\neq\emptyset$, and $\operatorname{dom}\varphi\subseteq\operatorname{int}\operatorname{dom}h$;* 5. a5

the first-order oracles of $f$, $g_{i}\ (i=1,\ldots,N)$, and $h$ are available.

Although, the problem (1.1) has a simple structure, it covers a broad range of optimization problems arising in signal and image processing, statistical and machine learning, control and system identification. Consequently, needless to say, there is a huge number of algorithmic studies around solving the optimization problems of the form (1.1). Among all of such methodologies, we are interested in the class of alternating minimization algorithms such as block coordinate descent [13, 16, 39, 45, 50, 51, 57, 58], block coordinate [29, 30, 38], and Gauss-Seidel methods [9, 17, 33], which assumes that all blocks are fixed except one and solves the corresponding auxiliary problem with respect to this block, update the latter block, and continue with the others. In particular, the proximal alternating minimization has received much attention in the last few years; see for example [4, 7, 8, 5, 6, 14]. Recently, the proximal alternating linearized minimization and its variation has been developed to handle (1.1); see for example [21, 48, 53].

Traditionally, the Lipschitz (Hölder) continuity of partial gradients of $f$ in (1.1) is a necessary tool for providing the convergence analysis of optimization algorithms; see, e.g., [21, 48]. It is, however, well-known that it is not the Lipschitz (Hölder) continuity of gradients playing a key role in such analysis, but one of its consequence: an upper estimation of $f$ including a Bregman distance called descent lemma ; cf. [11, 43]. This idea is central to convergence analysis of many optimization schemes requiring such an upper estimation; see, e.g., [2, 10, 11, 22, 55, 34, 35, 43]. In this paper, we propose a multi-block extension of the descent lemma given in [11, 43] and propose a Bregman proximal alternating linearized minimization (BPALM) algorithm and its adaptive version (A-BPALM) for (1.1).

1.1. Contribution

Our contribution is summarized as follows:

(Bregman proximal alternating linearized minimization) We introduce BPALM, a multi-block generalization of the proximal alternating linearized minimization (PALM) [21] using Bregman distances, and its adaptive version (A-BPALM). To do so, we extend the notion of relative smoothness [11, 43] to its multi-block counterpart to support a structured problem of the form (1.1). Owing to multi-block relative smoothness of $f$, unlike PALM, our algorithm does not need to know the local Lipschitz moduli of partial gradients ${\nabla}_{i}f$ ($i=1,\ldots,N$) and their lower and upper bounds, which are hard to provide in practice. 2. 2)

(Efficient framework for ONMF) Exploiting a suitable kernel for Bregman distance, it turns out that the objective of orthogonal nonnegative matrix factorization (ONMF) is multi-block relatively smooth, and the subproblems of our algorithms are solved in closed forms making them suitable for large-scale data analysis problems. To the best of our knowledge, BPALM and A-BPALM are the first algorithms with rigorous convergence theory for ONMF.

1.2. Related works

Closly related to our framework, there are two papers [40, 60]. However, we notice that [60] uses a sum separable kernel function which is just a special case of our multi-block kernel functions (see Example 2.2), and the paper only provides a limited convergence theory. Regarding [40], an algorithm (named B-PALM) proposed that is just a special case of our BPALM when $N=2$, $g_{1}=g_{2}=0$, and $f\in\mathcal{C}^{2}$, which restrict its applications. We stress that involving the block separable nonsmooth nonconvex functions $g_{i}$ and considering $N>2$ make our analysis different from those of [40].

1.3. Organization

This paper has four sections, besides this introductory section. In Section 2, we introduce the notion of multi-block relative smoothness, and verify the fundamental properties of Bregman proximal alternating linearized mapping. In Section 3, we introduce BPALM and A-BPALM and investigate their convergence analysis. In Section 4, we show that ONMF satisfies our assumptions, the related subproblems are solved in closed forms, and report our numerical results. Finally, Section 5 delivers some conclusions.

1.4. Notation

We denote by $\overline{\mathbb{R}}\coloneqq\mathbb{R}\cup{\mathopen{}\left\{\infty{}\mid{}{}\right\}\mathclose{}}$ the extended-real line. For the identity matrix $I_{n}$, we set $U_{i}\in\mathbb{R}^{n\times n_{i}}$ such that $I_{n}=(U_{1},\ldots,U_{N})\in\mathbb{R}^{n\times n}$. The open ball of radius $r\geq 0$ centered in $x\in\mathbb{R}^{p}$ is denoted as $\operatorname{B}(x;r)$. The set of cluster points of $(x^{k})_{k\in\mathbb{N}}$ is denoted as $\omega(x^{0})$. A function $f:\mathbb{R}^{p}\rightarrow\overline{\mathbb{R}}$ is proper if $f>-\infty$ and $f\not\equiv\infty$, in which case its domain is defined as the set $\operatorname{dom}f\coloneqq{\mathopen{}\left\{x\in\mathbb{R}^{p}{}\mid{}f(x)<\infty\right\}\mathclose{}}$. For $\alpha\in\mathbb{R}$, $[f\leq\alpha]{}\coloneqq{}{\mathopen{}\left\{x\in\mathbb{R}^{p}{}\mid{}f(x)\leq\alpha\right\}\mathclose{}}$ is the $\alpha$-(sub)level set of $f$; $[f\geq\alpha]$ and $[f=\alpha]$ are defined similarly. We say that $f$ is level bounded if $[f\leq\alpha]$ is bounded for all $\alpha\in\mathbb{R}$. A vector $v\in\partial f(x)$ is a subgradient of $f$ at $x$, and the set of all such vectors is called the subdifferential $\partial f(x)$ [52, Definition 8.3], i.e.

[TABLE]

2. Multi-block Bregman proximal alternating linearized mapping

We first establish the notion multi-block relative smoothness, which is an extension of the relative smoothness [11, 43] for problems of the form (1.1). We then introduce Bregman alternating linearized mapping and study some of its basic properties. For notation clarity, we will use bold lower-case letters (e.g., $\bm{x}$, $\bm{y}$, $\bm{z}$) for vectors in $\mathbb{R}^{\sum n_{i}}$ and use normal lower-case letters (e.g., $z$, $x_{i}$, $y_{i}$) for vectors in $\mathbb{R}^{n_{i}}$.

In order to extend the definition of Bregman distances for the multi-block problem (1.1), we first need to introduce the notion of multi-block kernel functions, which will coincide with the standard one (cf. [2, Definition 2.1]) if $N=1$.

Definition 2.1** (multi-block convexity and kernel function).**

Let $h:\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}$ be a proper and lsc function with $\operatorname{int}\operatorname{dom}h\neq\emptyset$ and such that $h\in\mathcal{C}^{1}(\operatorname{int}\operatorname{dom}h)$. For a fixed vector $\bm{x}\in\mathbb{R}^{n}$, we define the function $h_{\bm{x}}^{i}:\mathbb{R}^{n_{i}}\rightarrow\overline{\mathbb{R}}$ given by

[TABLE]

Then, we say that $h$ is

multi-block (strongly/strictly) convex* if the function $h_{\bm{x}}^{i}(\cdot)$ is (strongly/strictly) convex for all $\bm{x}\in\operatorname{dom}h$ and $i=1,\ldots,N$;* 2. 2)

multi-block locally strongly convex* around $\bm{x}^{\star}=(x_{1}^{\star},\ldots,x_{N}^{\star})$ if, for $i=1,\ldots,N$, there exists $\delta>0$ and $\sigma_{h}^{i}>0$ such that*

[TABLE] 3. 3)

a multi-block kernel function if $h$ is multi-block convex and $h_{\bm{x}}^{i}(\cdot)$ is $1$-coercive for all $\bm{x}\in\operatorname{dom}h$ and $i=1,\ldots,N$, i.e., $\lim_{\|z\|\to\infty}\tfrac{h_{\bm{x}}^{i}(z)}{\|z\|}=\infty$; 4. 4)

multi-block essentially smooth*, if for every sequence $(\bm{x}^{k})_{k\in\mathbb{N}}\subseteq\operatorname{int}\operatorname{dom}h$ converging to a boundary point of $\operatorname{dom}h$, we have $\|{\nabla}_{i}h(\bm{x}^{k})\|\to\infty$ for all $i=1,\ldots,N$;* 5. 5)

of multi block Legendre type if it is multi-block essentially smooth and multi-block strictly convex.

Example 2.2** (popular kernel functions).**

There are many kernel functions satisfying Definition 2.1. For example, for $N=1$, energy, Boltzmann-Shannon entropy, Fermi-Dirac entropy (cf. [12, Example 2.3]) and several examples in [43, Section 2]; and for $N=2$ see two examples in [40, Section 2]. Two important classes of multi-block kernels are sum separable kernels, i.e., $h(x_{1},\ldots,x_{N})=h_{1}(x_{1})+\ldots+h_{N}(x_{N})$, and product separable kernels, i.e., $h(x_{1},\ldots,x_{N})=h_{1}(x_{1})\times\ldots\times h_{N}(x_{N})$, see such a kernel for ONMF in Proposition 4.1. ∎

We now give the definition of Bregman distances (cf. [26]) for multi-block kernels.

Definition 2.3** (Bregman distance).**

For a kernel function $h$, the Bregman distance $\operatorname{D}_{h}:\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}$ is given by

[TABLE]

Fixing all blocks except the $i$-th one, the Bregman distance with respect to this block is given by

[TABLE]

which measures the proximity between $\bm{y}$ and $\bm{x}$ with respect to the $i$-th block of variables. Moreover, the kernel $h$ is multi-block convex if and only if $\operatorname{D}_{h}(\bm{x}+U_{i}(y_{i}-x_{i}),\bm{x})\geq 0$ for all $\bm{y}\in\operatorname{dom}h$ and $\bm{x}\in\operatorname{int}\operatorname{dom}h$ and $i=1,\ldots,N$. Note that if $h$ is multi-block strictly convex, then $\operatorname{D}_{h}(\bm{x}+U_{i}(y_{i}-x_{i}),\bm{x})=0$ ($i=1,\ldots,N$) if and only if $x_{i}=y_{i}$.

We are now in a position to present the notion of multi-block relative smoothness, which is the central tool for our analysis in Section 3.

Definition 2.4** (multi-block relative smoothness).**

Let $h:\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}$ be a multi-block kernel and let $f:\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}$ be a proper and lower semicontinuous function. If there exists $L_{i}\geq 0$ ($i=1,\ldots,N$) such that the functions $\phi_{\bm{x}}^{i}:\mathbb{R}^{n_{i}}\rightarrow\overline{\mathbb{R}}$ given by

[TABLE]

are convex for all $\bm{x}\in\operatorname{dom}h$ and $i=1,\ldots,N$, then, $f$ is called $(L_{1},\ldots,L_{N})$-smooth relative to $h$.

Note that if $N=1$, the multi-block relative smoothness is reduced to standard relative smoothness, which was introduced only recently in [11, 43]. In this case, if $f$ is $L$-Lipschitz continuous, then both $\nicefrac{{L}}{{2}}\|\cdot\|^{2}-f$ and $f-\nicefrac{{L}}{{2}}\|\cdot\|^{2}$ are convex, i.e., the relative smoothness of $f$ generalizes the notions of Lipschitz continuity using Bregman distances. If $N=2$, this definition will be reduced to the relative bi-smoothness given in [40] for $h,f\in\mathcal{C}^{2}$.

We next characterize the notion of multi-block relative smoothness.

Proposition 2.5** (characterization of multi-block relative smoothness).**

Let $h:\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}$ be a multi-block kernel and let $f:\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}$ be a proper lower semicontinuous function and $f\in\mathcal{C}^{1}$. Then, the following statements are equivalent:

1. (a)

$(L_{1},\ldots,L_{N})$-smooth relative to $h$; 2. (b)

for all $(\bm{x},\bm{y})\in\operatorname{int}\operatorname{dom}h\times\operatorname{int}\operatorname{dom}h$ and $i=1,\ldots,N$,

[TABLE] 3. (c)

for all $(\bm{x},\bm{y})\in\operatorname{int}\operatorname{dom}h\times\operatorname{int}\operatorname{dom}h$ and $i=1,\ldots,N$,

[TABLE] 4. (d)

if $f\in\mathcal{C}^{2}(\operatorname{int}\operatorname{dom}f)$ and $h\in\mathcal{C}^{2}(\mathrm{\bf int}\operatorname{dom}h)$ for all $\bm{x}\in\operatorname{int}\operatorname{dom}h$, then

[TABLE]

for $i=1,\ldots,N$.

{proof}

Fixing all the blocks except one of them, the results can be concluded in the same way as [43, Proposition 1.1].

2.1. Bregman proximal alternating linearized mapping

Recall that if $N=1$, for a kernel function $h:\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}$ and a proper lower semicontinuous function $g:\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}$, the Bregman proximal mapping is given by

[TABLE]

which is a generalization of the classical one using the Bregman distance (2.2) in place of the Euclidean distance; see, e.g., [27] and references therein. We note that

[TABLE]

which implies $\operatorname{dom}\operatorname{prox}_{\gamma g}^{h}\subset\mathbf{int}\operatorname{dom}h,\leavevmode\nobreak\ \operatorname{range}\operatorname{prox}_{\gamma g}^{h}\subset\operatorname{dom}g\cap\operatorname{dom}h$. The function $g$ is $h$-prox-bounded if there exists $\gamma>0$ such that $\min_{z}{\mathopen{}\left\{g(z)+\tfrac{1}{\gamma}\operatorname{D}_{h}(z,x){}\mid{}{}\right\}\mathclose{}}>-\infty$ for some $x\in\mathbb{R}^{n}$; cf. [2]. We next extend this definition to our multi-block setting.

Definition 2.6** (multi-block $h$-prox-boundedness).**

A function $g:\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}$ is multi-block $h$-prox-bounded if for each $i\in{\mathopen{}\left\{1,\ldots,N{}\mid{}{}\right\}\mathclose{}}$ there exists $\gamma_{i}>0$ and $\bm{x}\in\mathbb{R}^{n}$ such that

[TABLE]

The supremum of the set of all such $\gamma_{i}$ is the threshold $\gamma_{i,g}^{h}$ of the $h$-prox-boundedness, i.e.,

[TABLE]

For the problem (1.1), we have $g=\sum_{i=1}^{N}g_{i}$ leading to

[TABLE]

i.e., we therefore denote $\gamma_{g_{i}}^{h}:=\gamma_{i,g}^{h}$. If $g$ is multi-block $h$-prox-bounded for $\overline{\gamma}_{i}>0$, so is for all $\gamma_{i}\in{(0,\overline{\gamma}_{i})}$. We next present equivalent conditions to this notion.

Proposition 2.7** (characteristics of multi-block $h$-prox-boundedness).**

For a multi-block kernel function $h:\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}$ and proper and lsc functions $g_{i}:\mathbb{R}^{n_{i}}\rightarrow\overline{\mathbb{R}}$ with $i=1,\ldots,N$, the following statements are equivalent:

1. (a)

$g=\sum_{i=1}^{N}g_{i}$* is multi-block $h$-prox-bounded;* 2. (b)

for all $i=1,\ldots,N$ and $h_{\bm{x}}^{i}$ given in (2.1), $g_{i}+r_{i}h_{\bm{x}}^{i}$ is bounded below on $\mathbb{R}^{n_{i}}$ for some $r_{i}\in\mathbb{R}$; 3. (c)

for all $i=1,\ldots,N$, $\liminf_{\|z\|\to\infty}\nicefrac{{g_{i}(z)}}{{h_{\bm{x}}^{i}(z)}}>-\infty$.

{proof}

Suppose $g^{\nicefrac{{h}}{{\gamma_{i}}}}(\bm{x})>-\infty$ and let $r_{i}>\tfrac{1}{\gamma_{i}}$. Then, for all $i=1,\ldots,N$, it holds that

[TABLE]

Notice that $\tilde{g}_{i}$ is strictly convex and coercive, and as such is lower bounded. Conversely, suppose that $\alpha_{i}\coloneqq\inf g_{i}+r_{i}h_{\bm{x}}^{i}>-\infty$. Then, from (2.8), we obtain

[TABLE]

which is finite, owing to $1$-coercivity of $z\mapsto\tfrac{1-\gamma_{i}r_{i}}{\gamma_{i}}h_{\bm{x}}^{i}(z)-\tfrac{1}{\gamma_{i}}\langle{}{\nabla}_{i}h(\bm{x}){},{}z{}\rangle$.

Suppose that $\alpha_{i}\coloneqq\inf g_{i}+r_{i}h_{\bm{x}}^{i}>-\infty$. Since $h_{\bm{x}}^{i}(\cdot)$ is $1$-coercive, we have

[TABLE]

Conversely, suppose $\liminf_{\|z\|\to\infty}\nicefrac{{g_{i}(z)}}{{h_{\bm{x}}^{i}(z)}}>-\infty$. Then, there exists $\ell_{i},M_{i}\in\mathbb{R}$ such that $\nicefrac{{g_{i}(z)}}{{h_{\bm{x}}^{i}(z)}}\geq\ell_{i}$ whenever $\|z\|\geq M_{i}$. In particular

[TABLE]

where the last inequality follows from coercivity of $h_{\bm{x}}^{i}$. Since $\inf_{\|z\|\leq M_{i}}{g_{i}(z)+r_{i}h_{\bm{x}}^{i}(z)}>-\infty$ owing to lower semicontinuity, we conclude that $g_{i}+r_{i}h_{\bm{x}}^{i}$ is lower bounded on $\mathbb{R}^{n}$.

Let us now define the function $\operatorname{\mathcal{M}}_{\nicefrac{{h}}{{\gamma}}}:\mathbb{R}^{n}\times\mathbb{R}^{n}\rightarrow\overline{\mathbb{R}}$ as

[TABLE]

and the set-valued Bregman proximal alternating linearized mapping $\operatorname{T}_{\nicefrac{{h}}{{\gamma_{i}}}}^{i}:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n_{i}}$ as

[TABLE]

which reduces to the Bregman forward-backward splitting mapping if $N=1$; cf. [22, 2].

Remark 2.8** (majorization model).**

Note that invoking (b), the multi-block ($L_{1},\ldots,L_{N}$)-relative smoothness assumption of $f$ entails a majorization model

[TABLE]

for $\gamma_{i}\in(0,\nicefrac{{1}}{{L_{i}}})$. ∎

In the next lemma, we show that the cost function $\varphi$ is monotonically decreasing by minimizing the model (2.9) with respect to each block of variables.

Lemma 2.9** (Bregman proximal alternating inequality).**

Let the conditions in Assumption I hold, and let $\overline{z}\in\operatorname{T}_{\nicefrac{{h}}{{\gamma_{i}}}}^{i}(\bm{x})$ with $\gamma_{i}\in(0,\nicefrac{{1}}{{L_{i}}})$. Then,

[TABLE]

for all $i=1,\ldots,N$.

{proof}

For $i\in{\mathopen{}\left\{1,\ldots,N{}\mid{}{}\right\}\mathclose{}}$, (2.10) is simplified in the form

[TABLE]

Considering $\overline{z}\in\operatorname{T}_{\nicefrac{{h}}{{\gamma_{i}}}}^{i}(\bm{x})$, we have

[TABLE]

Since $f$ is $(L_{1},\ldots,L_{N})$-smooth relative to $h$, it follows from (b) for $\bm{x}$ and $y_{i}=\overline{z}$ that

[TABLE]

giving (2.11).

Recall that a function $\vartheta:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\overline{\mathbb{R}}$ with values $\vartheta(x,u)$ is level-bounded in $x$ locally uniformly in $u$ if for each $\bar{u}\in\mathbb{R}^{m}$ and $\alpha\in\mathbb{R}$ there is a neighborhood $\mathcal{U}$ of $\bar{u}$ along with a bounded set $B\subset\mathbb{R}^{n}$ such that $[\vartheta\leq\alpha]\subset B$ for all $u\in\mathcal{U}$, cf. [52]. Using this definition, the fundamental properties of the mapping $\operatorname{T}_{\nicefrac{{h}}{{\gamma_{i}}}}^{i}$ are investigated in the subsequent result.

Proposition 2.10** (properties of Bregman proximal alternating linearized mapping).**

Under conditions given in Assumption I for $i=1,\ldots,N$, the following statements are true:

$\operatorname{T}_{\nicefrac{{h}}{{\gamma_{i}}}}^{i}(\bm{x})$* is nonempty, compact, and outer semicontinuous (osc) for all $\bm{x}\in\operatorname{int}\operatorname{dom}h$;* 2. 2)

$\operatorname{dom}\operatorname{T}_{\nicefrac{{h}}{{\gamma_{i}}}}^{i}=\operatorname{int}\operatorname{dom}h$; 3. 3)

If $\overline{z}\in\operatorname{T}_{\nicefrac{{h}}{{\gamma_{i}}}}^{i}(\bm{x})$, then $\bm{x}+U_{i}(\overline{z}-x_{i})\subseteq\operatorname{dom}\varphi\subseteq\operatorname{int}\operatorname{dom}h$;

{proof}

For a fixed $\gamma_{i}^{0}\in(0,\gamma_{g_{i}}^{h})$ and a vector $\bm{x}\in\operatorname{int}\operatorname{dom}h$, let us define the function $\Phi_{i}:\mathbb{R}^{n_{i}}\times\mathbb{R}^{n}\times\mathbb{R}\rightarrow\overline{\mathbb{R}}$ given by

[TABLE]

Since $f$ and $g_{i}$ are proper and lsc, so is $\Phi_{i}$ on the set ${\mathopen{}\left\{(z,\bm{x},\gamma_{i}){}\mid{}\|z-x_{i}\|\leq\mu\gamma_{i},\leavevmode\nobreak\ 0\leq\gamma_{i}\leq\gamma_{i}^{0}\right\}\mathclose{}}$, for a constant $\mu>0$. We show that $\Phi_{i}$ is level-bounded in $z$ locally uniformly in $(\bm{x},\gamma)$. If it is not, then there exists $(\bm{x}^{k})_{k\in\mathbb{N}}\subset\operatorname{int}\operatorname{dom}h$, $(z^{k})_{k\in\mathbb{N}}$ with $\bm{x}^{k}+U_{i}(z^{k}-x_{i}^{k})\subset\operatorname{int}\operatorname{dom}h$, and $(\gamma_{i}^{k})_{k\in\mathbb{N}}\subset(0,\gamma_{i}^{0}]$ such that $\Phi_{i}(z^{k},\bm{x}^{k},\gamma_{i}^{k})\leq\beta<\infty$ with $(\bm{x}^{k},\gamma^{k})\to(\overline{\bm{x}},\overline{\gamma}_{i})$ and $\|z^{k}\|\to\infty$. This guarantees that, for sufficiently large $k$, $z^{k}\neq x_{i}^{k}$, i.e., $\gamma_{i}^{k}\in(0,\gamma_{0}]$ and

[TABLE]

Setting $\tilde{\gamma}_{i}\in(\gamma_{i}^{0},\gamma_{i}^{h})$, (b) ensures that there exists a constant $\tilde{\beta}\in\mathbb{R}$ such that

[TABLE]

Subtracting the last two inequalities, it holds that

[TABLE]

Expanding $\operatorname{D}_{h}(\bm{x}^{k}+U_{i}(z^{k}-x_{i}^{k}),\bm{x}^{k})$, dividing both sides by $\|z^{k}\|$, and taking limit from both sides of this inequality as $k\to\infty$, it can be deduced that

[TABLE]

This leads to the contradiction $+\infty\leq 0$, which implies that $\Phi_{i}$ is level-bounded. Therefore, all assumptions of the parametric minimization theorem [36, Theorem 2.2 and Corollary 2.2] are satisfied, i.e., Item 2). If $\overline{z}\in\operatorname{T}_{\nicefrac{{h}}{{\gamma_{i}}}}^{i}(\bm{x})$, then Lemma 2.9 implies that $\varphi(\bm{x}+U_{i}(\overline{z}-x_{i}))\leq\varphi(\bm{x})<\infty$ for $i=1,\ldots,N$, i.e., $\bm{x}+U_{i}(\overline{z}-x_{i})\subseteq\operatorname{dom}\varphi$, the second inclusion follows from a4.

Remark 2.11** (sum or product separable kernel).**

Let us observe the following.

If $h$ is an additive separable function, i.e., $h(x_{1},\ldots,x_{N})=h_{1}(x_{1})+\ldots+h_{N}(x_{N})$, then (2.10) can be written in the form

[TABLE] 2. 2)

If $h$ is product separable, i.e., $h(x_{1},\ldots,x_{N})=h_{1}(x_{1})\times\ldots\times h_{N}(x_{N})$, then

[TABLE]

where $\mu_{i}:=\nicefrac{{\gamma_{i}}}{{\eta_{\bm{x}}^{i}}}$ and $\eta_{\bm{x}}^{i}:=\prod_{j\neq i}h_{j}(x_{j})\neq 0$.

∎

3. Multi-block Bregman proximal alternating linearized minimization

We here introduce a multi-block proximal alternating linearized minimization algorithm and investigate its subsequential and global convergence, along with its convergence rate.

For a given point $\bm{x}^{k}=(x_{1}^{k},\ldots,x_{N}^{k})$, we set

[TABLE]

i.e., $\bm{x}^{k,0}=\bm{x}^{k}$ and $\bm{x}^{k,N}=\bm{x}^{k+1}$. Using this notation and (2.10), we next introduce the multi-block Bregman proximal alternating linearized minimization (BPALM) algorithm.

We note that each iteration of BPALM requires one call of the first-order oracle for the information needed in (3.1), and the iteration (3.1) are well-defined by Proposition 2.10. In addition, notice that if $N=1$, this algorithm reduces to the common (Bregman) proximal gradient (forward-backward) method [11, 15, 22]; if $N=2$, $h,f\in\mathcal{C}^{2}$, and $g_{1}=g_{2}=0$, then it reduces to B-PALM [40]; if $N=2$ and $h(\bm{x})=\tfrac{1}{2}(\|x_{1}\|^{2}+\|x_{2}\|^{2})$, it reduces to PALM [21]; if $h(\bm{x})=\tfrac{1}{2}\sum_{i=1}^{N}\|x_{i}\|^{2}$, then this algorithm is reduced to C-PALM [53].

We begin with showing some basic properties of the sequence generated by BPALM, involving a sufficient decrease condition.

Proposition 3.1** (sufficient decrease condition).**

Let Assumption I hold, and let $(\bm{x}^{k})_{k\in\mathbb{N}}$ be generated by BPALM. Then, the following statements are true:

the sequence $(\varphi(\bm{x}^{k}))_{k\in\mathbb{N}}$ is nonincreasing and

[TABLE]

where $\rho:=\min{\mathopen{}\left\{\tfrac{1-\gamma_{1}L_{1}}{\gamma_{1}},\ldots,\tfrac{1-\gamma_{N}L_{N}}{\gamma_{N}}{}\mid{}{}\right\}\mathclose{}}$; 2. 2)

we have

[TABLE]

i.e., $\lim_{k\to\infty}\operatorname{D}_{h}(\bm{x}^{k,i},\bm{x}^{k,i-1})=0$ for $i=1,\ldots,N$.

{proof}

Plugging $\overline{z}=x_{i}^{k,i}$ and $\bm{x}=\bm{x}^{k,i-1}$ into Lemma 2.9, it holds that

[TABLE]

Summing up both sides of (3.4) from $i=1$ to $N$, it follows that

[TABLE]

giving (3.2). Let us sum up both sides of (3.2) from $k=0$ to $q$:

[TABLE]

Taking the limit as $q\to+\infty$, (3.3) holds true. Together with $\operatorname{D}_{h}(\cdot,\cdot)\geq 0$, this proves the claim.

Let us consider the condition

[TABLE]

as a stopping criterion, for the accuracy parameter $\varepsilon>0$. Then, the first main consequence of Proposition 3.1 will provide us the iteration complexity of BPALM, which is the number of iterations needed for the stopping criterion (3.5) to be satisfied.

Corollary 3.2** (iteration complexity).**

Let Assumption I hold, and let $(\bm{x}^{k})_{k\in\mathbb{N}}$ be generated by BPALM, and let (3.5) be the stopping criterion. Then, BPALM will be terminated within $k\leq 1+\tfrac{\varphi(x^{0})-\inf\varphi}{\rho\varepsilon}$ iterations.

{proof}

Summing both sides of (3.2) over the first $\mathcal{K}\in\mathbb{N}$ iterations and telescoping the right hand side, it holds that

[TABLE]

Assuming that for all $(\mathcal{K}-1)$-th iterations the stopping criterion (3.5) is not satisfied, i.e., $\sum_{i=1}^{N}\operatorname{D}_{h}(\bm{x}^{k,i},\bm{x}^{k,i-1})>\varepsilon$, which leads to $\mathcal{K}\leq 1+\tfrac{\varphi(x^{0})-\inf\varphi}{\rho\varepsilon}$, giving the desired result.

In order to show the subsequential convergence of the sequence $(\bm{x}^{k})_{k\in\mathbb{N}}$ generated by BPALM, the next proposition will provide a lower bound for iterations gap $\|\bm{x}^{k+1}-\bm{x}^{k}\|$ using the subdifferential of $\partial\varphi(\bm{x}^{k+1})$.

Proposition 3.3** (subgradient lower bound for iterations gap).**

Let Assumption I hold, and let $(\bm{x}^{k})_{k\in\mathbb{N}}$ be generated by BPALM that we assume to be bounded. For a fixed $k\in\mathbb{N}$, we define

[TABLE]

Then, ${\mathopen{}\left(\mathcal{G}_{1}^{k+1},\ldots,\mathcal{G}_{N}^{k+1}\right)\mathclose{}}\in\partial\varphi(\bm{x}^{k+1})$ and

[TABLE]

with $\overline{c}:=\max{\mathopen{}\left\{\tfrac{\widetilde{L}+\gamma_{1}\widehat{L}}{\gamma_{1}},\ldots,\tfrac{\widetilde{L}+\gamma_{N-1}\widehat{L}}{\gamma{N-1}},\tfrac{\widetilde{L}_{N}+\gamma_{N}\widehat{L}_{N}}{\gamma_{N}}{}\mid{}{}\right\}\mathclose{}}$ in which $\widehat{L}$ and $\widetilde{L}>0$ are Lipschitz moduli of ${\nabla}_{i}f$, ${\nabla}_{i}h$ ($i=1,\ldots,N-1$) on bounded sets, and $\widehat{L}_{N}$ and $\widetilde{L}_{N}>0$ are Lipschitz moduli of ${\nabla}_{N}f$, ${\nabla}_{N}h$ on bounded sets, respectively.

{proof}

The optimality conditions for (3.1) ensures that there exists $q_{i}^{k,i}\in\partial g_{i}(x_{i}^{k,i})$ such that

[TABLE]

leading to

[TABLE]

On the other hand, owing to [5, Proposition 2.1], the subdifferential of $\varphi$ is given by

[TABLE]

i.e., for $\bm{x}=\bm{x}^{k+1}$,

[TABLE]

which means $(\mathcal{G}_{1}^{k+1},\ldots,\mathcal{G}_{N}^{k+1})\in\partial\varphi(\bm{x}^{k+1})$. It follows from the Lipschitz continuity of ${\nabla}_{i}f$, ${\nabla}_{i}h$ on bounded sets and the assumption of $(\bm{x}^{k})_{k\in\mathbb{N}}$ being bounded that there exist $\widehat{L}$, $\widehat{L}_{N}$, $\widetilde{L}>0$ and $\widetilde{L}_{N}>0$ such that

[TABLE]

for $i=1,\ldots,N-1$, and

[TABLE]

Invoking the last two inequalities, it can be concluded that

[TABLE]

as claimed.

Next, we proceed to derive the subsequential convergence of the sequence $(\bm{x}^{k})_{k\in\mathbb{N}}$ generated by BPALM: every cluster point of $(\bm{x}^{k})_{k\in\mathbb{N}}$ is a critical point of $\varphi$. Further, we explain some fundamental properties of the set of all cluster points $\omega(\bm{x}^{0})$ of the sequence $(\bm{x}^{k})_{k\in\mathbb{N}}$.

Theorem 3.4** (subsequential convergence and properties of $\omega(\bm{x}^{0})$).**

Let Assumption I hold, let the kernel $h$ be locally multi-block strongly convex, and let $(\bm{x}^{k})_{k\in\mathbb{N}}$ be generated by BPALM that we assume to be bounded. Then the following statements are true:

$\emptyset\neq\omega(\bm{x}^{0})\subset\mathbf{crit}\leavevmode\nobreak\ \varphi$; 2. 2)

$\lim_{k\to\infty}\operatorname{dist}{\mathopen{}\left(\bm{x}^{k},\omega(\bm{x}^{0})\right)\mathclose{}}=0$; 3. 3)

$\omega(\bm{x}^{0})$* is a nonempty, compact, and connected set;* 4. 4)

the objective function $\varphi$ is finite and constant on $\omega(\bm{x}^{0})$.

{proof}

For a limit point $\bm{x}^{\star}=(x_{1}^{\star},\ldots,x_{N}^{\star})$ of the sequence $(\bm{x}^{k})_{k\in\mathbb{N}}$, it follows from the boundedness of this sequence that there exists an infinite index set $\mathcal{J}\subset\mathbb{N}$ such that the subsequence $(\bm{x}^{k})_{k\in\mathcal{I}}$ converges to $\bm{x}^{\star}$ as $k\to\infty$. From the lower semicontinuity of $g_{i}$ ($i=1,\ldots,N$) and for $j\in\mathcal{J}$, it can be deduced that

[TABLE]

By (3.1), we get

[TABLE]

Using multi-block local strong convexity of $h$ around $\bm{x}^{\star}$ and invoking Item 2), there exist a neighborhood $\mathbf{B}(x_{i}^{\star},\varepsilon_{i}^{\star})$ for $\varepsilon_{i}^{\star}>0$, $\sigma_{i}^{\star}>0$, and $k_{i}^{0}\in\mathbb{N}$ such that for $k\geq k_{i}^{0}$ and $k\in\mathcal{J}$

[TABLE]

This indicates that the distance between two successive iterations goes to zero for large enough $k$. Since the sequence $(\bm{x}^{k})_{k\in\mathbb{N}}$ is bounded, ${\nabla}f$ and $h$ are continuous, substituting $k=k_{j}-1$ for $j\in\mathcal{J}$, taking the limit from both sides of the last inequality as $k\to\infty$, and (3.10), we come to

[TABLE]

and consequently,

[TABLE]

Further, Item 2) and Proposition 3.3 ensure ${\mathopen{}\left(\mathcal{G}_{1}^{k+1},\ldots,\mathcal{G}_{N}^{k+1}\right)\mathclose{}}\in\partial\varphi(\bm{x}^{k+1})$ and

[TABLE]

i.e., $\lim_{k\to\infty}{\mathopen{}\left(\mathcal{G}_{1}^{k+1},\ldots,\mathcal{G}_{N}^{k+1}\right)\mathclose{}}=(0_{n_{1}},\ldots,0_{n_{N}})$. Since the subdifferential mapping $\partial\varphi$ is closed, we have $(0_{n_{1}},\ldots,0_{n_{N}})\in\partial\varphi(x_{1}^{\star},\ldots,x_{N}^{\star})$, giving Item 1).

Item 2) is a direct consequence of Item 1), and Item 3) and Item 4) can be proved in the same way as [21, Lemma 5(iii)-(iv)].

3.1. Global convergence under Kurdyka-Łojasiewicz inequality

This section is devoted to the global convergence of BPALM under Kurdyka-Łojasiewicz inequality.

Definition 3.5** (KL property).**

A proper and lsc function $\varphi:\mathbb{R}^{n_{1}}\times\ldots\times\mathbb{R}^{n_{N}}\rightarrow\overline{\mathbb{R}}$ has the Kurdyka-Łojasiewicz property (KL property) at $\bm{x}^{\star}\in\operatorname{dom}\varphi$ if there exist a concave desingularizing function $\psi:[0,\eta]\to{[0,+\infty[}$ (for some $\eta>0$) and neighborhood $\operatorname{B}(\bm{x}^{\star};\varepsilon)$ with $\varepsilon>0$, such that

$\psi(0)=0$; 2. 2)

$\psi$* is of class $\mathcal{C}^{1}$ with $\psi>0$ on $(0,\eta)$;* 3. 3)

for all $\bm{x}\in\operatorname{B}(\bm{x}^{\star};\varepsilon)$ such that $\varphi(\bm{x}^{\star})<\varphi(\bm{x})<\varphi(\bm{x}^{\star})+\eta$ it holds that

[TABLE]

The set of all functions satisfying these conditions is denoted by $\Psi_{\eta}$.

The first inequality of this type is given in the seminal work of Łojasiewicz [41, 42] for analytic functions, which we nowadays call Łojasiewicz’s gradient inequality. Later, Kurdyka [37] showed that this in equality is valid for $\mathcal{C}^{1}$ functions whose graph belong to an $o$-minimal structure (see its definition in [59]). The first extensions of the KL property to nonsmooth functions was given by Bolte et al. [19, 18, 20].

The following two facts constitutes the crucial steps toward the establishment of the global convergence of the sequence generated by BPALM.

Fact 3.6** (uniformized KL property).**

[21, Lemma 6]** Let $\Omega$ be a compact set and $\zeta:\mathbb{R}^{d}\rightarrow\overline{\mathbb{R}}$ be a proper and lower semicontinuous functions. Assume that $\zeta$ is constant on $\Omega$ and satisfies the KL property at each point of $\Omega$. Then, there exists a $\varepsilon>0$, $\eta>0$, and $\psi\in\Psi_{\eta}$ such that for $\overline{u}$ and all $u$ in the intersection

[TABLE]

we have

[TABLE]

Fact 3.7**.**

[23*, Lemma 2.3]**

Let $(a_{k})_{k\in\mathbb{N}}$ and $(b_{k})_{k\in\mathbb{N}}$ be the sequences in $[0,+\infty)$ such that $\sum_{k=1}^{\infty}b_{k}<\infty$ and $a_{k+1}=\alpha a_{k}+b_{k}$ for all $k\in\mathbb{N}$ in which $\alpha<1$. Then, $\sum_{k=1}^{\infty}a_{k}<\infty$.*

Our subsequent main result indicates that the sequence $(\bm{x}^{k})_{k\in\mathbb{N}}$ generated by BPALM converges to a critical point $x^{\star}$ of $\varphi$ if it satisfies the KL property; cf. Definition 3.5.

Theorem 3.8** (global convergence).**

Let Assumption I hold, let the kernels $h$ be multi-block globally strongly convex with modulus $\sigma_{i}$ ($i=1,\ldots,N$), and let $(\bm{x}^{k})_{k\in\mathbb{N}}$ be generated by BPALM that we assume to be bounded. If $\psi$ is a KL function, then the following statements are true:

The sequence $(\bm{x}^{k})_{k\in\mathbb{N}}$ has finite length, i.e.,

[TABLE] 2. 2)

The sequence $(\bm{x}^{k})_{k\in\mathbb{N}}$ converges to a stationary point $\bm{x}^{\star}$ of $\varphi$.

{proof}

Let us define the sequence $(\mathcal{S}_{k})_{k\in\mathbb{N}}$ given by $\mathcal{S}_{k}:=\varphi(\bm{x}^{k})-\varphi^{\star}$, which is decreasing by Item 1), i.e., $(\mathcal{S}_{k})_{k\in\mathbb{N}}\to 0$. We now consider two cases: (i) there exists $\overline{k}\in\mathbb{N}$ such that $\mathcal{S}_{k}=0$; (ii) $\mathcal{S}_{k}>0$ for all $k\geq 1$.

In Case (i), invoking Item 1) implies that $\varphi(\bm{x}^{k})=\varphi^{\star}$ for all $k\geq\overline{k}$. It follows from Item 2) and multi-block strong convexity of $h$ that

[TABLE]

implying $\bm{x}^{k+1}=\bm{x}^{k}$ for all $k\geq\overline{k}$, which leads to Item 1).

In Case (ii), it holds that $\varphi(\bm{x}^{k})>\varphi^{\star}$ for all $k\geq 1$. From Item 3), the set of limit points $\omega(\bm{x}^{0})$ of $(\bm{x}^{k})_{k\in\mathbb{N}}$ is nonempty and compact and $\varphi$ is finite and constant on $\omega(\bm{x}^{0})$ due to Item 4). Moreover, the sequence $(\varphi(\bm{x}^{k}))_{k\in\mathbb{N}}$ is decreasing (Item 1)), i.e., for $\eta>0$, there exists a $k_{1}\in\mathbb{N}$ such that $\varphi^{\star}<\varphi(\bm{x}^{k})<\varphi^{\star}+\eta$ for all $k\geq k_{1}$. For $\varepsilon>0$, Item 2) implies that there exists $k_{2}\in\mathbb{N}$ such that $\operatorname{dist}(\bm{x}^{k},\omega(\bm{x}^{0}))<\varepsilon$ for $k\geq k_{2}$. Setting $k_{0}:=\max{\mathopen{}\left\{k_{1},k_{2}{}\mid{}{}\right\}\mathclose{}}$ and according to Fact 3.6, there exist $\varepsilon,\eta>0$ and a desingularization function $\psi$ such that for any element in

[TABLE]

the following inequality holds:

[TABLE]

Let us define $\Delta_{k}:=\psi(\varphi(\bm{x}^{k})-\varphi^{\star})=\psi(\mathcal{S}_{k})$. Then, it follows from the concavity of $\psi$ and Proposition 3.3 that

[TABLE]

with $\widehat{c}:=\nicefrac{{\overline{c}}}{{(\rho\min{\mathopen{}\left\{\sigma_{1},\ldots,\sigma_{N}{}\mid{}{}\right\}\mathclose{}})}}$. Using the arithmetic and quadratic means inequality, and applying the arithmetic and geometric means inequality, it can be concluded that

[TABLE]

We now define the sequences $(a_{k})_{k\in\mathbb{N}}$ and $(b_{k})_{k\in\mathbb{N}}$ as

[TABLE]

where $\sum_{i=1}^{\infty}b_{k}=\tfrac{\widehat{c}N}{2}\sum_{i=1}^{\infty}{\mathopen{}\left(\Delta_{i}-\Delta_{i+1}\right)\mathclose{}}=\Delta_{1}-\Delta_{\infty}=\Delta_{1}<\infty$. According to Fact 3.7, we infer $\sum_{k=1}^{\infty}a_{k}<\infty$, which proves Item 1).

By (3.12), the sequence $(\bm{x}^{k})_{k\in\mathbb{N}}$ is a Cauchy sequence, i.e., it converges to a stationary point $\bm{x}^{\star}$, giving the desired result.

Remark 3.9**.**

In Theorem 3.4 and Theorem 3.8, we implicitly assume that the sequence $(\bm{x}^{k})_{k\in\mathbb{N}}$ is bounded. This assumption is typical in convergence analysis of proximal-type algorithms for solving general non-convex non-smooth composite optimization problem, see e.g., [5, 22]. Proposition 3.1 shows that $\varphi(\bm{x}^{k})$ is non-increasing; hence, it is upper bounded by $\varphi(\bm{x}^{0})$. Therefore, the sequence $(\bm{x}^{k})_{k\in\mathbb{N}}$ would be bounded if $f(\cdot)$ has bounded level sets and $\sum_{i=1}^{n}g_{i}(x_{i})$ is bounded below. ∎

3.2. Convergence rate under Łojasiewicz-type inequality

We now investigate the convergence rate of the sequence generated by BPALM under KL inequality of Łojasiewicz type at $x^{\star}$ ($\psi(s):=\frac{\kappa}{1-\theta}s^{1-\theta}$ with $\theta\in[0,1)$), i.e., there exists $\varepsilon>0$ such that

[TABLE]

The following fact plays a key role in studying the convergence rate of the sequence generated by BPALM, where its proof can be found in [3, Lemma 1] and [24, Lemma 12].

Fact 3.10** (convergence rate of a sequence with positive elements).**

Let $(s_{k})_{k\in\mathbb{N}}$ be a sequence in $\mathbb{R}_{+}$ and let $\alpha$ and $\beta$ be some positive constants. Suppose that $s_{k}\to 0$ and that the sequence satisfies $s_{k}^{\alpha}\leq\beta(s_{k}-s_{k+1})$ for all $k$ sufficiently large. Then, the following assertions hold:

If $\alpha=0$, the sequences $(s_{k})_{k\in\mathbb{N}}$ converges to [math] in a finite number of steps; 2. 2)

If $\alpha\in(0,1]$, the sequences $(s_{k})_{k\in\mathbb{N}}$ converges linearly to [math] with rate $1-\nicefrac{{1}}{{\beta}}$, i.e., there exist $\lambda>0$ and $\tau\in[0,1)$ such that

[TABLE] 3. 3)

If $\alpha>1$, there exists $\mu>0$ such that for all $k$ sufficiently large

[TABLE]

We next derive the convergence rates of the sequences $(\bm{x}^{k})_{k\in\mathbb{N}}$ and $(\varphi(\bm{x}^{k}))_{k\in\mathbb{N}}$ under an additional assumption that the function $\varphi$ satisfies the KL inequality of Łojasiewicz type.

Theorem 3.11** (convergence rate).**

Let Assumption I hold, let the kernel $h$ be multi-block globally strongly convex with modulus $\sigma_{1},\ldots,\sigma_{N}$, and let the sequence $(\bm{x}^{k})_{k\in\mathbb{N}}$ generated by BPALM converges to $\bm{x}^{\star}$. If $\psi$ satisfies KL inequality of Łojasiewicz type (3.15), then the following assertions hold:

If $\theta=0$, then the sequences $(\bm{x}^{k})_{k\in\mathbb{N}}$ and $(\varphi(\bm{x}^{k}))_{k\in\mathbb{N}}$ converge in a finite number of steps to $\bm{x}^{\star}$ and $\varphi(\bm{x}^{\star})$, respectively; 2. 2)

If $\theta\in(0,\nicefrac{{1}}{{2}}]$, then there exist $\lambda_{1}>0$, $\mu_{1}>0$, $\tau,\overline{\tau}\in[0,1)$, and $\overline{k}\in\mathbb{N}$ such that

[TABLE] 3. 3)

If $\theta\in(\nicefrac{{1}}{{2}},1)$, then there exist $\lambda_{2}>0$, $\mu_{2}>0$, and $\overline{k}\in\mathbb{N}$ such that

[TABLE]

{proof}

The proof has two key parts.

In the first part, we show that there exists $\overline{k}\in\mathbb{N}$ such that for all $k\geq\overline{k}$ the following inequalities hold for $i=1,\ldots,N$:

[TABLE]

Let $\varepsilon>0$ be as described in (3.15) and $x^{k}\in\operatorname{B}(x^{\star};\varepsilon)$ for all $k\geq\tilde{k}$ and $\tilde{k}\in\mathbb{N}$. By the definitions of $a_{k}$ and $b_{k}$ in (3.14) and using (3.13), we get $a_{k+1}\leq\tfrac{1}{2}a_{k}+b_{k}$ for all $k\geq\tilde{k}$. Since $(\varphi)_{k\in\mathbb{N}}$ is nonincreasing,

[TABLE]

Together with the arithmetic and quadratic means inequality, $\psi(\mathcal{S}_{k})\leq\psi(\mathcal{S}_{k-1})$, and Item 1), this lead to

[TABLE]

On the other hand, for $i=1,\ldots,N$, we have

[TABLE]

This inequality, together with (3.17), yields

[TABLE]

leading to

[TABLE]

where $c:=\sqrt{\tfrac{2N}{\rho}}\max{\mathopen{}\left\{\nicefrac{{1}}{{\sqrt{\sigma_{1}}}},\ldots,\nicefrac{{1}}{{\sqrt{\sigma_{N}}}}{}\mid{}{}\right\}\mathclose{}}+\widehat{c}$ and $\psi(s):=\frac{\kappa}{1-\theta}s^{1-\theta}$. Let us consider the nonlinear equation

[TABLE]

which has a solution at $\mathcal{S}_{k-1}={\mathopen{}\left(\nicefrac{{\kappa}}{{1-\theta}}\right)\mathclose{}}^{\tfrac{2}{2\theta-2}}$. For $\hat{k}\in\mathbb{N}$ and $k\geq\hat{k}$, we assume that (3.18) holds and

[TABLE]

We now consider two cases: (a) $\theta\in[0,\nicefrac{{1}}{{2}}]$; (b) $\theta\in(\nicefrac{{1}}{{2}},1]$. In Case (a), if $\theta\in[0,\nicefrac{{1}}{{2}})$, then $\psi(\mathcal{S}_{k-1})\leq\sqrt{\mathcal{S}_{k-1}}$. If $\theta=\nicefrac{{1}}{{2}}$, then $\psi(\mathcal{S}_{k-1})=\tfrac{\kappa}{1-\theta}\sqrt{\mathcal{S}_{k-1}}$, i.e., $\max{\mathopen{}\left\{\sqrt{\mathcal{S}_{k-1}},\psi(\mathcal{S}_{k-1}){}\mid{}{}\right\}\mathclose{}}=\max{\mathopen{}\left\{1,\tfrac{\kappa}{1-\theta}{}\mid{}{}\right\}\mathclose{}}\sqrt{\mathcal{S}_{k-1}}$. Therefore, it holds that $\max{\mathopen{}\left\{\sqrt{\mathcal{S}_{k-1}},\psi(\mathcal{S}_{k-1}){}\mid{}{}\right\}\mathclose{}}\leq\max{\mathopen{}\left\{1,\tfrac{\kappa}{1-\theta}{}\mid{}{}\right\}\mathclose{}}\sqrt{\mathcal{S}_{k-1}}$. In Case (b), we have that

[TABLE]

i.e., $\max{\mathopen{}\left\{\sqrt{\mathcal{S}_{k-1}},\psi(\mathcal{S}_{k-1}){}\mid{}{}\right\}\mathclose{}}=\tfrac{\kappa}{1-\theta}\mathcal{S}_{k-1}^{1-\theta}$. Then, it follows from (3.18) that (3.16) holds for all $k\geq\overline{k}:=\max{\mathopen{}\left\{\tilde{k},\hat{k}{}\mid{}{}\right\}\mathclose{}}$.

In the second part of the proof, we will show the assertions in the statement of the theorem. For $(\mathcal{G}_{i}^{k},\ldots,\mathcal{G}_{N}^{k})\in\partial\varphi(\bm{x}^{k})$ as defined in Proposition 3.3, by Item 1), we infer

[TABLE]

with $\widehat{c}:=\frac{\rho}{2N\overline{c}^{2}\kappa^{2}}\min{\mathopen{}\left\{\sigma_{1},\ldots,\sigma_{N}{}\mid{}{}\right\}\mathclose{}}$ and for all $k\geq\overline{k}$. Hence, all assumptions of Fact 3.10 hold with $\alpha=2\theta$. Therefore, our results follows from this fact and (3.16).

3.3. Adaptive BPALM

The tightness of the $i$-th block upper estimation of the function $f$ given in (b) is dependent on the parameter $L_{i}>0$; however, in general, this parameter is a global information and it might not be tight locally, i.e., one may find a $L_{i}\geq\overline{L}_{i}(\bm{x})\geq 0$ such that

[TABLE]

for all $\bm{y}\in\operatorname{B}(\bm{x};\varepsilon_{1})$ with a small enough $\varepsilon_{1}>0$. Consequently, the majorization model described by $\operatorname{\mathcal{M}}_{\nicefrac{{h}}{{\gamma}}}$ may not be tight enough, which will consequently lead to smaller stepsizes $\gamma_{i}\in(0,\nicefrac{{1}}{{L_{i}}})$. In this case and in the case that $L_{1},\ldots,L_{N}$ are not available, one can retrieve them adaptively by applying a backtracking linesearch starting from a lower estimates; see, e.g., [1, 2, 44, 46, 56].

Putting together the above discussions, we propose an adaptive version of BPALM using a backtracking linesearch; see Algorithm 2.

We next provide an upper bound on the total number of calls of oracle after $k$ iterations of A-BPALM and those needed to satisfy (3.5).

Proposition 3.12** (worst-case oracle calls).**

Let $(\bm{x}^{k})_{k\in\mathbb{N}}$ be generated by A-BPALM. Then,

after at most $\tfrac{1}{\ln\nu_{1}}{\mathopen{}\left(\ln(\nu_{1}L_{i})-\ln\overline{L}_{i}^{0}\right)\mathclose{}}$ iterations the linesearch (Lines 4 to 7 of A-BPALM) is terminated; 2. 2)

the number of oracle call after $k$ full cycle $\mathcal{N}_{k}$ is bounded by

[TABLE] 3. 3)

the worst-case number of oracle calls to satisfy (3.5) is given by

[TABLE]

with $\overline{\rho}:=\min{\mathopen{}\left\{\nicefrac{{(1-\gamma_{1}^{0}\overline{L}_{1}^{0})}}{{\gamma_{1}^{0}}},\ldots,\nicefrac{{(1-\gamma_{N}^{0}\overline{L}_{N}^{0})}}{{\gamma_{N}^{0}}}{}\mid{}{}\right\}\mathclose{}}$.

{proof}

According to 7 and 10 of A-BPALM, we have $\overline{L}_{i}^{k+1}=\nu_{1}^{p_{i}^{k}}\overline{L}_{i}^{k}$, i.e.,

[TABLE]

giving Item 1). Hence, the total number of calls of oracle after $k$ iterations is given by

[TABLE]

giving Item 2).

Following the proof of Proposition 3.1 and since the sequence $(\nicefrac{{1-\gamma_{1}^{k}\overline{L}_{1}^{k}}}{{\gamma_{1}^{k}}})_{k\in\mathbb{N}}$ is increasing with respect to $k$, it is easy to see that

[TABLE]

On the other hand, 7 implies that $\nicefrac{{(1-\gamma_{i}^{k+1}\overline{L}_{i}^{k+1})}}{{\gamma_{i}^{k+1}}}\geq\nicefrac{{(1-\gamma_{i}^{0}\overline{L}_{i}^{0})}}{{\gamma_{i}^{0}}}$, $i=1,\ldots,N$, leading to

[TABLE]

Following the proof of Corollary 3.2, we have that BPALM will be terminated within $k\leq 1+\tfrac{\varphi(x^{0})-\inf\varphi}{\overline{\rho}\varepsilon}$ iterations. Together with Item 2), this implies that Item 3) is true.

Choosing appropriate constants $\overline{L}_{1}^{0},\ldots,\overline{L}_{N}^{0}$, Item 1) roughly speaking says that on average each full cycle of A-BPALM needs at most $2N$ oracle calls. Furthermore, in light of (3.19), Proposition 3.1 holds true by replacing $\rho$ with $\overline{\rho}$. Considering this replacement, all the results of Proposition 3.3, Theorem 3.4, Theorem 3.8, and Theorem 3.11 remain valid for A-BPALM.

Remark 3.13** (A-BPALM variant).**

We here notice that one may change Line 5 of A-BPALM as “set $\overline{L}_{i}^{k+1}=\nu_{1}^{p}\overline{L}_{i}^{0},\gamma_{i}^{k+1}=\nicefrac{{\gamma_{i}^{0}}}{{\nu_{1}^{p}}},p=p+1;$", which always start the backtracking procedure from $\overline{L}_{i}^{0}$ and $\gamma_{i}^{0}$. It is easy to see the results of Proposition 3.12 are still valid for this variant of A-BPALM. ∎

4. Application to orthogonal nonnegative matrix factorization

A natural way of analyzing large data sets is finding an effective way to represent them using dimensionality reduction methodologies. Nonnegative matrix factorization (NMF) is one such technique that has received much attention in the last few years; see, e.g., [28, 31, 32] and the references therein. In order to extract hidden and important features from data, NMF decomposes the data matrix into two factor matrices (usually much smaller than the original data matrix) by imposing componentwise nonnegativity and (possibly) sparsity constraints on these factor matrices. More precisely, let the data matrix be $X=[x_{1},x_{2},\ldots,x_{n}]\in\mathbb{R}_{+}^{m\times n}$ where each $x_{i}$ represents some data point. NMF seeks a decomposition of $X$ into a nonnegative $n\times r$ basis matrix $U=[u_{1},u_{2},\ldots,u_{r}]\in\mathbb{R}_{+}^{m\times r}$ and a nonnegative $r\times n$ coefficient matrix $V=[v_{1},v_{2},\ldots,v_{r}]^{T}\in\mathbb{R}_{+}^{r\times n}$ such that

[TABLE]

where $\mathbb{R}_{+}^{m\times n}$ is the set of $m\times n$ element-wise nonnegative matrices. Extensive research has been carried out on variants of NMF, and most studies in this area have focused on algorithmic developments, but with very limited convergence theory. This motivates us to study the application of BPALM and A-BPALM to a variant of NMF, namely orthogonal NMF (ONMF).

4.1. Orthogonal nonnegative matrix factorization

Besides the decomposition (4.1), the orthogonal nonnegative matrix factorization (ONMF) involves an additional orthogonality constraint $VV^{T}=I_{r}$ leading to the constrained optimization problem

[TABLE]

where $I_{r}\in\mathbb{R}^{r\times r}$ is the identity matrix. By imposing the matrix $V$ to be orthogonal (as well as nonnegative), ONMF imposes that each data points is only associated with one basis vector hence ONMF is closely related to clustering problems; see [49] and the references therein. Since the projection onto the set $C:={\mathopen{}\left\{(U,V)\in\mathbb{R}^{m\times r}\times\mathbb{R}^{r\times n}\mid U\geq 0,\leavevmode\nobreak\ V\geq 0,\leavevmode\nobreak\ VV^{T}=I_{r}{}\mid{}{}\right\}\mathclose{}}$ is costly, we here consider the penalized formulation

[TABLE]

for the penalty parameter $\lambda>0$. Introducing a product separable kernel, we next show that the objective function (4.3) is multi-block relatively smooth.

Proposition 4.1** (multi-block relative smoothness of ONMF objective).**

Let the function $h:\mathbb{R}^{m\times r}\times\mathbb{R}^{r\times n}\rightarrow\overline{\mathbb{R}}$ be a kernel given by

[TABLE]

Then the function $f:\mathbb{R}^{m\times r}\times\mathbb{R}^{r\times n}\rightarrow\overline{\mathbb{R}}$ given by $f(U,V):=\tfrac{1}{2}\|X-UV\|_{F}^{2}+\tfrac{\lambda}{2}\|I_{r}-VV^{T}\|_{F}^{2}$ is $(L_{1},L_{2})$-smooth relative to $h$ with

[TABLE]

{proof}

Using partial derivatives ${\nabla}_{U}f(U,V)=UVV^{T}-XV^{T}$, ${\nabla}_{UU}^{2}f(U,V)Z=ZVV^{T}$, and the Cauchy Schwarz inequality, it can be concluded that $\langle{}Z{},{}{\nabla}_{UU}^{2}f(U,V)Z{}\rangle\leq\|V\|_{F}^{2}\|Z\|_{F}^{2}$. On the other hand, ${\nabla}_{U}h(U,V)=\beta_{1}{\mathopen{}\left(\tfrac{\alpha_{2}}{4}\|V\|_{F}^{4}+\tfrac{\beta_{2}}{2}\|V\|_{F}^{2}+1\right)\mathclose{}}U$ and

[TABLE]

Together with (4.5), this yields

[TABLE]

which implies $L_{1}{\nabla}_{UU}^{2}h(U,V)-{\nabla}_{UU}^{2}f(U,V)\succeq 0$.

From ${\nabla}_{V}f(U,\cdot)(V)=U^{T}UV-U^{T}X+2\lambda(VV^{T}V-V)$ and the definition of directional derivative, we obtain

[TABLE]

This, $\langle{}Y_{1}{},{}Y_{2}{}\rangle:=\operatorname{tr}(Y_{1}^{T}Y_{2})$, basic properties of the trace, the Cauchy-Schwarz inequality, and the submultiplicative property of the Frobenius norm imply

[TABLE]

Plugging ${\nabla}_{V}h(U,V)={\mathopen{}\left(\tfrac{\beta_{1}}{2}\|U\|_{F}^{2}+1\right)\mathclose{}}{\mathopen{}\left(\alpha_{2}\|V\|_{F}^{2}+\beta_{2}\right)\mathclose{}}V$ into the directional derivative definition, we come to

[TABLE]

implying

[TABLE]

Hence, it follows from (4.5) that

[TABLE]

i.e., $L_{2}{\nabla}_{VV}^{2}h(U,V)-{\nabla}_{VV}^{2}f(U,V)\succeq 0$, as claimed.

The unconstrained version of the ONMF problem (4.2) is given by

[TABLE]

where $\delta_{U\geq 0}$ and $\delta_{V\geq 0}$ are the indicator functions of the sets $C_{1}:={\mathopen{}\left\{U\in\mathbb{R}^{m\times r}\mid U\geq 0{}\mid{}{}\right\}\mathclose{}}$ and $C_{2}:={\mathopen{}\left\{V\in\mathbb{R}^{r\times n}\mid V\geq 0{}\mid{}{}\right\}\mathclose{}}$, respectively. Comparing to (1.1), the next setting is recognized

[TABLE]

in which both $g_{1}(U)$ and $g_{2}(V)$ are nonsmooth and convex, and $f(U,V)$ is $(L_{1},L_{2})$-smooth relative to $h$ given in (4.4); cf. Proposition 4.1. For given $U^{k}$ and $V^{k}$, applying BPALM and A-BPALM to (4.6), $U^{k+1}$ and $V^{k+1}$ should be computed efficiently, which we study next.

Theorem 4.2** (closed-form solutions of the subproblem (3.1) for ONMF).**

Let $h_{1}:\mathbb{R}^{m\times r}\rightarrow\overline{\mathbb{R}}$ and $h_{2}:\mathbb{R}^{r\times n}\rightarrow\overline{\mathbb{R}}$ be the kernel functions given by

[TABLE]

i.e., $h(U,V)=h_{1}(U)h_{2}(V)$. For given $U^{k}$ and $V^{k}$, the problem (4.6), and the subproblem (3.1), the following assertions hold:

For $\eta_{1}=\tfrac{\alpha_{2}}{4}\|V^{k}\|_{F}^{4}+\tfrac{\beta_{2}}{2}\|V^{k}\|_{F}^{2}+1$ and $\mu_{1}:=\nicefrac{{\gamma_{1}}}{{(\beta_{1}\eta_{1})}}$, the iteration $U^{k+1}$ is given by

[TABLE] 2. 2)

For $\eta_{2}=\tfrac{\beta_{1}}{2}\|U^{k+1}\|_{F}^{2}+1$ and $\mu_{2}:=\nicefrac{{\gamma_{2}}}{{\eta_{2}}}$, the iteration $V^{k+1}$ is given by

[TABLE]

with ${\nabla}_{V}f(U^{k+1},V^{k})=(U^{k+1})^{T}U^{k+1}V^{k}-(U^{k+1})^{T}X+2\lambda(V^{k}(V^{k})^{T}V^{k}-V^{k})$ and

[TABLE]

where $\tau_{1}:=-\nicefrac{{\beta_{2}^{2}}}{{3}}$ and \tau_{2}:=\nicefrac{{{\mathopen{}\left(-2\beta_{2}^{3}-27\alpha_{2}\big{\|}\max{\mathopen{}\left\{(\alpha_{2}\|V^{k}\|_{F}^{2}+\beta_{2})V^{k}-\mu_{2}{\nabla}_{V}f(U^{k+1},V^{k}),0{}\mid{}{}\right\}\mathclose{}}\big{\|}_{F}\right)\mathclose{}}}}{{27}}.

{proof}

Setting $g_{1}:=\delta_{U\geq 0}$ and $f(U,V)=\tfrac{1}{2}\|X-UV\|_{F}^{2}+\tfrac{\lambda}{2}\|I_{r}-VV^{T}\|_{F}^{2}$, it follows from (2.12) that

[TABLE]

with ${\nabla}_{U}f(U^{k},V^{k})=U^{k}V^{k}(V^{k})^{T}-X(V^{k})^{T}$, giving (4.7).

By setting $g_{2}:=\delta_{V\geq 0}$ and invoking (2.12), we infer

[TABLE]

Let us consider the normal cone $\mathcal{N}_{V\geq 0}(V^{k+1})={\mathopen{}\left\{P\in\mathbb{R}^{r\times n}\mid V^{k+1}\odot P=0,\ P\leq 0{}\mid{}{}\right\}\mathclose{}}$ (see [54, Corollary 3.5]), where $V\odot P$ denotes the Hadamard products given pointwise by $(V\odot P)_{ij}:=V_{ij}P_{ij}$ for $i\in{1,\ldots,r}$ and $j\in{1,\ldots,n}$. The first-order optimality conditions for the latter identity leads to $G^{k}-(\alpha_{2}\|V^{k+1}\|_{F}^{2}+\beta_{2})V^{k+1}\in\mathcal{N}_{V\geq 0}(V^{k+1})$ with $G^{k}:={\nabla}h_{2}(V^{k})-\mu_{2}{\nabla}_{V}f(U^{k+1},V^{k})$. Let us consider two cases: (i) $G_{ij}\leq 0$; (ii) $G_{ij}>0$. In Case (i), $P_{ij}=G_{ij}^{k}-(\alpha_{2}\|V^{k+1}\|_{F}^{2}+\beta_{2})V_{ij}^{k+1}\leq 0$, i.e., $V_{ij}^{k+1}=0$. In Case (ii), if $V_{ij}^{k+1}=0$, then $P_{ij}=G_{ij}^{k}>0$, which contradicts $P\leq 0$, i.e., $G_{ij}^{k}-(\alpha_{2}\|V^{k+1}\|_{F}^{2}+\beta_{2})V_{ij}^{k+1}=0$. Combining both cases, we come to the equation

[TABLE]

i.e., there exists $t_{k}\in\mathbb{R}$ such that $t_{k}V^{k+1}=\operatorname{Proj}_{V\geq 0}(G^{k})$ that eventually lead to

[TABLE]

which is a Cardano equation and its solution is given by (4.9).

4.2. Preliminary numerical experiment

In this section, we report preliminary numerical experiments with BPALM and two variants of A-BPALM, namely,

A-BPALM1: the algorithm A-BPALM; 2. 2)

A-BPALM2: the variant of A-BPALM as described in Remark 3.13.

Since the unconstrained ONMF problem (4.6) involves the penalty term $\tfrac{\lambda}{2}\|I_{r}-VV^{T}\|_{F}^{2}$, we also consider a “continuation" variant of these algorithm that starts from some $\lambda>0$, run one of the above-mentioned algorithms until some stopping criterion holds and save its best point, and then it increases the penalty parameter and run the algorithm with the starting point as the best point of the last call, and it continues the procedure until we stop the algorithm. We refer to this procedure as continuation, which we will describe next in more details.

In our Implementation all the codes were written in MATLAB (publicly available at https://github.com/MasoudAhoo/BPALM) and runs were performed on a MacBook Pro with 2,8 GHz Intel Core i7 CPU and 16 GB RAM. On the basis of our preliminary experiments, we here set $\alpha_{2}=\beta_{2}=\beta_{1}=1$ to provide the relative smoothness constants as described in (4.5), and the related step-sizes are computed by $\gamma_{i}=\nicefrac{{1}}{{L_{i}}}-\epsilon$, when $\epsilon$ is set as the machine precision. For A-BPALM1 and A-BPALM1, we set $\nu=2$, and we also set $\overline{L}_{i}^{0}=0.01L_{i}$ for A-BPALM1 and $\overline{L}_{i}^{0}=0.1L_{i}$ for A-BPALM2. For the continuation version, we set $c=\nicefrac{{3}}{{2}}$.

We first report the experiment on a synthetic data set with $(m,n,r)=(200,2000,10)$. Our synthetic data set is generated as follows. We use the MATLAB command $\mathsf{rand}$ to generate random nonnegative matrices $U\in\mathbb{R}_{+}^{m\times r}$ and $R\in\mathbb{R}_{+}^{m\times n}$, then we generate a random orthogonal nonnegative matrix $V\in\mathbb{R}_{+}^{r\times n}$. Next, we set $X=UV$ to obtain the $m$-by-$n$ orthogonal decomposable matrix $X$, and finally add 5% of noise by $X=X+0.05\frac{\|X\|_{F}}{\|R\|_{F}}R$. Now, we use SVD-based initialization for providing starting points for our algorithms, see [25]. We here run our algorithms with both fixed penalty parameter and with the continuation scheme. For fixed penalty versions, we set $\lambda=10$, and for continuation versions we started with $\lambda=10$ and stopped the inner algorithms every 3 seconds and increased $\lambda$ by factor $c=\nicefrac{{3}}{{2}}$. We stopped the algorithms after 15 seconds of the running time.

The results of our implementation are illustrated in Figure 1. In this figure, Subfigure (a) stands for fixed penalty versions while Subfigure (b) stands for continuation versions. Hence, on Subfigure (b), the penalty $\lambda$ is progressively increased. We make two observations: (i) In both cases A-BPALM1 and A-BPALM2 outperform BPALM while A-BPALM1 is the best among them; (ii) The continuation schemes perform much better than the fixed penalty versions, especially for A-BPALM1. In fact, although the curve on Subfigure (b) corresponds to a larger value of $\lambda$, A-BPALM1 achieves a much lower function value; namely around 20 on Subfigure (a) vs. 0.2 on Subfigure (b). The reason is that increasing $\lambda$ leads to a better solution where the factor $V$ is closer to orthogonality hence closer to the ground truth.

We next report the performance of our algorithms on the Hubble telescope data set which is taken from [47]. Since the continuation versions of our algorithms perform better, we here only apply the continuation versions of BPALM, A-BPALM1, and A-BPALM2. We use the SVD-based initialization as in [49]. In this problem, each row of the matrix $X$ is a vectorized image of the Hubble telescope at a given wavelength for a total of $m=100$ wavelengths. Each image contains $n=128\times 128$ pixels. Since each pixel in the image contains mostly a single material, it makes sense to use ONMF to cluster the pixel according to the material they contain (see Figure 2 for an illustration). For this application problem, we report the final relative fidelity and orthogonal errors, i.e.,

[TABLE]

with respect to several initial values for the penalty parameter $\lambda$ in the continuation procedure Algorithm 3. The results of our implementations are reported in Table 1 and the final outputs of the algorithms, along with the ground true Hubble image, are illustrated in Figure 2.

From Table 1, we observe that A-BPALM1 attains the better error than BPALM and A-BPALM2 in the sense of both the relative fidelity and orthogonal errors. Further, we observe that the orthogonal errors $O_{error}$ produced by the algorithms are decreasing by increasing the initial penalty parameter $\lambda$. From Figure 2, A-BPALM1 provides slightly better quality image compared to BPALM and A-BPALM2 (look for example at the first basis image).

5. Conclusion

We have analysed two new alternating linearized minimization algorithms called BPALM and A-BPALM for solving the popular nonconvex nonsmooth optimization problem (1.1). Convergence analysis including the subsequential convergence, the global convergence and the convergence rate of the proposed algorithms is studied under the framework of multi-block relative smoothness and multi-block kernel functions. We emphasize that, to the best of our knowledge, BPALM and A-BPALM are the first algorithms with rigorous convergence guarantee for solving ONMF in the literature. We employ BPALM and A-BPALM to solve the orthogonal nonnegative matrix factorization problem. Some preliminary numerical tests are provided to illustrate the performance of our algorithms. A comprehensive numerical experiments with several data sets and comparison with state-of-the-art algorithms are out of the scope of the current paper, which we aim for future work.