(62) $\displaystyle\beta_{t}$ $\displaystyle=\frac{b_{t}}{B_{t+1}},$ (63)
$\displaystyle y_{t}$ $\displaystyle=(1-\beta_{t})x_{t}+\beta_{t}v_{t}.$ (64)

######  Lemma 1.

From [59, Lemma 4]. The Bregman divergence of the function $\|x\|^{i}$
satisfies, for $i\geq 2$,

$\|x\|^{i}-\|y\|^{i}-\nabla(\|y\|^{i})(x-y)\geq\frac{1}{2^{i-2}}\|x-y\|^{i}.$

######  Proposition 10.

The function $\phi_{t}$ is lower-bounded by

$\phi_{t}\geq\underbrace{\phi_{t}(v_{t})}_{=\phi_{t}^{\star}}+\frac{\lambda_{t}^{(1)}}{2}\|x-v_{t}\|^{2}+\frac{\lambda_{t}^{(2)}}{12}\|x-v_{t}\|^{3}$
(65)

where $v_{t}=\operatorname*{{arg\,min}}_{x}\phi_{t}(x)$.

###### Proof.

The first order condition on $\phi_{t}$ reads,

$\ell^{\prime}_{t}+\nabla\left(\frac{\lambda_{t}^{(1)}}{2}\|v_{t}-x_{0}\|^{2}+\frac{\lambda_{t}^{(2)}}{6}\|v_{t}-x_{0}\|^{3}\right)=0.$

Multiplying both sides by $(x-v_{t})$ gives

$\ell^{\prime}_{t}(x-v_{t})+\nabla\left(\frac{\lambda_{t}^{(1)}}{2}\|v_{t}-x_{0}\|^{2}+\frac{\lambda_{t}^{(2)}}{6}\|v_{t}-x_{0}\|^{3}\right)(x-v_{t})=0.$

Note that, since $\ell_{t}$ is an affine function,
$\ell^{\prime}_{t}(x-v_{t})=\ell_{t}(x)-\ell_{t}(v_{t})$. Hence,

$\ell_{t}(x)-\ell_{t}(v_{t})+\nabla\left(\frac{\lambda_{t}^{(1)}}{2}\|v_{t}-x_{0}\|^{2}+\frac{\lambda_{t}^{(2)}}{6}\|v_{t}-x_{0}\|^{3}\right)(x-v_{t})=0.$

Finally, adding
$\frac{\lambda_{t}^{(1)}}{2}\|x-x_{0}\|^{2}+\frac{\lambda_{t}^{(2)}}{6}\|x-x_{0}\|^{3}$
on both sides and after reorganizing the terms,

$\textstyle\phi_{t}(x)=\ell_{t}(v_{t})+\frac{\lambda_{t}^{(1)}}{2}\|x-x_{0}\|^{2}+\frac{\lambda_{t}^{(2)}}{6}\|x-x_{0}\|^{3}-\nabla\left(\frac{\lambda_{t}^{(1)}}{2}\|v_{t}-x_{0}\|^{2}+\frac{\lambda_{t}^{(2)}}{6}\|v_{t}-x_{0}\|^{3}\right)(x-v_{t}).$
(66)

From lemma 1 with $x=x-x_{0}$, $y=v_{t}-x_{0}$, and after reorganizing the
terms,

$\|x-x_{0}\|^{i}-\nabla(\|v_{t}-x_{0}\|^{i})(x-v_{t})\geq\frac{1}{2^{i-2}}\|x-v_{t}\|^{i}+\|v_{t}-x_{0}\|^{i}.$

Therefore, using the previous inequality with $i=2$ and $i=3$, eq. 66 becomes

$\phi_{t}(x)\geq\ell_{t}(v_{t})+\frac{\lambda_{t}^{(1)}}{2}\|v_{t}-x_{0}\|^{2}+\frac{\lambda_{t}^{(2)}}{6}\|v_{t}-x_{0}\|^{3}+\frac{\lambda_{t}^{(2)}}{2}\|v_{t}-x\|^{2}+\frac{\lambda_{t}^{(3)}}{12}\|v_{t}-x\|^{3}$

By definition of $\phi_{t}^{\star}=\phi_{t}(v_{t})$,

$\phi_{t}(x)\geq\phi_{t}^{\star}+\frac{\lambda_{t}^{(1)}}{2}\|v_{t}-x\|^{2}+\frac{\lambda_{t}^{(2)}}{12}\|v_{t}-x\|^{3}.$

∎

######  Proposition 11.

Let

$\gamma=\frac{\kappa_{D}}{\|D\|}\left(\frac{3}{2}\|\varepsilon\|+2\frac{\|(I-P)G\|}{M}\right).$

Then, under the assumptions of proposition 4 the condition

$\frac{\|f(x_{+})\|^{2}}{M\left(\gamma+\|D\alpha\|\right)}\leq-\nabla
f(x)^{T}D\alpha$

is guaranteed as long as $M\geq 2L$.

###### Proof.

The starting point is eq. 50 combined with proposition 4:

$\displaystyle\left\|\left(\frac{1}{2}D_{\dagger}\Gamma
D_{\dagger}^{T}+\frac{M\|D\alpha\|}{2}\right)D\alpha+\nabla f(x_{+})\right\|$
$\displaystyle\leq\frac{L}{2}\|D\alpha\|^{2}+\|D\alpha\|\left(\frac{L}{2}\frac{\|\varepsilon\|}{\|D\|}\kappa_{D}+\frac{\|(I-P)\nabla^{2}f(x)D\alpha\|}{\|D\alpha\|}\right)$
$\displaystyle\leq\frac{L}{2}\|D\alpha\|^{2}+\|D\alpha\|\left(\frac{L}{2}\frac{\|\varepsilon\|}{\|D\|}\kappa_{D}+\left(\|(I-P)G\|+L\|\varepsilon\|\right)\frac{\kappa_{D}}{\|D\|}\right)$
$\displaystyle\leq\frac{L}{2}\|D\alpha\|^{2}+\|D\alpha\|\left(\frac{3L}{2}\frac{\|\varepsilon\|}{\|D\|}\kappa_{D}+|(I-P)G\|\frac{\kappa_{D}}{\|D\|}\right)$

To simplify, let $\Gamma=MD\gamma D^{T}$. Hence,

$\left\|M\left(\frac{\|D\alpha\|+\gamma}{2}\right)D\alpha+\nabla
f(x_{+})\right\|\leq\frac{L}{2}\|D\alpha\|^{2}+\|D\alpha\|\left(\frac{3L}{2}\frac{\|\varepsilon\|}{\|D\|}\kappa_{D}+\|(I-P)G\|\frac{\kappa_{D}}{\|D\|}\right)$

Elevating to the square this inequality gives

$\displaystyle\left(M\left(\frac{\gamma+\|D\alpha\|}{2}\right)\right)^{2}\|D\alpha\|^{2}+\|\nabla
f(x_{+})\|^{2}+2\left(M\left(\frac{\gamma+\|D\alpha\|}{2}\right)\right)\nabla
f(x_{+})^{T}D\alpha$ $\displaystyle\leq\;\;$
$\displaystyle\|D\alpha\|^{2}\left(\frac{L}{2}\|D\alpha\|+\frac{L}{2}\frac{\|\varepsilon\|}{\|D\|}\kappa_{D}+\frac{\|(I-P)\nabla^{2}f(x)D\alpha\|}{\|D\alpha\|}\right)^{2}.$

The desired result holds if the following condition is satisfied,

$\displaystyle\left(M\left(\frac{\gamma+\|D\alpha\|}{2}\right)\right)^{2}\|D\alpha\|^{2}\geq\|D\alpha\|^{2}\left(\frac{L}{2}\|D\alpha\|+\frac{3L}{2}\frac{\|\varepsilon\|}{\|D\|}\kappa_{D}+\frac{\|(I-P)G\|\kappa_{D}}{\|D\|}\right)^{2}.$

After simplification of the squares,

$M\frac{\gamma+\|D\alpha\|}{2}\geq\frac{L}{2}\|D\alpha\|+\frac{3L}{2}\frac{\|\varepsilon\|}{\|D\|}\kappa_{D}+\frac{\|(I-P)G\|\kappa_{D}}{\|D\|}.$

Replacing $\gamma$ by its value gives

$M\frac{\|D\alpha\|+\frac{\kappa_{D}}{\|D\|}\left(\frac{3}{2}\|\varepsilon\|+2\frac{\|(I-P)G\|}{M}\right)}{2}\geq\frac{L}{2}\|D\alpha\|+\frac{3L}{2}\frac{\|\varepsilon\|}{\|D\|}\kappa_{D}+\frac{\|(I-P)G\|\kappa_{D}}{\|D\|}.$

The condition is simplified into

$(M-L)\frac{\|D\alpha\|}{2}+(M-2L)\frac{3}{2}\frac{\|\varepsilon\|\kappa_{D}}{\|D\|}\geq
0.$

This condition is implied by $M\geq 2L$. ∎

######  Proposition 12.

Under the same assumptions as proposition 7, if $M\geq 2L$, and if

$\gamma=\frac{\kappa_{D}}{\|D\|}\left(\frac{3}{2}\|\varepsilon\|+2\frac{\|(I-P)G\|}{M}\right)\leq\frac{(\sqrt{3}-1)\|D\alpha\|}{4},$

then

$\frac{2}{3^{3/4}}\frac{\|\nabla f(x_{+})\|^{3/2}}{\sqrt{M}}\leq-\nabla
f(x_{+})^{T}D\alpha.$

###### Proof.

The starting point is eq. 50,

$\displaystyle\left\|M\frac{\|D\alpha\|}{2}D\alpha+\nabla
f(x_{+})\right\|\leq\frac{L}{2}\|D\alpha\|^{2}+\|D\alpha\|\left(\frac{L}{2}\frac{\kappa_{D}\|\varepsilon\|}{\|D\|}+\frac{M\gamma}{2}+\frac{\left\|(I-P)\nabla^{2}f(x)D\alpha\right\|}{\|D\alpha\|}\right)$

Therefore, to obtain

$\displaystyle\left\|M\frac{\|D\alpha\|}{2}D\alpha+\nabla f(x_{+})\right\|\leq
M\left(\frac{\|D\alpha\|}{4}+\gamma\right)\|D\alpha\|,$

The following is sufficient,

$M\left(\frac{\|D\alpha\|}{4}+\gamma\right)\|D\alpha\|\geq\frac{L}{2}\|D\alpha\|^{2}+\|D\alpha\|\left(\frac{L}{2}\frac{\kappa_{D}\|\varepsilon\|}{\|D\|}+\frac{M\gamma}{2}+\frac{\left\|(I-P)\nabla^{2}f(x)D\alpha\right\|}{\|D\alpha\|}\right).$

Using proposition 4, the condition can be strengthened into

$\displaystyle\frac{M}{2}\left(\frac{\|D\alpha\|+\gamma}{2}\right)\|D\alpha\|$
$\displaystyle\geq\frac{L}{2}\|D\alpha\|^{2}+\|D\alpha\|\left(\frac{L}{2}\frac{\kappa_{D}\|\varepsilon\|}{\|D\|}+\frac{M\gamma}{2}+\left(\|(I-P)G\|+L\|\varepsilon\|\right)\frac{\kappa_{D}}{\|D\|}\right)$
$\displaystyle=\frac{L}{2}\|D\alpha\|^{2}+\|D\alpha\|\left(\frac{3L}{2}\frac{\kappa_{D}\|\varepsilon\|}{\|D\|}+\frac{M\gamma}{2}+\|(I-P)G\|\frac{\kappa_{D}}{\|D\|}\right)$

Defining

$\frac{\gamma}{2}=\left(\frac{3}{4}\frac{\kappa_{D}\|\varepsilon\|}{\|D\|}+\frac{\|(I-P)G\|\frac{\kappa_{D}}{\|D\|}}{M}\right)$

simplifies the condition into

$M\left(\frac{\|D\alpha\|}{4}+\gamma\right)\|D\alpha\|\geq\frac{L}{2}\|D\alpha\|^{2}+\|D\alpha\|\left(M\gamma+\frac{3(L-\frac{M}{2})}{2}\frac{\kappa_{D}\|\varepsilon\|}{\|D\|}\right)$

which is satisfied when $M>2L$. Now, assume that

$\gamma\leq\frac{(\sqrt{3}-1)\|D\alpha\|}{4}.$

Then,

$\left\|M\frac{\|D\alpha\|}{2}D\alpha+\nabla
f(x_{+})\right\|\leq\sqrt{3}\frac{M\|D\alpha\|^{2}}{4}.$

Elevating both sides to the square gives

$\|\nabla
f(x_{+})\|^{2}+\frac{3M^{2}\|D\alpha\|^{4}}{16}\leq-M\|D\alpha\|\nabla
f(x_{+})^{T}D\alpha$

Writing $r=\|D\alpha\|$,

$\frac{\|\nabla f(x_{+})\|^{2}}{Mr}+\frac{3Mr^{3}}{16}\leq-\nabla
f(x_{+})^{T}D\alpha.$

Using

$\frac{c_{1}}{r}+c_{2}r^{3}\geq
4c_{2}^{1/4}\left(\frac{c_{1}}{3}\right)^{3/4},$

the inequality becomes

$\displaystyle-\nabla f(x_{+})^{T}D\alpha$
$\displaystyle\geq\frac{M^{1/4}}{2}\frac{\|\nabla
f(x_{+})\|^{3/2}}{M^{3/4}}\frac{4}{3^{3/4}}$
$\displaystyle=\frac{2}{3^{3/4}}\frac{\|\nabla f(x_{+})\|^{3/2}}{\sqrt{M}}.$

∎

######  Proposition 13 (Termination of algorithm 6).

Let $f$ satisfies Assumption 1. Assume that Requirements 1b, 3 and 2 holds.
Then, once $M\geq 2L$, algorithm 6 terminates with ExitFlag equals to either
SmallStep or LargeStep. Moreover, if $M_{0}\leq L$, then the algorithm
terminates with $M\leq 4L$. Moreover, if the algorithm terminates with
ExitFlag equals to SmallStep, then

$\|D\alpha\|\leq\frac{4\gamma_{M}}{\sqrt{3}-1},\quad\gamma_{M}=\frac{\kappa_{D}}{\|D\|}\left(\frac{3}{2}\|\varepsilon\|+2\frac{\|(I-P)G\|}{M}\right).$

###### Proof.

Let

$\gamma_{M}=\frac{\kappa_{D}}{\|D\|}\left(\frac{3}{2}\|\varepsilon\|+2\frac{\|(I-P)G\|}{M}\right).$

Assume that $M\geq 2L$. If $\gamma_{M}\leq\frac{(\sqrt{3}-1)\|D\alpha\|}{4}$,
then, by proposition 12, the following condition is satisfied:

$\frac{2}{3^{3/4}}\frac{\|\nabla f(x_{+})\|^{3/2}}{\sqrt{M}}\leq-\nabla
f(x_{+})^{T}D\alpha.$

In this case the algorithm terminates with ExitFlag = LargeStep. In any case,
by proposition 11, the following conditions is always satisfied when $M\geq
2L$:

$\displaystyle\frac{\|f(x_{+})\|^{2}}{M\left(\gamma+\|D\alpha\|\right)}$
$\displaystyle\leq-\nabla f(x)^{T}D\alpha.$

Then, if $\gamma_{M}\geq\frac{(\sqrt{3}-1)\|D\alpha\|}{4}$, the algorithm
terminates with ExitFlag = SmallStep (otherwise the algorithm would have been
terminated with ExitFlag = LargeStep).

Since the algorithm doubles $M$ until one of the two condition is satisfied,
in the worst case, $M=4L$. ∎

######  Proposition 14.

If $\lambda_{t}^{(1)}$ and $\lambda_{t}^{(2)}$ satisfy

$\lambda^{(1)}_{t}\geq\frac{b_{t+1}^{2}}{B_{t}}M_{t+1}\left(\gamma_{t}+\|D_{t}\alpha_{t}\|\right),\quad\lambda^{(2)}_{t}\geq\frac{4}{\sqrt{3}}\frac{b_{t+1}^{3}}{B^{2}_{t}}M_{t+1},$

where
$\gamma_{t}=\frac{\kappa_{D_{t}}}{\|D_{t}\|}\left(\frac{3}{2}\|\varepsilon_{t}\|+2\frac{\|(I-P_{t})G_{t}\|}{M_{t+1}}\right).$
Then, the function $\phi$ satisfies

$B_{t}f(x_{t})\leq\phi_{t}(x),\qquad\phi_{t}(x)\leq
B_{t}f(x)+\frac{\lambda_{t}^{(1)}+\tilde{\lambda}^{(1)}}{2}\|x-x_{0}\|^{2}+\frac{\lambda_{t}^{(2)}+\tilde{\lambda}^{(2)}}{6}\|x-x_{0}\|^{3},$

where

$\tilde{\lambda}^{(1)}=\|\nabla f(x_{0})-P_{0}\nabla
f(x_{0})P_{0}\|+\delta\left(\frac{L\kappa+M_{1}\kappa^{2}}{2}\right),\quad\tilde{\lambda}^{(2)}=M_{1}+L.$

###### Proof.

The result is proven by recursion. At $t=1$, the condition
$B_{t}f(x_{t})\leq\phi_{t}(x)$ is obviously satisfied since

$f(x_{1})\leq\min_{v}\phi_{1}(v)=f(x_{1}).$

On the other hand, by proposition 8,

$\displaystyle f(x_{1})$
$\displaystyle\leq\min_{x}f(x)+\frac{\tilde{\lambda}^{(2)}}{6}\|x-x_{0}\|^{3}+\frac{\tilde{\lambda}^{(1)}}{2}\|x-x_{0}\|^{2}$
$\displaystyle\leq
f(x)+\frac{\tilde{\lambda}^{(2)}}{6}\|x-x_{0}\|^{3}+\frac{\tilde{\lambda}^{(1)}}{2}\|x-x_{0}\|^{2}.$

Therefore, the second condition holds by definition of $\phi$,

$\displaystyle\phi_{t}$
$\displaystyle=f(x_{1})+\frac{\lambda_{t}^{(1)}}{2}\|x-x_{0}\|^{2}+\frac{\lambda_{t}^{(2)}}{6}\|x-x_{0}\|^{3}$
$\displaystyle\leq\frac{\lambda_{1}^{(1)}+\tilde{\lambda}^{(1)}}{2}\|x-x_{0}\|^{2}+\frac{\lambda_{1}^{(2)}+\tilde{\lambda}^{(2)}}{6}\|x-x_{0}\|^{3}.$

Now, assume $t>1$, and $B_{t}f(x_{t})\leq\phi_{t}(x)$. Hence,

$\displaystyle\min_{x}\phi_{t+1}(x)$ $\displaystyle=$
$\displaystyle\min_{x}\ell_{t}(x)+b_{t}\left[f(x_{t+1})+\nabla
f(x_{t+1})(x-x_{t+1})\right]+\frac{\lambda_{t+1}^{(1)}}{2}\|x-x_{0}\|^{2}+\frac{\lambda_{t+1}^{(2)}}{6}\|x-x_{0}\|^{3}$
$\displaystyle=$
$\displaystyle\min_{x}\phi_{t}(x)+b_{t}\left[f(x_{t+1})+\nabla
f(x_{t+1})(x-x_{t+1})\right]$
$\displaystyle\quad+\frac{\lambda_{t+1}^{(1)}-\lambda_{t}^{(1)}}{2}\|x-x_{0}\|^{2}+\frac{\lambda_{t+1}^{(2)}-\lambda_{t}^{(2)}}{6}\|x-x_{0}\|^{3}$
$\displaystyle\geq$
$\displaystyle\min_{x}\phi_{t}(x)+b_{t}\left[f(x_{t+1})+\nabla
f(x_{t+1})(x-x_{t+1})\right]$
$\displaystyle\overset{\eqref{eq:lower_bound_phi_star}}{\geq}$
$\displaystyle\min_{x}\phi_{t}^{\star}+\frac{\lambda_{t}^{(1)}}{2}\|x-v_{t}\|^{2}+\frac{\lambda_{t}^{(2)}}{12}\|x-v_{t}\|^{3}+b_{t}\left[f(x_{t+1})+\nabla
f(x_{t+1})(x-x_{t+1})\right]$ $\displaystyle\geq$
$\displaystyle\min_{x}B_{t}f(x_{t})+\frac{\lambda_{t}^{(1)}}{2}\|x-v_{t}\|^{2}+\frac{\lambda_{t}^{(2)}}{12}\|x-v_{t}\|^{3}+b_{t}\left[f(x_{t+1})+\nabla
f(x_{t+1})(x-x_{t+1})\right]$
$\displaystyle\overset{\text{A.}\ref{assump:convex}}{\geq}$
$\displaystyle\min_{x}B_{t}f(x_{t+1})+\nabla
f(x_{t+1})(x_{t}-x_{t+1})+b_{t}\left[f(x_{t+1})+\nabla
f(x_{t+1})(x-x_{t+1})\right]$
$\displaystyle\quad+\frac{\lambda_{t}^{(1)}}{2}\|x-v_{t}\|^{2}+\frac{\lambda_{t}^{(2)}}{12}\|x-v_{t}\|^{3}$
$\displaystyle=$ $\displaystyle\min_{x}B_{t+1}f(x_{t+1})+\nabla
f(x_{t+1})(B_{t}x_{t}+b_{t}x-B_{t+1}x_{t+1})+\frac{\lambda_{t}^{(1)}}{2}\|x-v_{t}\|^{2}+\frac{\lambda_{t}^{(2)}}{12}\|x-v_{t}\|^{3}$
$\displaystyle\overset{\eqref{eq:def_y}}{=}$
$\displaystyle\min_{x}B_{t+1}f(x_{t+1})+B_{t+1}\nabla
f(x_{t+1})(y_{t}-x_{t+1})$ $\displaystyle\quad+b_{t}\nabla
f(x_{t+1})(x-v_{t})+\frac{\lambda_{t}^{(1)}}{2}\|x-v_{t}\|^{2}+\frac{\lambda_{t}^{(2)}}{12}\|x-v_{t}\|^{3}$

The inequality is satisfied if either

$\displaystyle\textbf{(a)}\quad 0$ $\displaystyle\leq B_{t+1}\nabla
f(x_{t+1})(y_{t}-x_{t+1})+b_{t}\nabla
f(x_{t+1})(x-v_{t})+\frac{\lambda_{t}^{(2)}}{12}\|x-v_{t}\|^{3},\;\;\text{or}$
$\displaystyle\textbf{({b})}\quad 0$ $\displaystyle\leq B_{t+1}\nabla
f(x_{t+1})(y_{t}-x_{t+1})+b_{t}\nabla
f(x_{t+1})(x-v_{t})+\frac{\lambda_{t}^{(1)}}{2}\|x-v_{t}\|^{2}.$

It remains now to find sufficient condition such that one of the previous
inequalities hold.

Define $x_{t+1}$ to be the output of algorithm 6 starting from $y_{t}$, hence
$y_{t}-x_{t+1}=-D_{t}\alpha_{t}$. The algorithm guarantees that

$\displaystyle\textbf{(a)}\quad-\nabla f(x_{t+1})^{T}D_{t}\alpha_{t}$
$\displaystyle\geq\frac{2}{3^{3/4}}\frac{\|\nabla
f(x_{t+1})\|^{3/2}}{\sqrt{M_{t+1}}}\quad\text{and}\;\;\quad\text{or}$ (67)
$\displaystyle\textbf{(b)}\quad-\nabla f(x_{t+1})^{T}D_{t}\alpha_{t}$
$\displaystyle\geq\frac{\|f(x_{t+1})\|^{2}}{M_{t+1}\left(\gamma_{t}+\|D_{t}\alpha_{t}\|\right)}$
(68)

Combining the expressions (a) and (b) leads to the following sufficient
conditions:

$\displaystyle 0$ $\displaystyle\leq B_{t+1}\frac{2}{3^{3/4}}\frac{\|\nabla
f(x_{t+1})\|^{3/2}}{\sqrt{M_{t+1}}}+b_{t}\nabla
f(x_{t+1})(x-v_{t})+\frac{\lambda_{t}^{(2)}}{12}\|x-v_{t}\|^{3},$ (69)
$\displaystyle\quad 0$ $\displaystyle\leq
B_{t+1}\frac{\|f(x_{t+1})\|^{2}}{M_{t+1}\left(\gamma_{t}+\|D_{t}\alpha_{t}\|\right)}+b_{t}\nabla
f(x_{t+1})(x-v_{t})+\frac{\lambda_{t}^{(1)}}{2}\|x-v_{t}\|^{2}.$ (70)

##### Case 1: equation (69).

Starting from the first order condition of the minimum of eq. 69 over $x$,

$b_{t}\nabla f(x_{t+1})+\frac{\lambda_{t}^{(2)}}{4}\|x-v_{t}\|(x-v_{t})=0.$
(71)

Multiplying eq. 71 by $(x-v_{t})$ gives

$b_{t}\nabla f(x_{t+1})(x-v_{t})=-\frac{\lambda_{t}^{(2)}}{4}\|x-v_{t}\|^{3}$

Hence, when $x$ satisfies eq. 71,

$b_{t}\nabla
f(x_{t+1})(x-v_{t})+\frac{\lambda_{t}^{(2)}}{12}\|x-v_{t}\|^{3}=-\frac{\lambda_{t}^{(2)}}{6}\|x-v_{t}\|^{3}.$
(72)

Going back to eq. 71, after isolating $x-v_{t}$,

$\displaystyle(x-v_{t})=-\frac{4b_{t}}{\lambda_{t}^{(2)}}\nabla
f(x_{t+1})\frac{1}{\|x-v_{t}\|}$

Therefore, after taking the norm and changing the power,

$\displaystyle\|x-v_{t}\|^{3}$
$\displaystyle=\left(\frac{4b_{t}}{\lambda_{t}^{(2)}}\|\nabla
f(x_{t+1})\|\right)^{3/2},$
$\displaystyle\Leftrightarrow\frac{\lambda^{(2)}_{t}}{6}\|x-v_{t}\|^{3}$
$\displaystyle=\frac{\lambda^{(2)}_{t}}{6}\left(\frac{4b_{t}}{\lambda_{t}^{(2)}}\|\nabla
f(x_{t+1})\|\right)^{3/2}$
$\displaystyle=\frac{4}{3\sqrt{\lambda^{(2)}_{t}}}\left(b_{t}\|\nabla
f(x_{t+1})\|\right)^{3/2}.$

After using eq. 72 and injecting the minimal value makes the condition eq. 69
stronger:

$0\leq B_{t+1}\frac{2}{3^{3/4}}\frac{\|\nabla
f(x_{t+1})\|^{3/2}}{\sqrt{M_{t+1}}}-\frac{4}{3\sqrt{\lambda^{(2)}_{t}}}\left(b_{t}\|\nabla
f(x_{t+1})\|\right)^{3/2}.$

Hence, if $\lambda^{(2)}_{t}$ satisfies

$B_{t+1}\frac{2}{3^{3/4}\sqrt{M_{t+1}}}\geq\frac{4}{3\sqrt{\lambda^{(2)}_{t}}}b_{t}^{(3/2)}\quad\Leftrightarrow\quad\lambda^{(2)}_{t}\geq\frac{4}{\sqrt{3}}\frac{b_{t}^{3}}{B^{2}_{t+1}}M_{t+1},$
(73)

then eq. 69 is satisfied.

##### Case 2: equation (70).

Starting from the first order condition of the minimum of eq. 70 over $x$,

$b_{t+1}\nabla f(x_{t+1})+\lambda_{t}^{(1)}(x-v_{t}).$ (74)

Hence,

$(x-v_{t})=-\frac{b_{t}\nabla f(x_{t+1})}{\lambda_{t}^{(1)}}.$

Injecting the value back in eq. 70 gives

$B_{t+1}\frac{\|f(x_{t+1})\|^{2}}{M\left(\gamma_{t}+\|D_{t}\alpha_{t}\|\right)}-b^{2}_{t}\frac{\|\nabla
f(x_{t+1})\|^{2}}{\lambda_{t}^{(1)}}+\frac{1}{2}b^{2}_{t}\frac{\|\nabla
f(x_{t+1})\|^{2}}{\lambda_{t}^{(1)}}.$

Therefore, if the following condition holds,

$\frac{B_{t+1}}{2M_{t+1}\left(\gamma_{t}+\|D_{t}\alpha_{t}\|\right)}\geq\frac{b^{2}_{t}}{\lambda_{t}^{(1)}}\quad\Leftrightarrow\quad\lambda^{(1)}_{t}\geq\frac{b_{t}^{2}}{2B_{t+1}}M_{t+1}\left(\gamma_{t}+\|D_{t}\alpha_{t}\|\right),$

then eq. 70 is satisfied.

∎

######  Proposition 15.

Let $f$ satisfies Assumption 1. Then, under Requirements 3, 2 and 1b,
$\lambda_{t}^{(1)}$ and $\lambda_{t}^{(2)}$ in algorithm 7 are bounded by

$\displaystyle\lambda_{t}^{(1)}$ $\displaystyle\leq
30\cdot\frac{b_{t+1}^{2}}{B_{t}}\kappa_{D}\left(\delta\max\\{4L,M_{0}\\}+\max_{i=0\ldots
t}\|(I-P_{i})\nabla f(x_{i})P_{i})\|\right)$ (75)
$\displaystyle\lambda_{t}^{(2)}$
$\displaystyle\leq\frac{L}{2}\delta+\max_{i=0\ldots t}\|(I-P_{i})\nabla
f(x_{i})P_{i})\|.$ (76)

###### Proof.

Since algorithm 7 doubles $\lambda_{t}^{(1)}$, $\lambda_{t}^{(2)}$ until
$\phi_{t}^{\star}\geq f(x_{t+1})$, then by proposition 14, both
$\lambda_{t}^{(1)}$, $\lambda_{t}^{(2)}$ achieves at most

$\lambda^{(1)}_{t}\leq
2\cdot\frac{b_{t+1}^{2}}{B_{t}}M_{t+1}\left(\gamma_{t}+\|D_{t}\alpha_{t}\|\right),\quad\lambda^{(2)}_{t}\leq
2\cdot\frac{4}{\sqrt{3}}\frac{b_{t+1}^{3}}{B^{2}_{t}}M_{t+1}.$

There are three cases to distinguish:

  1. 1.

The algorithm finishes with ExitFlag = LargeStep,

  2. 2.

The algorithm finishes with ExitFlag = SmallStep.

##### Case 1.

In this case, $\lambda^{(2)}_{t+1}$ may be updated. By proposition proposition
13, $M_{t}\leq 4L$ (unless $M_{0}\geq 4L)$. Hence, $\lambda^{(2)}_{t}$ is
bounded by

$\lambda^{(2)}_{t}\leq
2\cdot\frac{4}{\sqrt{3}}\frac{b_{t+1}^{3}}{B^{2}_{t}}\max\\{M_{0},4L\\}\leq
5\frac{b_{t+1}^{3}}{B^{2}_{t}}\max\\{M_{0},4L\\}.$

##### Case 2.

In this case, $\lambda^{(1)}_{t+1}$ may be updated. By proposition 13, and by
Requirements 3 and 2,

$\displaystyle M_{t+1}\left(\gamma_{t}+\|D_{t}\alpha_{t}\|\right)$
$\displaystyle\leq\frac{\sqrt{3}+1}{\sqrt{3}-1}M_{t+1}\gamma_{t}$
$\displaystyle=\frac{\sqrt{3}+1}{\sqrt{3}-1}\frac{\kappa_{D_{t}}}{\|D_{t}\|}\left(\frac{3}{2}\|\varepsilon_{t}\|M_{t+1}+2\|(I-P_{t})G_{t}\|\right),$
$\displaystyle\leq\frac{\sqrt{3}+1}{\sqrt{3}-1}\left(\frac{3}{2}\delta\kappa_{D}\max\\{4L,M_{0}\\}+2\kappa_{D}\frac{\|(I-P_{t})G_{t}\|}{\|D_{t}\|}\right).$

In addition, by Theorem 1 and Requirement 2,

$\displaystyle\frac{\|(I-P_{t})G_{t}\|}{\|D_{t}\|}$
$\displaystyle\leq\frac{\|(I-P_{t})(G_{t}-\nabla
f(x_{t})D_{t})\|+\|(I-P_{t})\nabla f(x_{t})D_{t})\|}{\|D_{t}\|}$
$\displaystyle\leq\frac{\frac{L}{2}\|\varepsilon_{t}\|+\|(I-P_{t})\nabla
f(x_{t})D_{t})\|}{\|D_{t}\|},$
$\displaystyle=\frac{\frac{L}{2}\|\varepsilon_{t}\|+\|(I-P_{t})\nabla
f(x_{t})P_{t}D_{t})\|}{\|D_{t}\|},$
$\displaystyle\leq\frac{L}{2}\delta+\max_{i=0\ldots t}\|(I-P_{i})\nabla
f(x_{i})P_{i})\|.$

Hence,

$\displaystyle M_{t+1}\left(\gamma_{t}+\|D_{t}\alpha_{t}\|\right)$
$\displaystyle\leq\frac{\sqrt{3}+1}{\sqrt{3}-1}\left(\frac{3}{2}\delta\kappa_{D}\max\\{4L,M_{0}\\}+2\kappa_{D}\left(\frac{L}{2}\delta+\max_{i=0\ldots
t}\|(I-P_{i})\nabla f(x_{i})P_{i})\|\right)\right),$
$\displaystyle\leq\frac{\sqrt{3}+1}{\sqrt{3}-1}\left(\frac{7}{4}\delta\kappa_{D}\max\\{4L,M_{0}\\}+2\kappa_{D}\max_{i=0\ldots
t}\|(I-P_{i})\nabla f(x_{i})P_{i})\|\right).$ $\displaystyle\leq
7.5\kappa_{D}\left(\delta\max\\{4L,M_{0}\\}+\max_{i=0\ldots
t}\|(I-P_{i})\nabla f(x_{i})P_{i})\|\right).$

Therefore,

$\displaystyle\lambda^{(1)}_{t}$ $\displaystyle\leq
2\cdot\frac{b_{t+1}^{2}}{B_{t}}M_{t+1}\left(\gamma_{t}+\|D_{t}\alpha_{t}\|\right)$
$\displaystyle\leq
30\cdot\frac{b_{t+1}^{2}}{B_{t}}\kappa_{D}\left(\delta\max\\{4L,M_{0}\\}+\max_{i=0\ldots
t}\|(I-P_{i})\nabla f(x_{i})P_{i})\|\right)$

∎

###  G.5 Missing proofs from Sections 2 and 3

See 1

###### Proof.

Using Cauchy-Schwartz with eq. 2 gives that, for all $v$,

$v^{T}\left(\nabla f(y)-\nabla
f(z)-\nabla^{2}f(z)(y-z)\right)\leq\frac{L\|v\|}{2}\|y-z\|^{2}.$

Let $v=v_{i}$, $y=y_{i}$, and $z=z_{i}$. By the definition of $Y,\,Z,\,D,\,G$
in eq. 6,

$v_{i}^{T}\left(g_{i}-\nabla^{2}f(z_{i})d_{i}\right)\leq\frac{L\|v_{i}\|}{2}\|d_{i}\|^{2}.$

Introducing $\nabla^{2}f(x)$ gives

$\displaystyle v_{i}^{T}\left(g_{i}-\nabla^{2}f(z_{i})d_{i}\right)$
$\displaystyle=v_{i}^{T}\left(g_{i}-\nabla^{2}f(x)d_{i}\right)+v_{i}^{T}(\nabla^{2}f(z_{i})-\nabla^{2}f(x))d_{i}.$

Since the Hessian is $L$-Lipchitz-continuous Assumption 1,
$(\nabla^{2}f(z_{i})-\nabla^{2}f(x))d_{i}\leq L\|d_{i}\|\|z_{i}-x\|$.
Therefore, by the definition of $\varepsilon_{i}$,

$\displaystyle v_{i}^{T}\left(g_{i}-\nabla^{2}f(x)d_{i}\right)$
$\displaystyle\leq\frac{L\|v_{i}\|\varepsilon_{i}}{2}.$ (77)

Let $v_{i}=\operatorname{sign}(\alpha_{i})w$. Summing all inequalities
multiplied by $|\alpha_{i}|$ gives the first desired result:

$\displaystyle w^{T}\left(G-\nabla^{2}f(x)D\right)\alpha$
$\displaystyle\leq\frac{L\|w\|\sum_{i=1}^{N}\varepsilon_{i}|\alpha_{i}|}{2}.$

The second result is rather straightforward, since eq. 77 with $v_{i}=w$ gives

$w^{T}\left(g_{i}-\nabla^{2}f(x)d_{i}\right)\leq\frac{L\|w\|\varepsilon_{i}}{2}.$

Therefore,

$\sqrt{\sum_{i=1}^{N}\left(w^{T}\left(g_{i}-\nabla^{2}f(x)d_{i}\right)\right)^{2}}\leq\|w\|\sqrt{\sum_{i=1}^{N}\left\|g_{i}-\nabla^{2}f(x)d_{i}\right\|^{2}}\leq\|w\|\sqrt{\sum_{i=1}^{N}L\varepsilon_{i}^{2}}\leq\frac{L\|w\|\|\varepsilon\|}{2}.$

∎

See 2

###### Proof.

The inequality eq. Type II bound is a direct consequence of eq. 2 (with
$y=x_{+}$, $z=x$) combined with eq. 9,

$\displaystyle\|\nabla f(x_{+})-\nabla
f(x)-\nabla^{2}f(x)D\alpha\|\leq\frac{L}{2}\|D\alpha\|^{2}$
$\displaystyle\Leftrightarrow\;\;$ $\displaystyle w^{T}\left(\nabla
f(x_{+})-\nabla
f(x)-\nabla^{2}f(x)D\alpha\right)\leq\frac{L\|w\|}{2}\|D\alpha\|^{2}$
$\displaystyle\Leftrightarrow\;\;$ $\displaystyle w^{T}\nabla
f(x_{+})\leq\frac{L\|w\|}{2}\|D\alpha\|^{2}+w^{T}\left(\nabla
f(x)+\nabla^{2}f(x)D\alpha\right)$ $\displaystyle\Leftrightarrow\;\;$
$\displaystyle w^{T}\nabla
f(x_{+})\overset{\lx@cref{creftype~refnum}{eq:bound_hessian_scalar}}{\leq}\frac{L\|w\|}{2}\left(\|D\alpha\|^{2}+\sum_{i=1}^{N}|\alpha_{i}|\varepsilon_{i}\right)+w^{T}\left(\nabla
f(x)+G\alpha\right)$ $\displaystyle\Leftrightarrow\;\;$ $\displaystyle
w^{T}\nabla
f(x_{+})\leq\|w\|\left(\frac{L}{2}\left(\|D\alpha\|^{2}+\sum_{i=1}^{N}|\alpha_{i}|\varepsilon_{i}\right)+\|\nabla
f(x)+G\alpha\|\right)$

Setting $w=\nabla f(x_{+})$ gives eq. Type II bound.

The inequality eq. Type I bound instead comes from eq. 3 combined with eq. 9.
Indeed,

$\displaystyle f(x_{+})$ $\displaystyle\leq f(x)+\nabla
f(x)D\alpha+\frac{1}{2}(D\alpha)^{T}\nabla^{2}f(x)(D\alpha)+\frac{L}{6}\|D\alpha\|^{3}$
$\displaystyle\overset{\lx@cref{creftype~refnum}{eq:bound_hessian_dalpha}}{\leq}f(x)+\nabla
f(x)D\alpha+\frac{1}{2}\left((D\alpha)^{T}G\alpha+\frac{L\|D\alpha\|}{2}\sum_{i=1}^{N}|\alpha_{i}|\varepsilon_{i}\right)+\frac{L}{6}\|D\alpha\|^{3}$

It remains to use the followings bounds:

$\displaystyle\sum_{i=1}^{N}|\alpha_{i}|\varepsilon_{i}=\alpha^{T}(\operatorname{sign}(\alpha)\odot\varepsilon)$
$\displaystyle\leq\|\alpha\|\|\varepsilon\|,$ $\displaystyle\|D\alpha\|$
$\displaystyle\leq\|D\|\|\alpha\|.$

All together,

$f(x_{+})\leq f(x)+\nabla
f(x)D\alpha+\frac{1}{2}(D\alpha)^{T}G\alpha+\frac{L}{4}\|\alpha\|^{2}\|D\|\|\varepsilon\|+\frac{L}{6}\|D\alpha\|^{3}$

Finally, since $(D\alpha)^{T}G\alpha$ is a quadratic form, only the symmetric
counterpart of $D^{T}G$ counts. That means,
$(D\alpha)^{T}G\alpha=\alpha^{T}\frac{D^{T}G+G^{T}D}{2}\alpha$. Hence, writing
$H=\frac{D^{T}G+G^{T}D}{2}+\mathrm{I}\frac{L}{2}\|D\|\|\varepsilon\|$ gives
the desired result,

$f(x_{+})\leq f(x)+\nabla
f(x)D\alpha+\frac{\alpha^{T}H\alpha}{2}+\frac{L}{6}\|D\alpha\|^{3}.$

∎

See 3

###### Proof.

The proof is done by recursion.

Since $D_{1}$ is simply a vector of norm one, $D_{1}^{T}D_{1}=1$. Moreover,
$D_{1}=\frac{\nabla f(x_{0})}{\|\nabla f(x_{0})\|}$. Hence, is obviously span
$\nabla f(x_{0})$.

Assume that $D_{t-1}$ is orthogonal. Then, potentially removing one column
does not affect its orthogonality. Therefore, $\tilde{D}_{t-1}$ is orthogonal.
Now, consider the vector

$\displaystyle v$
$\displaystyle=\frac{(I-\tilde{D}_{t-1}(\tilde{D}_{t-1}^{T}\tilde{D}_{t-1})^{-1}\tilde{D}_{t-1}^{T})\nabla
f(x_{t})}{\|(I-\tilde{D}_{t-1}(\tilde{D}_{t-1}^{T}\tilde{D}_{t-1})^{-1}\tilde{D}_{t-1}^{T})\nabla
f(x_{t})\|}$ $\displaystyle=\frac{(I-\tilde{D}_{t-1}\tilde{D}_{t-1}^{T})\nabla
f(x_{t})}{\|(I-\tilde{D}_{t-1}\tilde{D}_{t-1}^{T})\nabla f(x_{t})\|},$

where the second equality is obtained using the orthogonality of
$\tilde{D}_{t-1}$. This corresponds to a normalized orthogonal projection of
the vector $\nabla f(x_{t})$ onto the orthogonal columns span of
$\tilde{D}_{t-1}$. Since

$\tilde{D}_{t-1}(I-\tilde{D}_{t-1}\tilde{D}_{t-1}^{T})\nabla f(x_{t})=0,$

the matrix $D_{t}=[\tilde{D}_{t-1}v]$ is orthogonal:

$D_{t}^{T}D_{t}=\begin{bmatrix}I&0\\\ 0&1\end{bmatrix}.$

Finally, the matrix $D_{t}$ indeed spans $\nabla f(x_{t})$, since $P_{t}\nabla
f(x_{t})=\nabla f(x_{t})$:

$\displaystyle P_{t}\nabla f(x)=D_{t}D_{t}^{T}\nabla f(x_{t})$
$\displaystyle=\tilde{D}_{t-1}\tilde{D}_{t-1}^{T}\nabla f(x_{t})+vv^{T}\nabla
f(x_{t})$ $\displaystyle=\tilde{D}_{t-1}\tilde{D}_{t-1}^{T}\nabla
f(x_{t})+\frac{(I-\tilde{D}_{t-1}\tilde{D}_{t-1}^{T})\nabla f(x_{t})\nabla
f(x_{t})^{T}(I-\tilde{D}_{t-1}\tilde{D}_{t-1}^{T})}{\|(I-\tilde{D}_{t-1}\tilde{D}_{t-1}^{T})\nabla
f(x_{t})\|^{2}}\nabla f(x_{t})$
$\displaystyle=\tilde{D}_{t-1}\tilde{D}_{t-1}^{T}\nabla
f(x_{t})+(I-\tilde{D}_{t-1}\tilde{D}_{t-1}^{T})\nabla f(x_{t})\frac{\nabla
f(x_{t})^{T}(I-\tilde{D}_{t-1}\tilde{D}_{t-1}^{T})^{2}\nabla
f(x_{t})}{\|(I-\tilde{D}_{t-1}\tilde{D}_{t-1}^{T})\nabla f(x_{t})\|^{2}}$
$\displaystyle=\tilde{D}_{t-1}\tilde{D}_{t-1}^{T}\nabla
f(x_{t})+(I-\tilde{D}_{t-1}\tilde{D}_{t-1}^{T})\nabla f(x_{t})$
$\displaystyle=\nabla f(x_{t}).$

∎

See 4

###### Proof.

Using eq. 40, at each iteration, after the while loop, the first-order
condition of the subroutine algorithm 4 reads

$\displaystyle D_{t}^{T}\nabla
f(x_{t})+H_{t}\alpha_{t+1}+\frac{M_{t+1}}{2}D_{t}^{T}D_{t}\alpha_{t+1}\|D_{t}\alpha_{t+1}\|=0.$
(78)

The subscript $t$ is dropped for clarity. After multiplying by $\alpha$,

$\nabla f(x_{t})^{T}D\alpha+\alpha^{T}H\alpha+\frac{M}{2}\|D\alpha\|^{3}=0.$

In addition, multiplying both times by $\alpha$ the second-order condition eq.
41 gives

$\alpha^{T}H\alpha\geq-\frac{M}{2}\|D\alpha\|^{3}.$

which gives, after replacing it in eq. 78,

$\displaystyle\nabla
f(x_{t})^{T}D\alpha\leq-\frac{M}{2}\|D\alpha\|^{3}+\frac{M}{2}\|D\alpha\|^{3}=0.$
(79)

Injecting eqs. 78 and 79 into the while condition of algorithm 4 gives the
desired result:

$\displaystyle f(x_{+})$ $\displaystyle\leq f(x)+\nabla
f(x)^{T}D\alpha+\frac{1}{2}\alpha^{T}H\alpha+\frac{M\|D\alpha\|^{3}}{6},$ (80)
$\displaystyle=f(x)-\frac{1}{2}\nabla
f(x)^{T}D\alpha-\frac{M\|D\alpha\|^{3}}{12}$ $\displaystyle\leq
f(x)-\frac{M\|D\alpha\|^{3}}{12}.$

Where eq. 80 is guaranteed if $M>L$. Therefore, in the worst case, $M<2L$.
Finally, after $t$ iterations, the number of total gradient calls is bounded
by $2t+\log_{2}\left(\frac{M_{0}}{L}\right)$ as shown in [61]. ∎

See 5

###### Proof.

The starting inequality is eq. 51:

$\|\nabla
f(x_{+})\|\leq\frac{L+M}{2}\|D\alpha\|^{2}+\|D\alpha\|\left(\frac{\|\varepsilon\|}{\|D\|}\left(\frac{L+M\kappa_{D}}{2}\right)\kappa_{D}+\|(I-P)\nabla^{2}f(x)P\|\right).$

The result is obtained by decomposing the inequality using a maximum,

$\displaystyle\|\nabla f(x_{+})\|$
$\displaystyle\leq\max\left\\{(L+M)\|D\alpha\|^{2}\;;\;2\|D\alpha\|\left(\frac{\|\varepsilon\|}{\|D\|}\left(\frac{L+M\kappa_{D}}{2}\right)\kappa_{D}+\|(I-P)\nabla^{2}f(x)P\|\right)\right\\}.$

In the first case,

$\|D\alpha\|\geq\sqrt{\frac{\|\nabla f(x_{+})\|}{L+M}},$ (81)

while in the second case,

$\|D\alpha\|\geq\frac{\|\nabla
f(x_{+})\|}{\frac{\|\varepsilon\|}{\|D\|}\left(\frac{L+M\kappa_{D}}{2}\right)\kappa_{D}+\|(I-P)\nabla^{2}f(x)P\|}.$

Let $C_{t}$ be defined as

$C_{t}=\frac{\|\varepsilon_{t}\|}{\|D_{t}\|}\left(\frac{L+M_{t+1}\kappa_{D_{t}}}{2}\right)\kappa_{D_{t}}+\|(I-P_{t})\nabla^{2}f(x_{t})P_{t}\|.$

Then, using Requirements 2 and 3, and since $M<2L$ by Theorem 4,

$C_{t}\leq C=\delta
L\left(\frac{1+2\kappa}{2}\right)\kappa+\max_{t}\|(I-P_{t})\nabla^{2}f(x_{t})P_{t}\|$

Therefore,

$\|D\alpha\|\geq\frac{\|\nabla f(x_{+})\|}{C}.$ (82)

At each iteration $t$, combining eqs. 81 and 82 into Theorem 4 gives

$f(x_{t})-f(x_{t+1})\geq\frac{M_{t+1}}{12}\|\underbrace{x_{t+1}-x_{t}}_{=D_{t}\alpha_{t}}\|^{3}\geq\frac{M_{t+1}}{12}\min\left\\{\left(\frac{\|\nabla
f(x_{+})\|}{L+M_{t+1}}\right)^{3/2}\;;\;\left(\frac{\|\nabla
f(x_{+})\|}{C}\right)^{3}\right\\}$

Therefore,

$\displaystyle f(x_{0})-f^{\star}$ $\displaystyle\geq f(x_{0})-f(x_{t})$
$\displaystyle=\sum_{i=0}^{t-1}f(x_{i})-f(x_{i+1})$
$\displaystyle\geq\sum_{i=0}^{t-1}\left(\frac{M_{i+1}}{12}\|x_{i+1}-x_{i}\|^{3}\right)$
$\displaystyle\geq\sum_{i=0}^{t-1}\min_{t}\frac{M_{i+1}}{12}\left\\{\left(\frac{\|\nabla
f(x_{i+1})\|}{L+M_{i+1}}\right)^{3/2}\;;\;\left(\frac{\|\nabla
f(x_{i+1})\|}{C}\right)^{3}\right\\}$ $\displaystyle\geq
t\min_{i\in[0,t-1]}\frac{M_{i+1}}{12}\min\left\\{\left(\frac{\|\nabla
f(x_{i+1})\|}{L+M_{i+1}}\right)^{3/2}\;;\;\left(\frac{\|\nabla
f(x_{i+1})\|}{C}\right)^{3}\right\\}$ $\displaystyle\geq
t\frac{M_{\min}}{12}\min\left\\{\min_{i\in[1,t]}\left(\frac{\|\nabla
f(x_{i})\|}{3L}\right)^{3/2}\;;\;\min_{i\in[1,t]}\left(\frac{\|\nabla
f(x_{i})\|}{C}\right)^{3}\right\\}$

After analyzing separately each case of the minimum, either

$\left(\frac{\min\limits_{i\in[1,t]}\|\nabla f(x_{i})\|}{3L}\right)^{3/2}\leq
12\frac{f(x_{0})-f^{\star}}{tM_{\min}}\quad\text{or}\quad\left(\frac{\min\limits_{i\in[1,t]}\|\nabla
f(x_{t+1})\|}{C}\right)^{3}\leq 12\frac{f(x_{0})-f^{\star}}{tM_{\min}}.$

It remains to simplify to obtain the desired result,

$\min_{i=1\ldots t}\|\nabla
f(x_{i})\|\leq\max\left\\{\frac{3L}{t^{2/3}}\left(12\frac{f(x_{0})-f^{\star}}{M_{\min}}\right)^{2/3}\;;\;\left(\frac{C}{t^{1/3}}\right)\left(12\frac{f(x_{0})-f^{\star}}{M_{\min}}\right)^{1/3}\right\\}.$

∎

See 6

###### Proof.

Starting from the inequality in proposition 8,

$\displaystyle f(x_{t+1})\leq
f(y)+\frac{M_{t+1}+L}{6}\|y-x_{t}\|^{3}+\frac{\|y-x_{t}\|^{2}}{2}C_{2}^{(t)},$

where

$C_{2}^{(t)}=\|\nabla^{2}f(x_{t})-P_{t}\nabla^{2}f(x_{t})P_{t}\|+\delta\frac{L\kappa+M_{t+1}\kappa^{2}}{2},$

and setting $y=(1-\beta_{t})x_{t}+\beta_{t}x^{\star}$ and
$f(x^{\star})=f^{\star}$ gives

$\displaystyle f(x_{t+1})-f^{\star}\leq
f((1-\beta_{t})x_{t}+\beta_{t}x^{\star})-f^{\star}+\frac{M_{t+1}+L}{6}\beta_{t}^{3}\|x_{t}-x^{\star}\|^{3}+\frac{\beta_{t}^{2}\|x_{t}-x^{\star}\|^{2}}{2}C_{2}^{(t)}.$

Because the function is star-convex,

$\displaystyle
f(x_{t+1})-f^{\star}\leq(1-\beta_{t})(f(x_{t})-f^{\star})+\frac{M_{t+1}+L}{6}\beta_{t}^{3}\|x_{t}-x^{\star}\|^{3}+\frac{\beta_{t}^{2}\|x_{t}-x^{\star}\|^{2}}{2}C_{2}^{(t)}.$

Since algorithm 4 ensure a decrease in the function value, the iterate $x_{t}$
satisfies

$x_{t}\in\\{x:f(x\leq f(x_{0}))\\},$

and therefore, $\|x_{t}-x^{\star}\|\leq R$ by Assumption 2. In addition,
$M<2L$ by Theorem 4. The inequality now becomes

$\displaystyle(f(x_{t+1})-f^{\star})\leq(1-\beta_{t})(f(x_{t})-f^{\star})+\beta_{t}^{3}\frac{LR^{3}}{2}+\beta_{t}^{2}\frac{R^{2}C_{2}^{(t)}}{2}.$
(83)

Finally, since $M<2L$, the scalar $C_{2}^{t}$ is bounded over time by $C_{2}$:

$C_{2}^{(t)}\leq C_{2}\stackrel{{\scriptstyle\text{def}}}{{=}}\delta
L\frac{\kappa+2\kappa^{2}}{2}+\max_{t}\|\nabla^{2}f(x_{t})-P_{t}\nabla^{2}f(x_{t})P_{t}\|.$

Now, let

  * •

$B_{t}=\frac{t(t+1)(t+2)}{6}$,

  * •

$b_{t}:B_{t}=B_{t-1}+b_{t}$, hence $b_{t}=\frac{t(t+1)}{2}$, and

  * •

$\beta_{t}=\frac{b_{t+1}}{B_{t+1}}$.

Therefore, for $t\geq 1$,

$1=\frac{B_{t}}{B_{t}}=\frac{B_{t-1}}{B_{t}}+\frac{b_{t}}{B_{t}}=\frac{B_{t-1}}{B_{t}}+\beta_{t-1}\quad\Rightarrow\quad
1-\beta_{t-1}=\frac{B_{t-1}}{B_{t}}.$

Injecting those relations in eq. 83 gives

$\displaystyle(f(x_{t+1})-f^{\star})\leq\frac{B_{t}}{B_{t+1}}(f(x_{t})-f^{\star})+\left(\frac{b_{t+1}}{B_{t+1}}\right)^{3}\frac{LR^{3}}{2}+\left(\frac{b_{t+1}}{B_{t+1}}\right)^{2}\frac{R^{2}C_{2}}{2},$

hence the recursion

$\displaystyle B_{t+1}(f(x_{t+1})-f^{\star})$ $\displaystyle\leq
B_{t}(f(x_{t})-f^{\star})+\frac{b^{3}_{t+1}}{B^{2}_{t+1}}\frac{LR^{3}}{2}+\frac{b^{2}_{t+1}}{B_{t+1}}\frac{R^{2}C_{2}}{2}$
$\displaystyle\leq
B_{0}(f(x_{t})-f^{\star})+\sum_{i=0}^{t}\frac{b^{3}_{i+1}}{B^{2}_{i+1}}\frac{LR^{3}}{2}+\sum_{i=0}^{t}\frac{b^{2}_{i+1}}{B_{i+1}}\frac{R^{2}C_{2}}{2}.$
$\displaystyle(f(x_{t+1})-f^{\star})\leq\frac{B_{0}}{B_{t+1}}(f(x_{t})-f^{\star})+\frac{\sum_{i=0}^{t}\frac{b^{3}_{i+1}}{B^{2}_{i+1}}}{B_{t+1}}\frac{LR^{3}}{2}+\frac{\sum_{i=0}^{t}\frac{b^{2}_{i+1}}{B_{i+1}}}{B_{t+1}}\frac{R^{2}C_{2}}{2}.$

Therefore, the rate reads By the definition of $b_{t}$ and $B_{t}$,

$\displaystyle\frac{b^{3}_{i+1}}{B^{2}_{i+1}}$
$\displaystyle=\frac{36}{8}\frac{(i+1)^{3}(i+2)^{3}}{(i+1)^{2}(i+2)^{2}(i+3)^{2}}=\frac{9}{2}\frac{(i+1)(i+2)}{(i+3)^{2}}\leq\frac{9}{2},$
$\displaystyle\frac{b^{2}_{i+1}}{B_{i+1}}$
$\displaystyle=\frac{6}{4}\frac{(i+1)^{2}(i+2)^{2}}{(i+1)(i+2)(i+3)}=\frac{3}{2}\frac{(i+2)}{(i+3)}(i+1)\leq\frac{3}{2}(i+1).$

Hence,

$\displaystyle\frac{\sum_{i=0}^{t}\frac{b^{3}_{i+1}}{B^{2}_{i+1}}}{B_{t+1}}$
$\displaystyle\leq\frac{\frac{9}{2}(t+1)}{\frac{(t+1)(t+2)(t+3)}{6}}\leq\frac{27}{(t+2)(t+3)},$
$\displaystyle\frac{\sum_{i=0}^{t}\frac{b^{2}_{i+1}}{B_{i+1}}}{B_{t+1}}$
$\displaystyle\leq\frac{\sum_{i=0}^{t}\frac{3}{2}(i+1)}{\frac{(t+1)(t+2)(t+3)}{6}}=\frac{\frac{3}{4}(t+2)(t+1)}{\frac{(t+1)(t+2)(t+3)}{6}}=\frac{9}{2(t+3)}.$

Shifting from $t+1$ tp $t$ gives the desired result,

$(f(x_{t})-f^{\star})\leq
6\frac{f(x_{t})-f^{\star}}{t(t+1)(t+2)}+\frac{1}{(t+1)(t+2)}\frac{L(3R)^{3}}{2}+\frac{1}{t+2}\frac{C_{2}(3R)^{2}}{4}.$

∎

See 7

###### Proof.

The proof technique is similar to [44]. Starting from proposition 9 with
$x=x_{t}$,

$\displaystyle\mathbb{E}f(x_{t+1})\leq$
$\displaystyle\left(1-\frac{N}{d}\right)f(x_{t})+\frac{N}{d}f(y)+\frac{N}{d}\frac{(M_{t+1}+L)}{6}\|y-x_{t}\|^{3}$
$\displaystyle+\frac{N}{d}\frac{\|y-x_{t}\|^{2}}{2}\left(\delta\frac{L\kappa+M_{t+1}\kappa^{2}}{2}+\frac{(d-N)}{d}\|\nabla^{2}f(x_{t})\|\right),$

where the expectation is taken with $D_{0},\ldots,D_{t-1}$ fixed. Using the
inequality $M_{t+1}\leq 2L$ gives

$\displaystyle\mathbb{E}f(x_{t+1})\leq$
$\displaystyle\left(1-\frac{N}{d}\right)f(x_{t})+\frac{N}{d}\left(f(y)+\frac{\|y-x_{t}\|^{2}}{2}C_{3}+\frac{L}{2}\|y-x_{t}\|^{3}\right)$

where

$C_{3}\stackrel{{\scriptstyle\text{def}}}{{=}}\left(\delta
L\frac{\kappa+2\kappa^{2}}{2}+\frac{(d-N)}{d}\max_{i\in[0,t]}\|\nabla^{2}f(x_{i})\|\right).$

Let $y=\beta_{t}x^{\star}+(1-\beta_{t})x_{t}$, $\beta_{t}\in[0,1]$. After
using Assumption 4 and Assumption 2,

$\displaystyle\mathbb{E}f(x_{t+1})$
$\displaystyle\leq\left(1-\frac{N}{d}\right)f(x_{t})+\frac{N}{d}\left(f\Big{(}\beta_{t}x^{\star}+(1-\beta_{t})x_{t}\big{)}+\beta_{t}^{2}\frac{C_{3}R^{2}}{2}+\beta_{t}^{3}\frac{LR^{3}}{2}\right)$
$\displaystyle\leq\left(1-\frac{N}{d}\right)f(x_{t})+\frac{N}{d}\left(\beta_{t}f(x^{\star})+(1-\beta_{t})f(x_{t})+\beta_{t}^{2}\frac{C_{3}R^{2}}{2}+\beta_{t}^{3}\frac{LR^{3}}{2}\right)$
$\displaystyle=\left(1-\frac{N}{d}\right)f(x_{t})+\frac{N}{d}\left(\beta_{t}f(x^{\star})+(1-\beta_{t})f(x_{t})+\beta_{t}^{2}\frac{C_{3}R^{2}}{2}+\beta_{t}^{3}\frac{LR^{3}}{2}\right),$
$\displaystyle=\left(1-\beta_{t}\frac{N}{d}\right)f(x_{t})+\frac{N}{d}\left(\beta_{t}f(x^{\star})+\beta_{t}^{2}\frac{C_{3}R^{2}}{2}+\beta_{t}^{3}\frac{LR^{3}}{2}\right).$

Hence, the recursion

$\left(\mathbb{E}f(x_{t+1})-f^{\star}\right)\leq\left(1-\beta_{t}\frac{N}{d}\right)(f(x_{t})-f^{\star})+\frac{N}{d}\left(\beta_{t}^{2}\frac{C_{3}R^{2}}{2}+\beta_{t}^{3}\frac{LR^{3}}{2}\right).$

Now, define

$\displaystyle b_{t}$ $\displaystyle=t^{2},$ $\displaystyle B_{t}$
$\displaystyle=B_{0}+\sum_{i=0}^{t}b_{i},\;\;B_{0}=\frac{4}{3}\left(\frac{d}{N}\right)^{3}$
$\displaystyle\beta_{t}$
$\displaystyle=\frac{d}{N}\frac{b_{t+1}}{B_{t+1}}\;\;\Rightarrow\;1-\frac{N}{d}\beta_{t}=\frac{B_{t}}{B_{t+1}}.$

Replacing those relations in the recursion gives

$\displaystyle B_{t+1}\left(\mathbb{E}f(x_{t+1})-f^{\star}\right)$
$\displaystyle\leq$ $\displaystyle
B_{t}(f(x_{t})-f^{\star})+\frac{N}{dB_{t+1}}\left(\left(\frac{d}{N}\frac{b_{t+1}}{B_{t+1}}\right)^{2}\frac{C_{3}R^{2}}{2}+\left(\frac{d}{N}\frac{b_{t+1}}{B_{t+1}}\right)^{3}\frac{LR^{3}}{2}\right)$
$\displaystyle=$ $\displaystyle
B_{t}(f(x_{t})-f^{\star})+\frac{d}{N}\frac{b^{2}_{t+1}}{B_{t+1}}\frac{C_{3}R^{2}}{2}+\frac{d^{2}}{N^{2}}\frac{b^{3}_{t+1}}{B^{2}_{t+1}}\frac{LR^{3}}{2}$

Expanding the inequality gives

$\displaystyle B_{t+1}\left(\mathbb{E}f(x_{t+1})-f^{\star}\right)\leq
B_{0}(f(x_{0})-f^{\star})+\frac{d}{N}\sum_{t=0}^{t+1}\frac{b^{2}_{i+1}}{B_{i+1}}\frac{C_{3}R^{2}}{2}+\frac{d^{2}}{N^{2}}\sum_{t=0}^{t+1}\frac{b^{3}_{i+1}}{B^{2}_{i+1}}\frac{LR^{3}}{2}$

Since

$\displaystyle B_{t}$ $\displaystyle=B_{0}+\sum_{i=1}^{t}\geq
B_{0}+\int_{0}^{t}x^{2}\mathrm{d}x=B_{0}+\frac{t^{3}}{3}$
$\displaystyle\sum_{i=0}^{t}\frac{b_{t}^{2}}{B_{t}}$
$\displaystyle\leq\sum_{i=0}^{t}\frac{i^{4}}{B_{0}+i^{3}/3}\leq 3t^{2},$
$\displaystyle\sum_{i=0}^{t}\frac{b_{t}^{3}}{B^{2}_{t}}$
$\displaystyle\leq\sum_{i=0}^{t}\frac{i^{6}}{(B_{0}+i^{3}/3)^{2}}\leq 9t,$

the bound becomes

$\displaystyle B_{t+1}\left(\mathbb{E}f(x_{t+1})-f^{\star}\right)\leq
B_{0}(f(x_{0})-f^{\star})+\frac{d}{N}3t^{2}\frac{C_{3}R^{2}}{2}+\frac{d^{2}}{N^{2}}9t\frac{LR^{3}}{2}$

Dividing both sides by $B_{t+1}$ gives

$\displaystyle\mathbb{E}f(x_{t+1})-f^{\star}\leq\frac{B_{0}}{B_{0}+\frac{(t+1)^{3}}{3}}(f(x_{0})-f^{\star})+\frac{d}{N}\frac{3(t+1)^{2}}{B_{0}+\frac{(t+1)^{3}}{3}}\frac{C_{3}R^{2}}{2}+\frac{d^{2}}{N^{2}}\frac{9(t+1)}{B_{0}+\frac{(t+1)^{3}}{3}}\frac{LR^{3}}{2}.$

After the following simplifications,

$\displaystyle\frac{B_{0}}{B_{0}+(t+1)^{3}/3}$
$\displaystyle=\frac{1}{1+\frac{(t+1)^{3}}{3B_{0}}}=\frac{1}{1+\frac{1}{4}\left(\frac{N}{d}(t+1)\right)^{3}},$
$\displaystyle\frac{3(t+1)^{2}}{B_{0}+(t+1)^{3}/3}$
$\displaystyle=\frac{3}{B_{0}}\frac{(t+1)^{3}}{1+\frac{(t+1)^{3}}{3B_{0}}}\frac{1}{t+1}\leq\frac{3}{B_{0}}3B_{0}\frac{1}{t+1}=\frac{9}{t+1},$
$\displaystyle\frac{9(t+1)}{B_{0}+\frac{(t+1)^{3}}{3}}$
$\displaystyle=\frac{9}{B_{0}}\frac{(t+1)^{3}}{\frac{(t+1)^{3}}{3B_{0}}}\frac{1}{(t+1)^{2}}\leq\frac{9}{B_{0}}3B_{0}\frac{1}{(t+1)^{2}}=\frac{27}{(t+1)^{2}},$

the inequality finally becomes (after shifting from $t+1$ to $t$),

$\displaystyle\mathbb{E}f(x_{t})-f^{\star}\leq\frac{1}{1+\frac{1}{4}\left[\frac{N}{d}t\right]^{3}}(f(x_{0})-f^{\star})+\frac{1}{\left[\frac{N}{d}t\right]^{2}}\frac{L(3R)^{3}}{2}+\frac{1}{\left[\frac{N}{d}t\right]}\frac{C_{3}(3R)^{2}}{2}.$

∎

See 8

###### Proof.

By construction of $\phi_{t}(x)$, from proposition 14 and Assumption 2,

$\displaystyle B_{t}f(x_{t})$ $\displaystyle\leq\min_{x}\phi_{t}(x)$ (84)
$\displaystyle\leq\phi_{t}(x^{\star})$ (85) $\displaystyle\leq
B_{t}f(x^{\star})+\frac{\lambda_{t}^{(1)}+\tilde{\lambda}^{(1)}}{2}\|x^{\star}-x_{0}\|^{2}+\frac{\lambda_{t}^{(2)}+\tilde{\lambda}^{(2)}}{6}\|x^{\star}-x_{0}\|^{3}$
(86) $\displaystyle\leq
B_{t}f(x^{\star})+\frac{\lambda_{t}^{(1)}+\tilde{\lambda}^{(1)}}{2}R^{2}+\frac{\lambda_{t}^{(2)}+\tilde{\lambda}^{(2)}}{6}R^{3}$
(87) $\displaystyle\Rightarrow f(x_{t})-f^{\star}$
$\displaystyle\leq\frac{\lambda_{t}^{(1)}+\tilde{\lambda}^{(1)}}{2B_{t}}R^{2}+\frac{\lambda_{t}^{(2)}+\tilde{\lambda}^{(2)}}{6B_{t}}R^{3}.$
(88)

By proposition 15, the following bounds holds:

$\displaystyle\lambda_{t}^{(1)}$ $\displaystyle\leq
30\cdot\frac{b_{t+1}^{2}}{B_{t}}\kappa_{D}\left(\delta\max\\{4L,M_{0}\\}+\max_{i=0\ldots
t}\|(I-P_{i})\nabla f(x_{i})P_{i})\|\right),$ $\displaystyle\lambda_{t}^{(2)}$
$\displaystyle\leq 5\frac{b_{t+1}^{3}}{B^{2}_{t}}\max\\{M_{0},4L\\}.$

Since $\frac{b_{t+1}}{B_{t}}=\frac{3}{(t+3)}$,

$\displaystyle\frac{b_{t+1}^{3}}{B^{3}_{t}}=\frac{3^{3}}{(t+3)^{3}},\qquad\frac{b_{t+1}^{2}}{B^{2}_{t}}=\frac{3^{2}}{(t+3)^{2}}.$
(89)

Therefore,

$\displaystyle f(x_{t})-f^{\star}\leq$ $\displaystyle
30\cdot\kappa_{D}\left(\delta\max\\{4L,M_{0}\\}+\max_{i=0\ldots
t}\|(I-P_{i})\nabla f(x_{i})P_{i})\|\right)\frac{(3R)^{2}}{(t+3)^{2}}$
$\displaystyle+5\max\\{M_{0},4L\\}\left(\frac{3R}{t+3}\right)^{3}$
$\displaystyle+\frac{\frac{\tilde{\lambda}^{(1)}R^{2}}{2}+\frac{\tilde{\lambda}^{(2)}R^{3}}{6}}{(t+1)^{3}}.$

∎


# Continuous max-flow augmentation of self-supervised few-shot learning on
SPECT left ventricles

###### Abstract

Single-Photon Emission Computed Tomography (SPECT) left ventricular assessment
protocols are important for detecting ischemia in high-risk patients. To
quantitatively measure myocardial function, clinicians depend on commercially
available solutions to segment and reorient the left ventricle (LV) for
evaluation. Based on large normal datasets, the segmentation performance and
the high price of these solutions can hinder the availability of reliable and
precise localization of the LV delineation. To overcome the aforementioned
shortcomings this paper aims to give a recipe for diagnostic centers as well
as for clinics to automatically segment the myocardium based on small and low-
quality labels on reconstructed SPECT, complete field-of-view (FOV) volumes. A
combination of Continuous Max-Flow (CMF) with prior shape information is
developed to augment the 3D U-Net self-supervised learning (SSL) approach on
various geometries of SPECT apparatus. Experimental results on the acquired
dataset have shown a 5-10% increase in quantitative metrics based on the
previous State-of-the-Art (SOTA) solutions, suggesting a good plausible way to
tackle the few-shot SSL problem on high-noise SPECT cardiac datasets.

Index Terms—  Cardiac SPECT, Myocardial Perfusion Imaging, 3D U-net, CNN,
Augmentation

##  1 Introduction

Myocardial Perfusion Imaging (MPI) is aimed at getting results on hypertrophic
cardiomyopathy and infarction with high specificity but at rather low
sensitivity [1]. This functional imaging technique with various
radiopharmaceuticals is aimed at expressing metabolic tissue-level processes.
As a side effect, implicitly purports the myocardial blood flow to discover
ischemia in patients, since the administration of the pharmaceutical is
injected intravenously. The MPI protocol is standardized by the European
Association of Nuclear Medicine (EANM) and the American Society of Nuclear
Cardiology (ASNC) by specifying the pharmaceutical administration as well as
the stress or rest testing in different scenarios. However, with the wide
range of imaging companies and equipment of SPECT cameras, assuring image
quality remains a challenging task up to this day. Combining this high
variability of acquisition with lowered doses in MPI, image processing becomes
increasingly hard, especially segmentation.

To address and solve the aforementioned problems of myocardium segmentation in
SPECT MPI, different machine learning algorithms have been developed [2, 3].
Specifically convolutional neural networks (CNN) are one of the best
candidates to segment the left ventricle [4]. Supervised training was
conducted on cropped cardiac images with a 3D V-Net architecture by [5]
utilizing binary cross entropy (BCE) with Dice similarity, resulting in a
compound loss function. By taking the already great results further, Zhang et
al. [6] proposed a complete FOV-based spatial transformer network (STN) with
3D CNN to solve the segmentation and reorientation problem as well. Alongside
the 3D CNN, another 2D U-Net-based network was proposed in [7], which is a
good improvement over a dynamic programming (DP), and comparable performance
was shown to the QGS software package [8]. A different approach combining STN
and CNN was used on whole image segmentation and reorientation, rendering
great improvement over the latter approaches. And finally, the current SOTA
results are conducted by [9], combining 3D V-Net with shape priors on left
ventricles, and in [10] multi-scale 3D U-Net with STN segmentation is used to
segment complete FOV volumes with high performance.

The literature is rich in the application and development of various learning
strategies and architectures for the segmentation problem in complete FOV –
whole reconstructed volume – segmentation, however, it becomes quite narrow,
when one approaches the problem with small patient data at hand, namely with
few-shot learning (FSL) or unsupervision. In this study small amount of
labeled dataset ($<$10 patients) is considered with various conditions, where
the predictions are tested on different imaging geometries and dosages on
multiple patients.

##  2 Methods

In this paper, to segment the left ventricle on myocardial perfusion imaging
SPECT, the complete reconstructed volume – whole image – is considered. The 3D
volume of the scan is denoted $\Omega(x)$ as a closed domain. As the model of
segmentation is thought of as continuous flow in three dimensions [11, 12] the
source and sink terminal vertices are denoted as $s$ and $t$. At every
$x\in\Omega$, $p(x)$ is the passing, $p_{s}(x)$ the directed flow from source
to the current position, $p_{t}(x)$ is the directed sink flow from $x$ to $t$.

###  2.1 3D U-Net SSL with augmentation

The main backbone of the segmentation is 3D U-Net architecture [13] with the
following setup. The encoder part of the network consisted of three layers
with Instance Normalization (IN) and PReLU activation functions. The various
sized inputs were upsampled to $128^{3}$ sizes and were fed to the encoder
part in $B\times C\times D\times H\times W$ order, B is batch, C is channel, H
is height and W is the width. It should be noted that during the training of
the encoder part of the network, the decoder was replaced by two dense layers
with PReLU activation functions with Batch Normalization (BN) layer between
with a dropout of $1/2$.

With very small labeled data at hand, the algorithm used heavy augmentation to
reach the desired performance. For both stages of the self-supervised learning
and supervised fine-tuning the augmentation pipeline consisted of the
following processes, random blur, Gaussian noise, random affine, and elastic
transformations with large parameter variations. These steps were applied in
distinct modules in the pipeline and were obfuscated based on the best
probabilities. However, after further testing, it was shown, that the network
needed more augmentation to perform well. As in [14], to improve segmentation
performance, the network had to solve jigsaw puzzles. The complete abdomen
reconstructed 3D SPECT volumes were disassembled into tiles and randomly
shuffled given as a task to solve for the network. The get the best out of
this Representation Learning (RL) technique a large number of permutations,
1000 were used to perform the self-supervised pretext task. From the
precomputed permutations the network had to choose the ones with maximal
Hamming distances [14].

After the training of the encoder, it was loaded into a new network and the
decoder was optimized by supervised-fine tuning on the 10 patient labeled
dataset. The model utilized AdamW [15], with LogCosh loss function and a
learning rate of 1e-3 with decay by a factor of $0.01$.

###  2.2 Continuous Max-Flow enhancement

The CNN 3D U-Net SSL [16] will generate a prediction map, the probability of
each pixel belonging the the myocardium. This map is denoted by $\lambda(x),\
x\in\Omega$, where $\lambda(x)\in[0,1]$ and $\lambda\in BV(\Omega)$, where
$BV$ denotes the space of functions with bounded variation. This renders the
original Mumford-Shah model [17] convex based on the Truncation lemma. The
maximum flow problem is the dual of the minimal cut problem, which can be
formalized as the following

$\begin{split}L_{c}(p_{s},p_{t},p,\lambda)&=\int_{\Omega}p_{s}dx\\\
&+\int_{\Omega}\lambda(\nabla\cdot p-p_{s}+p_{t})dx\\\
&-\frac{c}{2}\left\lVert\nabla\cdot p-p_{s}-p_{t}\right\rVert,\end{split}$ (1)

where $c>0$ and the solution is using the method by [18] and [19] as

$p^{k+1}=\Pi_{\alpha}\left(p^{k}+c\nabla(\nabla\cdot p^{k}-F^{k})\right),$ (2)

where $\Pi_{\alpha}$ is a projection onto the convex set $C_{\alpha}=\\{q\ |\
\left\lVert q\right\rVert\leq\alpha\\}$.

During the investigations on the test set, the Signal-to-Noise-Ratio rendered
so low on different 201TI Chloride and some 99mTc MIBI studies that the
optimization got stuck in a local optima. To increase the performance of the
scheme, a shape prior is introduced in the CMF method, with the aid of density
estimation techniques.

###  2.3 Shape prior incorporation with density estimation

Applying the 3D U-Net SSL approach on relatively normal patients showed good
results in the evaluation. However, on hypoperfused myocardium, the approach
was suboptimal as it was measured. Even increasing the dataset during the
self-supervised pretext task the network was unable to segment dilated hearts.

To overcome these challenges, the Continuous Max-Flow enhancement incorporated
prior information about the cardiac model. Since SPECT-based functional
imaging varies a lot even in the same patient, with identical pharmacons at
different times of the day, the analytical approach of the cardiac model is
challenging. To relax the task of cardiac modeling of functional imaging,
statistical evidence can help to build prior information about the different
left ventricles.

A similar approach is taken as in [20], where the outlined training data is
encoded as a Gaussian distribution in higher dimensions. Given the negative
logarithm of the probability, the Mahalanobis-type distance-based shape loss
can be incorporated. Let $z$ be points in the higher dimensional space, then

$-\log{P(z)}=\frac{1}{2}(z-\mu)^{t}\Sigma^{-1}_{\perp}\frac{1}{2}(z-\mu),$ (3)

where $\mu$ is the mean shape and $\Sigma^{-1}_{\perp}$ is the inverse of the
regularized covariance matrix as in [20]. This distance can act as a cost
function in a probabilistic setting.

The problem with the latter approach is the assumption of the Gaussian
approximation of the training shapes, which renders it obsolete for SPECT
segmentation. To handle the non-gaussian nature and the complicated shape
deformations of the training surfaces, the proposed model utilizes kernel
space techniques [21], to give a good density estimation on prior variation.

##  3 Materials

###  3.1 Data

Simulated data consisted of mathematical phantoms[22], whereas patient data
consisted of 102 patient cases with different myocardial states in the
perfusion protocol (non-gated) were evaluated. We created 3 different groups
for various tests for our automation. (A) Perfusion studies with 53 patients,
23 females, and 30 males with an average age of 68.2 years, a body weight of
75.23, and a height of 164 cm. (B) Attenuation corrected perfusion studies
consisted of 51 patients, 35 females, 16 males with an average age of 63.92
years, a body weight of 79.31 kgs and height of 166.62 cm. (C) Motion artifact
perfusion studies of 36 patients, 19 female, 17 male with an average age of
53.97 years, the weight of 81 kgs, and height of 165.3 cm.

###  3.2 Acquisition

Patient data acquisition was performed on different versions of Mediso
AnyScan, AnyScan Trio SPECT, and SPECT/CT systems. The CT Image scanning
parameters: SPIRAL acquisition type; CT Dose Index Volume 0.22; tube voltage
150 kVp; and tube current 30-50 mAs. Myocardial Perfusion SPECT (MPS) studies
were performed at 4 participating sites with various isotopes 99mTc MIBI,
99mTc Tetrofosmin and 201TI Chloride. Acquisitions were performed with both 64
and 128 matrices, with 64 projections on a two-headed camera and 96
projections on a triple-headed camera. The angular interval was 180 on both
dual and triple-headed cameras. Studies with attenuation maps were
reconstructed with the TeraTomo method. Samples without attenuation correction
(AC) utilized OSEM reconstruction method [23] to create the complete FOV
volume stack.

##  4 Results

Averaged performance metrics
---
|  Precision |  Recall |  IoU |  Dice score
Proposed model |  0.7034 |  0.7286 |  0.5423 |  0.7065
3D U-Net SSL |  0.5558 |  0.8551 |  0.5079 |  0.6737
Zhu et. al |  0.0393 |  0.8456 |  0.0385 |  0.0725
Table 1: Averaged segmentation results on the patient dataset. Our approach is
capable of outperforming [9] and 3D U-Net SSL [16] as well on most of the
metrics.

###  4.1 Quantitative evaluation

The segmentation accuracy was evaluated on the following metrics. Accuracy =
$\frac{TP+TN}{TP+TN+FP+FN}$, Precision = $\frac{TP}{TP+FP}$, Recall =
$\frac{TP}{TP+FN}$, Intersection over Union (IoU) = $\frac{TP}{TP+FN+FP}$,
DICE coefficient = $\frac{2TP}{2TP+FP+FN}$, where T stands for true, F false,
N is negative, P is positive.

###  4.2 Implementation.

The code111https://github.com/JacksonFurrier/isbi_2024_code.git is written in
PyTorch [24] and models were trained on a single Nvidia A100 GPU [25] and a
1080Ti GPU as well. The annotated data consisted of 10 patients and the
complete dataset contained 102 patients reconstructed complete FOV volumes.
The implementation uses the TorchIO [2], Scikit-learn [26] package, and
extensively the GeomLoss [27] library for the Sinkhorn loss and
differentiation of cost functions.

###  4.3 Experiments

As seen in Table 1 the combined method of 3D SSL U-Net with CMF Shape prior
enhancement reached the best results in segmenting the left ventricular volume
of the complete-FOV reconstructed images. As the metric shows, the recall
score is relatively higher, which lead to the conjecture that the method is
over-predicting cardiac wall thickness. This problem can be two-fold, one
possibility that the training data, along with the model of the left ventricle
assumes thicker walls or the reconstruction method from kernel space lacks
precision. This will need further tests and investigations, with combined
oracle testing with MRI.

Knowing that the 3D SSL U-Net performance was deteriorating in hypoperfused
cardiac patients, the tests consisted of evaluation with ill-conditioned
patients. As seen on 1 the CMF enhancement outperformed the original 3D SSL
U-Net approach, however, as it was previously mentioned the “thickening” of
the cardiac wall persists.

(a)

(b)

Fig. 1: Enhancement of the 3D SSL U-Net [16] approach with CMF Shape prior
segmentation. The first column depicts the SPECT volume of the cardiac area,
the second column shows the segmentation with 3D SSL U-Net, and the third
column shows the results with CMF shape prior enhancement. Two patients with
different cardiac conditions, 1(a) with stable perfusion defect and 1(b)
patient showing signs of inferior perfusion defect.

##  5 Conclusion

In this study, a different approach is taken in “augmenting” the 3D U-Net SSL
technique for segmenting normal and heavily hypoperfused SPECT left
ventricles. The contribution lies in the combination of 3D SSL U-Net with a
Shape prior-based Continuous Max-Flow optimization. The results show that
further investigation of shape prior information involvement in deep learning
segmentation tasks could be a viable approach to overcome the issues of a few
labeled datasets.

##  6 Compliance with ethical standards

  * •

“This study was performed in line with the principles of the Declaration of
Helsinki. Approval was granted by the Ethics Committee of Semmelweis
University”

##  7 Acknowledgments

The research was supported by the project No. 2019-1.3.1-KK-2019-00011
financed by the National Research, Development and Innovation Fund of Hungary
under the Establishment of Competence Centers, Development of Research
Infrastructure Programme funding scheme. Further, on behalf of the MedisoLab
project, we are grateful for the possibility of using ELKH Cloud (see Héder et
al. 2022 [25]; https://science-cloud.hu/) which helped us achieve the results
published in this paper.

## References

  * [1] M Elizabeth Oates and Vincent Sorrell,  Myocardial Perfusion Imaging - Beyond the Left Ventricle,  Springer International Publishing, 2017.
  * [2] Ernest V Garcia,  “Artificial intelligence in nuclear cardiology: Preparing for the fifth industrial revolution,”  Journal of Nuclear Cardiology, vol. 28, pp. 1199–1202, 2021.
  * [3] Ioannis D Apostolopoulos, Nikolaos I Papandrianos, Anna Feleki, Serafeim Moustakidis, and Elpiniki I Papageorgiou,  “Deep learning-enhanced nuclear medicine spect imaging applied to cardiac studies,”  EJNMMI physics, vol. 10, no. 1, pp. 6, 2023.
  * [4] Waseem Hijazi and Robert JH Miller,  “Deep learning to automate spect mpi myocardial reorientation,”  Journal of Nuclear Cardiology, pp. 1–2, 2023.
  * [5] Tonghe Wang, Yang Lei, Haipeng Tang, Zhuo He, Richard Castillo, Cheng Wang, Dianfu Li, Kristin Higgins, Tian Liu, Walter J Curran, et al.,  “A learning-based automatic segmentation and quantification method on left ventricle in gated myocardial perfusion spect imaging: A feasibility study,”  Journal of Nuclear Cardiology, vol. 27, pp. 976–987, 2020.
  * [6] Duo Zhang, P Hendrik Pretorius, Kaixian Lin, Weibing Miao, Jingsong Li, Michael A King, and Wentao Zhu,  “A novel deep-learning–based approach for automatic reorientation of 3d cardiac spect images,”  European Journal of Nuclear Medicine and Molecular Imaging, vol. 48, pp. 3457–3468, 2021.
  * [7] Haixing Wen, Qiuyue Wei, Jin-Long Huang, Shih-Chuan Tsai, Chi-Yen Wang, Kuo-Feng Chiang, Yu Deng, Xiongtao Cui, Rui Gao, Weihua Zhou, et al.,  “Analysis on spect myocardial perfusion imaging with a tool derived from dynamic programming to deep learning,”  Optik, vol. 240, pp. 166842, 2021.
  * [8] Guido Germano, Paul B. Kavanagh, Piotr J. Slomka, Serge D. Van Kriekinge, Geoff Pollard, and Daniel S. Berman,  “Quantitation in gated perfusion SPECT imaging: The Cedars-Sinai approach,”  Journal of Nuclear Cardiology, vol. 14, no. 4 SPEC. ISS., pp. 433–454, 2007.
  * [9] Fubao Zhu, Longxi Li, Jinyu Zhao, Chen Zhao, Shaojie Tang, Jiaofen Nan, Yanting Li, Zhongqiang Zhao, Jianzhou Shi, Zenghong Chen, Chuang Han, Zhixin Jiang, and Weihua Zhou,  “A new method incorporating deep learning with shape priors for left ventricular segmentation in myocardial perfusion spect images,”  Computers in Biology and Medicine, vol. 160, pp. 106954, 2023.
  * [10] Yangfan Ni, Duo Zhang, Gege Ma, Lijun Lu, Zhongke Huang, and Wentao Zhu,  “A multi-scale spatial transformer u-net for simultaneously automatic reorientation and segmentation of 3d nuclear cardiac images,”  arXiv preprint arXiv:2310.10095, 2023.
  * [11] Gilbert Strang,  “Maximal flow through a domain,”  Mathematical Programming, vol. 26, pp. 123–143, 1983.
  * [12] Gilbert Strang,  “Maximum flows and minimum cuts in the plane,”  Advances in Applied Mathematics and Global Optimization: In Honor of Gilbert Strang, pp. 1–11, 2009.
  * [13] Özgün Çiçek, Ahmed Abdulkadir, Soeren S Lienkamp, Thomas Brox, and Olaf Ronneberger,  “3d u-net: learning dense volumetric segmentation from sparse annotation,”  in Medical Image Computing and Computer-Assisted Intervention–MICCAI 2016: 19th International Conference, Athens, Greece, October 17-21, 2016, Proceedings, Part II 19. Springer, 2016, pp. 424–432.
  * [14] Mehdi Noroozi and Paolo Favaro,  “Unsupervised learning of visual representations by solving jigsaw puzzles,”  in Computer Vision–ECCV 2016: 14th European Conference, Amsterdam, The Netherlands, October 11-14, 2016, Proceedings, Part VI. Springer, 2016, pp. 69–84.
  * [15] Ilya Loshchilov and Frank Hutter,  “Decoupled weight decay regularization,”  arXiv preprint arXiv:1711.05101, 2017.
  * [16] Adam Istvan Szucs et al.,  “Self-supervised segmentation of myocardial perfusion imaging spect left ventricles,”  in 2023 the 10th International Conference on Bioinformatics Research and Applications (ICBRA) (ICBRA 2023), September 22-24, 2023, Barcelona, Spain. ACM, 2023.
  * [17] David Bryant Mumford and Jayant Shah,  “Optimal approximations by piecewise smooth functions and associated variational problems,”  Communications on pure and applied mathematics, 1989.
  * [18] Jing Yuan, Egil Bae, and Xue-Cheng Tai,  “A study on continuous max-flow and min-cut approaches,”  in 2010 ieee computer society conference on computer vision and pattern recognition. IEEE, 2010, pp. 2217–2224.
  * [19] Antonin Chambolle,  “An algorithm for total variation minimization and applications,”  Journal of Mathematical imaging and vision, vol. 20, pp. 89–97, 2004.
  * [20] Daniel Cremers, Florian Tischhäuser, Joachim Weickert, and Christoph Schnörr,  “Diffusion snakes: Introducing statistical shape knowledge into the mumford-shah functional,”  International journal of computer vision, vol. 50, pp. 295–313, 2002.
  * [21] Daniel Cremers, Timo Kohlberger, and Christoph Schnörr,  “Shape statistics in kernel space for variational image segmentation,”  Pattern Recognition, vol. 36, no. 9, pp. 1929–1943, 2003.
  * [22] W. P. Segars, G. Sturgeon, S. Mendonca, Jason Grimes, and B. M. W. Tsui,  “4d xcat phantom for multimodality imaging research,”  Medical Physics, vol. 37, no. 9, pp. 4902–4915, 2010.
  * [23] H Malcolm Hudson and Richard S Larkin,  “Accelerated image reconstruction using ordered subsets of projection data,”  IEEE transactions on medical imaging, vol. 13, no. 4, pp. 601–609, 1994.
  * [24] Adam Paszke and et al.,  “Pytorch: An imperative style, high-performance deep learning library,”  pp. 8024–8035. Curran Associates, Inc., 2019.
  * [25] Mihály Héder, Ernő Rigó, Dorottya Medgyesi, Róbert Lovas, Szabolcs Tenczer, Attila Farkas, Márk Benjámin Emődi, József Kadlecsik, and Péter Kacsuk,  “The past, present and future of the elkh cloud,”  INFORMÁCIÓS TÁRSADALOM: TÁRSADALOMTUDOMÁNYI FOLYÓIRAT, vol. 22, no. 2, pp. 128–137, 2022.
  * [26] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay,  “Scikit-learn: Machine learning in Python,”  Journal of Machine Learning Research, vol. 12, pp. 2825–2830, 2011\.
  * [27] Jean Feydy, Thibault Séjourné, François-Xavier Vialard, Shun-ichi Amari, Alain Trouve, and Gabriel Peyré,  “Interpolating between optimal transport and mmd using sinkhorn divergences,”  in The 22nd International Conference on Artificial Intelligence and Statistics, 2019, pp. 2681–2690.


# Exploring medium mass nuclei using effective chiral nucleon-nucleon
interactions

Sota Yoshida<EMAIL_ADDRESS>Institute for Promotion of Higher
Academic Education, Utsunomiya University, Mine, Utsunomiya, 321-8505, Japan

###### Abstract

In the past decades, it has been shown that the three-body force is necessary
for a fully microscopic understanding of nuclear many-body systems, and thus
efficient schemes for storing and utilizing three-body matrix elements have
been developed. However, three-body matrix elements still demand a lot of
memory and/or storage and it readily exceeds the limitation of the state-of-
the-art supercomputers. We explore the possibility of mimicking the effects of
three-body forces in medium-mass region, especially lower $pf$ shell, by
density-dependent effective nucleon-nucleon interactions combining with a
modern family of chiral nucleon-nucleon interaction and an uncertainty
quantification on the coupling constants. We found that some effective
nucleon-nucleon potentials give decent description of ground state energies,
but also systematic underestimation in charge radii and deficiencies relevant
to the $pf$-shell region.

††preprint: …

## Introduction.

Chiral effective field theory (EFT) [1, 2, 3, 4] provides a systematic and
hierarchical description of nucleon-nucleon (NN) force and three-nucleon force
(3NF), and has become a standard approach for microscopic nuclear structure
and reaction studies. The potentials are parametrized in terms of low-energy
constants (LECs), and the LECs at NN-sector are determined to reproduce e.g.
scattering data.

While many applications of chiral potentials so far were usually based on the
NN force up to next-to-next-to-next-to-leading order (N3LO) plus the 3NF up to
next-to-next-to-leading order (NNLO), the NN part (and sometimes 3NF too) in
structure and reaction calculations is gradually being replaced by higher
order one such as next-to-next-to-next-to-next-to-leading order (N4LO) [5],
see e.g. [6, 7, 8] for recent applications.

In this work, we explore the capability of effective NN interactions,
consisting of N4LO NN force and density-dependent effective NN force from 3NF
at NNLO [9, 10, *KohnoErratum] without introducing explicit 3NFs at all.
Although it is a drastic attempt, this should make it easier, in terms of
computational burden, to explore the structure of various nuclei across the
nuclear chart using microscopic many-body methods such as the post-Hartree-
Fock methods and configuration interaction methods. The main focus of this
paper is on a specific intermediate-mass region, in particular a lower
$pf$-shell region, as the testing ground. We found that our effective NN
interactions are qualitatively well, and still need some improvements for more
quantitative discussions. We have also introduce a newly developed Julia
package covering chiral potentials and various many-body methods for nuclear
structure calculations. The package enables us more flexible studies such as
uncertainty quantification in medium-mass regions.

## Chiral potentials.

For the NN part, we use the so-called Entem-Machleidt-Nosyk (EMN) interaction
[5] up to N4LO with cutoffs $\Lambda=500$ MeV and $\Lambda_{\mathrm{SFR}}=700$
MeV, where $\Lambda_{\mathrm{SFR}}$ is a cutoff for the spectral function
regularization [12] for two-pion exchange (TPE) contributions. Detailed
expressions of the TPE contributions are described in Ref. [13]. Then, we
soften the NN interaction with the similarity renormalization group (SRG) in
the momentum space to get better convergence of binding energies with respect
to the model space size, which is typically specified by $e_{\mathrm{max}}$
truncation for harmonic oscillator quanta $e=2n+l\leq e_{\mathrm{max}}$.

Besides, we attempt to approximate both induced and genuine 3NFs by a density-
dependent NN force. We refer to a density-dependent NN force from three-body
forces as "2n3n" to distinguish it from NN and induced/genuine 3NF. Since the
detailed derivation of the 2n3n contribution at NNLO can be found in Refs. [9,
10, *KohnoErratum], we do not repeat the derivation here. Note that this
approximation and the matrix elements have been partly used in our previous
study [14] to investigate the impact of the choice of three-body forces on the
neutron driplines around $Z=9$–$12$.

As discussed in Ref. [5] and references therein, the TPE contribution in the
3NF at N3LO and N4LO have essentially the same mathematical structure at NNLO
and thereby can be approximated in terms of effective LECs at NNLO, i.e.,
$c_{1},c_{3},c_{4}$. Therefore, we regard the LECs for TPE contribution in the
3NF as independent parameters, and refer to them as
$\tilde{c}_{1},\tilde{c}_{3},\tilde{c}_{4}$ to distinguish from the ones that
appear in NN sector. The other LECs for 3NF appear at NNLO, i.e., the one-pion
exchange term $c_{D}$ and contact term $c_{E}$. The 2n3n terms are also
dependent on the Fermi momentum $k_{F}$, and its dependence on target nuclei
is nontrivial. For this reason, we fix it as $k_{F}=1.35$ $\mathrm{fm}^{-1}$
(corresponding to an empirical normal nuclear density $\rho=0.166$
$\mathrm{fm}^{-3}$), and restrict our considerations to the nuclei having
similar mass numbers and binding energy per nucleon, 40,48,52,54Ca and 56Ni,
and use them to fix the LECs.

Throughout this study, we work in the two-body space only, expecting both
induced 3NF via SRG and genuine 3NF can be mimicked by the 2n3n approximation.
If it works quantitatively well, this can be one solution to the problem of
computational costs for three-body matrix elements. If not, it indicates that
the structure of 3NF is much richer than that of a density dependent NN
interaction at NNLO and one needs to consider sub-leading contributions or
three-body matrix elements rather explicitly. One of the aims of this work is
to consider and clarify the capability of such effective NN potentials in
medium mass regions using a modern chiral NN potential up to N4LO.

To summarize, we work with both NN and 2n3n contributions up to N4LO (in an
effective viewpoints for the latter), and the free parameters in our potential
to be calibrated below are the five LECs for the 2n3n terms,
$\vec{c}=\\{\tilde{c}_{1},\tilde{c}_{3},\tilde{c}_{4},c_{D},c_{E}\\}$.

## Many-body methods.

In this work, we adopted a couple of many-body methods to calculate medium-
mass nuclei, and all the methods below are originally implemented in
NuclearToolkit.jl [15]. One is the so-called Hartree-Fock many-body
perturbation theory (HF-MBPT). This gives faster evaluations for g.s.
properties of medium-mass nuclei and thus will be used to calibrate LECs for
the 2n3n. With the HF-MBPT, one can evaluate perturbative corrections to the
estimation of observables at the HF level, known as Møller-Plesset
perturbation theory in the field of chemistry. More details can be found in
Refs. [16, 17, 18] on which our implementation is based.

Another one is the in-medium similarity renormalization group (IMSRG) [19].
The calculation is performed with the Magnus formulation of IMSRG [20] and the
arctangent generator [21]. We truncate the IMSRG flow equation up to the
normal-ordered two-body level at each step, which is the so-called IMSRG(2)
truncation. The comparisons with HF-MBPT and IMSRG(2) will be discussed later.

In addition to these, the valence-space IMSRG (VS-IMSRG) [22] is applied to
derive shell-model effective interactions for a target model space. The
generator to decouple between valence space and core/outer space is again the
arctangent generator. More precisely, we employ a more generalized one with
the denominator delta $\Delta$ introduced in Ref. [23] to decouple a multi-
shell valence space. However, $\Delta$ is now set as 0. Valence shell-model
calculations are performed using both NuclearToolkit.jl and KSHELL [24,
*KSHELL2]. The latter is more appropriate for larger systems.

## NuclearToolkit.jl.

We developed a new Julia package for nuclear structure calculations, named
NuclearToolkit.jl [15], covering chiral EFT potentials, HF-MBPT, IMSRG/VS-
IMSRG, and the valence shell model. The shell-model part of this package had
been developed as an independent code named ShellModel.jl, and used in the
previous study [26] to demonstrate the shell-model emulator with the
eigenvector continuation. The package is neither a wrapper nor a Julia
translation of existing codes such as EM potential [1], imsrg [27], KSHELL
[24, *KSHELL2], etc., while these pioneering works provide us with a certain
guideline to develop NuclearToolkit.jl. The package is designed as an open
source software for reproducibility and transparency of upcoming works and for
further acceleration of international research collaborations in the
community. It is automated by GitHub Actions to generate the document from
docstring on top of each function, deploy the document, and run the test code
to avoid destructive changes.

With the help of features of Julia language [28] (high readability,
performance, portability, package manager, etc.), one can use the package on a
laptop or remote environment without root permissions in the almost same way.
The author would be pleased if this package could promoted open science, and
served as a basis for exciting studies by readers.

## Calibration of LECs for the 2n3n contributions.

Next, we explain our procedures to calibrate the five LECs for the 2n3n terms.
The aim of this work is to explore the capability of a density-dependent NN
interaction to spare discrepancy between the NN-only results and the nature
without introducing explicit 3NFs. Parameter calibration and uncertainty
quantification in the following will not be exhaustive to examine, e.g.,
dependences on the priors and adopted data, but enough for the purpose. Our
strategy can be divided into the following three steps and more technical
details are described later: (i). Defining a prior distribution to encode the
so-called naturalness of LECs and a likelihood function for target nuclei.
(ii). Carrying out Markov chain Monte Carlo (MCMC) samplings over the five-
dimensional space with relatively smaller model spaces $e_{\mathrm{max}}=6,8$,
$\lambda_{\mathrm{SRG}}=2.0$ fm-1, and the oscillator parameter
$\hbar\omega=20$ MeV. (iii). Evaluating an optimal value for
$e_{\mathrm{max}}\geq 12$ based on the maximum a posteriori (MAP) value and
correlation seen in the results with smaller $e_{\mathrm{max}}=6,8$.

For step (i), it is expected that the g.s. energies of the five nuclei listed
above are not enough to constrain all the degrees of freedom among the five
LECs, while carrying out a sampling with many observables is intractable.
Therefore, we adopted independent Gaussian priors for the five LECs to softly
encode the naturalness and the deviations of
$\tilde{c}_{1},\tilde{c}_{3},\tilde{c}_{4}$ from $c_{1},c_{3},c_{4}$ are
imposed to be smaller. The log-prior can be written as

$\displaystyle\log
Pr(\vec{c})=-\sum^{5}_{i=1}\frac{(c_{i}-\mu_{i})^{2}}{2\sigma^{2}_{i}}.$ (1)

The mean values $\\{\mu_{i}\\}$ are $c_{1},c_{3},c_{4}$ of EMN potential [5]
and zero for $c_{D}$ and $c_{E}$. The widths of Gaussian $\\{\sigma_{i}\\}$
are 1.0 except for $\tilde{c}_{1}$. From analyses on $c_{1}$ at NN-sector [5,
13], we anticipated the $c_{1}$ having a smaller uncertainty, so the width is
set $\sigma=0.5$. We use g.s. energies of 40,48,52,54Ca and 56Ni as data to
constrain the LECs, and the log-likelihood $L(\vec{c})$ is given with g.s.
energies per nucleon

$\displaystyle
L(\vec{c})=-\frac{1}{2}\sum^{5}_{n=1}\left(\frac{E^{\mathrm{exp.}}_{n}-E^{\mathrm{th.}}_{n}(\vec{c})}{A_{n}}\right)^{2},$
(2)

where $E^{\mathrm{exp.}}$ is the experimental data, and $E^{\mathrm{th.}}$ is
the calculated result with given LECs, $\vec{c}$.

In step (ii), we adopted the so-called affine invariant MCMC method [29] as
the MCMC sampler and HF-MBPT as the many-body method to calculate g.s.
energies. The affine invariant MCMC method has several advantages: samples
tend to show much shorter autocorrelation time compared to ordinary methods
and one essentially does not need fine-tuning of the parameters in the sampler
such as proposal distribution for random-walk Metropolis-Hastings, leap-flog
parameters for Hamiltonian Monte Carlo, temperatures and their intervals for
replica exchange methods, etc. Since the sampler is an ensemble method, one
can easily benefit from multi-node parallelizations. We set the number of
walkers, i.e. the number of nodes to be used, as 64 and generated $\geq
20,000$ samples after certain burn-in steps to confirm convergence and
reproducibility by different runs.

Figure 1: MCMC samples over five LECs
($\tilde{c}_{1},\tilde{c}_{3},\tilde{c}_{4},c_{D}$, $c_{E}$). The lower
(upper) non-diagonal panels show correlation plots of samples drawn from
posteriors with $e_{\mathrm{max}}=6(8)$ truncation, and the diagonal panels
show histograms, where orange (hatched green) ones correspond to
$e_{\mathrm{max}}=6(8)$ results. The dotted lines show the prior distribution.
The cross (blue) and star (red) symbols are mean and MAP, respectively.

Regarding (iii), we plot the MCMC samples over five-dimensional LECs space in
FIG. 1. The panels in the lower and upper triangle, respectively, show density
plots of the MCMC samples for $e_{\mathrm{max}}=6$ and $8$ and covariance
between the two LECs. The cross and star symbols show the mean values and the
maximum a posteriori (MAP) values, respectively. The MAP values are now
evaluated by the parameters having the largest posterior ($\propto$ prior
$\times$ likelihood) among the MCMC samples. The diagonal panels show the
histograms of MCMC samples for $e_{\mathrm{max}}=6(8)$ results and the MAP
values. Under our settings (2n3n, prior, data, etc.), $\tilde{c}_{1}$ shows
almost no correlations with the others, whereas $c_{E}$ show rather modest
correlations with $\tilde{c}_{3},\tilde{c}_{4},c_{D}$. Note that, from our
analysis,
$\tilde{c}_{1}=-1.0,\tilde{c}_{3}=-0.53,\tilde{c}_{4}=-1.45,c_{D}=1.0,c_{E}=-0.25$
gives a similar size of repulsion in terms of the g.s. energies of the target
nuclei. To take account of model space convergence, which is typically seen at
$e_{\mathrm{max}}\sim 12$-$14$ for 40Ca with $\lambda_{\mathrm{SRG}}\leq 2.0$
fm-1, optimal values for the lower $pf$-shell region at $e_{\mathrm{max}}=8$
need to be modified so as to give a bit more repulsive results to the
experimental values. Indeed, the MAP values at $e_{\mathrm{max}}=8$ give more
repulsive results than those by the MAP at $e_{\mathrm{max}}=6$. Since the
$\tilde{c}_{1}$ is now uncorrelated to the others, we impose a role to gain
such repulsion on $\tilde{c}_{1}$. It should be noted that this procedure can
be valid only when adopting the current approximation and g.s. energies of
40,48,52,54Ca and 56Ni.

Table 1: The parameter sets to be used in this study. The MAP* is constructed from maximum a posteriori (MAP) at $e_{\mathrm{max}}=8$ of the MCMC samples, and $\tilde{c}_{1}$ is slightly shifted to the larger value to gain a repulsion. The other ones will be used in IMSRG/VS-IMSRG calculation to see some possible difference in the nuclear structure. LECs |  | $\tilde{c}_{1}$ | $\tilde{c}_{3}$ | $\tilde{c}_{4}$ | $c_{D}$ | $c_{E}$
---|---|---|---|---|---|---
MAP* |  | -1.100 | -4.753 | 4.089 | -0.277 | 0.820
A |  | -1.100 | -3.736 | 4.514 | -0.122 | 0.344
B |  | -1.600 | -6.205 | 3.326 | -0.901 | 1.520
C |  | -0.950 | -3.415 | 4.594 | -1.400 | -0.182
Figure 2: Results with MAP* and convergence pattern of g.s. energies and
charge radii with respect to model space size and oscillator parameter. The
open and filled symbols correspond to HF-MBPT and IMSRG(2) results,
respectively. The model space sizes, $e_{\mathrm{max}}=6,8,10,12$, are shown
by different symbols.

In this way, we constructed a tentative MAP value for larger
$e_{\mathrm{max}}$, which will be called MAP* in the followings. In addition
to the MAP*, we will also consider three-different parameters generated from
the posterior distributions whose $\tilde{c}_{1}$ value is slightly changed
for the reason above. The aim to consider these parameters is to see the
effect on the nuclear structure by some degrees of freedom, which may not be
well constrained by only the g.s. energies. We list them in Table 1.

## Ground state properties with HF-MBPT/IMSRG(2).

We discuss ground state properties, energy and charge radius, of sub-shell
closed nuclei, 40,48,52,54Ca and 56Ni. In FIG. 2, the convergence pattern of
g.s. energies and charge radii with respect to the $e_{\mathrm{max}}$
truncation and the oscillator parameter $\hbar\omega$ are shown. The open and
filled symbols correspond to HF-MBPT and IMSRG(2) results, respectively. Note
that the HF-MBPT is the 3rd (2nd) order for g.s. energies (charge radii), and
more concrete expressions for expectation values for scaler operators can be
found in, e.g., [18]. As a whole, the results are nearly converged at
$e_{\mathrm{max}}=12$ and show rather flat behavior around $\hbar\omega=20$
MeV, so we will use interactions with $e_{\mathrm{max}}=12,\hbar\omega=20$ MeV
in the followings.

Table 2: The IMSRG(2) results of g.s. energies using the four parameters with
$\hbar\omega=20$ MeV, $e_{\mathrm{max}}=12$.

$E_{\mathrm{g.s.}}$ (MeV) |  |  40Ca |  48Ca |  52Ca |  54Ca |  56Ni
---|---|---|---|---|---|---
MAP* |  | -351.855 | -414.203 | -440.373 | -449.032 | -476.644
A |  | -349.161 | -414.154 | -443.224 | -453.427 | -476.886
B |  | -351.703 | -409.548 | -430.634 | -436.610 | -471.810
C |  | -344.348 | -409.784 | -439.942 | -450.746 | -472.148
Exp. |  | -342.052 | -416.001 | -438.328 | -445.365 | -483.996

Table 3: The IMSRG(2) results of charge radii using the four parameters with
$\hbar\omega=20$ MeV, $e_{\mathrm{max}}=12$.

$R_{\mathrm{ch}}$ (fm) |  |  40Ca |  48Ca |  52Ca |  54Ca |  56Ni
---|---|---|---|---|---|---
MAP* |  | 2.8796 | 2.8582 | 2.8689 | 2.8823 | 3.0851
A |  | 2.8631 | 2.8350 | 2.8421 | 2.8539 | 3.0587
B |  | 2.9068 | 2.8934 | 2.9100 | 2.9263 | 3.1209
C |  | 2.8608 | 2.8300 | 2.8353 | 2.8464 | 3.0516
Exp. |  | 3.4776 | 3.4771 | - | - | -

In Tables 3 and 3, we summarized IMSRG(2) results for ground state properties
of the target nuclei. The g.s. energies by our interactions are well, but
systematic underestimation is seen in charge radii. This is consistent with
the trade-off between binding energies and charge radii seen in previous
studies, e.g., [30], but underestimation by our 2n3n interactions are more
severe and almost irrelevant to the choice of LECs.

Figure 3: Charge radii of the target nuclei. The NN-only results with
$e_{\mathrm{max}}$ and different SRG resolution scales, $\lambda=2.0$
$\mathrm{fm}^{-1}$ and $\lambda=\infty$ (unevolved) are show by circles and
triangles, respectively. The closed violin plots show uncertainties of 2n3n
contributions, which are evaluated by the MCMC samples ($e_{\mathrm{max}}=6$)
shown FIG. 1, and the open violin plots show a naive estimates on charge radii
by NN (bare) + 2n3n. See the main text for more details.

Now let us look at possible origins of the underestimation from an uncertainty
analysis. We perform HF-MBPT(2) evaluation of charge radii for 10,000 MCMC
samples with $e_{\mathrm{max}}=6$ and $\hbar\omega=20$ MeV, and the results
are drawn as the violin plot (filled one) in FIG. 3. As shown in FIG. 2,
charge radii for 40,48Ca are almost converged at $e_{\mathrm{max}}=6$ and
close to IMSRG(2) results, so only the HF-MBPT results with a smaller model
space is enough to estimate the possible 2n3n contributions on charge radii.
In addition to this, we show the NN-only results with SRG evolved interaction
($\lambda_{\mathrm{SRG}}=2.0$ fm-1,$e_{\mathrm{max}}=12$, circle symbols) and
unevolved one ($e_{\mathrm{max}}=12$, triangle symbols). From these, we found
that there is no room for the 2n3n contribution, giving decent binding
energies, to explain the observed charge radii.

The open violin plots show naive estimates of results under the bare NN plus
2n3n; the open violins are drawn merely by shifting the closed violins by
offsets,
$R_{\mathrm{ch}}(\mathrm{e12},\mathrm{bare})-R_{\mathrm{ch}}(\mathrm{e12},\lambda=2.0\
\mathrm{fm}^{-1})$. This offset roughly shows the size of induced many-body
forces, which are now ignored for brevity. If one employs the harder NN
interaction, the size of 2n3n needed to explain binding energies must be
smaller, and thereby the 2n3n contributions on charge radii are to be smaller
too. In this sense, the shifted values (open violin plots) in FIG. 2 may give
us upper values by effective NN potentials.

We also note that, in the current approximation, all the contributions are
truncated up to the two-body level by definition and there is no explicit
normal ordered one-body (NO1B) contribution from both induced and genuine 3NF.
This can deteriorate the consistent description of binding energies and charge
radii.

Now, we explore the extensibility of the 2n3n approximation over the nuclear
chart. As already mentioned, density-dependent interactions have a dependence
on Fermi momentum $k_{F}$ and we fixed it as the momentum from an empirical
normal nuclear density. Performing HF-MBPT calculation for 4He, 16O, 90Zr, and
132Sn, we naively estimate the validity of 2n3n interaction with a common
$k_{F}$ in different mass regions. For lighter nuclei, 4He and 16O, the
binding energies are underestimated by about $10\%$, whereas the binding
energies are overestimated for heavier ones, 90Zr and 132Sn, by about $10\%$.
Regarding the slower convergence for heavier nuclei, fully converged results
with respect to the model space size will be worse for 90Zr and 132Sn. We also
note that the charge radii are systematically underestimated as in the Ca and
Ni regions, probably due to the same reasons as above. From this observation,
it is not straightforward to use the 2n3n approximation with a common $k_{F}$
in different mass ranges. One possible improvement is to introduce a
phenomenological parameter to scale the 2n3n effect as a function of the mass
number and to introduce a correction of the isospin asymmetry [9], but this is
left as a future study in a different context.

## Valence space results of lower $pf$-shell nuclei.

Figure 4:  Systematics of first 2+ excitation energies of even calcium
isotopes calculated by the VS-ISMRG. The star and square symbols show the
results with TNO or ENO, and the others are using effective interactions
derived with the 40Ca reference.

Next, we discuss the shell-model results for some lower $pf$-shell nuclei
using the effective interactions derived by the VS-IMSRG from our effective NN
potentials. We focus on the first 2+ energies for Ca isotopes, and low-lying
spectra for 56Ni.

In FIG. 4, systematics of the first 2+ excited states of even calcium isotopes
are shown for various different effective interactions. The model space for
the effective interactions is the $pf$ shell on top of the 40Ca core, and the
reference state is 40Ca for our interactions, MAP*, A, B, and C. As partly
shown by open star symbols, the dependence on the reference state, which is
examined using target normal ordering (TNO)/ensemble normal ordering (ENO)
[22], is minor for calcium isotopes. One can see that most of the interactions
give too low excitation energies especially for non-magic Ca isotopes. The
$N=28$ ($N=34$) gap is underestimated (overestimated) as a whole; only the
exception is the "Stroberg2017" result at $N=28$. The interaction
"Stroberg2017" is derived by the VS-IMSRG approach with consistently SRG
evolved EM NN interaction at N3LO plus 3NF at NNLO, and the interaction files
are available on the author’s GitHub repository [31].

From the trend seen in FIG. 4, the systematic underestimation of the first 2+
energies for non-magic Ca isotopes indicates that pairing components in
$0f7/2$ and $0f5/2$ in the derived interactions are weak in general. Let us
look at the results with an unevolved N4LO NN interaction, the pentagon
symbols in the figure. While the differences in the 3NF part are reflected as
differences in 2+ energies, the results are essentially the same as those by
the unevolved NN potential. In this sense, one of the origins for weak pairing
components is the NN interaction itself, although it is nontrivial to identify
the origins of such trends and results are consequences of non-perturbative
IMSRG/VS-IMSRG flow under given potentials.

Figure 5: Low-lying states of 56Ni by the VS-IMSRG approach. The effective
interactions are the same as in FIG. 4.

Next, we show the shell-model results of 56Ni in FIG. 5. All the effective
interactions are derived for the $pf$ shell on top of the 40Ca core using TNO
and diagonalized within the full $pf$-shell space. As can be expected from the
results of 48Ca, the $N(Z)=28$ gap is underestimated by our effective
interactions, leading to severe underestimation of the first 2+ energy.
Indeed, all the VS-IMSRG results with our interaction giving Ex.$(2+)\sim 500$
keV show that the occupation numbers of proton/neutron $0f7/2$ orbit are about
$4$-$5$ for the $0^{+}_{1}$ and $2^{+}_{1}$ states. The VS-IMSRG result by
Ref. [31] gives a better description of the spectrum. The effective
interaction gives a rather prominent effective $N(Z)=28$ gap as partly shown
48Ca results above. This is partly due to the large gap between $0f7/2$ and
$0p3/2$, $\sim 5-6$ MeV, in the single particle energies.

It should be noted that one can find the VS-IMSRG results with the so-called
"EM 1.8/2.0" interaction [32], giving the 2+ systematics of Ca and Ni isotopes
well described [33] and slightly overestimated $N(Z)=28$ gap. The possible
deficiency of pairing components and $N(Z)=28$ gap in NN potential discussed
above would be compensated by the genuine 3NF, which may lead to various
successful descriptions by the EM1.8/2.0 interaction. The systematic
overestimation of 2+ energies of Ni isotopes by the VS-IMSRG with the
EM1.8/2.0 is attributed to higher order contributions of IMSRG, and then the
IMSRG(3) truncation has begun to be considered with lighter nuclei [34].

The $N(Z)=28$ gap in the VS-IMSRG approach is certainly related to the many-
body machinery, but it is also likely to be related to the absence of an
explicit NO1B term from 3NF in effective NN force. It is desirable to consider
this point from more diverse perspectives.

## Summary and outlook.

In this work, we explored the capability of chiral effective nucleon-nucleon
(NN) potentials without introducing explicit three-nucleon force (3NF). We
employed a modern chiral NN interaction up to N4LO [13, 5] and a density-
dependent interaction [9, 10, *KohnoErratum], and calibrated the low-energy
constants (LECs) for the density-dependent NN interaction by MCMC samplings.
We found the binding energies of some lower $pf$-shell nuclei are well
explained by our effective NN potentials, whereas systematic underestimates
were seen in charge radii and that is almost irrelevant to the choice of the
LECs. From an uncertainty analysis on the LECs and the results with a bare NN
interaction, we discussed some possible origins of such underestimations. From
the VS-IMSRG results, significant underestimation is seen in the pairing of
the $0f$ orbits and the $N(Z)=28$ gap. The deficiency in the current effective
interactions may be caused by a combination of multiple factors, including the
NN potential itself, density-dependent approximation of 3NF, free-space SRG
and the treatment of induced many-body forces, many-body methods, etc. For
example, if we extend the model space to $pf$-$sdg$ shell in the VS-IMSRG
decoupling phase, the 2+ energies of Ca isotopes and 56Ni improve by a few
hundred keV. The same trend can be seen in previous work using MBPT [35] and
when one employs the extended Kuo-Krenciglowa method [36] to derive
interactions in the $pf$-$sdg$ space and diagonalizes it within the $pf$
shell. This indicates that the inclusion of the $0g9/2$ orbit or others in the
valence space can be important to derive effective interactions in this
region. The $pf$-shell (and beyond) nuclei would give us good testing grounds
for these compound problems.

It should be also noted that sub-leading density-dependent NN forces in higher
order have been already derived in Refs. [37, 38, 39], and applications of
these are left as future works along this study.

It is expected that this work and the package, NuclearToolkit.jl, will enable
more detailed and multi-angle analyses on chiral potentials, many-body
methods, and structure of medium-mass nuclei.

## Acknowledgments.

The author thanks Noritaka Shimizu, Takayuki Miyagi, and Tokuro Fukui for
discussions at the early stage of the development of NuclearToolkit.jl, and
Michio Kohno for discussions on the density-dependent NN force. This work was
supported by JSPS KAKENHI (Grants No. 22K14030) and partially by the
computational resource of Fujitsu PRIMERGY CX400M1/CX2550M5 (Oakbridge-CX) at
the Information Technology Center, The University of Tokyo.

## References

  * Machleidt and Entem [2011] R. Machleidt and D. Entem, Phys. Rept. 503, 1 (2011).
  * Machleidt and Sammarruca [2016] R. Machleidt and F. Sammarruca, Physica Scripta 91, 083007 (2016).
  * Epelbaum  _et al._ [2009] E. Epelbaum, H.-W. Hammer, and U.-G. Meißner, Rev. Mod. Phys. 81, 1773 (2009).
  * Epelbaum  _et al._ [2020] E. Epelbaum, H. Krebs, and P. Reinert, Frontiers in Physics 8 (2020), 10.3389/fphy.2020.00098.
  * Entem  _et al._ [2017] D. R. Entem, R. Machleidt, and Y. Nosyk, Phys. Rev. C 96, 024004 (2017).
  * Hüther  _et al._ [2020] T. Hüther, K. Vobig, K. Hebeler, R. Machleidt, and R. Roth, Physics Letters B 808, 135651 (2020).
  * Kravvaris  _et al._ [2020a] K. Kravvaris, K. R. Quinlan, S. Quaglioni, K. A. Wendt, and P. Navrátil, Phys. Rev. C 102, 024616 (2020a).
  * Kravvaris  _et al._ [2020b] K. Kravvaris, K. R. Quinlan, S. Quaglioni, K. A. Wendt, and P. Navrátil, Phys. Rev. C 102, 024616 (2020b).
  * Holt  _et al._ [2010] J. W. Holt, N. Kaiser, and W. Weise, Phys. Rev. C 81, 024002 (2010).
  * Kohno [2013] M. Kohno, Phys. Rev. C 88, 064005 (2013).
  * Kohno [2017] M. Kohno, Phys. Rev. C 96, 059903 (2017).
  * Epelbaum  _et al._ [2004] E. Epelbaum, W. Glöckle, and U.-G. Meißner, The European Physical Journal A - Hadrons and Nuclei 19, 125 (2004).
  * Entem  _et al._ [2015] D. R. Entem, N. Kaiser, R. Machleidt, and Y. Nosyk, Phys. Rev. C 91, 014002 (2015).
  * Tsunoda  _et al._ [2020] N. Tsunoda, T. Otsuka, K. Takayanagi, N. Shimizu, T. Suzuki, Y. Utsuno, S. Yoshida, and H. Ueno, Nature 587, 66 (2020).
  * [15] S. Yoshida, https://github.com/SotaYoshida/NuclearToolkit.jl.
  * Shavitt and Bartlett [2009] I. Shavitt and R. J. Bartlett,  _Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory_, Cambridge Molecular Science (Cambridge University Press, 2009).
  * Tichai  _et al._ [2020] A. Tichai, R. Roth, and T. Duguet, Frontiers in Physics 8 (2020), 10.3389/fphy.2020.00164, and Doctoral thesis by A. Tichai @TU Darmstadt, 2017.
  * Miyagi  _et al._ [2022] T. Miyagi, S. R. Stroberg, P. Navrátil, K. Hebeler, and J. D. Holt, Phys. Rev. C 105, 014302 (2022).
  * Tsukiyama  _et al._ [2011] K. Tsukiyama, S. K. Bogner, and A. Schwenk, Phys. Rev. Lett. 106, 222502 (2011).
  * Morris  _et al._ [2015] T. D. Morris, N. M. Parzuchowski, and S. K. Bogner, Phys. Rev. C 92, 034331 (2015).
  * Hergert  _et al._ [2016] H. Hergert, S. Bogner, T. Morris, A. Schwenk, and K. Tsukiyama, Physics Reports 621, 165 (2016), memorial Volume in Honor of Gerald E. Brown.
  * Stroberg  _et al._ [2019] S. R. Stroberg, H. Hergert, S. K. Bogner, and J. D. Holt, Annual Review of Nuclear and Particle Science 69, 307 (2019).
  * Miyagi  _et al._ [2020] T. Miyagi, S. R. Stroberg, J. D. Holt, and N. Shimizu, Phys. Rev. C 102, 034320 (2020).
  * Shimizu [2013] N. Shimizu, “Nuclear shell-model code for massive parallel computation, "kshell",” (2013), arXiv:1310.5431 .
  * Shimizu  _et al._ [2019] N. Shimizu, T. Mizusaki, Y. Utsuno, and Y. Tsunoda, Computer Physics Communications 244, 372 (2019).
  * Yoshida and Shimizu [2022] S. Yoshida and N. Shimizu, Progress of Theoretical and Experimental Physics 2022 (2022), 10.1093/ptep/ptac057, 053D02.
  * [27] S. R. Stroberg, https://github.com/ragnarstroberg/imsrg.
  * Bezanson  _et al._ [2017] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, SIAM Review 59, 65 (2017).
  * Goodman and Weare [2010] J. Goodman and J. Weare, Communications in Applied Mathematics and Computational Science 5, 65 (2010).
  * Hoppe  _et al._ [2019] J. Hoppe, C. Drischler, K. Hebeler, A. Schwenk, and J. Simonis, Phys. Rev. C 100, 024318 (2019).
  * Stroberg  _et al._ [2017] S. R. Stroberg, A. Calci, H. Hergert, J. D. Holt, S. K. Bogner, R. Roth, and A. Schwenk, Phys. Rev. Lett. 118, 032502 (2017), interaction files are available through https://github.com/ragnarstroberg/PRL118032502.
  * Hebeler  _et al._ [2011] K. Hebeler, S. K. Bogner, R. J. Furnstahl, A. Nogga, and A. Schwenk, Phys. Rev. C 83, 031301 (2011).
  * Simonis  _et al._ [2017] J. Simonis, S. R. Stroberg, K. Hebeler, J. D. Holt, and A. Schwenk, Phys. Rev. C 96, 014303 (2017).
  * Heinz  _et al._ [2021] M. Heinz, A. Tichai, J. Hoppe, K. Hebeler, and A. Schwenk, Phys. Rev. C 103, 044318 (2021).
  * Holt  _et al._ [2014] J. D. Holt, J. Menéndez, J. Simonis, and A. Schwenk, Phys. Rev. C 90, 024312 (2014).
  * Tsunoda  _et al._ [2017] N. Tsunoda, T. Otsuka, N. Shimizu, M. Hjorth-Jensen, K. Takayanagi, and T. Suzuki, Phys. Rev. C 95, 021304 (2017), note that the interaction for the $pf$-$sdg$ shell is unpublished.
  * Kaiser and Niessner [2018] N. Kaiser and V. Niessner, Phys. Rev. C 98, 054002 (2018).
  * Kaiser and Singh [2019] N. Kaiser and B. Singh, Phys. Rev. C 100, 014002 (2019).
  * Kaiser [2020] N. Kaiser, Phys. Rev. C 101, 014001 (2020).