A majorization–minimization algorithm for split feasibility problems

Abstract

The classical multi-set split feasibility problem seeks a point in the intersection of finitely many closed convex domain constraints, whose image under a linear mapping also lies in the intersection of finitely many closed convex range constraints. Split feasibility generalizes important inverse problems including convex feasibility, linear complementarity, and regression with constraint sets. When a feasible point does not exist, solution methods that proceed by minimizing a proximity function can be used to obtain optimal approximate solutions to the problem. We present an extension of the proximity function approach that generalizes the linear split feasibility problem to allow for non-linear mappings. Our algorithm is based on the principle of majorization–minimization, is amenable to quasi-Newton acceleration, and comes complete with convergence guarantees under mild assumptions. Furthermore, we show that the Euclidean norm appearing in the proximity function of the non-linear split feasibility problem can be replaced by arbitrary Bregman divergences. We explore several examples illustrating the merits of non-linear formulations over the linear case, with a focus on optimization for intensity-modulated radiation therapy.

Appendix

1.1 Proof of Proposition 3

Our proof of Proposition 3 relies on showing the algorithm map $\psi$ that carries out the Armijo backtracking line search is closed at all non-stationary points of f. Recall that a point-to-set map A is closed at a point $\varvec{x}$, if for any sequence $\varvec{x}_j \rightarrow \varvec{x}$ such that $\varvec{y}_j \in A(\varvec{x}_j)$ converges to $\varvec{y}$, it follows that $\varvec{y}\in A(\varvec{x})$. Define the following point-to-set map

$$\begin{aligned} S(\varvec{x}, \varvec{v}) = \left\{ \varvec{y}: \varvec{y}= \varvec{x}+ \sigma ^k \varvec{v}\text { for }k \in \mathbb {Z}_+\text { such that }\, f(\varvec{y}) \le f(\varvec{x}) + \alpha \sigma ^{k} df(\varvec{x}) \varvec{v}\right\} , \end{aligned}$$

and let G denote the map

$$\begin{aligned} G(\varvec{x}) = \begin{pmatrix} \varvec{x}\\ \nabla f(\varvec{x}) \end{pmatrix}. \end{aligned}$$

The map G is continuous since f is Lipschitz differentiable. Let $A = S \circ G(\varvec{x})$. By Corollary 2 in Chapter 7.7 of [49], the composite mapping A is closed at $\varvec{x}$, if the mapping S is closed at $G(\varvec{x})$. We will state and prove a slightly more general result.

Proposition 4

The mapping S is closed at all points $(\varvec{x}, \varvec{v})$ provided that $\varvec{v}\not = \varvec{0}$.

Proof

Consider sequences $\{ \varvec{x}_j \}$, $\{ \varvec{v}_j \}$ such that $\varvec{x}_j \rightarrow \varvec{x}, \varvec{v}_j \rightarrow \varvec{v}$, and let $\varvec{y}_k \in S(\varvec{x}_{k}, \varvec{v}_k)$ with $\varvec{y}_k \rightarrow \varvec{y}$. For every j, there is some $k_j$ such that $\varvec{y}_j = \varvec{x}_j + \sigma ^{k_j} \varvec{v}_j$, where

$$\begin{aligned} \sigma ^{k_j} = \log \left( \frac{ ||\varvec{y}_j - \varvec{x}_j ||}{||\varvec{v}_j ||} \right) . \end{aligned}$$

Taking limits on both side above yields

$$\begin{aligned} \sigma ^{\bar{k}} = \frac{ ||\varvec{y}- \varvec{x}||}{||\varvec{v}||}, \text { where }\bar{k} = \lim _{j \rightarrow \infty } k_j. \end{aligned}$$

Because $\varvec{y}_k \in S(\varvec{x}_{k}, \varvec{v}_k)$, it holds for each j that

$$\begin{aligned} f(\varvec{y}_j) \le f(\varvec{x}_j) + \alpha \sigma ^{k_j} df(\varvec{x}_j ) \varvec{v}_j . \end{aligned}$$

(22)

Since $\{k_j\}$ is a sequence of integers, it assumes only finitely many distinct values before converging to the constant sequence $\{ \bar{k} \}$; let $k^\dagger$ denote the maximum of these values. Then replacing $k_j$ by $k^\dagger$ and letting $j \rightarrow \infty$ in (22), together with continuity of f and $df(\varvec{x})$, we have that

$$\begin{aligned} f(\varvec{y}) \le f(\varvec{x}) + \alpha \sigma ^{k^\dagger } df(\varvec{x}) \varvec{v}. \end{aligned}$$

That is, $\varvec{y}\in S(\varvec{x}, \varvec{v})$, proving that S is closed at $(\varvec{x}, \varvec{v})$. $\square$

Proposition 4 informs us that the map A is closed at all non-stationary points of f. We are now ready to prove Proposition 3.

Proof

Fix an initial iterate $\varvec{x}_0$. Note that the set $\mathcal {L}_f(\varvec{x}_0) \equiv \{ \varvec{x}: f(\varvec{x}) \le f(\varvec{x}_0)\}$ is compact since f is coercive and the modified quasi-Newton method generates monotonically decreasing values. Since $\mathcal {L}_f(\varvec{x}_0)$ is compact, the sequence $\{\varvec{x}_{k}\}$ has a convergent subsequence whose limit is in $\mathcal {L}_f(\varvec{x}_0)$; denote this as $\varvec{x}_{k_l} \rightarrow \varvec{x}_\star$ as $l \rightarrow \infty$. Our goal is to show that $\varvec{x}_\star$ must be a stationary point of f. To the contrary, suppose that $\varvec{x}_\star$ is not a stationary point of f.

Let $\varvec{y}_{k_l} = \psi (\varvec{x}_{k_l}) \in A(\varvec{x}_{k_l})$. Note that $\varvec{y}_{k_l} \in \mathcal {L}_f(\varvec{x}_0)$ and therefore the sequence $\{\varvec{y}_{k_l}\}$ has a convergent subsequence $\{\varvec{y}_{k_{l_j}}\}$. Denote this sequence’s limit $\varvec{y}_\star$. Note that the map A is closed at $\varvec{x}_\star$, since $\varvec{x}_\star$ is not a stationary point of f. Therefore, by the definition of closed maps, we have that $\varvec{y}_\star \in A(\varvec{x}_\star )$ and consequently

$$\begin{aligned} f(\varvec{y}_\star ) \le f(\varvec{x}_\star ) + \alpha \sigma ^k df(\varvec{x}_\star ) \varvec{H}(\varvec{x}_\star )^{-1}\nabla f(\varvec{x}_\star ) < f(\varvec{x}_\star ), \end{aligned}$$

(23)

for some positive integer k. On the other hand, since f is Lipschitz-differentiable, it is continuous; therefore $\lim f(\varvec{x}_{k_l}) = f(\varvec{x}_\star )$. Moreover, for all $k_l$ we have that

$$\begin{aligned} f(\varvec{x}_\star ) \le f(\varvec{x}_{k_{l+1}}) \le f(\varvec{x}_{k_l + 1}) = f(\varphi (\varvec{x}_{k_l})) \le f(\psi \circ \psi (\varvec{x}_{k_l})) \le f(\psi (\varvec{x}_{k_l})). \end{aligned}$$

By continuity of f, we have that $f(\varvec{x}_\star ) \le f(\varvec{y}_\star )$, contradicting the inequality established in (23). We conclude that $\varvec{x}_\star$ must be a stationary point of f. $\square$

1.2 $\beta$-divergence

The $\beta$-divergence is defined

$$\begin{aligned} D(\varvec{x},\varvec{y}) = {\left\{ \begin{array}{ll} \displaystyle \sum _{j}\dfrac{1}{\beta (\beta -1)} x_j^\beta + \sum _{j}\dfrac{1}{\beta } y_j^{\beta } - \dfrac{1}{\beta -1} x_j y_j^{\beta -1} &{} \beta \in \mathbb {R}\setminus \{0, 1 \} \\ \sum _j x_j \log \left( \dfrac{x_j}{y_j}\right) - x_j + y_j &{} \beta = 1 \\ \sum _j \dfrac{x_j}{y_j} - \log \left( \dfrac{x_j}{y_j}\right) - 1 &{} \beta = 0 \end{array}\right. } \end{aligned}$$

We see that the $\beta$-divergence corresponds to the Kullback–Leibler and Itakura–Saito divergences when $\beta =1,0$ respectively. Below we discuss the projection onto a hyperplane for the case of $\beta \in \mathbb {R}\setminus \{0, 1 \}$.

The function $\phi (x) = \dfrac{1}{\beta (\beta -1)} \varvec{x}^\beta$ generates the $\beta$-divergence

$$\begin{aligned} {D(\varvec{x},\varvec{y}) = \sum _{j}\dfrac{1}{\beta (\beta -1)} x_j^\beta + \sum _{j}\dfrac{1}{\beta } y_j^{\beta } - \dfrac{1}{\beta -1} x_j y_j^{\beta -1} .} \end{aligned}$$

Its gradient is

$$\begin{aligned} {\nabla _\beta \phi (\varvec{x}) = \dfrac{1}{\beta -1}\varvec{x}^{\beta -1}, } \end{aligned}$$

and recall the Fenchel conjugate of $\phi$ is given by

$$\begin{aligned} {\phi ^\star (\varvec{z}) = \sup _{\varvec{x}} \bigg ( \langle \varvec{z}, \varvec{x}\rangle - \phi (\varvec{x}) \bigg ) = \sup _{\varvec{x}} \bigg ( \langle \varvec{z}, \varvec{x}\rangle - \sum _j \dfrac{1}{\beta (\beta -1)} x_j^\beta \bigg ).} \end{aligned}$$

Defining $h(\varvec{x}) = \langle \varvec{z}, \varvec{x}\rangle - \sum _j \dfrac{1}{\beta (\beta -1)} x_j^\beta$, and differentiating with respect to $x_j$:

$$\begin{aligned} \nabla _{x_j} h&= z_j - \dfrac{1}{\beta -1}x_j^{\beta -1} = 0 \\ x_j&= (\beta -1)^{\frac{1}{\beta -1}} z_j^{\frac{1}{\beta -1}} . \end{aligned}$$

Plugging into $\phi ^\star (\varvec{x})$,

$$\begin{aligned} \phi ^\star (\varvec{z})&= \sum _j z_j (\beta -1)^{\frac{1}{\beta -1}} z_j^{\frac{1}{\beta -1}} - \sum _j \dfrac{1}{\beta (\beta -1)} \left( (\beta -1)^{\frac{1}{\beta -1}} z_j^{\frac{1}{\beta -1}}\right) ^\beta \\&=(\beta -1)^{\frac{1}{\beta -1}} \sum _j z_j^{\frac{\beta }{\beta -1}} - \dfrac{ (\beta -1)^{\frac{1}{\beta -1}} }{\beta } \sum _j z_j^{\frac{\beta }{\beta -1}} \\&=(\beta -1)^{\frac{1}{\beta -1}} \bigg ( 1 - \dfrac{1 }{\beta } \bigg ) \sum _j z_j^{\frac{\beta }{\beta -1}} . \end{aligned}$$

Finally, differentiating the Fenchel conjugate yields

$$\begin{aligned} \nabla \phi ^\star (\varvec{z}) =(\beta -1)^{\frac{1}{\beta -1}} \bigg ( 1 - \dfrac{1 }{\beta } \bigg ) \frac{\beta }{\beta -1} \varvec{z}^{\frac{1}{\beta -1}} = (\beta -1)^{\frac{1}{\beta -1}} \varvec{z}^{\frac{1}{\beta -1} } . \end{aligned}$$

Thus, the projection of $\varvec{x}$ onto a hyperplane is given by

$$\begin{aligned} {\mathcal {P}^\phi _{Q(a,c)} (\varvec{x}) =(\beta -1)^{\frac{1}{\beta -1}} \bigg ( \dfrac{1}{\beta -1} \varvec{x}^{\beta -1} - \widetilde{\gamma } a\bigg ) ^{\frac{1}{\beta -1}},} \end{aligned}$$

where $\widetilde{\gamma } = \text {arg}\min _{\gamma } (\beta -1)^{\frac{1}{\beta -1}} \bigg ( 1 - \dfrac{1 }{\beta } \bigg ) \sum _j \bigg ( \dfrac{1}{\beta -1} x_j^{\beta -1} - \gamma a_j\bigg )^{\frac{\beta }{\beta -1}} + \gamma c .$

1.3 Per-iteration complexity of IMRT study

We detail the calculations behind our per iteration complexity remarks. Note that in the IMRT dataset considered, $p \ll n < m$.

1.3.1 Linear

How many flops are required to compute a single MM update given by Eq. (8) in the linear case?

$$\begin{aligned} \nabla f(\varvec{x}_{k})= & {} \sum _i \varvec{v}_i\big (\varvec{x}_{k} - \mathcal {P}_{C_i}(\varvec{x}_{k})\big ) + \sum _j w_j \varvec{A}^t\big (\varvec{A}\varvec{x}_{k} - \mathcal {P}_{Q_j}(\varvec{A}\varvec{x}_{k})\big ) \\ \varvec{x}_{k+1}= & {} \varvec{x}_{k} - \big [v\varvec{I}+ w\varvec{A}^t\varvec{A}\big ]^{-1}\nabla f(\varvec{x}_{k}). \end{aligned}$$

We first tally the flops required to compute the gradient $\nabla f(\varvec{x}_{k})$. In the IMRT example, the first sum $\sum _i v_i(\varvec{x}_{k} - \mathcal {P}_{C_i}(\varvec{x}_{k}))$ requires $\mathcal {O}(n)$ flops. The matrix-vector product $\varvec{z}_k = \varvec{A}\varvec{x}_{k}$ requires $\mathcal {O}(mn)$ flops. The sum $\varvec{y}_k = \sum _j w_j(\varvec{z}_k - \mathcal {P}_{Q_j}(\varvec{z}_k))$ requires $\mathcal {O}(m)$ flops. The matrix-vector product $\varvec{A}^t\varvec{y}_k$ requires $\mathcal {O}(mn)$ flops. Adding the two sums requires $\mathcal {O}(n)$ flops. In total, the gradient requires $\mathcal {O}(mn)$ flops.

We next tally the flops required to compute the MM update. Forming the matrix $v\varvec{I}+ w\varvec{A}^t\varvec{A}$ requires $\mathcal {O}(mn^2)$ flops. Computing its Cholesky factorization requires $\mathcal {O}(n^3)$ flops. We only need to compute the factorization once and can cache it. Subsequent iterations will require $\mathcal {O}(n^2)$ flops. Thus, the amount of work required to compute a single MM update for the linear formulation is $\mathcal {O}(mn)$.

Note that the exact update in the linear case given by

$$\begin{aligned} \varvec{x}_{k+1} = (v\varvec{I}+ w\varvec{A}^t\varvec{A})^{-1} \left[ \sum _i \varvec{v}_i \mathcal {P}_{C_i}(\varvec{x}_{k}) + \varvec{A}^t \sum _j w_j \varvec{A}\mathcal {P}_{Q_j}(\varvec{A}\varvec{x}_{k}) \right] \end{aligned}$$

has the same complexity as the update considered above.

1.3.2 Non-linear

We next consider the number of flops required for an MM update in the non-linear case:

$$\begin{aligned} \nabla f(\varvec{x}_{k})= & {} \sum _i v_i\big (\varvec{x}_{k} - \mathcal {P}_{C_i}(\varvec{x}_{k})\big ) + \nabla h(\varvec{x}_{k})\sum _j w_j \big (h(\varvec{x}_{k}) - \mathcal {P}_{Q_j}(h(\varvec{x}_{k}))\big ) \\ \varvec{x}_{k+1}= & {} \varvec{x}_{k} - \eta _{\varvec{x}_{k}}\big [v\varvec{I}+ w\nabla h(\varvec{x}_{k})dh(\varvec{x}_{k})\big ]^{-1}\nabla f(\varvec{x}_{k}) \\= & {} \varvec{x}_{k} - \eta _{\varvec{x}_{k}}\left[ \frac{1}{v}\varvec{I}- \frac{w}{v^2}\nabla h(\varvec{x}_{k}) \left( \varvec{I}+ \frac{w}{v}dh(\varvec{x}_{k})\nabla h(\varvec{x}_{k}) \right) ^{-1} dh(\varvec{x}_{k}) \right] \nabla f(\varvec{x}_{k}) \\= & {} \varvec{x}_{k} - \eta _{\varvec{x}_{k}}\left[ \frac{1}{v}\nabla f(\varvec{x}_{k}) - \frac{w}{v}\nabla h(\varvec{x}_{k}) \big (v\varvec{I}+ { {wdh}}(\varvec{x}_{k})\nabla h(\varvec{x}_{k}) \big )^{-1} dh(\varvec{x}_{k})\nabla f(\varvec{x}_{k}) \right] . \end{aligned}$$

We first tally the flops required to compute the gradient $\nabla f(\varvec{x}_{k})$. The first sum $\sum _i v_i(\varvec{x}_{k} - \mathcal {P}_{C_i}(\varvec{x}_{k}))$ requires $\mathcal {O}(n)$ flops. Computing $\varvec{A}\varvec{x}_{k}$ requires $\mathcal {O}(mn)$ flops. Computing the sum $\sum _j w_j (h(\varvec{x}_{k}) - \mathcal {P}_{Q_j}(h(\varvec{x}_{k})))$ requires $\mathcal {O}(m)$ flops. Computing $\nabla h(\varvec{x}_{k})$ requires $\mathcal {O}(mn)$ flops, and computing its product with the sum term $\nabla h(\varvec{x}_{k})\sum _j w_j (h(\varvec{x}_{k}) - \mathcal {P}_{Q_j}(h(\varvec{x}_{k})))$ requires $\mathcal {O}(mp)$ flops. In total, the gradient requires $\mathcal {O}(mn)$ flops. As in the linear case, the dominant cost are matrix-vector products involving the matrix $\varvec{A}$.

We next tally the flops required to compute the MM update. Forming the matrix $v\varvec{I}+ w\nabla h(\varvec{x}_{k})dh(\varvec{x}_{k})$ requires $\mathcal {O}(np^2)$ flops. Computing its Cholesky factorization requires $\mathcal {O}(p^3)$ flops. Computing $dh(\varvec{x}_{k})\nabla f(\varvec{x}_{k})$ requires $\mathcal {O}(np)$ flops. Computing $\left( v\varvec{I}+ wdh(\varvec{x}_{k})\nabla h(\varvec{x}_{k}) \right) ^{-1} dh(\varvec{x}_{k})\nabla f(\varvec{x}_{k})$ requires $\mathcal {O}(p^2)$ flops. The product $\nabla h(\varvec{x}_{k})\left( v\varvec{I}+ wdh(\varvec{x}_{k})\nabla h(\varvec{x}_{k}) \right) ^{-1} dh(\varvec{x}_{k})\nabla f(\varvec{x}_{k})$ requires $\mathcal {O}(np)$ flops. Computing a candidate update requires $\mathcal {O}(n)$ flops. An objective function evaluation requires $\mathcal {O}(mn)$ flops. Thus, including the line search, the amount of work required to compute a single MM update for the linear formulation is $\mathcal {O}(\max \{mn,np^2\})$. When $p^2 < m$, we conclude that the computational work for a single MM update for the non-linear case is comparable to that of the linear case. In practice, reducing the problem size via a non-linear formulation may additionally reduce the number of MM updates.