=    min_i = 0, ⋯k-1 1/k ∑_j = i+1^k q^j-i/k =  min_i = 0, ⋯k-1 q^1/k/k + 1/k ∑_j = i+2^k q^j-i/k ≥q^1/k/k
and we have:
1-q/k q ∑_i = 1^k q^i/k/1+(1/q-1)q^i/k    =    1-q/k ∑_i = 1^k q^i/k/q+(1-q)q^i/k

≥    1-q/k ∑_i = 1^k q^i/k/q^i/k+(1-q)q^i/k

≥    1-q/k ∑_i = 1^k 1/1+(1-q) = 1-q/kk/2-q = 1-q/2-q
Let $k = -\log(q)$, we have

sup_Ψ∈ inf_F ∈q R(Ψ, F) ≥inf_F ∈q R(Ψ_q, k, F) ≥min(-e^-1/log(q), 1-q/2-q)
For $q \in ]0, 0.4[$, it's easy to show that $\frac{-e^{-1}}{\log(q)} \leq 0.5 \leq  \frac{1-q}{2-q}$, therefore:
sup_Ψ∈ inf_F ∈q R(Ψ, F) ≥-e^-1/log(q) = c_1/log(1/q)  c_1 = e^-1

Let $q \in ]0, 1[$ and $K \in \mathbb{N}^*$. Define $\varepsilon = q^{\frac{1}{K}} \in ]q, 1[$ and $a_k = w \varepsilon^k \in [q, 1[$ for $k = 1 \cdots K$. Define the following familie of distributions $F_k \in \Fb_0 (w,\{q\})$:
_k(p)    =
1      if  p ≤a_k

1/1+(1/q-1)p-a_k/w-a_k      if  p ≤w

0      if  p >w.

For all $1 \le k \le K$, the revenue curve associated with $\bF_{k}(p)$ is monotone non-increasing on $[a_k,w]$. Furthermore, the revenue curve is convex in $[a_k,w]$ as:
(p _k(p))^'    =
1      if  p ≤a_k

w(1-a_k/wq)/(1+(1/q-1)p-a_k/w-a_k)^2      if  p ≤w

0      if  p >w.

G(p)    =    1/K ∑_i=1^K  p F̅_i(p)/(F_i).
Note that $G(p)$ is convex on each interval $[a_{k},a_{k+1}]$ for for $k=1,...,K-1$. Hence, the sup  on an interval must be attained at one of the extreme points of the interval. Therefore, using Yao's principle <cit.>, we have:
sup_Ψ∈ inf_F ∈_0 (w,{q}) R(Ψ, F) ≤sup_p ≥0 G(p)    =    1/K  max_1 ≤k ≤K  ∑_i=1^K  p F̅_i(p)/(F_i)

=    1/K  max_1 ≤k ≤K  { ∑_i=1^k  a_k /a_i + ∑_i=k+1^K  a_k /a_i  1/1+(1/q-1)a_k-a_i/w-a_i }

=    1/K  max_1 ≤k ≤K  { ∑_i=1^k ε^k-i + ∑_i=k+1^K ε^k-i 1-ε^i/1- ε^i+1/q(ε^k-ε^i)-ε^k+ε^i }

=    1/K  max_1 ≤k ≤K  { ∑_i=0^k-1 ε^i + ∑_i=k+1^K ε^k-i 1-ε^i/1-ε^k+1/q(ε^k-ε^i) }

=    1/K  max_1 ≤k ≤K  { ∑_i=0^k-1 ε^i + ∑_i=k+1^K  1-ε^i/ε^i-k-ε^i+1/q(ε^i-ε^2i-k) }

=    1/K  max_1 ≤k ≤K  { ∑_i=0^k-1 ε^i + ∑_i=k+1^K  1-ε^i/ε^i(ε^-k-1)+1/qε^i(1-ε^i-k) }

≤    1/K  max_1 ≤k ≤K  { 1-ε^k/1-ε + ∑_i=k+1^K  ε^-i-1/(ε^-k-1)+1/q(1-ε) }

≤    1/K  max_1 ≤k ≤K  { 1-ε^k/1-ε + ∑_i=k+1^K  ε^-i/(ε^-k-1)+1/q(1-ε) }

≤    1/K  max_1 ≤k ≤K  { 1-ε^k/1-ε +  1/1/ε-11/ε^K+1-1/ε^k+1/(ε^-k-1)+1/q(1-ε) }

≤    1/K (1-ε)  max_1 ≤k ≤K  { 1-ε^k +  1/ε^K-1/ε^k/(ε^-k-1)+1/q(1-ε) }

≤    1/K (1-ε)  max_1 ≤k ≤K  { 1-ε^k +  1/q/(ε^-k-1)+1/q(1-ε) }

≤    1/K (1-ε)  max_1 ≤k ≤K  { 1-ε^k +  1/q(ε^-k-1)+(1-ε) }

≤    1/K (1-ε)  ( 1-q +  1/1-ε) = 1/K (1-q^1/K)  ( 1-q +  1/1-q^1/K)

By choosing $K = \log(\frac{1}{q})$, we get:
sup_Ψ∈ inf_F ∈_0 (w,{q}) R(Ψ, F) ≤1/log(1/q) (1-e^-1)  ( 1-q +  1/1-e^-1) ≤c_2/log1/q  c_2 =  1/1-e^-1 ( 1 +  1/1-e^-1)
Therefore $\sup_{\Psi \in \MechSet} \sup_{\Psi \in \MechSet} \inf _{F \in \Fb_0 (w,\{q\})} R(\Psi, F) \to 0$ as $q \to 0$.

Proof of <Ref>.

This proof is divided into two steps. In the first step, we will show the lower bound by analyzing the performance of a specific mechanism. Then in a second step, we will derive the upper through the analysis of a family of hard cases when $q$ is close to $1$.

Through the proof we will assume that $q \ge 3/4$ since we are interested in the limit when  when $q$ is close to $1$.

Step 1: Lower bound
Let  us define the following measure parameterized by $a, b \geq 0$:
d Ψ(u) =
a       u = 1

b (u _0(u, 0, 1, 1,q ))'/u _0(u, 0, 1, 1,q )       u > 1.

Note that $d \Psi(u) \ge 0$  since  the revenue function $u \bH_{0}(u, 0, 1, 1,q )$ is increasing in $[1, \infty[.$ Let us determine the condition on  the parameters $a$ and $b$ so that $\Psi(u) $ is a distribution. For that we the following
∫_0^∞ dΨ(u) = 1,
this implies that
a + b log(1 _0(1, 0, 1, 1,q )/lim_u →∞ u _0(u, 0, 1, 1,q )) = 1,
since $\bH_{0}(1, 0, 1, 1,q ) = 1 $ and $\lim_{u \to \infty} u \bH_{0}(u, 0, 1, 1,q ) = 1-q$, then we get that
a + b log(1/1-q) = 1.
Hence the relation between $a$ and $b$ is as follows
b  = 1-a/log(1/1-q)    a ∈[0,1].

Using <ref>, we have

inf_F ∈q R(Ψ, F)    ≥    inf_F ∈_0(w,{q}) R(Ψ, F)
=    min{ inf_x ∈[q, 1]  1/x [ ∫_0^x  u dΨ̂(u) + ∫_x^1   u _0(u, x,1, 1,q ) dΨ̂(u) ],

inf_x ∈[1, +∞[ 1/x  _0(x , 0,1, 1,q ) ∫_0^x  u _0(u, 0,1, 1,q )  dΨ̂_q, k(u) }

We will analyze each term separately depending if $x \le 1$ or $x \ge 1.$

Case 1: $x \le 1$
We have
1/x [ ∫_0^x  u dΨ(u) + ∫_x^1   u _0(u, x,1, 1,q ) dΨ(u) ] (a)= aq/x ≥aq,
where the last inequality is due to the fact that $x \le 1$ and $(a)$ is due to $d\Psi(u) = 0$ for $u<1$, $d\Psi(1) = a$ and $1 \bH_{0}(1, x,1, 1,q ) = q$.

Hence we conclude that

inf_x ∈[q, 1]  1/x [ ∫_0^x  u dΨ̂(u) + ∫_x^1   u _0(u, x,1, 1,q ) dΨ̂(u) ] ≥aq.

Case 2: $x \ge 1$
Let us now analyze the second term, we have
∫_0^x  u _0(u, 0,1, 1,q )  dΨ̂_q, k(u)     =     ∫_0^1  u _0(u, 0,1, 1,q )  dΨ̂_q, k(u) +   ∫_1^x  u _0(u, 0,1, 1,q )  dΨ̂_q, k(u)

(a)=    aq + b ∫_1^x  u _0(u, 0,1, 1,q ) (u _0(u, 0, 1, 1,q ))'/u _0(u, 0, 1, 1,q ) du

=    aq + b  [u _0(u, 0, 1, 1,q )]_1^x

≥    aq + b (x _0(x, 0, 1, 1,q ) - q ),
$(a)$ is due to $d\Psi(u) = 0$ for $u<1$, $d\Psi(1) = a$ and $1 \bH_{0}(1, x,1, 1,q ) = q$.  Hence we conclude that
1/x _0(x, 0, 1, 1,q ) ∫_0^x  u _0(u, 0,1, 1,q )  dΨ̂_q, k(u)
≥     qa-b/x _0(x, 0, 1, 1,q )+ b

≥    (a-b)(1-q) + b,
where the last inequality we used the fact that $x \to x  \bH_{0}(x , 0,1, w,q )$ is monotone in $[1, \infty[$ and that $\lim_{\infty} \bH_{0}(x , 0,1, w,q ) = (1-q)/q.$

Thus we conclude that

inf_x ∈[1, +∞[ 1/x _0(x, 0, 1, 1,q ) ∫_0^x  u _0(u, 0,1, 1,q )  dΨ̂_q, k(u)
≥    (a-b)(1-q) + b,

By combining (<ref>), (<ref>) and (<ref>) we get that
inf_F ∈q R(Ψ, F) ≥minaq, (a-b)(1-q) + b.

Now let us set a = q/2(q-1/2)log(1/1-q) +q.

Note that $a\in [0,1]$ as $q \in [3/4, 1],$ this also leads to the fact that
(a-b)(1-q) + b = aq = q^2/2(q-1/2)log(1/1-q) +q(a)
≥ 9/4 1/2(3/4-1/2)log(1/1-q) +1    =    9/2 log(1/1-q) + 4
(b)≥    9/6 1/log(1/1-q) ,
where the (a) inequality stems from the fact $q \in [3/4, 1],$ and (b) we have $\log\left(\frac{1}{1-q}\right) \ge 1.$

Hence we we get that
inf_F ∈q R(Ψ, F) ≥9/6 1/log(1/1-q).
This conclude the lower bound.

Step 2: Upper bound
To show the upper bound we will introduce a family of worst cases.  We consider the family of distributions $(F_{r, 1, q})_{r \geq 1} \in \Fb_0 (w,\{q\})$ and the following weight distribution:

dλ(r)    =
0      if  r ≤1

1/log(1/1-q)(r _0(r | (0,1), (1,q) ))'/(r _0(r | (0,1), (1,q) ))      if  r ≥1

One can verify that  $\int_0^{\infty} d\lambda(r) =1$ and that  $ d\lambda(r) \ge 0.$ Let us define
G(p)    =    ∫_0^∞  p F_r, 1, q(p)/(F_r, 1, q) dλ(r)

Therefore, using Yao's principle <cit.>, we have:

sup_Ψ∈ inf_F ∈_0 (w,{q}) R(Ψ, F)    ≤    sup_p ≥0 G(p)
= sup_p ≥0 ∫_0^∞  p F_r, 1, q(p)/(F_r, 1, q) dλ(r).

Now let us analyze the intergral term. We have
∫_0^∞  p F_r, 1, q(p)/(F_r, 1, q) dλ(r)
=    ∫_p^∞  p F_r, 1, q(p)/(r _0(r | (0,1), (1,q) )) 1/log(1/1-q)(r _0(r | (0,1), (1,q) ))'/(r _0(r | (0,1), (1,q) )) dr

=    1/log(1/1-q)   p _0(p | (0,1), (1,q) ) ∫_p^∞ (r _0(r | (0,1), (1,q) ))'/(r _0(r | (0,1), (1,q) ))^2 dr

=    1/log(1/1-q)   p _0(p | (0,1), (1,q)) (1/p _0(p | (0,1), (1,q) )  - lim_p →∞ 1/p _0(p | (0,1), (1,q) ) )

=    1/log(1/1-q) ( 1- (1/q-1 ) p _0(p | (0,1), (1,q) ) )

≤    1/log(1/1-q) ( 1- (1/q-1 ) 1 _0(1 | (0,1), (1,q) ) ),
where the last inequality is due to the fact that the revenue curve $p \to p \bH_{0}(p | (0,1), (1,q) )$ non-decreasing in $[1, \infty[$

Hence we conclude that for any $p \ge 0$, we have
∫_0^∞  p F_r, 1, q(p)/(F_r, 1, q) dλ(r)
≤    1/log(1/1-q) ( 1- (1/q-1 ) 1 _0(1 | (0,1), (1,q) ) ) ≤1/log(1/1-q),

By using, the last inequality with (<ref>), we conclude the result.

OLD PROOF

Let $q \geq \frac{1}{2}$. Let $a, b \geq 0$ and define the following mechanism:
d Ψ(u) =
a       u = 1

b (u _0(u, 0, 1, 1,q ))'/u _0(u, 0, 1, 1,q )       u > 1

is $d \Psi(u)$ always positive?
$\Psi$ is a mechanism, therefore
∫d Ψ(u) = 1        a + b log(1 _0(1, 0, 1, 1,q )/lim_u →∞ u _0(u, 0, 1, 1,q )) = 1

a + b log(1/1-q) = 1

b  = 1-a/log(1/1-q)    a ∈[0,1]
Using <ref>, we have:
inf_F ∈q R(Ψ, F) ≥inf_F ∈_0(w,{q}) R(Ψ, F)     =    min{ inf_x ∈[q, 1]  1/x [ ∫_0^x  u dΨ̂(u) + ∫_x^1   u _0(u, x,1, 1,q ) dΨ̂(u) ],

inf_x ∈[1, +∞[ 1/x  _0(x , 0,1, w,q ) ∫_0^x  u _0(u, 0,1, 1,q )  dΨ̂_q, k(u) }

=    min{ inf_x ∈[q, 1]  aq/x, inf_x ∈[1, +∞[ 1/x  _0(x , 0,1, w,q ) aq
+b ∫_1^x  u _0(u, 0,1, 1,q ) (u _0(u, 0, 1, 1,q ))'/u _0(u, 0, 1, 1,q ) dΨ̂_q, k(u) }

=    min{aq, inf_x ∈[1, +∞[ q a-b/x  _0(x , 0,1, w,q )  + b }

=    min{aq, q a-b/1  _0(1 , 0,1, w,q )  + b,lim_x →∞ q a-b/x  _0(x , 0,1, w,q )  + b}

=    min{aq, a, (a-b)(1-q) } = min{aq, (a-b)(1-q) }

=    min{aq, a(1-q(1+1/log(1/1-q))) + q/log(1/1-q) }
We choose $a \in [0,1]$, such that $aq = a(1-q(1+\frac{1}{\log\left(\frac{1}{1-q}\right)})) + \frac{q}{\log\left(\frac{1}{1-q}\right)}$, i.e:
a = q/2(q-1/2)log(1/1-q) +q ∈[0,1]
sup_Ψ∈ inf_F ∈q R(Ψ, F) ≥inf_F ∈q R(Ψ, F)    ≥    q^2/2(q-1/2)log(1/1-q) +q

≥    1/log(1/1-q) q^2/2(q-1/2) + 1/log(1/1-q)

≥    1/log(1/1-q) 1/4(2(1-1/2) + 1/log(1/1-1/2))

=    1/log(1/1-q) 1/4(1 + 1/log(2)) = c_3/log(1/1-q)    c_3 := 1/4(1 + 1/log(2))

UPPER BOUND
Let $q \in ]1/2, 1[$. We consider the family of distributions $(F_{r, 1, q})_{r \geq 1} \in \Fb_0 (w,\{q\})$ and the following weight distribution:

dλ(r)    =
0      if  r ≤1

1/log(1/1-q)(r _α(r | (0,1), (1,q) ))'/(r _α(r | (0,1), (1,q) ))      if  r ≥1

We can easily verify that $\int_0^{\infty} d\lambda(r)$ IS $ d\lambda(r) \ge 0$? . Define
G(p)    =    ∫_0^∞  p F_r, 1, q(p)/(F_r, 1, q) dλ(r)

Therefore, using Yao's principle <cit.>, we have:
sup_Ψ∈ inf_F ∈_0 (w,{q}) R(Ψ, F)    ≤    sup_p ≥0 G(p)

=    sup_p ≥0 ∫_0^∞  p F_r, 1, q(p)/(F_r, 1, q) dλ(r)

=    sup_p ≥0 ∫_p^∞  p F_r, 1, q(p)/(r _α(r | (0,1), (1,q) )) 1/log(1/1-q)(r _α(r | (0,1), (1,q) ))'/(r _α(r | (0,1), (1,q) )) dr

=    1/log(1/1-q) sup_p ≥0  p _α(p | (0,1), (1,q) ) ∫_p^∞ (r _α(r | (0,1), (1,q) ))'/(r _α(r | (0,1), (1,q) ))^2 dr

=    1/log(1/1-q) sup_p ≥0  p _α(p | (0,1), (1,q) ) (1/p _α(p | (0,1), (1,q) )

- lim_p →∞ 1/p _α(p | (0,1), (1,q) ))

=    1/log(1/1-q) sup_p ≥0  (1- (1/q-1) p _α(p | (0,1), (1,q) ))

=    1/log(1/1-q) (1- (1/q-1  ) 1 _α(1 | (0,1), (1,q) ) ) = 1/log(1/1-q) (1- (1/q-1) q)

=    q/log(1/1-q) ≤c_4/log(1/1-q)     c_4 = 1

§ UPPER BOUND LINEAR PROGRAM AND IMPLEMENTATION PARAMETERS

In this section, we show that one can obtain an upper bound on the maxmin ratio $\Rb (\MechSet, \aDistSet{[\ql, \qh]})$ by solving a linear program and provide the parameters used to obtain the results depicted in the main text.

Upper-bound LP:

For $\alpha \in [0,1]$, $N > 1$, $b \geq \pin$ and $\eta \in ]0, 1[$, consider the following finite sequence of prices  $\IN = (a_i)_{0 \leq i \leq 2N+1}$ in $[\rl[\ql] , \min(b, \rhone[\qh])]$:
a_i    =
[] + i/N (  1-η-[] )       0 ≤i ≤N.

+ i-N-1/N ( min(b, [])-1 )       N+1 ≤i ≤2N+1.

Fix a mechanism $\Psi \in \MechSet$. Using <ref>, let $\hat{\Psi}(u) = \Psi(u w)$ for $u \geq 0$ and $p_j := \int_{a_j}^{a_{j+1}} d\hat{\Psi}(u)$ for $0 \leq j \leq 2N$ and $p_{2N+1} := \int_{a_{2N+1}}^{\rhone[\qh]} d\hat{\Psi}(u)$. Then we have:
inf_F ∈[, ] R(Ψ, F)     =    min{ inf_x ∈[1/+1, 1]  1/x [ ∫_[] ^x  u dΨ̂(u) + ∫_x^1   u _α(u, x,1, 1,) dΨ̂(u) ],

inf_x ∈[1, 1/αq] 1/x  _α(x , 0,1, w,q ) ∫_[] ^x  u _α(u, 0,1, 1,)  dΨ̂(u) }

=    min{ min_i = 0 ⋯Ninf_x ∈[a_i, a_i+1[  1/x [ ∫_[] ^1   H^left_α, q(x, u) dΨ̂(u)  ],

min_i = N+1 ⋯2Ninf_x ∈[a_i, a_i+1[ 1/ H^right_α, q(x) ∫_[] ^x   H^right_α, q(u) dΨ̂(u),

inf_x ∈[a_2N+1, [][ 1/ H^right_α, q(x) ∫_[] ^x   H^right_α, q(u) dΨ̂(u),
Such that, for $(x,u) \in [\rl[\ql] , 1]^2$ and $t \in [1, \min(b, \rhone[\qh])]$
H^left_α, q(x, u)     =
u      [] ≤u≤x

u _α(u, x,1, 1,)       x ≤u≤1.

H^right_α, q(u)     =    u _α(u, 0,1, 1,)
Fix any $u \in [\rl[\ql] , 1]$, $x \to  H^{left}_{\alpha, q}(x, u) = u \Gad{\Gainv{q} \frac{u-x}{1-x}}$ is increasing by composition as $\Gad{x}$ and $x \to \frac{u-x}{1-x}$ are non-increasing. In turn, the revenue function $ H^{right}_{\alpha, q}$ is increasing in $[1, \min(b, \rhone[\qh])]$  as $\rhone[\qh]$ is the optimal reserve price. Using these properties, we have:
inf_F ∈[, ] R(Ψ, F) ≤min{ min_i = 0 ⋯N  1/a_i [ ∫_[] ^1   H^left_α, q(a_i+1, u) dΨ̂(u)  ],

min_i = N+1 ⋯2N 1/ H^right_α, q(a_i) ∫_[] ^a_i   H^right_α, q(u) dΨ̂(u) }

inf_x ∈[a_2N+1, [][ 1/ H^right_α, q(x) ∫_[] ^x   H^right_α, q(u) dΨ̂(u)

=    min{ min_i = 0 ⋯N 1/a_i ∑_j=0^i  ∫_a_j ^a_j+1 u dΨ̂(u)  + ∑_j=i+1^N ∫_a_j ^a_j+1 H^left_α, q(a_i+1, u) dΨ̂(u),

min_i = N+1 ⋯2N 1/a_i _α(a_i | (0, 1) ,(1 ,)) ∑_j=0^i-1  ∫_a_j ^a_j+1 u _α(u | (0, 1) ,(1 ,))  dΨ̂(u),

inf_x ∈[a_2N+1, [][ 1/ H^right_α, q(x) ∑_j=0^2N  ∫_a_j ^a_j+1 u _α(u | (0, 1) ,(1 ,))  dΨ̂(u)

+ ∫_a_2N+1 ^[] u _α(u | (0, 1) ,(1 ,))  dΨ̂(u)

≤    min{ min_i = 0 ⋯N 1/a_i ∑_j=0^i  a_j+1∫_a_j ^a_j+1 dΨ̂(u)  + ∑_j=i+1^N H^left_α, q(a_i+1, a_j) ∫_a_j ^a_j+1 dΨ̂(u),

min_i = N+1 ⋯2N 1/a_i _α(a_i | (0, 1) ,(1 ,)) ∑_j=0^i-1 a_j+1 _α(a_j+1 | (0, 1) ,(1 ,)) ∫_a_j ^a_j+1   dΨ̂(u),

inf_x ∈[a_2N+1, [][ 1/ H^right_α, q(a_2N+1) ∑_j=0^2N a_j+1 _α(a_j+1 | (0, 1) ,(1 ,)) ∫_a_j ^a_j+1   dΨ̂(u)

+ [] _α([] | (0, 1) ,(1 ,)) ∫_a_2N+1 ^[] dΨ̂(u)

=    min{ min_i = 0 ⋯N 1/a_i ∑_j=0^i  a_j+1p_j  + ∑_j=i+1^N H^left_α, q(a_i+1, a_j) p_j,

min_i = N+1 ⋯2N 1/a_i _α(a_i | (0, 1) ,(1 ,)) ∑_j=0^i-1 a_j+1 _α(a_j+1 | (0, 1) ,(1 ,)) p_j,

1/ H^right_α, q(a_2N+1) ∑_j=0^2N a_j+1 _α(a_j+1 | (0, 1) ,(1 ,)) p_j + [] _α([] | (0, 1) ,(1 ,)) p_2N+1

:=    c_U(α, , , , p)
Therefore, we have:
inf_F ∈[, ] R(Ψ, F)    ≤    max_p c_U(α, , , , p) = _U_α, , ,
_U_α, , ,     =   maxc

1/a_i [∑_j=0^i  a_j+1p_j  + ∑_j=i+1^N  a_j _α(a_j | (a_i+1, 1) ,(1 ,)) p_j] ≥c   i=0,...N,

1/a_i _α(a_i | (0, 1) ,(1 ,)) ∑_j=0^i-1 a_j+1 _α(a_j+1 | (0, 1) ,(1 ,)) p_j  ≥c   i=N+1,...2N-1,

1/a_2N+1 _α(a_2N+1 | (0, 1) ,(1 ,)) [ ∑_j=0^2N a_j+1 _α(a_j+1 | (0, 1) ,(1 ,)) p_j
+ [] _α([] | (0, 1) ,(1 ,))] p_2N+1 ≥c

∑_j=0^2N p_j ≤1

p_i ≥0   i=0,...2N-1,

By taking the sup over the mechanisms $\Psi$, we get:
(, [, ])    =     sup_Ψ∈ inf_F ∈[, ] R(Ψ, F) ≤_U_α, , ,
Therefore, we obtain
(, [, ]) ≤_U_α, , ,

Implementation parameters:
For all reported values in the main text, we use the following sequence in the LPs
a_i    =
[] + i/N (  1-η-[] )       0 ≤i ≤N.

+ i-N-1/N ( min(b, [])-1 )       N+1 ≤i ≤2N+1.

With $N = 2500, \eta = 10^{-5}, b = 250$.

§ ILLUSTRATION OF NEAR OPTIMAL MECHANISMS FOR SECTION <REF>

[ xmin=0, xmax=250,
ymin=0, ymax=1,
restrict y to domain=0:1.0,
minor tick num=1,
axis line style=->,,
xlabel= $p$,
ylabel= $\Psi(p)$,
width=7cm, legend pos=north west,
height=7cm, cycle list name=black white]
[blue, very thick ]  table [ col sep=comma ] optmech_001.csv;
[red, very thick]  table [ col sep=comma ] optmech_025.csv;
[black!60!green, very thick ]  table [ col sep=comma ] optmech_050.csv;
[black!100, very thick ]  table [ col sep=comma ] optmech_075.csv;

*regular distributions

[ xmin=0, xmax=4,
ymin=0, ymax=1,
restrict y to domain=0:1.0,
minor tick num=1,
axis line style=->,
xlabel= $p$,
ylabel= $\Psi(p)$,
width=7cm, legend pos=south east,
height=7cm, cycle list name=black white]
[blue, very thick ]  table [ col sep=comma ] optmech_101.csv;
[red, very thick ]  table [ col sep=comma ] optmech_125.csv;
[black!60!green, very thick ]  table [ col sep=comma ] optmech_150.csv;
[black!100, very thick ]  table [ col sep=comma ] optmech_175.csv;

*mhr distributions

Illustration of Near Optimal Mechanisms. The figure depicts near optimal pricing distributions for $\pin = 1$, $q \in \{0.01, 0.25, 0.5, 0.75\}$. The left panel corresponds to regular distributions (plotted using a log scale) and the right panel to  mhr distributions (on a regular scale).


∎

11institutetext: Lionel Tondji 22institutetext: Institute for Analysis and
Algebra, TU Braunschweig, 38092 Braunschweig, Germany,
22email<EMAIL_ADDRESS>33institutetext: Ion Necoara
44institutetext: Automatic Control and Systems Engineering Department,
National University of Science and Technology Politehnica Bucharest, 060042
Bucharest, Romania. 55institutetext: Gheorghe Mihoc-Caius Iacob Institute of
Mathematical Statistics and Applied Mathematics of the Romanian Academy,
050711 Bucharest, Romania,
55email<EMAIL_ADDRESS>66institutetext: Dirk A. Lorenz 77institutetext:
Institute for Analysis and Algebra, TU Braunschweig, 38092 Braunschweig,
Germany. 88institutetext: Center for Industrial Mathematics, University of
Bremen, Postfach 330440, 28334 Bremen, Germany,
88email<EMAIL_ADDRESS>
# Acceleration and restart for the randomized Bregman-Kaczmarz method††thanks:
This work has received funding from: the European Union’s Framework Programme
for Research and Innovation Horizon 2020 (2014-2020) under the Marie
Sklodowska-Curie Grant Agreement No. 861137; UEFISCDI PN-III-P4-PCE-2021-0720,
under project L2O-MOC, no 70/2022.

Lionel Tondji  Ion Necoara  Dirk A. Lorenz

###### Abstract

Optimizing strongly convex functions subject to linear constraints is a
fundamental problem with numerous applications. In this work, we propose a
block (accelerated) randomized Bregman-Kaczmarz method that only uses a block
of constraints in each iteration to tackle this problem. We consider a dual
formulation of this problem in order to deal in an efficient way with the
linear constraints. Using convex tools, we show that the corresponding dual
function satisfies the Polyak-Lojasiewicz (PL) property, provided that the
primal objective function is strongly convex and verifies additionally some
other mild assumptions. However, adapting the existing theory on coordinate
descent methods to our dual formulation can only give us sublinear convergence
results in the dual space. In order to obtain convergence results in some
criterion corresponding to the primal (original) problem, we transfer our
algorithm to the primal space, which combined with the PL property allows us
to get linear convergence rates. More specifically, we provide a theoretical
analysis of the convergence of our proposed method under different assumptions
on the objective and demonstrate in the numerical experiments its superior
efficiency and speed up compared to existing methods for the same problem.

###### Keywords:

Linear systems, sparse solutions, randomized Bregman-Kaczmarz method,
acceleration, restart, convergence rates.

###### MSC:

65F10, 68W20, 90C25

##  1 Introduction

In this work, we consider the fundamental problem of approximating solutions
of large-scale linear systems of the form:

$\mathbf{A}x=b$ (1.1)

with matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$ and right hand side
$b\in\mathbb{R}^{m}$. We consider the case that the full matrix is not
accessible simultaneously, but that one can only work with a single or block
of rows of the system (1.1) at a time. Given the possibility of multiple
solutions of (1.1), we set out to find the unique solution, characterized by
the function $f$, i.e.,

$\hat{f}\overset{\text{def}}{=}\min_{x\in\mathbb{R}^{n}}f(x)\quad\text{subject
to}\quad\mathbf{A}x=b,$ (1.2)

with a general strongly convex function $f$. However, we will not assume
smoothness (not even differentiability) of $f$. One possible example is
$f(x)=\lambda\cdot\|x\|_{1}+\frac{1}{2}\|x\|^{2}_{2}$ and it is known that
this function favors sparse solutions for appropriate choices of $\lambda>0$,
see cai2009linearized ; LSW14 ; LWSM14 ; LS19 ; yin2010analysis . Similarly,
by dividing the components of $x$ into $M$ groups, $x=(x_{(1)},\dots,x_{(M)})$
with $x_{(j)}\in\mathbb{R}^{n_{j}}$, the function
$f(x)=\lambda\cdot\sum_{j=1}^{M}\|x_{(j)}\|_{2}+\frac{1}{2}\cdot\|x\|_{2}^{2}$
favours group sparsity. We assume throughout the paper that both $m$ and $n$
are large and that the system is consistent. The $i$-th row of $\mathbf{A}$ is
denoted by $a_{i}^{\top}$ and we assume $a_{i}\neq 0$ for all
$i\in[m]:=\\{1,\dots,m\\}$. Since $f$ is strongly convex, problem (1.2) has a
unique solution $\hat{x}$. In applications, it is usually sufficient to find a
point which is not too far from $\hat{x}$. In particular, one chooses the
error tolerance $\varepsilon>0$ and aims to find a point $x$ satisfying
$\|x-\hat{x}\|^{2}_{2}\leq\varepsilon$. Since our method is of a stochastic
nature, the iterates $x$ are random vectors. Hence, our goal will be to
compute approximate solutions of (1.2) that fulfills
$\mathbb{E}[\|x-\hat{x}\|^{2}_{2}]\leq\varepsilon$, where $\mathbb{E}[\cdot]$
denotes the expectation w.r.t. randomness of the algorithm.

###  1.1 Related work

The linear system (1.1) may be so large that full matrix operations are very
expensive or even infeasible. Then, it appears desirable to use iterative
algorithms with low computational cost and storage per iteration that produce
good approximate solutions of (1.2) after relatively few iterations. The
Kaczmarz method and its randomized variants gower2019adaptive ;
gower2015randomized ; necoara2019faster ; SV09  are used to compute the
minimum $\ell_{2}$-norm solutions of consistent linear systems. In each
iteration $k$, a row vector $a_{i}^{\top}$ of $\mathbf{A}$ is chosen at random
from the system (1.1) and the current iterate $x^{k}$ is projected onto the
solution space of that equation to obtain $x^{k+1}$. Note that this update
rule requires only low cost per iteration and storage of order
$\mathcal{O}(n)$. Recently, a new variant of the Randomized Kaczmarz (RK)
method namely the Randomized Sparse Kaczmarz method (RSK) LWSM14 ; LS19  with
almost the same low cost and storage requirements has shown good performance
in approximating sparse solutions of large consistent linear systems. The
papers LSW14 ; LS19  analyze this method by interpreting it as a sequential,
randomized Bregman projection method (where the Bregman projection is done
with respect to the function $f$), while P15  connects it with the coordinate
descent method via duality. Variations of RSK including block/averaging
variants can be found in du2020randomized ; necoara2019faster ;
needell2014paved ; P15 ; tondji2022faster ; tondji2023adaptive , momentum
variants lorenz2023minimal , sampling scheme tondji2021linear  and adaptations
to least squares problems are given in Du19 ; schopfer2022extended ;
zouzias2013randomized . In those variants, one usually needs to have access to
more than one row of the matrix $\mathbf{A}$ at the cost of increasing the
memory. Coordinate descent methods are often considered in the context of
minimizing composite convex objective functions fercoq2015accelerated ;
necoara2016 ; nesterov2012efficiency . At each iteration, they update only one
coordinate of the vector of variables using partial derivatives and proximal
operators rather than the whole gradient and the full proximal operator,
respectively. A randomized accelerated coordinate descent method was proposed
in nesterov2012efficiency  for smooth convex functions and later it was
extended to composite convex functions in e.g., fercoq2015accelerated ;
necoara2016 . Note that accelerated methods transform the coordinate descent,
for which the optimality gap decreases as $\mathcal{O}(1/k)$, into an
algorithm with optimal $\mathcal{O}(1/k^{2})$ complexity under mild additional
computational cost fercoq2015accelerated ; nesterov2012efficiency . By
exploiting the connection between coordinate descent methods and Kaczmarz
methods, an accelerated method was proposed in 10178390  to solve problem
(1.2) with $\mathcal{O}(1/k^{2})$ convergence rate. The present study builds
upon the initial findings in 10178390  and expands upon them with additional
restart algorithm, new theoretical analysis for the primal and dual problem,
new linear convergence results and new experiments.

###  1.2 Contribution

To the best of our knowledge, accelerated Bregman-Kaczmarz methods for solving
linearly constrained optimization problems with strongly convex (possibly
nonsmooth) objective function have not been considered before in the
literature. In this work, we provide such methods for solving this type of
problems. Hence, this study has the following main contributions:

  1. (a)

Beyond interpreting the Bregman-Kaczmarz method as a dual coordinate descent,
we propose a primal accelerated coordinate gradient algorithm using only
blocks of rows of the matrix defining the constraints. We refer to it as Block
Accelerated Randomized Bregman-Kaczmarz method (ARBK) with more details in
Algorithm 2. We have also extended the method and propose a restart scheme for
Algorithm 2 that exhibits faster convergence.

  2. (b)

Using tools from convex theory, we show that the corresponding dual function
satisfies the Polyak-Lojasiewicz (PL) property, provided that the primal
objective function is strongly convex and verifies additionally some other
mild assumptions. By exploiting the connection between primal and dual
updates, we obtain convergence as a byproduct and convergence rates under
different assumptions on the objective function which have not been available
so far. In particular, we prove that our restart accelerated method leads to
faster (i.e., linear rate) convergence than its standard counterpart.

  3. (c)

We also validate the superior efficiency of our approach empirically and we
provide an implementation of our algorithm in Python.

###  1.3 Outline

The remainder of the paper is organized as follows. Section 2 provides
notations and a brief overview on convexity and Bregman distances. In Section
3 we state our method and give an interpretation of it in the dual space.
Section 4 provides convergence guarantees for our proposed method. In Section
5, numerical experiments demonstrate the effectiveness of our method and
provide insight regarding its behavior and its hyper-parameters. Finally,
Section 6 draws some conclusions.

##  2 Notation and basic notions

For a positive integer $n$ we denote
$[n]\overset{\text{def}}{=}\\{1,2,\dots,n\\}$. Given a symmetric positive
definite matrix $\mathbf{B},$ we denote the induced inner product by

$\langle x,y\rangle_{\mathbf{B}}\overset{\text{def}}{=}\langle
x,\mathbf{B}y\rangle=\sum\limits_{i,j\in[n]}x_{(i)}\mathbf{B}_{ij}y_{(j)},\quad
x,y\in\mathbb{R}^{n},$

and its induced norm by
$\|.\|_{\mathbf{B}}^{2}\overset{\text{def}}{=}\langle\cdot,\cdot\rangle_{\mathbf{B}}$.
We use the short-hand notation $\|.\|_{2}$ to mean $\|.\|_{\mathbf{I}_{n}}$.
For an $n\times m$ real matrix $\mathbf{A}$ we denote by
$\mathcal{R}(\mathbf{A}),\|\mathbf{A}\|_{F}$ and $a_{i}^{\top}$ its range
space, its Frobenius norm and its $i$-th row, respectively. By $e$ we denote
the Euler’s number and by $\lceil x\rceil$ we denote the smallest integer
greater than or equal to $x$ with $x\in\mathbb{R}$. We use
$\textbf{Diag}(d_{1},d_{2},\dots,d_{M})$ to denote the diagonal matrix with
$d_{1},d_{2},\dots,d_{M}$ on the diagonal. For a random vector $x_{i}$ that
depends on a random index $i$ (where $i\in[n]$ is chosen with probability
$p_{i}$) we denote
$\mathbb{E}[x_{i}]\overset{\text{def}}{=}\sum_{i\in[n]}p_{i}x_{i}$ and we will
just write $\mathbb{E}[x_{i}]$ when the probability distribution is clear from
the context. Given a vector $x\in\mathbb{R}^{n}$, we define the soft shrinkage
operator, which acts componentwise on a vector $x$ as

$(S_{\lambda}(x))_{j}=\max\\{|x_{j}|-\lambda,0\\}\cdot\operatorname{sign}(x_{j})\quad\forall
j\in[n].$ (2.1)

We will assume that the vector $y\in\mathbb{R}^{m}$ has the following block
partition $y=(y_{(1)},y_{(2)},\dots,y_{(M)})$, composed of $M$ blocks, where
$y_{(i)}\in\mathbb{R}^{m_{i}}$ and $\sum_{i\in[M]}m_{i}=m$. Using the matrices
$\mathbf{U}_{i}\in\mathbb{R}^{m\times m_{i}},i\in[M]$ which partition the
identity as
$\mathbf{I}_{m}=(\mathbf{U}_{1},\mathbf{U}_{2},\dots,\mathbf{U}_{M})\in\mathbb{R}^{m\times
m}$, we get

$y_{(i)}=\mathbf{U}_{i}^{\top}y\in\mathbb{R}^{m_{i}},\;\;\forall
y\in\mathbb{R}^{m},\;\;\forall i\in[M]\quad\text{and}\quad
y=\sum_{i\in[M]}\mathbf{U}_{i}y_{(i)}.$

Similarly, for a given differentiable function
$g:\mathbb{R}^{m}\to\mathbb{R}$, the vector of partial derivatives
corresponding to the variables in $y_{(i)}$ is
$\nabla_{(i)}g(y)=\mathbf{U}^{\top}_{i}\nabla g(y)\in\mathbb{R}^{m_{i}}$. A
similar partition is applied for $x\in\mathbb{R}^{n}$ and with
$\mathbf{V}_{i}\in\mathbb{R}^{n\times n_{i}}$ where $\sum_{i\in[M]}n_{i}=n$ we
write
$\mathbf{I}_{n}=(\mathbf{V}_{1},\mathbf{V}_{2},\dots,\mathbf{V}_{M})\in\mathbb{R}^{n\times
n}$. Hence, we get

$x_{(i)}=\mathbf{V}_{i}^{\top}x\in\mathbb{R}^{n_{i}},\forall
x\in\mathbb{R}^{n},\forall i\in[M]\quad\text{and}\quad
x=\sum_{i\in[M]}\mathbf{V}_{i}x_{(i)}.$

Now, we collect some basic notions on convexity and the Bregman distance, our
exposition mainly follows from nesterov2003introductory ; Roc:70 .

###  2.1 Basic Notions

Let $f:\mathbb{R}^{n}\to\mathbb{R}$ be convex (note that we assume that $f$ is
finite everywhere, hence it is also continuous). The _subdifferential_ of $f$
at any $x\in\mathbb{R}^{n}$ is defined by

$\partial f(x)\overset{\text{def}}{=}\\{x^{*}\in\mathbb{R}^{n}|f(y)\geq
f(x)+\langle x^{*},y-x\rangle\;\;\forall y\in\mathbb{R}^{n}\\},$

which is nonempty, compact and convex set. The function
$f:\mathbb{R}^{n}\to\mathbb{R}$ is said to be _$\sigma$ -strongly convex_ with
$\sigma>0$, if for all $x,y\in\mathbb{R}^{n}$ and subgradients
$x^{*}\in\partial f(x)$ we have

$f(y)\geq f(x)+\langle
x^{*},y-x\rangle+\tfrac{\sigma}{2}\cdot\|y-x\|_{2}^{2}\,.$

If $f$ is $\sigma$-strongly convex, then $f$ is coercive, i.e.

$\lim_{\|x\|_{2}\to\infty}f(x)=\infty\,,$

The Fenchel conjugate of an $\sigma$-strongly convex function $f$, denoted
$f^{*}:\mathbb{R}^{n}\to\mathbb{R}$, given by

$f^{*}(x^{*})\overset{\text{def}}{=}\sup_{y\in\mathbb{R}^{n}}\langle
x^{*}\,,\,y\rangle-f(y),$

is also convex, finite everywhere and coercive. Additionally, $f^{*}$ is
differentiable with a _Lipschitz-continuous gradient_ with constant
$L_{f^{*}}=\frac{1}{\sigma}$, i.e., for all $x^{*},y^{*}\in\mathbb{R}^{n}$ we
have

$\|\nabla f^{*}(x^{*})-\nabla f^{*}(y^{*})\|_{2}\leq
L_{f^{*}}\cdot\|x^{*}-y^{*}\|_{2}\,,$

which implies the estimate

$\displaystyle f^{*}(y^{*})\leq f^{*}(x^{*})\\!+\langle\nabla
f^{*}(x^{*})\,,\,y^{*}-x^{*}\rangle+\tfrac{L_{f^{*}}}{2}\|x^{*}-y^{*}\|_{2}^{2}.$
(2.2)

The function $f^{*}$ is said to have block-coordinate Lipschitz continuous
gradient if there exists $L_{f^{*},i}>0$ such that

$\displaystyle\|\nabla_{(i)}f^{*}(x^{*}+V_{i}h_{(i)})-\nabla_{(i)}f^{*}(x^{*})\|\leq
L_{f^{*},i}\cdot\|h_{(i)}\|,$ (2.3)

for all $x^{*}\in\mathbb{R}^{n},h_{(i)}\in\mathbb{R}^{n_{i}},i\in[M].$

###### Definition 1

The _Bregman distance_ $D_{f}^{x^{*}}(x,y)$ between $x,y\in\mathbb{R}^{n}$
with respect to a convex function $f$ and a subgradient $x^{*}\in\partial
f(x)$ is defined as

$D_{f}^{x^{*}}(x,y)\overset{\text{def}}{=}f(y)-f(x)-\langle
x^{*}\,,\,y-x\rangle\,.$

Fenchel’s equality states that $f(x)+f^{*}(x^{*})=\langle x\,,\,x^{*}\rangle$
if $x^{*}\in\partial f(x)$ and implies that the Bregman distance can be
written as

$D_{f}^{x^{*}}(x,y)=f^{*}(x^{*})-\langle x^{*}\,,\,y\rangle+f(y)\,.$

######  Example 1 ((LS19, , Example 2.3))

The objective function

$f(x)\overset{\text{def}}{=}\lambda\cdot\|x\|_{1}+\tfrac{1}{2}\cdot\|x\|_{2}^{2}$
(2.4)

is strongly convex with constant $\sigma=1$ and its conjugate function can be
computed with the soft shrinkage operator from (2.1)

$f^{*}(x^{*})=\tfrac{1}{2}\cdot\|S_{\lambda}(x^{*})\|_{2}^{2},\quad\mbox{with}\quad\nabla
f^{*}(x^{*})=S_{\lambda}(x^{*})\,.$

Its Fenchel conjugate $f^{*}$ has componentwise Lipschitz gradient with
constants $L_{f^{*},i}=1$ and for any $x^{*}=x+\lambda\cdot s\in\partial
f(x)$, with $s\in\partial(\|x\|_{1})$, we have

$D_{f}^{x^{*}}(x,y)=\frac{1}{2}\cdot\|x-y\|_{2}^{2}+\lambda\cdot(\|y\|_{1}-\langle
s\,,\,y\rangle)\,.$

In particular, $D_{f}^{x^{*}}(x,y)=\frac{1}{2}\|x-y\|_{2}^{2}$ for
$\lambda=0.$

The following inequality is crucial for the convergence analysis of the
randomized algorithms. It immediately follows from the definition of the
Bregman distance and the assumption of strong convexity of $f$, see LSW14 .
For $x,y\in\mathbb{R}^{n}$ and $x^{*}\in\partial f(x)$ we have

$\displaystyle D_{f}^{x^{*}}(x,y)\geq\frac{\sigma}{2}\|x-y\|_{2}^{2}.$ (2.5)

###  2.2 Calmness, Linear regularity and Linear growth conditions

In this section, we are going to define notions such as calmness, linear
regularity, and linear growth schopfer2022extended , which are necessary
conditions for error bounds to hold. Let $B_{2}$ denote the closed unit ball
of the $\ell_{2}$-norm.

###### Definition 2

The (set-valued) subdifferential mapping $\partial
f:\mathbb{R}^{n}\rightrightarrows\mathbb{R}^{n}$ is _calm_ at $\hat{x}$ if
there are constants $\varepsilon,L>0$ such that

$\partial f(x)\subset\partial f(\hat{x})+L\cdot\|x-\hat{x}\|_{2}\cdot
B_{2}\quad\mbox{for any $x$ with}\quad\|x-\hat{x}\|_{2}\leq\varepsilon\,.$
(2.6)

Note that calmness is a local growth condition akin to Lipschitz-continuity of
a gradient mapping, but for fixed $\hat{x}$. Of course, any Lipschitz-
continuous gradient mapping is calm everywhere.

###### Example 2

  * (a)

The subdifferential mapping of any convex piecewise linear or quadratic
function is calm everywhere. In particular, this holds for the function
defined in eq (2.4).

  * (b)

The subdifferential mapping of
$f(x)=\lambda\cdot\|x\|_{2}+\tfrac{1}{2}\cdot\|x\|_{2}^{2}$ is calm
everywhere. For more examples on calmness, we refer the reader to
schopfer2022extended .

###### Definition 3

Let $\partial f(\hat{x})\cap\mathcal{R}(A^{\top})\neq\emptyset.$ Then, the
collection $\\{\partial f(\hat{x}),\mathcal{R}(\mathbf{A}^{\top})\\}$ is
_linearly regular_ , if there is a constant $\nu>0$ such that for all
$x^{*}\in\mathbb{R}^{n}$ we have

$\operatorname{dist}(x^{*},\partial
f(\hat{x})\cap\mathcal{R}(\mathbf{A}^{\top}))\leq\nu\cdot\bigg{(}\operatorname{dist}(x^{*},\partial
f(\hat{x}))+\operatorname{dist}(x^{*},\mathcal{R}(\mathbf{A}^{\top}))\bigg{)}.$
(2.7)

The collection $\\{\partial f(\hat{x}),\mathcal{R}(\mathbf{A}^{\top})\\}$ is
always linearly regular if, for example, $\partial f(\hat{x})$ is polyhedral
(which holds e.g., for piecewise linear/quadratic functions $f$).

###### Definition 4

We say the subdifferential mapping of $f$ _grows at most linearly_ if there
exist $\rho_{1},\rho_{2}\geq 0$ such that for all $x\in\mathbb{R}^{n}$ and
$x^{*}\in\partial f(x)$ we have

$\|x^{*}\|_{2}\leq\rho_{1}\cdot\|x\|_{2}+\rho_{2}\,.$ (2.8)

Any Lipschitz-continuous gradient mapping grows at most linearly. Furthermore,
the subdifferential mappings of the functions in Example 2 grow at most
linearly. Next, we give the definition of the Polyak-Lojasiewicz (PL)
inequality which is a lower bound on the size of the gradient as the value of
$f$ increases necoara2019linear .

###### Definition 5

A differentiable function satisfies PL inequality if there exists some $\mu>0$
such that

$f(x)-\hat{f}\leq\frac{1}{2\mu}\|\nabla f(x)\|^{2}_{2},\quad\forall
x\,\in\text{dom}(f).$ (2.9)

Any strongly-convex function satisfies the PL inequality. However, there are
many important functions satisfying the PL inequality which are not strongly
convex (see e.g., Theorem 4.3 below).

##  3 Coordinate descent and Bregman-Kaczmarz methods

Note that by a proper scalling of $f$ we can always assume the strong
convexity constant $\sigma=1$. Therefore, in the sequel we consider 1-strongly
convex objective function $f$. To connect the Kaczmarz method to coordinate
descent algorithm, we consider the dual problem to problem (1.2). Using the
conjugate of $f$ we can write the dual objective as infimum of the Lagrangian
$\mathcal{L}(x,y)\overset{\text{def}}{=}f(x)-y^{\top}(\mathbf{A}x-b)$ as

$\displaystyle\inf_{x}\mathcal{L}(x,y)$
$\displaystyle=\inf_{x}(f(x)-y^{\top}(\mathbf{A}x-b))=b^{\top}y-f^{*}(\mathbf{A}^{\top}y).$

Hence, the dual problem of (1.2) is given by:

$\hat{\Psi}\overset{\text{def}}{=}\min_{y\in\mathbb{R}^{m}}\left[\Psi(y):=f^{*}(\mathbf{A}^{\top}y)-b^{\top}y\right].$
(3.1)

The primal-dual optimality system

$\displaystyle\mathbf{A}\hat{x}=b,\quad\hat{x}=\nabla
f^{*}(\mathbf{A}^{\top}\hat{y}),$ (3.2)

coupled the primal and dual optimal solutions $\hat{x}$ and $\hat{y}$ via
$\hat{x}=\nabla f^{*}(\mathbf{A}^{\top}\hat{y})$ and we also see that the
optimal solution set $\mathcal{Y}^{*}$ of (3.1) is nonempty. Further, the dual
objective $\Psi$ is unconstrained, differentiable and its gradient is given by
the following expression

$\nabla\Psi(y)=\mathbf{A}\nabla f^{*}(\mathbf{A}^{\top}y)-b.$

Moreover, the gradient $\nabla\Psi$ of the dual function is Lipschitz and
block-wise Lipschitz continuous w.r.t. the Euclidean norm $\|\cdot\|_{2}$,
with constants necoara2022linear :

$L_{\Psi}\overset{\text{def}}{=}\|\mathbf{A}\|^{2}_{2}\ \text{ and }\
L_{\Psi,i}\overset{\text{def}}{=}L_{i}\overset{\text{def}}{=}\|\mathbf{A}_{(i)}\|^{2},$
(3.3)

respectively. The following lemma relates the Bregman distance to the optimal
primal solution to the distance to optimality of the dual problem.

###### Lemma 1

Let $b\in\mathcal{R}(\mathbf{A})$ and $(y^{k},x^{k})$ be two sequences such
that $x^{k}=\nabla f^{*}(\mathbf{A}^{\top}y^{k})$. Then, it holds

$\displaystyle D_{f}^{d^{k}}(x^{k},\hat{x})=\Psi(y^{k})-\hat{\Psi}.$ (3.4)

###### Proof

By Definition 1 we have that for $d^{k}=\mathbf{A}^{\top}y^{k}$ the following
is valid:

$\displaystyle D_{f}^{d^{k}}(x^{k},\hat{x})$
$\displaystyle=f^{*}(d^{k})+f(\hat{x})-\langle d^{k},\hat{x}\rangle$
$\displaystyle=f^{*}(\mathbf{A}^{\top}y^{k})-\langle
b,y^{k}\rangle+f(\hat{x})$ $\displaystyle=\Psi(y^{k})+\hat{f},$

and since $\hat{\Psi}=-\max-\Psi=-\hat{f}$ (by strong duality), the assertion
follows.

The randomized block coordinate descent update applied to minimize $\Psi$ from
(3.1) reads, see necoara2016 ; nesterov2012efficiency :

$\displaystyle y^{k+1}$
$\displaystyle=y^{k}-\frac{1}{\|\mathbf{A}_{(i)}\|^{2}}\mathbf{U}_{i}\bigg{(}\mathbf{A}_{(i)}\nabla
f^{*}(\mathbf{A}^{\top}y^{k})-b_{(i)}\bigg{)}$
$\displaystyle=y^{k}-\frac{1}{\|\mathbf{A}_{(i)}\|^{2}}\mathbf{U}_{i}\bigg{(}\mathbf{A}_{(i)}x^{k}-b_{(i)}\bigg{)},$
(3.5)

where $i=i_{k}\in[M]$ is a block row index chosen randomly. For the sampling
probabilities of the blocks of rows we will assume the following model which
depends on a parameter $\alpha\in[0,1]$:

$p_{\alpha,M}(i)\overset{\text{def}}{=}p_{i}=L_{i}^{\alpha}\cdot\bigg{[}\sum\limits_{j=1}^{M}L_{j}^{\alpha}\bigg{]}^{-1},$
(3.6)

with $L_{i}$ from (3.3). The case $\alpha=0$ corresponds to uniform
probabilities, and the case $\alpha=1$ corresponds to the case where the
$i$-th block is sampled with probability proportional to the block Lischitz
constant. Multiplying both sides in (3.5) by $\mathbf{A}^{\top}$, we have

$\mathbf{A}^{\top}y^{k+1}=\mathbf{A}^{\top}y^{k}-\frac{1}{\|\mathbf{A}_{(i)}\|^{2}}\mathbf{A}^{\top}\mathbf{U}_{i}\bigg{(}\mathbf{A}_{(i)}x^{k}-b_{(i)}\bigg{)}.$

Using the notations

$\displaystyle d^{k}=\mathbf{A}^{\top}y^{k},\quad x^{k}=\nabla f^{*}(d^{k}),$
(3.7)

and the fact that $\mathbf{A}_{(i)}=\mathbf{U}_{i}^{\top}\mathbf{A}$ the block
dual coordinate descent iterates transfer in the primal space to block
Bregman-Kaczmarz iterates as shown in Algorithm 1. Note that in Algorithm 1,
setting $f(x)=\tfrac{1}{2}\|x\|^{2}_{2}$ gives us the block Randomized
Kaczmarz (RK) method gower2019adaptive ; necoara2019faster ; SV09 , while
$f(x)=\lambda\|x\|_{1}+\tfrac{1}{2}\|x\|^{2}_{2}$ gives us a block variant of
the Randomized Sparse Kaczmarz (RSK) method LWSM14 ; LS19 .

Algorithm 1 Block Bregman-Kaczmarz method (BK)

1: choose $d^{0}\in\mathbb{R}^{n}$ and set $x^{0}=\nabla f^{*}(d^{0})$. Choose
the number of block $M$

2: Output: (approximate) solution of
$\min\limits_{x\in\mathbb{R}^{n}}f(x)\quad\text{s.t.}\quad\mathbf{A}x=b$.

3: initialize $k=0$

4: repeat

5: choose a block row index $i=i_{k}\in[M]$ randomly

6: update
$d^{k+1}=d^{k}-\frac{1}{\|\mathbf{A}_{(i)}\|^{2}}\mathbf{A}^{\top}_{(i)}\bigg{(}\mathbf{A}_{(i)}x^{k}-b_{(i)}\bigg{)}$

7: update $x^{k+1}=\nabla f^{*}(d^{k+1})$

8: increment $k=k+1$

9: until a stopping criterion is satisfied

10: return $x^{k+1}$

A more general version of randomized coordinate descent, namely APPROX has
been proposed in fercoq2015accelerated  to accelerate proximal coordinate
descent methods for the minimization of composite functions. The APPROX update
applied to solve our dual problem (3.1) is given by:

$\displaystyle v^{k}$ $\displaystyle=(1-\theta_{k})y^{k}+\theta_{k}z^{k}$
(3.8) $\displaystyle z^{k+1}_{(i)}$
$\displaystyle=z^{k}_{(i)}-\frac{1}{\theta_{k}M\|\mathbf{A}_{(i)}\|^{2}}\cdot\bigg{(}\mathbf{A}_{(i)}\nabla
f^{*}(\mathbf{A}^{\top}v^{k})-b_{(i)}\bigg{)}$ $\displaystyle y^{k+1}$
$\displaystyle=v^{k}+M\theta_{k}(z^{k+1}-z^{k})$ $\displaystyle\theta_{k+1}$
$\displaystyle=\frac{\sqrt{\theta_{k}^{4}+4\theta_{k}^{2}}-\theta_{k}^{2}}{2}.$

In the updates (3.8), $M$ and $\theta_{k}$ represent respectively the number
of blocks and a sequence used for interpolation. The number of blocks is fixed
by the user and the sequence $\theta_{k}$ is initialised with
$\theta_{0}=\frac{1}{M}$. The second equation in (3.8) can be translated to
full vector operation as

$z^{k+1}=z^{k}-\frac{1}{\theta_{k}M\|\mathbf{A}_{(i)}\|^{2}}\mathbf{U}_{i}\cdot\bigg{(}\mathbf{A}_{(i)}\nabla
f^{*}(\mathbf{A}^{\top}v^{k})-b_{(i)}\bigg{)}.$

We transfer this algorithm to the primal space using the relation
$x^{k}=\nabla f^{*}(\mathbf{A}^{\top}y^{k})$ and this results in a primal
block accelerated randomized Bregman Kaczmarz scheme, see Algorithm 2.

Algorithm 2 Block Accelerated Randomized Bregman Kaczmarz method (ARBK)

1: Choose $y^{0}\in\mathbb{R}^{m}$ and set $z^{0}=y^{0}$. Choose the number of
block $M$ and set $\theta_{0}=\frac{1}{M}$, $b\in\mathbb{R}^{m}$,
$\mathbf{A}\in\mathbb{R}^{m\times n}$.

2: Output: (approximate) solution $x^{k+1}=\nabla
f^{*}(\mathbf{A}^{\top}y^{k+1})$ of
$\min\limits_{x\in\mathbb{R}^{n}:\mathbf{A}x=b}f(x)$

3: initialize $k=0$

4: repeat

5: update $v^{k}=(1-\theta_{k})y^{k}+\theta_{k}z^{k}$

6: choose a block row index $i=i_{k}\in[M]$ randomly

7: update
$z^{k+1}=z^{k}-\frac{1}{M\theta_{k}}\frac{\mathbf{U}_{i}\cdot\big{(}\mathbf{A}_{(i)}\nabla
f^{*}(\mathbf{A}^{\top}v^{k})-b_{(i)}\big{)}}{\|\mathbf{A}_{(i)}\|^{2}}$

8: update $y^{k+1}=v^{k}+M\theta_{k}(z^{k+1}-z^{k})$

9: update
$\theta_{k+1}=\dfrac{\sqrt{\theta_{k}^{4}+4\theta_{k}^{2}}-\theta_{k}^{2}}{2}$

10: increment $k=k+1$

11: until a stopping criterion is satisfied

12: return $y^{k+1}$

###### Remark 1

We have written the APPROX (3.8) and ARBK (Algorithm 2) in a unified framework
to emphasize their similarities. However, the two methods are used for two
different purposes. The APPROX method (3.8) outputs $y^{k}$, an approximate
solution of the dual problem (3.1), whereas Algorithm 2 returns $x^{k}$, an
approximate solution of our primal problem (1.2).

We also use the following relation of the sequences, that is, the iterates of
APPROX (3.8) and Algorithm 2 satisfy for all $k\geq 1$,

$d^{k}=\mathbf{A}^{\top}y^{k},\quad x^{k}=\nabla f^{*}(d^{k}).$

Note that Algorithm 1 can be recovered from Algorithm 2 by setting
$\theta_{k}=\theta_{0}$ for all $k\geq 0$.

##  4 Convergence Analysis

In this section we first review a basic convergence result for APPROX from
(3.8) that is valid in the dual space, and which will be used later to build
convergence results for the proposed method Algorithm 2 in the primal space.
We first recall the following properties on the sequence
$\\{\theta_{k}\\}_{k\geq 0}$.

######  Lemma 2 ((tseng2008accelerated, , Note 2))

The sequence $\theta_{k}$ defined by $\theta_{0}\leq 1$ and
$\theta_{k+1}=\frac{\sqrt{\theta_{k}^{4}+4\theta_{k}^{2}}-\theta_{k}^{2}}{2}$
satisfies

$\displaystyle\frac{(2-\theta_{0})}{k+(2-\theta_{0})/\theta_{0}}$
$\displaystyle\leq\theta_{k}\leq\frac{2}{k+2/\theta_{0}},$ (4.1)
$\displaystyle\frac{1-\theta_{k+1}}{\theta^{2}_{k+1}}$
$\displaystyle=\frac{1}{\theta^{2}_{k}},\quad\forall k=0,1,\dots$
$\displaystyle\theta_{k+1}$ $\displaystyle\leq\theta_{k},\quad\forall
k=0,1,\dots$

The next lemma gives a recurrence relation for the dual sequences generated by
Algorithm 2.

######  Lemma 3 ((fercoq2015accelerated, , see equation (52)))

Consider the linearly constrained optimization problem (1.2) with $1$-strongly
convex $f:\mathbb{R}^{n}\to\mathbb{R}$. Let $y^{k},z^{k}$ be the sequences
generated by Algorithm 2, $\theta_{0}=\frac{1}{M}$ and any solution
$\hat{y}\in\mathcal{Y}^{*}$. Then, it holds:

$\displaystyle\frac{1}{\theta^{2}_{k-1}}\mathbb{E}[\Psi(y^{k})-\hat{\Psi}]+\frac{1}{2\theta^{2}_{0}}\mathbb{E}[\|z^{k}-\hat{y}\|^{2}_{\mathbf{B}}]\leq\frac{1-\theta_{0}}{\theta^{2}_{0}}(\Psi(y^{0})-\hat{\Psi})+\frac{1}{2\theta^{2}_{0}}\|y^{0}-\hat{y}\|^{2}_{\mathbf{B}},$
(4.2)

with
${\mathbf{B}}=\textbf{Diag}(\|\mathbf{A}_{(1)}\|^{2},\|\mathbf{A}_{(2)}\|^{2},\dots,\|\mathbf{A}_{(M)}\|^{2})$.

For our proposed methods we are interested in showing convergence results for
the primal iterates $x^{k}$ given in (3.7) and not for the dual iterates
$y^{k}$. By showing that the APPROX and ARBK methods are equivalent, and using
recent theoretical results fercoq2015accelerated ; necoara2016 ;
nesterov2012efficiency ; 10178390  and Lemma 1, we obtain the following
convergence results for ARBK (and consequently also for BK) as a byproduct. We
denote (similarly to nesterov2012efficiency ) the size of the initial level
set in the standard Euclidean norm by

$R_{\alpha}(y^{0})=\text{max}_{y}\big{\\{}\text{min}_{\hat{y}\in\mathcal{Y}^{*}}\|y-\hat{y}\|_{\alpha}:\Psi(y)\leq\Psi(y^{0})\big{\\}},$
(4.3)

where

$\|y\|_{\alpha}^{2}\overset{\text{def}}{=}\sum\limits_{i=1}^{M}L_{i}^{\alpha}\|y_{(i)}\|^{2}_{2},\quad\quad\alpha\in[0,1]$
(4.4)

Moreover, we also denote

$\bar{L}_{\alpha}\overset{\text{def}}{=}\frac{1}{M}\sum\limits_{j=1}^{M}L_{j}^{\alpha}\
\text{ and }\ \bar{L}\overset{\text{def}}{=}\bar{L}_{1}$ (4.5)

###  4.1 Sublinear convergence results

In this section, we are going to show sublinear convergence results for both
Algorithm 1 and 2.

###### Theorem 4.1

Consider the linearly constrained optimization problem (1.2) with $1$-strongly
convex $f:\mathbb{R}^{n}\to\mathbb{R}$. Let $x^{k},d^{k}$ be the sequences
generated by BK i.e. Algorithm 1. At every iteration $k$ of the BK method, a
block index $i=i_{k}\in[M]$ is sampled with probability $p_{\alpha,M}(i)$ from
(3.6) and $\alpha\in[0,1]$. Then, it holds that

$\displaystyle\frac{1}{2}\mathbb{E}\big{[}\|x^{k}-\hat{x}\|^{2}_{2}\big{]}\leq\mathbb{E}\big{[}D_{f}^{d^{k}}(x^{k},\hat{x})\big{]}\leq\frac{2M\bar{L}_{\alpha}}{k+4}R^{2}_{1-\alpha}(y^{0}).$
(4.6)

###### Proof

Let $y^{k}$ be the sequence generated by iteration (3.5). It holds from
(nesterov2012efficiency, , Theorem 1):

$\mathbb{E}\bigg{[}\Psi(y^{k})-\hat{\Psi}\bigg{]}\leq\frac{2}{k+4}\cdot\bigg{[}\sum\limits_{j=1}^{M}L_{j}^{\alpha}\bigg{]}\cdot
R^{2}_{1-\alpha}(y^{0})\quad\forall k\geq 0.$

Furthermore, from Lemma 1, it holds
$D_{f}^{d^{k}}(x^{k},\hat{x})=\Psi(y^{k})-\hat{\Psi}$ and the first estimate
follows from the fact that $f$ is $1$-strongly convex, hence satisfies eq.
(2.5).

###### Remark 2

For the two extremal variants of the choice of probabilities in (3.6) we get
in (4.6) the following results:

  1. (a)

Uniform probabilities (i.e., $\alpha=0$):

$\displaystyle\mathbb{E}\big{[}\|x^{k}-\hat{x}\|^{2}_{2}\big{]}\leq\frac{4M}{k+4}R^{2}_{1}(y^{0}).$

  2. (b)

Probabilities proportional to $L_{i}$ (i.e., $\alpha=1$):

$\displaystyle\mathbb{E}\big{[}\|x^{k}-\hat{x}\|^{2}_{2}\big{]}\leq\frac{4M\bar{L}}{k+4}R^{2}_{0}(y^{0}).$

The next theorem shows convergence results for our accelerated method,
Algorithm 2.

###### Theorem 4.2

Consider the linearly constrained optimization problem (1.2) with $1$-strongly
convex $f:\mathbb{R}^{n}\to\mathbb{R}$. Let $x^{k},d^{k}$ be the sequences
generated by ARBK (Algorithm 2) with uniform probabilities,
$\theta_{0}=\tfrac{1}{M}$ and any $\hat{y}\in\mathcal{Y}^{*}$. Then, it holds:

$\displaystyle\frac{1}{2}\mathbb{E}\big{[}\|x^{k}-\hat{x}\|^{2}_{2}\big{]}\leq\mathbb{E}\big{[}D_{f}^{d^{k}}(x^{k},\hat{x})\big{]}\leq\frac{4M^{2}}{(k-1+2M)^{2}}C_{0},$
(4.7)

where

$C_{0}=\bigg{(}1-\frac{1}{M}\bigg{)}D_{f}^{d^{0}}(x^{0},\hat{x})+\frac{1}{2}\|y^{0}-\hat{y}\|^{2}_{\mathbf{B}}$

with $\mathbf{B}$ from Lemma 3.

###### Proof

From Lemma 3 and Lemma 1, it holds:

$\displaystyle\frac{1}{\theta^{2}_{k-1}}\mathbb{E}\big{[}D_{f}^{d^{k}}(x^{k},\hat{x})\big{]}$
$\displaystyle\leq\frac{1-\theta_{0}}{\theta^{2}_{0}}D_{f}^{d^{0}}(x^{0},\hat{x})+\frac{1}{2\theta^{2}_{0}}\|y^{0}-\hat{y}\|^{2}_{\mathbf{B}},$
$\displaystyle\mathbb{E}\big{[}D_{f}^{d^{k}}(x^{k},\hat{x})\big{]}$
$\displaystyle\leq\frac{\theta^{2}_{k-1}}{\theta^{2}_{0}}\bigg{(}(1-\theta_{0})D_{f}^{d^{0}}(x^{0},\hat{x})+\frac{1}{2}\|y^{0}-\hat{y}\|^{2}_{\mathbf{B}}\bigg{)},$
$\displaystyle\mathbb{E}\big{[}D_{f}^{d^{k}}(x^{k},\hat{x})\big{]}$
$\displaystyle\leq\frac{4M^{2}}{(k-1+2M)^{2}}\bigg{(}(1-\theta_{0})D_{f}^{d^{0}}(x^{0},\hat{x})+\frac{1}{2}\|y^{0}-\hat{y}\|^{2}_{\mathbf{B}}\bigg{)},$

where we used the first inequality in Lemma 2, i.e.,
$\frac{\theta^{2}_{k-1}}{\theta^{2}_{0}}\leq\frac{4}{(k-1+2M)^{2}}\cdot\frac{1}{\frac{1}{M^{2}}}=\frac{4M^{2}}{(k-1+2M)^{2}}$.
The first inequality in Theorem 4.2 follows from Eq (2.5).

Theorem 4.2 shows that the iterates $x^{k}$ of Algorithm 2 (ARBK) converge at
the rate $\mathcal{O}(1/k^{2})$, thus accelerating its standard counterpart,
Algorithm 1 (BK), which has $\mathcal{O}(1/k)$ rate of convergence (cf.
Theorem 4.1). To the best of our knowledge, accelerated Kaczmarz variants have
not yet been proposed for problem (1.2) and the convergence guarantees for
ARBK algorithm from Theorem 4.2 are new. Moreover, the convergence rates from
both theorems (Theorem 4.1 and Theorem 4.2) depend on the number of blocks,
$M$, proving that it is beneficial to consider blocks of the matrix
$\mathbf{A}$ in our previous algorithms. Under some regularity assumptions on
the objective function $f$, see Section 2, we can even show linear convergence
for the BK and ARBK methods, respectively.

###  4.2 Linear convergence results

Consider the feasible, convex and linearly constrained optimization problem
(1.2). To obtain convergence rates for the solution algorithm, we will
estimate the Bregman distance of the iterates to the solution $\hat{x}$ by
error bounds of the form
$D_{f}^{x^{*}}(x,\hat{x})\leq\theta\cdot\|\mathbf{A}x-b\|^{2}_{2}$. We know
that such error bounds always hold if $f$ has a Lipschitz-continuous gradient
necoara2022linear . But they also hold under weaker conditions. For example,
the following inequality holds for problem (1.2) with objective function
defined by (2.4):

$D_{f}^{x^{*}}(x,\hat{x})\leq\theta(\hat{x})\cdot\|\mathbf{A}x-b\|^{2}_{2},$

we refer the reader to schopfer2022extended  for more details on the constant
$\theta(\hat{x})$. To clarify the assumptions under which such error bounds
hold for more general objective functions, notions such as calmness, linear
regularity, and linear growth schopfer2022extended  have been defined in
Section 2.2. The next result shows that the dual function satisfies some PL
inequality.

######  Theorem 4.3 ((schopfer2022extended, , Theorem 3.9))

Consider the linearly constrained optimization problem (1.2) with strongly
convex objective function $f:\mathbb{R}^{n}\to\mathbb{R}$. If Assumptions 2
and 4 hold at the unique solution $\hat{x}$ of (1.2) and if the collection
$\\{\partial f(\hat{x}),\mathcal{R}(\mathbf{A}^{\top})\\}$ satisfies
Assumption 3, then $\Psi$ satisfies the PL inequality (2.9), i.e., there
exists $\gamma(\hat{x})>0$ such that for all $x\in\mathbb{R}^{n}$ and
$d\in\partial f(x)\cap\mathcal{R}(\mathbf{A}^{\top})$ we have the global error
bound :

$D_{f}^{d}(x,\hat{x})=\Psi(y)-\hat{\Psi}\leq\gamma(\hat{x})\cdot\|\nabla\Psi(y)\|^{2}_{2}=\gamma(\hat{x})\cdot\|\mathbf{A}x-b\|_{2}^{2}.$
(4.8)

The following Theorem gives a convergence result for the iterates generated by
Algorithm 1.

###### Theorem 4.4

Consider the linearly constrained optimization problem (1.2) with $1$-strongly
convex $f:\mathbb{R}^{n}\to\mathbb{R}$. Assume that $f$ is calm at the unique
solution $\hat{x}$ of (1.2) and that $\partial f$ grows at most linearly (cf.
Def. 4, and Def. 2, respectively), and the collection $\\{\partial
f(\hat{x}),\mathcal{R}(\mathbf{A}^{\top})\\}$ is linearly regular (cf. Def.
3). Let $x^{k}$ be the sequence generated by BK (Algorithm 1) with probability

$p_{\alpha,M}(i)=p_{i}=L_{i}^{\alpha}\cdot\bigg{[}\sum\limits_{j=1}^{M}L_{j}^{\alpha}\bigg{]}^{-1},\;\;\text{where}\;\;\alpha\in[0,1].$

Then, it holds

$\displaystyle\mathbb{E}[D_{f}^{d^{k+1}}(x^{k+1},\hat{x})]\leq\bigg{(}1-\frac{1}{2M\gamma(\hat{x})\bar{L}_{\alpha}\bar{L}^{1-\alpha}}\bigg{)}\mathbb{E}[D_{f}^{d^{k}}(x^{k},\hat{x})],$
(4.9)

and

$\displaystyle\mathbb{E}[\|x^{k}-\hat{x}\|_{2}^{2}]\leq
2\bigg{(}1-\frac{1}{2M\gamma(\hat{x})\bar{L}_{\alpha}\bar{L}^{1-\alpha}}\bigg{)}^{k}D_{f}^{d^{0}}(x^{0},\hat{x})$
(4.10)

###### Proof

Let $y^{k}$ be the sequence generated by iteration (3.5). From the fact that
the gradient of function $\Psi$ is coordinatewise Lipschitz continuous with
constants $L_{i}=\|\mathbf{A}_{(i)}\|^{2}$, we have:

$\displaystyle\Psi(y^{k+1})\leq\Psi(y^{k})-\frac{1}{2\|\mathbf{A}_{(i)}\|^{2}}\|\nabla_{(i)}\Psi(y^{k})\|^{2},$

and thanks to Theorem 4.3 we have the following relation

$\Psi(y)-\min\Psi\leq\gamma(\hat{x})\cdot\|\nabla\Psi(y)\|^{2}_{2}.$

Then with
$p_{M}(i)=p_{i}=L_{i}^{\alpha}\cdot\big{[}\sum\limits_{j=1}^{M}L_{j}^{\alpha}\big{]}^{-1},\;\;\text{where}\;\;\alpha\in[0,1]$,
we have:

$\displaystyle\mathbb{E}_{i_{k}}\bigg{[}\Psi(y^{k+1})-\hat{\Psi}\bigg{]}$
$\displaystyle\leq\Psi(y^{k})-\hat{\Psi}-\sum\limits_{i=1}^{M}\frac{p_{i}}{2L_{i}}\|\nabla_{(i)}\Psi(y^{k})\|^{2}_{2}$
$\displaystyle=\Psi(y^{k})-\hat{\Psi}-\frac{1}{2\sum\limits_{j=1}^{M}L_{j}^{\alpha}}\sum\limits_{i=1}^{M}\frac{1}{L_{i}^{1-\alpha}}\|\nabla_{(i)}\Psi(y^{k})\|^{2}_{2}$
$\displaystyle\leq\Psi(y^{k})-\hat{\Psi}-\frac{1}{2M\bar{L}_{\alpha}\bar{L}^{1-\alpha}}\cdot\|\nabla\Psi(y^{k})\|^{2}_{2},$
$\displaystyle\leq\bigg{(}1-\frac{1}{2M\gamma(\hat{x})\bar{L}_{\alpha}\bar{L}^{1-\alpha}}\bigg{)}(\Psi(y^{k})-\hat{\Psi}).$

Then, using the rule of total expectation gives us:

$\displaystyle\mathbb{E}\big{[}\Psi(y^{k+1})-\hat{\Psi}\big{]}\leq\bigg{(}1-\frac{1}{2M\gamma(\hat{x})\bar{L}_{\alpha}\bar{L}^{1-\alpha}}\bigg{)}\cdot\mathbb{E}\big{[}\Psi(y^{k})-\hat{\Psi}\big{]}.$

Furthermore, from Lemma 1, it holds
$D_{f}^{d^{k}}(x^{k},\hat{x})=\Psi(y^{k})-\hat{\Psi}$, which gives Equation
(4.9). Unrolling Eq (4.9) give us :

$\displaystyle\mathbb{E}[D_{f}^{d^{k}}(x^{k},\hat{x})]\leq\bigg{(}1-\frac{1}{2M\gamma(\hat{x})\bar{L}_{\alpha}\bar{L}^{1-\alpha}}\bigg{)}^{k}\cdot
D_{f}^{d^{0}}(x^{0},\hat{x}).$

The second estimate follows from the fact that $f$ is $1-$strongly convex,
hence satisfies $\frac{1}{2}\|x^{k}-\hat{x}\|^{2}_{2}\leq
D^{d^{k}}_{f}(x^{k},\hat{x})$.

###### Remark 3

The linear convergence rate from the previous theorem depends on the number of
blocks $M$ and shows that considering blocks is beneficial for BK (Algorithm
1). Indeed, Theorem 4.4 tells us that taking $M=m$ gives us the slowest
convergence rate and thus one should consider updating the blocks having a
size greater than 1. Note that Theorem 4.4 generalizes the convergence results
from P15 ; LS19 . In fact, to get the linear convergence rate in P15  the
authors needed the objective function to be smooth which is not the case in
our theorem. The linear convergence result in LS19  can be obtained from our
Theorem 4.4 by setting $\alpha=1,M=m$. In this case
$M\bar{L}=\|\mathbf{A}\|^{2}_{F}$, which corresponds to the result in (LS19, ,
Theorem 3.2).

A well-known drawback of accelerated first-order methods like Algorithm 2
types is their oscillating behavior that slows them down (see e.g Figure 2 and
3 below). As a first path towards solving these issues, restart schemes have
been proposed in the literature, which have been demonstrated to improve
convergence in practice by removing the oscillatory behavior necoara2019linear
. Practically speaking in a restart scheme, Algorithm 2 is stopped when a
certain criterion is met and then restarted using the last value provided by
the method as the new initial condition. In the next section, we present a
restart scheme for ARBK (Algorithm 2) that exhibits linear convergence for our
problem (1.2).

###  4.3 Restarted ARBK

It has been shown that under the quadratic growth property accelerated
gradient methods can have improved complexity bounds in combination with an
appropriate restarting procedure, see e.g. fercoq2020restarting ;
necoara2019linear ; necoara2022linear . Our goal is to design a restarted ARBK
akin to the restart strategies in fercoq2020restarting ; necoara2022linear
with similar properties based on the results in the previous sections. Having
defined a set of integers $\\{K_{0},K_{1},...\\}$, at which frequencies one
wishes to restart the method, we can write the restarted version of ARBK. In
the following version of ARBK in Algorithm 3 we restart with these
frequencies, but only update the dual variable if we have decreased the
objective during the last restart period.

Algorithm 3 Restarted-ARBK (RARBK)

1: Choose $y^{0}\in\mathbb{R}^{m}$ and set $\tilde{y}^{0}=y^{0}$.

2: choose restart periods $\\{K_{0},\dots,K_{r},\dots\\}\subseteq\mathbb{N}$.

3: initialize $r=0$

4: repeat

5: $\bar{y}^{r+1}=\text{ARBK}(f,\tilde{y}^{r},K_{r})$

6: update $\tilde{y}^{r+1}=\begin{cases}\bar{y}^{r+1},&\text{if
}\Psi(\bar{y}^{r+1})\leq\Psi(\tilde{y}^{r})\\\ \\\ \tilde{y}^{r},&\text{if
}\Psi(\bar{y}^{r+1})>\Psi(\tilde{y}^{r})\end{cases}$

7: increment $r=r+1$

8: until a stopping criterion is satisfied

9: return $\tilde{x}^{r+1}=\nabla f^{*}(\mathbf{A}^{\top}\tilde{y}^{r+1})$.

The next lemma gives an estimate of the progress we make in the primal
variable in one restart period.

######  Lemma 4 (Conditional restarting at $y^{k}$)

Let $(y^{k},x^{k})$ be the iterates of ARBK (Algorithm 2). Let

$\displaystyle\tilde{y}$ $\displaystyle=\begin{cases}y^{k},&\text{if
}\Psi(y^{k})\leq\Psi(y^{0})\\\ y^{0},&\text{if
}\Psi(y^{k})>\Psi(y^{0})\end{cases}$ (4.11)

and set $\tilde{d}=\mathbf{A}^{\top}\tilde{y}$, $\tilde{x}=\nabla
f^{*}(\tilde{d})$. Then, we have

$\mathbb{E}[D^{\tilde{d}}_{f}(\tilde{x},\hat{x})]\leq\theta^{2}_{k-1}\Bigg{(}\frac{1-\theta_{0}}{\theta^{2}_{0}}+\frac{L_{\text{max}}}{\gamma(\hat{x})\theta^{2}_{0}}\Bigg{)}\cdot
D^{d^{0}}_{f}(x^{0},\hat{x}),$ (4.12)

with
$L_{\text{max}}=\max\limits_{i\in[M]}\|\mathbf{A}_{(i)}\|^{2},\theta_{0}=\frac{1}{M}$.
Moreover, given $\zeta<1$, if

$k\geq\frac{2}{\theta_{0}}\Bigg{(}\sqrt{\frac{L_{\text{max}}+\gamma(\hat{x})}{\zeta\gamma(\hat{x})}}-1\Bigg{)}+1,$
(4.13)

then $\mathbb{E}[D^{\tilde{d}}_{f}(\tilde{x},\hat{x})]\leq\zeta\cdot
D^{d^{0}}_{f}(x^{0},\hat{x})$.

###### Proof

Since $\Psi$ satisfies now the PL inequality (2.9) it also satisfies the
quadratic growth condition, i.e.,

$\Psi(y^{k})-\hat{\Psi}\geq\frac{\gamma(\hat{x})}{2}\text{dist}(y^{k},\mathcal{Y}^{*})^{2}.$
(4.14)

From Lemma 3, the following holds for the iterates of ARBK:

$\displaystyle\mathbb{E}[\Psi(y^{k})-\hat{\Psi}]$
$\displaystyle\leq\theta^{2}_{k-1}\Bigg{(}\frac{1-\theta_{0}}{\theta^{2}_{0}}(\Psi(y^{0})-\hat{\Psi})+\frac{1}{2\theta^{2}_{0}}\text{dist}_{\mathbf{B}}(y^{0},\mathcal{Y}*)^{2}\Bigg{)}$
$\displaystyle\leq\theta^{2}_{k-1}\Bigg{(}\frac{1-\theta_{0}}{\theta^{2}_{0}}(\Psi(y^{0})-\hat{\Psi})+\frac{L_{\text{max}}}{2\theta^{2}_{0}}\text{dist}(y^{0},\mathcal{Y}*)^{2}\Bigg{)}$
$\displaystyle\leq\theta^{2}_{k-1}\Bigg{(}\frac{1-\theta_{0}}{\theta^{2}_{0}}+\frac{L_{\text{max}}}{\gamma(\hat{x})\theta^{2}_{0}}\Bigg{)}(\Psi(y^{0})-\hat{\Psi}),$

where $\text{dist}_{\mathbf{B}}(y^{0},\mathcal{Y}*)$ denotes the distance of
$y^{0}$ to the set $\mathcal{Y}*$ with respect to the norm
$\|\cdot\|_{\mathbf{B}}$ with
${\mathbf{B}}=\textbf{Diag}(\|\mathbf{A}_{(1)}\|^{2},\|\mathbf{A}_{(2)}\|^{2},\dots,\|\mathbf{A}_{(M)}\|^{2})$
and $L_{\text{max}}=\max_{i\in[M]}\|\mathbf{A}_{(i)}\|^{2}$. Furthermore from
Lemma 1, it holds $D_{f}^{d^{k}}(x^{k},\hat{x})=\Psi(y^{k})-\hat{\Psi}$ and
from condition (4.11), we have $D^{\tilde{d}}_{f}(\tilde{x},\hat{x})\leq
D^{d}_{f}(x^{k},\hat{x})$ which gives Equation (4.12). Condition (4.13) is
equivalent to:

$\frac{4}{(k-1+2M)^{2}}\Bigg{(}\frac{1-\theta_{0}}{\theta^{2}_{0}}+\frac{L_{\text{max}}}{\gamma(\hat{x})\theta^{2}_{0}}\Bigg{)}\leq\zeta.$

###### Corollary 1

Denote
$K(\zeta)=\bigg{\lceil}{2M\bigg{(}\sqrt{\frac{L_{\text{max}}+\gamma(\hat{x})}{\zeta\gamma(\hat{x})}}-1\bigg{)}+1\bigg{\rceil}}$.
If the restart periods $\\{K_{0},\dots,K_{r},\dots\\}$ are all equal to
$K(\zeta)$, then the iterates of Algorithm 3 satisfy:

$\mathbb{E}[D^{\tilde{d}^{r}}_{f}(\tilde{x}^{r},\hat{x})]\leq\zeta^{r}\cdot
D^{d^{0}}_{f}(x^{0},\hat{x}),$

where the iterate $\tilde{x}^{r}$ is obtained after restarting $r$ times ARBK
performing each time $K(\alpha)$ iterations.

###### Proof

By the definition of $\tilde{x}^{r}$, we know that for all $r$,
$D^{\tilde{d}^{r+1}}_{f}(\tilde{x}^{r+1},\hat{x})\leq\zeta\cdot
D^{\tilde{d}^{r}}_{f}(\tilde{x}^{r},\hat{x})$. Applying recursively Lemma 4
gives us the desired linear convergence.

From relation (2.5) and Corollary 1 we can also derive a linear rate in terms
of the expected quadratic distance of this primal sequence $\tilde{x}^{r}$ to
the optimal solution $\hat{x}$

$\mathbb{E}[\|\tilde{x}^{r}-\hat{x}\|_{2}^{2}]\leq 2\zeta^{r}\cdot
D^{d^{0}}_{f}(x^{0},\hat{x})$

The iterate $\tilde{x}^{r}$ is obtained after running $r$ times Algorithm 3 of
$K(\zeta)$ iterations, which suggests $\zeta=e^{-2}$, or equivalently a fixed
restart period (see also necoara2019linear ):

$K^{*}=\bigg{\lceil}{2eM\bigg{(}\sqrt{\frac{L_{\text{max}}+\gamma(\hat{x})}{\gamma(\hat{x})}}-1\bigg{)}+1\bigg{\rceil}}$
(4.15)

Even though the quantity $\gamma(\hat{x})$ is not known, Corollary 1 shows
that expected linear convergence can be obtained for restarted ARBK with a
long enough period $K$. In fact, we can easily see that those restarting
periods depend on the block size $M$ and hence it is beneficial to work with
blocks of matrix less than $m$. Further, we present a restarting variation of
ARBK as proposed in fercoq2020restarting , which does not require explicit
knowledge of $\gamma(\hat{x})$. Finally, we fix the length of the first epoch,
$K_{0}$, to

$K_{0}=\bigg{\lceil}{2eM\bigg{(}\sqrt{\frac{L_{\text{max}}+\bar{\gamma}(\hat{x})}{\bar{\gamma}(\hat{x})}}-1\bigg{)}+1\bigg{\rceil}}$

where $\bar{\gamma}(\hat{x})$ is an estimate of the unknown constant
$\gamma(\hat{x})$.

###### Theorem 4.5

Let the sequence $\\{K_{j}\\}_{j\geq 0}\subset\mathbb{N}$ satisfy:

  1. (a)

$K_{2^{j}-1}=2^{j}K_{0}$ for all $j\in\mathbb{N}$

  2. (b)

$|\\{0\leq r<2^{p}-1|\,\,K_{r}=2^{j}K_{0}\\}|=2^{p-1-j}$ for all $j\in[p]$.

Then, after $p$ epochs of Algorithm 3 we have the following linear rate in
expectation:

$\mathbb{E}[D^{\tilde{d}^{k}}_{f}(\tilde{x}^{k},\hat{x})]\leq\bigg{(}e^{-\frac{4}{(p+2)2^{i}K_{0}}}\bigg{)}^{k}\cdot
D^{d^{0}}_{f}(x^{0},\hat{x}),$

with $i=\lceil\max(0,\text{log}_{2}(K^{*}/K_{0}))\rceil$ and
$k=\sum\limits_{j=0}^{2^{p}-1}K_{j}$ is the total number of iterations.

###### Proof

Let us define
$\Delta_{r}=\mathbb{E}[D^{\tilde{d}^{r}}_{f}(\tilde{x}^{r},\hat{x})]$ and

$c_{i}(p)=|\\{l<2^{p}-1|\,\,K_{l}\geq
2^{i}K_{0}\\}|+1=1+\sum\limits_{k=i}^{p-1}2^{p-1-k}=2^{p-i}$

Note that $c_{i}(p)$ represents the number of restarts such that $K_{l}\geq
K^{*}$, i.e. the number of restarts for which Lemma 4 applies. When
$K_{l}<K^{*}$, we invoke the guaranteed Bregman distance decrease enforced in
Algorithm 3. Then it follows from Lemma 4 that

$\Delta_{2^{p}-1}\leq e^{-2c_{i}(p)}\cdot\Delta_{0}.$ (4.16)

Moreover, we have :

$\displaystyle k=\sum\limits_{j=0}^{2^{p}-1}K_{j}$
$\displaystyle=\sum\limits_{j=0}^{p-1}|\\{0\leq
r<2^{p}-1|\,\,K_{r}=2^{j}K_{0}\\}|\cdot 2^{j}K_{0}+K_{2^{p}-1}$
$\displaystyle=\sum\limits_{j=0}^{p-1}2^{p-1-j}\cdot 2^{j}K_{0}+K_{2^{p}-1}$
$\displaystyle=(p+2)2^{p-1}K_{0}$

and the relation (4.16) implies:

$\Delta_{k}\leq\bigg{(}e^{-2\frac{c_{i}(p)}{\sum_{j=0}^{2^{p}-1}K_{j}}}\bigg{)}^{k}\Delta_{0}=\bigg{(}e^{-\frac{2^{p-i+1}}{(p+2)2^{p-1}K_{0}}}\bigg{)}^{k}\Delta_{0}$

which confirms the above result.

###### Remark 4

  1. (a)

The linear convergence rate stated in Theorem 4.4 clearly implies the
following estimate on the total number of iterations required by BK to obtain
an $\varepsilon$-suboptimal solution in expectation:

$\mathcal{O}\bigg{(}2M\gamma(\hat{x})\bar{L}_{\alpha}\bar{L}^{1-\alpha}\text{log}\bigg{(}\frac{1}{\varepsilon}\bigg{)}\bigg{)}.$

  2. (b)

It follows from Theorem 4.5 that an upper bound on the total number of
iterations performed by the Restarted-ARBK scheme to attain an
$\varepsilon$-suboptimal solution in expectation is given by:

$\mathcal{O}\bigg{(}M\sqrt{\frac{L_{\text{max}}}{\gamma(\hat{x})}}\text{log}\bigg{(}\frac{1}{\varepsilon}\bigg{)}\text{log}_{2}\bigg{(}\text{log}\bigg{(}\frac{f(\hat{x})}{\varepsilon}\bigg{)}\sqrt{\frac{\bar{\gamma}(\hat{x})}{\gamma(\hat{x})}}\bigg{)}\bigg{)}.$

  3. (c)

If we assume that $f(\hat{x})\leq 1$ and $\gamma(\hat{x})$ known, then the
Restarted-ARBK is better than BK when
$\text{log}_{2}(\text{log}(1/\varepsilon))\leq
2\gamma\sqrt{L_{\text{max}}\gamma}$. For instance, let’s suppose that the rows
of the matrix are normalized and set $M=m$, then $L_{\text{max}}=1$. If
$\gamma=10^{1/6}$, then Restarted-ARBK has a better worst-case complexity than
ARBK for all accuracies $\varepsilon>10^{-6}$, $\gamma=10^{1/3}$ for
$\varepsilon>10^{-35}$ and $\gamma=10^{2/3}$ for $\varepsilon>10^{-455391}$.

##  5 Experiments

In this section, we study the computational behaviour of the proposed
algorithms, BK (Algorithm 1) ARBK (Algorithm 2), and RARBK (Algorithm 3), for
the application of finding sparse solutions of linear systems. The first part
shows a visual illustration on a simple $\mathbb{R}^{2}$ example.

###  5.1 Visualizing the acceleration mechanism

This section shows the graphical illustration of our acceleration mechanism.
Our goal is to shed more light on how the proposed algorithm works in
practice. For simplicity, we illustrate this by comparing BK and ARBK.

Figure 1: Graphical interpretation of the randomized Bregman Kaczmarz method
(BK) and the accelerated randomized Bregman Kaczmarz method (ARBK) using our
scheme in a simple example with only two hyperplanes $H_{i}=\\{x:\langle
a_{i},x\rangle=b_{i}\\}$, with $i=1,2$, and a unique solution $\hat{x}$.

In Fig. 1 we present a simple $\mathbb{R}^{2}$ illustration of the difference
between the workings of BK and ARBK with $f(x)=\frac{1}{2}\|x\|^{2}_{2}$. Our
goal is to show graphically how the addition of interpolation leads to
acceleration. In the example of Fig. 1, the performance of ARBK is similar to
the performance of BK until iterate $x_{3}$. After this point, the
interpolation parameter becomes more effective and the ARBK method
accelerates. This behavior appears also in our experiments in the next section
where we work with matrices with many rows. There we can notice that the
interpolation parameter seems to become more effective after the first $m+1$
iterations, where $m$ is the number of hyperplanes.

###  5.2 Numerical Experiments

We present several experiments to demonstrate the effectiveness of Algorithm 2
and 3 under various conditions111The code that produces the figures in this
paper is available at https://github.com/tondji/arbk.. The simulations were
performed in Python on an Intel Core i7 computer with 16GB RAM. For all the
experiments we consider $f(x)=\lambda\cdot\|x\|_{1}+\frac{1}{2}\|x\|^{2}_{2}$,
where $\lambda$ is the sparsity parameter.

####  5.2.1 Synthetic experiments

In the first part, synthetic data for the experiments is generated as follows:
all elements of the data matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$ are
chosen independent and identically distributed from the standard normal
distribution $\mathcal{N}(0,1)$. We constructed overdetermined, square, and
underdetermined linear systems. To construct sparse solutions
$\hat{x}\in\mathbb{R}^{n}$, we generate a random vector $y$ from the standard
normal distribution $\mathcal{N}(0,1)$ and we set
$\hat{x}=S_{\lambda}(\mathbf{A}^{\top}y)$, which comes from Eq (3.2) and the
corresponding right hand side is $b=\mathbf{A}\hat{x}\in\mathbb{R}^{m}$. For
each experiment, we run independent trials each starting with the initial
iterate $x_{0}=0$ and the number of blocks $M$ will always be specified. We
measure performance by plotting the relative residual error
$\|\mathbf{A}x-b\|_{2}/\|b\|_{2}$ and the relative error
$\|x_{k}-\hat{x}\|_{2}/\|\hat{x}\|_{2}$ against the number of epochs (number
of passing through all rows of $\mathbf{A}$). All the experiments are run for
a maximum number of epochs in case the tolerance tol is not reached. During
every experiment, the $i-$th row block of the matrix $\mathbf{A}$ is sampled
with probabilities according to Eq (3.6). We compared the four following
methods for $200*\max(m,n)$ epochs:

  1. (a)

BK: The randomized Bregman-Kaczmarz method, i.e., Algorithm 1, with
$\lambda=15$ and $M\in\\{m/4,m/2\\}$.

  2. (b)

NRBK: The Nesterov randomized Bregman-Kaczmarz method (nesterov2012efficiency,
, Method ACDM in Section 5), with $\lambda=15$ and $M\in\\{m/4,m/2\\}$.

  3. (c)

ARBK: The accelerated randomized Bregman-Kaczmarz method, i.e., Algorithm 2,
with $\lambda=15$ and $M\in\\{m/4,m/2\\}$.

  4. (d)

RARBK: The restarted accelerated randomized Bregman-Kaczmarz method, i.e.,
Algorithm 3, with $\lambda=15$, $M\in\\{m/4,m/2\\}$ with restarting periods
$K^{*}\in\\{165M,200M\\}$.

In Figure 2 and Figure 5, we report relative errors and residuals of the
Bregman Kaczmarz method (BK), the accelerated Bregman Kaczmarz method (ARBK),
it restarted version (RARBK) and the Nesterov acceleration (NRBK) with block
size equal to $M=m/4$ and $M=m/2$ for undetermined and squared matrices,
respectively. For the RARBK method, we chose as restart period $K^{*}=165M$ in
Figure 2 and $K^{*}=200M$ in Figure 5. We observe that the ARBK method
accelerates the BK method and the restarted version of ARBK, namely RARBK,
gives us the fastest convergence, thus accelerating the ARBK method. From
Table 1 and Table 4, we also see that the RARBK method needs the least
computational time to achieve a relative error of $10^{-6}$, followed by the
ARBK method (in fact, RARBK is about 4 times faster than BK).

In Figure 3 and Figure 4, we analyze the effects of the number of blocks $M$
in the relative error and residual for Algorithm 2 and Algorithm 3. We used an
underdetermined and consistent system, with $M\in\\{m,m/2,m/4,m/10,m/100\\}$ ,
$\lambda=15$ and we choose as restarting period $K^{*}=165M$. As the number of
blocks $M$ decreases, i.e., increasing the number of rows used in one
iteration, we see a corresponding decrease in the number of iterations needed
to reach a certain accuracy, however at some point decreasing $M$ does not
improve the method. Moreover, we observed that the RARBK compared to ARBK
always needs fewer iterations to reach the same accuracy. These observations
are in line with Table 2 and Table 3 in terms of computational times.

(a) Relative Residual

(b) Relative Error

Figure 2: A comparison of Bregman Kaczmarz (BK-blue), accelerated randomized Bregman Kaczmarz (ARBK-red), Nesterov Acceleration scheme (NRBK-green) and restarted version of ARBK (RARBK-black), $m=500,n=784$, sparsity $s=408$, $M=m/4$, $\lambda=15$, $\kappa(A)=8.98$, tol=$10^{-6}$, $K^{*}=165M$. | CPU time (s)
---|---
BK-125 | 46.58
ARBK-125 | 24.65
NRBK-125 | 53.16
RARBK-125 | 11.86
Table 1: CPU time until the relative error or residual falls below tol for
Figure 2.

(a) Relative Residual

(b) Relative Error

Figure 3: A comparison of the accelerated randomized Bregman Kaczmarz (ARBK) for different numbers of blocks $M$, with $m=500,n=784$, sparsity $s=408$, $\lambda=15$, $\kappa(A)=8.98$, tol=$10^{-6}$, $K^{*}=165M$. | CPU time (s)
---|---
ARBK-500 | 65.28
ARBK-250 | 51.82
ARBK-125 | 31.08
ARBK-50 | 11.30
ARBK-5 | 2.01
Table 2: CPU time until the relative error or residual falls below tol for
Figure 3.

(a) Relative Residual

(b) Relative Error

Figure 4: A comparison of restart accelerated randomized Bregman Kaczmarz (RARBK) for different numbers of blocks $M$, with $m=500,n=784$, sparsity $s=408$, $\lambda=15$, $\kappa(A)=8.98$, tol=$10^{-6}$, $K^{*}=165M$. | CPU time (s)
---|---
RARBK-500 | 41.94
RARBK-250 | 18.25
RARBK-125 | 7.87
RARBK-50 | 2.93
RARBK-5 | 0.53
Table 3: CPU time until the relative error or residual falls below tol for
Figure 4.

(a) Relative Residual

(b) Relative Error

Figure 5: A comparison of Bregman Kaczmarz (BK-blue), accelerated randomized Bregman Kaczmarz (ARBK-red), the Nesterov Acceleration scheme (NRBK-green) and the restarted version of ARBK (RARBK-black), $m=700,n=700$, sparsity $s=404$, $M=m/2$, $\lambda=15$, $\kappa(A)=1150.06$, tol=$10^{-6}$, $K^{*}=200M$. | CPU time (s)
---|---
BK-350 | 56.70
ARBK-350 | 35.21
NRBK-350 | 62.13
RARBK-350 | 27.38
Table 4: CPU time until the relative error or residual falls below tol for
Figure 5.

####  5.2.2 Computerized tomography

In this section, we consider the application of computerized tomography (CT),
where the Radon transform implementation from the Python library skimage was
used to build a system matrix for a parallel beam CT for a phantom of size
$N\times N$ with $N=50$ and with $60$ equispaced angles. Hence, the system
matrix $\mathbf{A}$ has size $3000\times 2500$, i.e., $m=3000$ and $n=2500$.
We interpret the projection for each angle as one block $\mathbf{A}_{(i)}$,
i.e., we have $M=60$ and each block has the size $m_{i}=50$. The ground truth
solution $\hat{x}$ is fairly sparse (see below) and we generated the exact
right-hand side simply as $b=\mathbf{A}\hat{x}$. We used different methods for
reconstruction, each run for $10$ epochs using the function
$f(x)=\lambda\cdot\|x\|_{1}+\frac{1}{2}\|x\|^{2}_{2}$:

  1. (a)

BK: Randomized Bregman-Kaczmarz method, i.e., Algorithm 1, with $\lambda=30$
and $M=60$.

  2. (b)

NRBK: Nesterov randomized Bregman-Kaczmarz method (nesterov2012efficiency, ,
Method ACDM in Section 5), with $\lambda=30$ and $M=60$.

  3. (c)

ARBK: Accelerated randomized Bregman-Kaczmarz method, i.e., Algorithm 2, with
$\lambda=30$ and $M=60$.

  4. (d)

RARBK: Restarted accelerated randomized Bregman-Kaczmarz method, i.e.,
Algorithm 3, with $\lambda=30$, $M=60$ and restarting periods $K^{*}=165M$.

In this part, BK algorithm represents our baseline comparison. It can be
applied by just choosing $\lambda$ by trial and error to adapt the sparsity of
the outcome. Figure 6 shows the ground truth data and the reconstructions by
BK, NRBK, ARBK and RARBK. As can be seen, all methods produce a clear
background due to the sparsity enforcing $f$. Both ARBK and RARBK methods lead
to better reconstruction. Table 5 shows the CPU time for the four methods in
order to reach a relative error or relative residual below tol. We notice
again that RARBK is the fastest method (about 2.5 times faster than BK).

Figure 6: Ground truth solution and reconstruction by the different methods in Figure 7. | CPU time (s)
---|---
BK-60 | 202.17
ARBK-60 | 165.35
NRBK-60 | 203.63
RARBK-60 | 84.84
Table 5: CPU time until the relative error or residual falls below tol for
Figure 7.

In Figure 7, we show the relative residuals and the relative reconstruction
errors for all methods. While the BK and NRBK methods bring the residual down
to some value, the accelerated method (ARBK) and its restarted version (RARBK
) are able to go way below that and it seems that they even go down further.
It is important to note that the choice of the restarting periods $K^{*}$ is
crucial for making the RARBK work properly and accelerating the ARBK method.
The parameter $K^{*}$ is difficult to estimate (see Eq (4.15)) as it is
related to the quadratic growth of the dual objective $\gamma(\hat{x})$.
Nevertheless a restarting strategy with long enough periods $K^{*}$ leads to
good results in most experiments.

(a) Relative Residual

(b) Relative Error

Figure 7: A comparison of Bregman Kaczmarz (BK-blue), accelerated randomized
Bregman Kaczmarz (ARBK-red), Nesterov Acceleration scheme (NRBK-green) and
restarted version of ARBK (RARBK-black), $m=3000,n=2500$, sparsity $s=912$,
$M=60$, $\lambda=30$, $\kappa(A)=5411.08$, tol=$10^{-5}$, $K=165M$.

##  6 Conclusions

In this work, we have proposed block (accelerated) randomized Bregman-Kaczmarz
methods for solving convex linearly constrained optimization problems. Our
methods are of Kaczmarz-type, i.e., they use a block of constraints in each
iteration. We have used a recently proposed acceleration variant for the
randomized coordinate descent and transferred it to the primal space. A
theoretical analysis of the convergence of the proposed nonaccelerated and
accelerated algorithms was given. Numerical experiments have shown that the
proposed methods are more efficient and faster than the existing methods for
solving the same problem.

## References

  * (1) Jian-Feng Cai, Stanley Osher, and Zuowei Shen.  Linearized Bregman iterations for compressed sensing.  Mathematics of Computation, 78(267):1515–1536, 2009.
  * (2) Kui Du.  Tight upper bounds for the convergence of the randomized extended Kaczmarz and Gauss–Seidel algorithms.  Numerical Linear Algebra with Applications, 26(3):e2233, 2019.
  * (3) Kui Du, Wu-Tao Si, and Xiao-Hui Sun.  Randomized extended average block Kaczmarz for solving least squares.  SIAM Journal on Scientific Computing, 42(6):A3541–A3559, 2020.
  * (4) Olivier Fercoq and Zheng Qu.  Restarting the accelerated coordinate descent method with a rough strong convexity estimate.  Computational Optimization and Applications, 75(1):63–91, 2020\.
  * (5) Olivier Fercoq and Peter Richtárik.  Accelerated, parallel, and proximal coordinate descent.  SIAM Journal on Optimization, 25(4):1997–2023, 2015.
  * (6) Robert M. Gower, Denali Molitor, Jacob Moorman, and Deanna Needell.  On adaptive sketch-and-project for solving linear systems.  SIAM Journal on Matrix Analysis and Applications, 42(2):954–989, 2021.
  * (7) Robert M Gower and Peter Richtárik.  Randomized iterative methods for linear systems.  SIAM Journal on Matrix Analysis and Applications, 36(4):1660–1690, 2015.
  * (8) D. A. Lorenz, F. Schöpfer, and S. Wenger.  The linearized Bregman method via split feasibility problems: Analysis and generalizations.  SIAM J. Imaging Sciences, 7(2):1237–1262, 2014.
  * (9) Dirk A Lorenz, Stephan Wenger, Frank Schöpfer, and Marcus Magnor.  A sparse Kaczmarz solver and a linearized Bregman method for online compressed sensing.  In 2014 IEEE international conference on image processing (ICIP), pages 1347–1351. IEEE, 2014.
  * (10) Dirk A Lorenz and Maximilian Winkler.  Minimal error momentum Bregman-Kaczmarz.  arXiv preprint arXiv:2307.15435, 2023.
  * (11) Ion Necoara.  Faster randomized block Kaczmarz algorithms.  SIAM Journal on Matrix Analysis and Applications, 40(4):1425–1452, 2019.
  * (12) Ion Necoara and Dragos Clipici.  Parallel random coordinate descent method for composite minimization: Convergence analysis and error bounds.  SIAM Journal on Optimization, 26(1):197–226, 2016.
  * (13) Ion Necoara and Olivier Fercoq.  Linear convergence of random dual coordinate descent on nonpolyhedral convex problems.  Mathematics of Operations Research, 2022.
  * (14) Ion Necoara, Yu Nesterov, and Francois Glineur.  Linear convergence of first order methods for non-strongly convex optimization.  Mathematical Programming, 175(1):69–107, 2019.
  * (15) Deanna Needell and Joel A Tropp.  Paved with good intentions: analysis of a randomized block Kaczmarz method.  Linear Algebra and its Applications, 441:199–221, 2014.
  * (16) Yu Nesterov.  Efficiency of coordinate descent methods on huge-scale optimization problems.  SIAM Journal on Optimization, 22(2):341–362, 2012.
  * (17) Yurii Nesterov.  Introductory lectures on convex optimization: A basic course, volume 87.  Springer Science & Business Media, 2003.
  * (18) Stefania Petra.  Randomized sparse block Kaczmarz as randomized dual block-coordinate descent.  Analele Stiintifice Ale Universitatii Ovidius Constanta-Seria Matematica, 23(3):129–149, 2015.
  * (19) Tyrrell Rockafellar and Roger J-B Wets.  Variational analysis, volume 317.  Springer Science & Business Media, 2009.
  * (20) Frank Schöpfer and Dirk A Lorenz.  Linear convergence of the randomized sparse Kaczmarz method.  Mathematical Programming, 173(1):509–536, 2019.
  * (21) Frank Schöpfer, Dirk A Lorenz, Lionel Tondji, and Maximilian Winkler.  Extended randomized Kaczmarz method for sparse least squares and impulsive noise problems.  Lineare Algebra and Applications, 652:132–154, 2022.
  * (22) Thomas Strohmer and Roman Vershynin.  A randomized Kaczmarz algorithm with exponential convergence.  Journal of Fourier Analysis and Applications, 15(2):262–278, 2009\.
  * (23) Lionel Tondji and Dirk A Lorenz.  Faster randomized block sparse Kaczmarz by averaging.  Numerical Algorithms, pages 1–35, 2022.
  * (24) Lionel Tondji, Dirk A. Lorenz, and Ion Necoara.  An accelerated randomized Bregman-Kaczmarz method for strongly convex linearly constraint optimization*.  In 2023 European Control Conference (ECC), pages 1–6, 2023.
  * (25) Lionel Tondji, Idriss Tondji, and Dirk A Lorenz.  Adaptive Bregman-Kaczmarz: An approach to solve linear inverse problems with independent noise exactly.  arXiv preprint arXiv:2309.06186, 2023.
  * (26) Lionel N Tondji, Maximilian Winkler, and Dirk A Lorenz.  Linear convergence of a weighted randomized sparse Kaczmarz method.  In Conference: International Conference on Computational Harmonic Analysis, 2021.
  * (27) Paul Tseng.  On accelerated proximal gradient methods for convex-concave optimization.  submitted to SIAM Journal on Optimization, 2(3), 2008.
  * (28) Wotao Yin.  Analysis and generalizations of the linearized Bregman method.  SIAM Journal on Imaging Sciences, 3(4):856–877, 2010.
  * (29) Anastasios Zouzias and Nikolaos M Freris.  Randomized extended Kaczmarz for solving least squares.  SIAM Journal on Matrix Analysis and Applications, 34(2):773–793, 2013.