A Random Block-Coordinate Douglas-Rachford Splitting Method with Low Computational Complexity for Binary Logistic Regression
Luis M. Briceno-Arias Giovanni Chierchia Emilie Chouzenoux Jean-Christophe Pesquet
math.OC cs.AI







Random Douglas-Rachford Splitting with Low Computational Complexity






L. M. Briceño-Arias Departamento de Matemática, Universidad Técnica Federico Santa María, Av Espanã 1681, Valparaíso, Chile.
G. Chierchia (corresponding author) and E. Chouzenoux Université Paris Est, LIGM, CNRS UMR 8049, ESIEE Paris, UPEM, Noisy-le-Grand, France.
E. Chouzenoux and J.-C. Pesquet Center for Visual Computing, INRIA Saclay, CentraleSupélec, University Paris-Saclay, Gif sur Yvette, France.





A Random Block-Coordinate Douglas-Rachford Splitting Method with Low Computational Complexity for Binary Logistic Regression
This work was partly supported by the the CNRS MASTODONS project under grant 2016TABASCO.

    
Luis M. Briceño-Arias
Giovanni Chierchia
Emilie Chouzenoux
Jean-Christophe Pesquet

    Received  2009 month day; accepted  2009  month day
========================================================================================================================================================================================================================




In this paper, we propose a new optimization algorithm for sparse logistic regression based on a stochastic version of the Douglas-Rachford splitting method. Our algorithm sweeps the training set by randomly selecting a mini-batch of data at each iteration, and it allows us to update the variables in a block coordinate manner.
Our approach leverages the proximity operator of the logistic loss, which is expressed with the 
generalized Lambert W function.

Experiments carried out on standard datasets demonstrate the efficiency of our approach w.r.t. stochastic gradient-like methods.





§ INTRODUCTION

Sparse classification algorithms have gained much popularity in the context of supervised learning, thanks to their ability to discard irrelevant features during the training stage. Such algorithms aim at learning a weighted linear combination of basis functions that fits the training data
, while encouraging as many weights as possible to be equal to zero. This amounts to solving an optimization problem that involves a loss function plus a sparse regularization term. Different types of classifiers arise by varying the loss function, the most popular being the hinge and the logistic losses <cit.>. 

In the context of supervised learning, sparse regularization traces back to the work of Bradley and Mangasarian <cit.>, who showed that the ℓ_1-norm can efficiently perform feature selection by shrinking small coefficients to zero. Other forms of regularization have also been studied, such as the ℓ_0-norm <cit.>, the ℓ_p-norm with p > 0 <cit.>, the ℓ_∞-norm <cit.>, and other nonconvex terms <cit.>. Mixed-norms have been investigated as well, due to their ability to impose a more structured form of sparsity <cit.>.

Many efficient learning algorithms exist in the case of quadratic regularization, by benefiting from the advantages brought by Lagrangian duality <cit.>. This is unfortunately not true for sparse regularization, because the dual formulation is as difficult to solve as the primal one. Consequently, sparse linear classifiers are usually trained through the direct resolution of the primal optimization problem. Among the possible approaches, one can resort to linear programming <cit.>, gradient-like methods <cit.>, proximal algorithms <cit.>, and other optimization techniques <cit.>. 

Nowadays, it is well known that random updates can significantly reduce the computational time when a quadratic regularization is used <cit.>. Therefore, a great deal of attention has been paid recently to stochastic approaches capable of handling a sparse regularization <cit.>. The list of investigated techniques includes block-coordinate descent strategies <cit.>, stochastic forward-backward iterations <cit.>, random Douglas-Rachford splitting methods <cit.>, random primal-dual proximal algorithms 

<cit.>, and stochastic majorization-minimization methods <cit.>. 



In this paper, we propose a random-sweeping block-coordinate Douglas-Rachford splitting method. In addition to the stochastic behavior, it presents three distinctive features with respect to related approaches <cit.>. Firstly, the matrix to be inverted at the initial step is block-diagonal, while in the concurrent approaches, it did not present any specific structure. The block diagonal property implies that the inversion step actually amounts to inverting a set of smaller size matrices. Secondly, the proposed algorithm can take advantage explicitly from a strong convexity property possibly fulfilled by some of the functions involved in the optimization problem. Finally, the dual variables appear explicitly in the proposed scheme, making it possible to use clever block-coordinate descent strategies <cit.>.

Moreover, the proposed algorithm appears to be well tailored to binary logistic regression with sparse regularization. In contrast to gradient-like methods, our approach deals with the logistic loss through its proximity operator. This results in an algorithm that is not tied up to the Lipschitz constant of the loss function, possibly leading to larger updates per iteration. In this regard, our second contribution is to show that the proximity operator of the binary logistic loss can be expressed in closed form using the generalized Lambert W function <cit.>. We also provide comparisons with state-of-the-art stochastic methods using benchmark datasets.

The paper is organized as follows. In Section <ref>, we derive the proposed Douglas-Rachford algorithm. In Section <ref>, we introduce sparse logistic regression, along with the proximal operator of the logistic loss. In Section <ref>, we evaluate our approach on standard datasets, and compare it to three sparse classification algorithms proposed in the literature <cit.>. Finally, conclusions are drawn in Section <ref>.




Notation: Γ_0() denotes the set of proper, lower semicontinuous, convex functions from a real Hilbert space  to ]-∞,+∞]. 

Let  ψ∈Γ_0(). For every ν∈, the subdifferential of ψ at ν is ∂ψ(ν) = {ξ∈ | (∀ζ∈) ⟨ζ-ν|ξ⟩ +ψ(ν)≤ψ(ζ)}, the proximity operator of ψ at ν is _ψ(ν)= *argmin_ξ∈ 1/2ξ-ν^2 + ψ(ν), and the conjugate of ψ is ψ^*= sup_ξ∈⟨ξ|·⟩ -ψ(ξ) in Γ_0(). 
The adjoint of a bounded linear operator A from  to a real Hilbert space  is denoted by A^*
Let (Ω,ℱ,ℙ) be the underlying probability space,
the σ-algebra generated by a family Φ of random variables is denoted by σ(Φ).



§ OPTIMIZATION METHOD

Throughout this section, _1,…,_B, _1,…,_L are separable real Hilbert spaces. In addition,
= _1 ⊕⋯⊕_B denotes the Hilbertian sum of 
_1,…,_B.
Any vector v∈ can thus be uniquely decomposed as (v_b)_1≤ b 
≤ B where, for every b∈{1,…,B},
v_b ∈_b. In the following, a similar notation will be used to denote vectors in any product space (in bold) and their components.

We will now aim at solving the following problem.


For every b∈{1,…,B} and for every ℓ∈{1,…,L}, let f_b ∈Γ_0(_b), let h_ℓ_ℓ→
be a differentiable convex function with β_ℓ-Lipschitz gradient, for some β_ℓ∈ 
]0,+∞[,

and let A_ℓ,b be a linear bounded operator from _b to _ℓ.  
The problem is to

    _w∈ ∑_b=1^B f_b(w_b) + ∑_ℓ=1^L 
    h_ℓ(∑_b=1^B A_ℓ,b w_b),

under the assumption that the set of solutions ℰ is nonempty.


In order to address Problem <ref>, we propose to employ the random-sweeping 
block-coordinate version of the Douglas-Rachford splitting method with stochastic errors 
described in Algorithm <ref>. Let us define

    (∀ b∈{1,…,B})   C_b = (+τ_b ∑_ℓ=1^L 
    γ_ℓ/1+γ_ℓρ_ℓ A_ℓ,b^*   A_ℓ,b)^-1_b→_b,

where (τ_b)_1≤ b ≤ B and (γ_ℓ)_1≤ℓ≤ 1 are the positive constants introduced in Algorithm <ref>.

The next result establishes the convergence of the proposed algorithm.

	For every b∈{1,…,B}, let w_b^[0], t_b^[0] and (a_b^[i])_i∈ 
	be 
	_b-valued random variables and, for every ℓ∈{1,…,L}, let 
	v_ℓ^[0] and s_ℓ^[0] be 
	_ℓ^B-valued random variables and let (d_ℓ^[i])_i∈ be
	_ℓ-valued random variables. In addition, let 
	(ε^[i])_i∈  be identically 
		distributed {0,1}^B+L∖{0}-valued random variables 
and in Algorithm <ref> assume that


  *   (∀ i∈)
σ(ε^[i]) and 
χ^[i]=σ(t^[0],…,t^[i],s^[0],
…,s^[i])
are independent;

  *  (∀ 
b∈{1,…,B}) ∑_i∈√((a_b^[i]^2|χ^[i])) < 
+∞; 

  *  
(∀ℓ∈{1,…,L}) ∑_i∈√((d_ℓ^[i]^2|χ^[i]))
 < +∞;

  *  (∀ b∈{1,…,B}) [ε_b^[0]=1] > 0 and (∀ℓ∈{1,…,L}) [ε_B+ℓ^[0]=1] > 0.

Then, the sequence (^[i])_i∈ generated by Algorithm <ref> 
converges weakly -a.s. to an ℰ-valued 
random variable.









Problem <ref> can be reformulated as minimizing
+ ∘A
where 

    →]-∞,+∞]↦∑_b=1^B f_b(w_b)
       A→↦(A_ℓ,1w_1,…,A_ℓ,Bw_B)_1≤ℓ≤ L
       →↦∑_ℓ=1^Lh_ℓ(Λ_ℓv_ℓ 
    )
    
    (∀ℓ∈{1,…,L})     Λ_ℓ_ℓ^B→_ℓv_ℓ↦∑_b=1^Bv_ℓ,b

and v=(v_ℓ)_1≤ℓ≤ L denotes a generic element of 
= _1^B⊕⋯⊕_L^B  with 
v_ℓ=(v_ℓ,b)_1≤ b≤ B∈_ℓ^B for every 
ℓ∈{1,…,L}.
Since dom()=, from <cit.>,
Problem <ref> is equivalent to 

    find  ∈  such that  0∈∂() + ^* 
    ∇(),

which, from <cit.> is also equivalent to

    find  (,) ∈×  such that 
      (0,0) ∈N(,)+ 
    S(,),

where
N (,)↦∂()×∂^*() is maximally monotone and S 
(,)↦(^*,-) is a skewed linear operator.
Note that ^*↦ (∑_ℓ=1^LA_ℓ,b^*v_ℓ,b)_1≤ b≤ B
and, from (<ref>), (<ref>) and 
<cit.>,
 ∂↦×_b=1^B∂ f_b(w_b) and 
∂^*↦×_ℓ=1^L∂ 
(h_ℓ∘Λ_ℓ)^*(v_ℓ). Since, for every 
ℓ∈{1,…,L}, 
h_ℓ∘Λ_ℓ is convex differentiable with 
a Bβ_ℓ-Lipschitzian gradient
∇(h_ℓ∘Λ_ℓ)=Λ_ℓ^*
∘∇ 
h_ℓ∘Λ_ℓ,
it follows from Baillon-Haddad theorem 
<cit.> that 
∇(h_ℓ∘Λ_ℓ) is (Bβ_ℓ)^-1-cocoercive and, 
hence,
∂ 
(h_ℓ∘Λ_ℓ)^*
=(∇(h_ℓ∘Λ_ℓ))^-1
is (Bβ_ℓ)^-1-strongly 
monotone. Therefore, for every ρ_ℓ∈[0,(Bβ_ℓ)^-1],  
(h_ℓ∘Λ_ℓ)^* is 
ρ_ℓ-strongly convex.
By defining 

    (∀ℓ∈{1,…,L})   φ_ℓ=(h_ℓ∘Λ_ℓ)^*-ρ_ℓ·^2/2,

it follows from <cit.> that, for every ℓ∈{1,…,L},
φ_ℓ∈Γ_0(_ℓ^B) and, hence,
∂φ_ℓ:=∂(h_ℓ∘Λ_ℓ)^*
-ρ_ℓ 
is maximally monotone.
Consequently,
Problem <ref> can be rewritten equivalently as

    find    (,) ∈×  such that 
      
    
    (∀ b∈{1,…,B})   0∈∂ f_b(w_b)+B_b(,)
    
    (∀ℓ∈{1,…,L})   
    0∈∂φ_ℓ(_ℓ)+B_ℓ(,),

which, for strictly positive constants (τ_b)_1≤ b≤ B and (γ_ℓ)_1≤ℓ≤ L, is equivalent to

    find    (,) ∈×  such that 
      
    
    (∀ b∈{1,…,B})   0∈τ_b ∂ f_b(w_b)+τ_b B_b(,)
    
    (∀ℓ∈{1,…,L})   
    0∈γ_ℓ∂φ_ℓ(_ℓ)+γ_ℓ B_ℓ(,),

where 

    B_b(,)↦∑_ℓ=1^LA_ℓ,b^*v_ℓ,b
    
    B_ℓ(,)↦-
    (A_ℓ,1w_1,…,A_ℓ,Bw_B)+ρ_ℓv_ℓ.

Since S 
(,)↦(^*,-) and D→v↦ 
(ρ_ℓv_ℓ)_1≤ℓ≤ L
are linear and monotone operators in × and , respectively,
the operator 
    B(,)↦
    (^*,-+Dv)=((B_b(,))_1≤ b≤ 
    	B,(B_ℓ(,))_1≤ℓ≤ 
    	L)
 is maximally monotone
in ×. Therefore, by defining the strongly positive diagonal linear operator 

    U×   →×
    
    (,)    ↦ (T,Γ),

where T↦ (τ_b w_b)_1≤ b≤ B and 
Γ↦ (γ_ℓ_ℓ)_1≤ℓ≤ L,
the operator 

    UB(,)↦
    	(T^*,-Γ+
    	ΓDv)=((τ_bB_b(,))_1≤ b≤ 
    		B,(γ_ℓB_ℓ(,))_1≤ℓ≤ 
    		L)

is maximally monotone
in (×, ·_U^-1), where

    (∀  (,) ∈×)    (,) _U^-1 = √(∑_b=1^B τ_b^-1w_b^2 + ∑_ℓ 
    = 1^L γ_ℓ^-1_ℓ^2).

Note that the renormed product space (×, 
·_U^-1) is the Hilbert sum 
_1⊕⋯_B⊕^B_1⊕⋯^B_L where, for every 
b∈{1,…,B} and ℓ∈{1,…,L}, _b and ^B_ℓ are endowed 
by the norm 
·_τ_b 
w_b↦w_b/√(τ_b)
and
·_γ_ℓv_ℓ↦v_ℓ/√(γ_ℓ), respectively. Therefore,
since τ_b∂ f_b and γ_ℓ∂φ_ℓ are 
maximally monotone in (_b,·_τ_b) and 
(^B_ℓ,·_γ_ℓ), respectively,
we conclude that (<ref>) is a particular case of the primal inclusion in 
<cit.>.

Now we  write Algorithm <ref> as a particular case of the random 
block-coordinate 
Douglas-Rachford splitting proposed in <cit.> applied to 
(<ref>) in (×, ·_U^-1).
Given 
(t,s)∈×, let 








(w,v)=J_UB(t,
s)=(+
UB)^-1(t,s).
It follows from (<ref>) that

    =  - T^* 
    = (+ ΓD)^-1( + Γ),

which leads to

    = ( + T^*  (+ ΓD)^-1Γ)^-1( - T^*  (+ 
    ΓD)^-1).

In order to derive an explicit formula for the matrix inversion in (<ref>), 
set =  - T^*  (+ 
ΓD)^-1.



We have =+T^*  (+ 
ΓD)^-1Γ
and, since (<ref>) and D is diagonal, we obtain

    T^*  (+ 
    ΓD)^-1Γ↦(τ_b∑_ℓ=1^L
    γ_ℓ A_ℓ,b^*A_ℓ,bw_b/1+γ_ℓρ_ℓ)_1≤ b≤ B,

and, hence,

    (∀ b∈{1,…,B})    
    w_b=(+τ_b∑_ℓ=1^L  
    γ_ℓ A_ℓ,b^*A_ℓ,b/1+γ_ℓρ_ℓ)^-1z_b
    =C_b z_b.

Therefore, (<ref>) can be written equivalently as

    (∀ b∈{1,…,B})    
    w_b=C_b(t_b^[i]-
    τ_b∑_ℓ=1^LA_ℓ,b^*s_ℓ,b/1+γ_ℓρ_ℓ)=
    C_b(t_b-τ_b u_b),

where 

    (∀ b∈{1,…,B})      u_b= ∑_ℓ=1^L A_ℓ,b^*   
    s_ℓ,b/1+γ_ℓρ_ℓ.

Moreover, from (<ref>), we deduce that

    (∀ℓ∈{1,…,L})  _ℓ=_ℓ+
    	γ_ℓ (A_ℓ,bw_b)_1≤ b≤ B/1+γ_ℓρ_ℓ,

and, hence, we have J_UB(t,
s)↦((Q_b(t,
s))_1≤ b≤ B,(Q_ℓ(t,
s))_1≤ℓ≤ L), where 

    (∀ b∈{1,…,B})   Q_b (t,)↦ C_b(t_b-
    τ_bu_b)
     (∀ℓ∈{1,…,L})   Q_ℓ (t,)↦_ℓ+
    	γ_ℓ (A_ℓ,bQ_b(t,))_1≤ b≤ B/1+γ_ℓρ_ℓ.

Now, it follows from <cit.> that, for every 
b∈{1,…,B} and ℓ∈{1,…,L}, J_τ_b∂ 
f_b=_τ_b f_b and
 J_γ_ℓ∂φ_ℓ=_γ_ℓφ_ℓ
and,
















for every ℓ∈{1,…,L} 
	and 
(_ℓ,_ℓ)∈_ℓ^B×_ℓ^B, we have 

    _ℓ = _γ_ℓφ_ℓ_ℓ⇔     _ℓ-_ℓ/γ_ℓ∈∂φ_ℓ(_ℓ)
    ⇔     _ℓ-_ℓ/γ_ℓ∈∂(h_ℓ∘Λ_ℓ)^*_ℓ 
    -ρ_ℓ_ℓ
    ⇔     _ℓ-(1-γ_ℓρ_ℓ)_ℓ∈γ_ℓ∂(h_ℓ∘Λ_ℓ)^*_ℓ .

Therefore, if γ_ℓρ_ℓ  <1 we have that (<ref>) is 
equivalent to

    _ℓ = _γ_ℓ/1-γ_ℓρ_ℓ  
    	(h_ℓ∘Λ_ℓ)^*(_ℓ/1-γ_ℓρ_ℓ)

and, from Moreau's decomposition formula <cit.>, 
we 
obtain

    _ℓ =  1/1-γ_ℓρ_ℓ(_ℓ - γ_ℓ_1-γ_ℓρ_ℓ/γ_ℓ 
    	(h_ℓ∘Λ_ℓ)(_ℓ/γ_ℓ) 
    ).

Noting that, for every ℓ∈{1,…,L}, 
Λ_ℓ∘Λ_ℓ^*=B,
from <cit.>, (<ref>) and 
(<ref>) we deduce that

    (∀ b∈{1,…,B})   r_ℓ,b=
    1/B(1-γ_ℓρ_ℓ)(∑_d=1^B 
    z_ℓ,d-γ_ℓ_B(1-γ_ℓρ_ℓ)/γ_ℓ 
    	h_ℓ(1/γ_ℓ∑_d=1^B 
    z_ℓ,d))

and, hence,

    _γ_ℓφ_ℓ_ℓ=1/B(1-γ_ℓρ_ℓ)(∑_d=1^B 
    z_ℓ,d-γ_ℓ_B(1-γ_ℓρ_ℓ)/γ_ℓ 
    	h_ℓ(1/γ_ℓ∑_d=1^B 
    z_ℓ,d))_1≤ b≤ B.

Therefore, by defining, for every i∈ and ℓ∈{1,…,L}, 
e_ℓ^[i]∈_ℓ^B 
via

    (∀ℓ∈{1,…,L})  e_ℓ^[i] =  
    (-γ_ℓ/B(1-γ_ℓρ_ℓ) d_ℓ^[i])_1≤ b≤ B,

we deduce that Algorithm <ref> can be written equivalently as


For   i = 0, 1, …


    ⌊   For b =1,…,B
       ⌊   w_b^[i+1]=w_b^[i]+ε_b^[i](Q_b(t^[i],^[i])-w_b^[i])
       t_b^[i+1]=t_b^[i] + ε_b^[i]μ^[i]( _τ_b 
    	f_b(2w_b^[i+1]-t_b^[i]) + 
    a_b^[i] - w_b^[i+1])
    .
       For ℓ =1,…,L
       ⌊   _ℓ^[i+1]=_ℓ^[i]+ε_B+ℓ^[i](Q_ℓ(t^[i],^[i])-_ℓ^[i])
       _ℓ^[i+1]=_ℓ^[i] + ε_B+ℓ^[i]μ^[i]( _γ_ℓφ_ℓ(2_ℓ^[i+1]-_ℓ^[i]) + 
    e_ℓ^[i] - _ℓ^[i+1]).
    .
    .

Defining, for every i∈, a^[i]=((a_b^[i])_1≤ b≤ 
B,(e_ℓ^[i])_1≤ℓ≤ L)∈×, we have 

    ∑_i∈√((a^[i]^2_U^-1|χ^[i]))   =
    ∑_i∈√(∑_b=1^Bτ_b^-1(a_b^[i]^2|χ^[i])+
    	∑_ℓ=1^Lγ_ℓ^-1(e_ℓ^[i]^2|χ^[i]))
       ≤∑_b=1^B τ_b^-1/2∑_i∈√((a_b^[i]^2|χ^[i]))+
    	∑_ℓ=1^Lγ_ℓ^-1/2∑_i∈√((e_ℓ^[i]^2|χ^[i]))
       <+∞,

where the last inequality follows from <ref>, <ref>,
(<ref>) and

    ∑_i∈√((e_ℓ^[i]^2|χ^[i]))
    =γ_ℓ/√(B)(1-γ_ℓρ_ℓ)∑_i∈√((d_ℓ^[i]^2|χ^[i]))<+∞.

Altogether, since operator J_UB is weakly 
sequentially continuous because it is continuous and linear, the result follows from
<cit.> when the error term in the computation of
J_UB 
is zero.






























  * In Proposition <ref>, 

the binary variables ε_b^[i] and ε_B+ℓ^[i] signal whether the variables t_b^[i] and _ℓ^[i] are activated or not at iteration i. Assumption <ref> guarantees that each of the latter variables is activated with a nonzero probability at each iteration.
In particular, it must be pointed out that the variables  p_ℓ^[i] and q_ℓ^[i] only need to be computed when  ε_B+ℓ^[i]=1.

  * Note that Algorithm <ref> may look similar to the stochastic approach proposed 
in  <cit.> (see also <cit.>, and 
<cit.> for deterministic variants). It exhibits however three key differences. Most 
importantly, the operator inversions performed at the initial step amount to inverting a set of 
positive definite self-adjoint operators defined on the spaces (_b)_1≤ b ≤ B. We 
will see in our application that this reduces to invert a set of small size symmetric  positive 
definite matrices. Another advantage is that the smoothness of the functions (h_ℓ)_1 
	≤ℓ≤ L is taken into account, and a last one is that the dual variables  appear 
explicitly in the iterations.

  * If, for every ℓ∈{1,…,L}, ρ_ℓ = 0 and B = 1, Algorithm <ref> simplifies to Algorithm <ref>, where unnecessary indices have been dropped and we have set

    (∀ i ∈)  u^[i] = -τ u^[i]
    
    (∀ℓ∈{1,…,L}  s_ℓ^[i]  = -τ s_ℓ^[i].

In this case,

    C = (+τ∑_ℓ=1^L γ_ℓ A_ℓ^*   A_ℓ)^-1.


















































































































































































































When (∀ℓ∈{1,…,L}) γ_ℓ = 1/τ, it turns out  this algorithm is exactly the same as the one resulting from a direct application of <cit.><cit.>.


  * The situation when, for a given ℓ, h_ℓ∈Γ_0(_ℓ) is not Lipschitz-differentiable can be seen as the limit case when β_ℓ→ +∞. It can then be shown that Algorithm <ref> remains valid by setting ρ_ℓ = 0. 







§ SPARSE LOGISTIC REGRESSION

The proposed algorithm can be applied in the context of binary linear classification. 







A binary linear classifier can be modeled as a function that predicts the output y ∈{-1,+1} associated to a given input ∈^N. This prediction is defined through a linear combination of the input components, yielding the decision variable

    d_(x) = ( ^⊤),

where ∈^N is the weight vector to be estimated.
In supervised learning, this weight vector is determined from a set of input-output pairs

    𝒮 = { (_ℓ,y_ℓ) ∈^N ×{-1,+1} | ℓ∈{1,…,L}},

which is called training set. More precisely, the learning task can be defined as the trade-off between fitting the training data and reducing the model complexity, leading to an optimization problem expressed as

    _∈^N() + ∑_ℓ=1^L h(y_ℓ _ℓ^⊤),

where ∈Γ_0(^N) is a regularization function and
h∈Γ_0() stands for the loss function. In the context of sparse learning, a 
popular choice for the regularization is the ℓ_1-norm.





Although many choices for the loss function are possible, we are primarily interested in the logistic loss, which is detailed in the next section.





 §.§ Logistic regression

Logistic regression aims at maximizing the posterior probability density function of the weights given the training data, here assumed to 
be a realization of statistically independent input-output random variables. This leads us to

    _∈^N
    	φ() ∏_ℓ=1^L π(y_ℓ|_ℓ,)θ_ℓ(_ℓ|),

where φ  is the weight prior probability density function and, for every ℓ∈{1,…,L}, 
θ_ℓ is the conditional data likelihood of the ℓ-th input knowing the weight values,
while π(y_ℓ|_ℓ,) is the conditional probability of the ℓ-th output knowing the ℓ-th input and the weights.
Let us model this conditional probability with the sigmoid function defined as

    π(y_ℓ|_ℓ,) = 1/1+exp(-y_ℓ_ℓ^⊤),

and assume that the inputs and the weights are statistically independent and
that φ is log-concave. 

Then, the negative-logarithm of the energy in (<ref>) yields an instance of Problem (<ref>) in which

    (∀ v∈)     h(v) = log(1+exp(-v))

and, for every ∈^N, () = -logφ(). 
(The term ∏_ℓ=1^L θ_ℓ(_ℓ|) can be discarded since  the inputs and the weights are assumed statistically independent.)

The function in (<ref>) is called logistic loss. For completeness, note 
that other loss functions, leading to different kinds of classifiers, are the hinge loss <cit.>

    (∀ v∈)     h^hinge(v) = (max{0,1-v})^ q

with q∈{1,2}, and the Huber loss <cit.>

    (∀ v∈)     h^huber(v) = 
    
    0    if v≥1
    
    -v    if v≤-1
    1/4(v-1)^2    otherwise.

These functions can be also handled by the proposed algorithm.



 §.§ Optimization algorithm

Let us blockwise decompose the weight variable ∈^N as

    ^⊤=[ w_1^⊤ …  w_B^⊤],

where, for every b∈{1,…,B},  w_b∈^N_b and
N_1,…,N_B are strictly positive integers such that
N=N_1+…+N_B. Let us also decompose the input vector 
as ^⊤=[ x_1^⊤ …  x_B^⊤]. Finally, let us assume that the regularization function
is block-separable, i.e. f=⊕_b=1^Bf_b,
where, for every b∈{1,…,B}, f_b ∈Γ_0(^N_b).
A typical example of such functions is given by

    (∀ b ∈{1,…,B})
         f_b = λ·_κ_b,

where λ∈ [0,+∞[ and ·_κ_b, κ_b∈ [1,+∞],
denotes the ℓ_κ_b-norm of ^N_b.
In particular, when for every b∈{1,…,B} κ_b = 1,  reduces to the standard
ℓ_1 regularizer, whereas setting κ_b ≡ 2
results in a potential promoting group sparsity <cit.>.

In the context described above, (<ref>) is a particular case of 
Problem <ref>
where, for every b ∈{1,…,B},  _b = ^N_b, 
for every ℓ∈{1,…,L},
_ℓ =, h_ℓ=h and A_ℓ,b=y_ℓx_ℓ,b^⊤.
 Note that since, for every ℓ∈{1,…,L}, y_ℓ^2=1, A_ℓ,b^*A_ℓ,b =  
 x_ℓ,bx_ℓ,b^⊤. Moreover,
h defined in (<ref>) is twice differentiable with

    (∀ v ∈)    
    h'(v)   =-exp(-v)/1+exp(-v), 
    
    h”(v)   =exp(-v)/(1+exp(-v))^2.

 Since h” is maximized in v=0, we have
sup_v∈|h”(v)|=1/4, which implies that h' is 1/4-Lipschitz continuous
and we have thus, for every ℓ∈{1,…,L}, β_ℓ = 1/4.

The problem can thus be solved with Algorithm <ref>, the convergence of which 
is guaranteed almost surely under the assumptions of Proposition <ref>.



 §.§ Proximity operator

An efficient computation of the proximity operators of functions (f_b)_1≤ b ≤ B and h plays a crucial role in the implementation of Algorithm <ref>. 
There exists an extensive literature on the computation of the proximity operators of functions like (<ref>) <cit.>.
In particular,
when κ_b = 1 (resp. κ_b = 2), this proximity operator reduces to a 
component-wise (resp. blockwise) soft-thresholding <cit.>.



Regarding the logistic loss in (<ref>), although some numerical methods exist <cit.>, to the best of our knowledge, 
no thorough investigation of the form of its proximity operator has been made. 
The next proposition will contribute to fill such a void.

The result relies on the generalized W-Lambert function  recently 
introduced in <cit.>, defined via 

    (∀v̅∈)(∀ v∈)(∀ r∈]0,+∞[)  v̅ 
    (exp(v̅) + r)= v   ⇔  v̅ = W_r(v).
 
When r∈ [exp(-2),+∞[, W_r is uniquely defined and strictly increasing, but when r ∈ ]0,exp(-2)[, there exist three branches for W_r.
We will retain the only one which can take nonnegative values  (denoted by W_r,0 in <cit.>) and is also strictly increasing.
This function can be efficiently evaluated through a Newton-based method devised by Mező et al. <cit.> and available on line.[https://sites.google.com/site/istvanmezo81/others]


	
Let γ∈]0,+∞[ and let h v↦log(1+exp(-v)). We have

    (∀ v∈)  _γ h (v) = v + W_exp(-v)(γexp(-v ) ).



Let v∈ and γ∈]0,+∞[. For every p∈, it follows from the 
definition of _γ h, (<ref>) and (<ref>) that

    p=_γ h(v)     ⇔    
    v-p=-γexp(-p)/1+exp(-p)=-γ/exp(p)+1
       ⇔    
    (p-v)(exp(p)+1)=γ
       ⇔    
    (p-v)(exp(p-v)+exp(-v))=γexp(-v)
       ⇔    
    p-v=W_exp(-v)(γexp(-v))

and the result follows.


From a numerical standpoint, it must be emphasized that the exponentiation in (<ref>) may be problematic, as it yields an arithmetic overflow when v tends to -∞. To overcome this issue, one can use the asymptotic equivalence[Hereafter, following Landau's notation,  we will write that F(v) = ϑ(G(v)),
where F→ and G→, if F(v)/G(v) → 0
as v → +∞ (or v → -∞).]
between the proximity operator of the logistic function and other more tractable functions.

Let γ∈]0,+∞[ and let h v↦log(1+exp(-v)). Then, as 
v → -∞,

    _γ h(v) 
    = v + γ(1 - exp(γ + v) + (1+γ)  exp( 2(γ + v))) 
    + ϑ(exp(2v)).
 


Define

    (∀ v∈)   φ(v) = W_exp(-v)(γexp(-v)).

According to Proposition <ref>, 

    φ=_γ h-.

It follows from  (<ref>) that 
ran  φ= ran(_γ h-)⊂]0,γ[.
Moreover, from <cit.> we deduce that φ=-γ_h^*/γ(·/γ), 
which is decreasing and continuous by virtue of
<cit.>. Therefore, lim_v→-∞φ(v) exists and
from (<ref>) we have

    (∀ v∈)   φ(v)exp(φ(v))=(γ-φ(v))exp(-v).

Since the left side of the equality above is bounded, we deduce that 
lim_v→-∞φ(v)=γ. Subsequently, we define

    (∀ v∈)   u(v)=φ(v)-γ  satisfying  lim_v→-∞u(v)=0.

Hence, (<ref>) can be rewritten as

    (γ+u(v))exp(γ)exp(u(v))=-u(v)exp(-v)

and by using the first order Taylor expansion around ξ=0, exp(ξ)=1+ξ+ϑ(ξ) 
and the fact that u(v)^2+(γ+u(v))ϑ(u(v))=ϑ(u(v)),
we obtain

    u(v)   =-exp(γ+v)(γ+u(v))(1+u(v)+ϑ(u(v)))
       =-exp(γ+v)(γ+(γ+1)u(v)+ϑ(u(v)))
       =-γexp(γ+v)-(γ+1)exp(γ+v)u(v)-exp(γ+v)ϑ(u(v)).

We deduce from this relation that

    u(v)=-γexp(γ+v)+exp(γ+v)ϑ(u(v))/1+(γ+1)exp(γ+v).

It follows that

    lim_v → -∞ u(v)exp(-v) =- γexp(γ),

which implies that exp(γ+v)ϑ(u(v))=ϑ(exp(2v))


and, from (<ref>) we obtain

    u(v)=-γexp(γ+v)/1+(γ+1)exp(γ+v)+ϑ(exp(2v)).

Combining (<ref>), (<ref>), and (<ref>) yields

    _γ h(v)   =v+γ(1-exp(γ+v)/1+(γ+1)exp(γ+v))+ϑ(exp(2v))
       =v+γ(1+γexp(γ+v)/1+(γ+1)exp(γ+v))+ϑ(exp(2v))
       =v+γ(1-exp(γ+v)+(1+γ)exp(2(γ+v)))
    +ϑ(exp(2v)),

where the last equality follows from the second order Taylor expansion around ξ = 0.






§ EXPERIMENTAL RESULTS

In order to assess the performance of Algorithm <ref>, we performed the training on standard datasets[http://www.csie.ntu.edu.tw/∼cjlin/libsvmtools/datasets/binary.html][http://yann.lecun.com/exdb/mnist] (see Table <ref>), and we compared it with the following approaches.


  * Stochastic Forward-Backward splitting (SFB) <cit.>

    w^[0]∈^N
       For i=0,1,…
       ⌊   Select ^[i]⊂{1,…,L}
    
    
    
        w^[i+1] = _γ_i f( w^[i] - γ_i ∑_ℓ∈^[i] y_ℓ x_ℓ h'(y_ℓ x_ℓ^⊤ w^[i]) )
    .

where (γ_i)_i∈ is a decreasing sequence of positive values.


  * 
Regularized Dual Averaging (RDA) <cit.>

    w^[0]∈^N,  z^[0] = 0
       For i=0,1,…
       ⌊   Select ^[i]⊂{1,…,L}
    
        z^[i+1] = z^[i] + ∑_ℓ∈^[i] y_ℓ x_ℓ h'(y_ℓ x_ℓ^⊤ w^[i])
    
    
    
        w^[i+1] = _γ_i f( - γ_i   z^[i+1])
    .

where (γ_i)_i∈ is a decreasing sequence of positive values.


  * Block-Coordinate Primal-Dual splitting (BCPD) <cit.>

    w^[0]∈^N,  v^[0]∈^L
    
       For i=0,1,…
       ⌊   Select ^[i]⊂{1,…,L}
    
       w^[i+1] = _τ f(w^[i] - τ u^[i])
    
        (∀ℓ∈^[i])   v^[i+1]_ℓ = _σ h^*(v^[i]_ℓ + σ y_ℓ x_ℓ^⊤(2w^[i+1] - w^[i]))
    
        (∀ℓ∉_i)   v^[i+1]_ℓ = v^[i]_ℓ
    
       u^[i+1] = u^[i] + ∑_ℓ∈^[i](v_ℓ^[i+1]-v_ℓ^[i])y_ℓ x_ℓ
    .

where τ>0 and σ>0 are such that τσ∑_ℓ=1^L x_ℓ x_ℓ^⊤≤ 1.


The algorithmic parameters are reported in Table <ref>. For all the algorithms, mini-batches of size 1000 were randomly selected using a uniform distribution, and the initial vector w^[0] was randomly drawn from the normal distribution with zero mean and unit variance. For the datasets with more than two classes (MNIST and RCV1), the “one-versus-all” approach is used <cit.>. All experiments were carried out with Matlab 2015a on an Intel i7 CPU at 3.40 GHz and 12 GB of RAM.





Table <ref> reports the classification performance achieved by the compared algorithms, which includes the classification errors computed on a (disjoint) test set, as well as the sparsity degree of the solution. For all the considered datasets, the regularization parameter λ was selected with a cross-validation procedure. The results show that the proposed algorithm finds a solution that yields the same accuracy as state-of-the-art methods, while being sparser than the ones produced by gradient-like methods (SFB and RDA). 

Figure <ref> reports the training performance versus time of the considered algorithms, which includes the criterion in (<ref>), and the distance to the solution w^[∞] obtained after many iterations for each compared method. The results indicate that the proposed approach converges faster to a smaller value of the objective criterion. This could be related to the implicit preconditioning present in Algorithm <ref> through the matrix Q. Another interesting feature of our algorithm is the free choice of parameters γ and μ_i. Conversely, in both SFB and RDA, the parameter γ_i (also referred to as learning rate) needs to be carefully selected by hand, causing such algorithms to slow down or even diverge if the learning rate is chosen too small or too high. 







§ CONCLUSION

In this paper, we have proposed a block-coordinate Douglas-Rachford algorithm for sparse logistic regression. In contrast to gradient-like methods, our approach relies on the proximity operator of the logistic loss, for which we derived a closed-form expression that can be efficiently implemented. Thanks to this feature, our approach removes restrictions on the choice of the algorithm parameters, unlike gradient-like methods, for which it is essential that the learning rate is carefully chosen. This is confirmed by our numerical results, which indicate that the training performance of the proposed algorithm compares favorably with state-of-the-art stochastic methods. 

IEEEtran



	
