Block-Coordinatewise Methods for Sparse Optimization with Nonsmooth Regularization

Block-Coordinatewise Methods for Sparse Optimization with Nonsmooth Regularization

Paul Tseng

Mathematics, University of Washington Seattle

Hong Kong Polytechnic University February 20, 2008

(Joint works with Sylvain Sardy (Univ. Geneva) and Sangwoon Yun (NUS)) Abstract

This is a talk given at UBC, April 2007.

Talk Outline

• Block-Coordinate Minimization

? Properties

? Application I: Basis Pursuit/Lasso

• Block-Coordinate Gradient Descent

? Properties

? Application II: Group Lasso for Logistic Regression

? Application III: Sparse Inverse Covariance Estimation

• Conclusions & Future Work

Block-Coordinate Minimization

min

x=(x₁,...,x_n)f (x) f : <ⁿ → < is cont. diff.

Given x ∈ <ⁿ, choose i ∈ {1, ..., n}. Update x^new = arg min

u|u_j=x_j ∀j6=i

f (u).

Repeat until “convergence”.

Gauss-Seidel: Choose i cyclically, 1, 2,..., n, 1, 2,...

Gauss-Southwell: Choose i with |_∂x^∂f

i(x)| maximum.

Example: min

x=(x₁,x₂)(x₁ + x₂)² +

4(x₁ − x₂)²

Properties

• This method extends to update a block of coordinates at each iteration.

• It is simple, and efficient for large “weakly coupled” problems (off-block-diagonals of ∇²f (x) not too large).

• If f is convex, then very cluster point of the x-sequence is a minimizer. ^{Zadeh ’70}

• If f is nonconvex, then G-Seidel can cycle ^{Powell ’73} but G-Southwell still converges.

• Can get stuck if f is nondifferentiable.

Example: min

x=(x₁,x₂)(x₁ + x₂)² +

4(x₁ − x₂)² + |x₁ + x₂|

But, if the nondifferentiable part is separable, then convergence is possible.

Example: min

x=(x₁,x₂)(x₁ + x₂)² + 1

4(x₁ − x₂)² + |x₁| + |x₂|

Application I: Basis Pursuit/Lasso

minx F_c(x) := kAx − bk²₂ + ckxk₁

Tibshirani ’96, Fu ’98 Osborne et al. ’98

Chen, Donoho, Saunders ’99 ...

A ∈ <^m×n, b ∈ <^m, c ≥ 0.

• Typically m ≥ 1000, n ≥ 8000, and A is dense. k · k₁ is nonsmooth.

6. .

_

• Can reformulate this as a convex QP and solve using an IP method. ^Chen,

Donoho, Saunders ’99

Assume the columns of A come from an overcomplete set of basis functions associated with a fast transform (e.g., wavelet packets).

BCM for BP

Given x, choose I ⊆ {1, ..., n} with |I| = m and {A_i}_i∈I orthog. Update x^new = arg min

u_i=x_i ∀i6∈I

F_c(u) ^{has closed}

← ^{form soln} Repeat until “convergence”.

Gauss-Southwell: Choose I to maximize min

v∈∂_xIF_c(x)

kvk₂.

• Finds I in O(n + m log m) opers. by algorithm of Coifman & Wickerhauser.

• x-sequence is bounded & each cluster point minimizes F_c. Sardy, Bruce, T ’00

Electronic surveillance

0 10 20 30 40 50 60

bp.cptable.spectrogram

50 55 60 65

0 10 20 30 40 50 60

ew.cptable.spectrogram

microseconds

KHz

-40 -30 -20 -10 0 10 20

m = 2¹¹ = 2048, c = 4, b: top image, A: local cosine transform, all but 4 levels

Method efficiency

2500 2000

1500

1000

500

1 5 10

















 









 


















 


          !" #$

%&

CPU time (log scale) CPU time (log scale)

Real part Imaginary part

Objective

1 5 10

1000 500 1500 2000 2500 3000

BCR

IP IP

BCR

Comparing CPU times of IP and BCM (S-Plus, Sun Ultra 1).

• IP method requires fewer iterations, but each iteration is expensive (many CG steps per iteration). No good preconditioner for CG is known.

• BCM method requires more iterations, but each iteration is cheap.

• Convergence of BCM depends crucially on

• differentiability of k · k²₂

• separability of k · k₁

• convexity of F_c (stationarity ⇒ global minimum)

Generalization to ML Estimation with `

-Regularization

minx −`(Ax; b) + c X

i∈J

|x_i|

` is log likelihood, {A_i}_i6∈J are lin. indep “coarse-scale Wavelets”, c ≥ 0

• −`(y; b) = ¹₂ky − bk²₂ Gaussian noise

• −`(y; b) =

m

X

i=1

(y_i − b_i ln y_i) (y_i ≥ 0) Poisson noise

Can solve this problem by adapting IP method. But IP method is slow (many CG steps per IP iteration) Antoniadis, Sardy, T ’04.

6. .

_

Adapt BCM method?

Optimization with Nonsmooth Regularization

minx F_c(x) := f (x) + cP (x)

f : <ⁿ → < is cont. diff. c ≥ 0.

P : <ⁿ → (−∞, ∞] is proper, convex, lsc, and block-separable, i.e., P (x) = P

I∈C P_I(x_I) (I ∈ C partition {1, ..., n}).

• P (x) = kxk₁ Basis Pursuit/Lasso

• P (x) = P

I∈C kx_Ik₂ group Lasso

• P (x) = n0 if l ≤ x ≤ u

∞ else bound constrained NLP

Idea

Block-Coord. Gradient Descent Method

For x ∈ domP, ∅ 6= I ⊆ {1, ..., n}, and H  0, let d_H(x; I) and q_H(x; I) be the optimal soln and obj. value of

min

d|d_i=0 ∀i6∈I

{∇f (x)^Td + 1

2d^THd + cP (x + d) − cP (x)} direc.

subprob

Facts:

• d_H(x; {1, ..., n}) = 0 ⇔ F_c⁰(x; d) ≥ 0 ∀d ∈ <ⁿ. stationarity

• H is diagonal, P is separable ⇒ d_H(x; I) = X

i∈I

d_H(x; i),

q_H(x; I) = X

i∈I

q_H(x; i). separab.

• H is diagonal, P is “simple” ⇒ d_H(x; I) has “closed form”.

Given x ∈ domP, choose I ⊆ {1, ..., n}, H  0.

Update x^new = x + αd_H(x; I) (α > 0) until “convergence.”

Gauss-Seidel: Choose I ∈ C cyclically.

Gauss-Southwell: Choose I with

q_D(x; I) ≤ υ q_D(x; {1, ..., n})

(0 < υ ≤ 1, D  0 is diagonal, e.g., D = I or D = diag(H)).

Inexact Armijo LS: α = largest element of {1, β, β², ...} satisfying F_c(x + αd) − F_c(x) ≤ 0.1 α q_H(x; I) (0 < β < 1)

Convergence properties ^{T, Yun ’06}:

(a) If λI  D, H  ¯λI (0 < λ ≤ ¯λ), then every cluster point of the x-sequence generated by BCGD method is a stationary point of F_c.

(b) If in addition P and f satisfy any of the following assumptions, then the x-sequence converges linearly in the root sense.

A1 f is strongly convex, ∇f is Lipschitz cont. on domP. A2 f is (nonconvex) quadratic. P is polyhedral.

A3 f (x) = g(Ax) + q^Tx, where A ∈ <^m×n, q ∈ <ⁿ, g is strongly convex, ∇g is Lipschitz cont. on <^m. P is polyhedral.

Note: BCGD has stronger global convergence property (and cheaper iteration) than BCM.

Application II: Group Lasso for Logistic Regression

minx f (x) + c X

I∈C

ω_Ikx_Ik₂

Yuan, Lin ’06 Kim3 _’06

Meier, van de Geer, B ¨uhlmann ’06 ...

c > 0, ω_I > 0.

f (x) = PN

j=1 log

1 + e^{aTj x}

− y_ja^T_j x (a_j ∈ <ⁿ, y_j ∈ {0, 1})

• f is convex, cont. diff. k · k₂ is convex, nonsmooth. In prediction of short DNA motifs, n > 1000, N > 11, 000.

• BCM-GSeidel has been used Yuan, Lin ’06, but each iteration is expensive. Every cluster point of the x-sequence is a minimizer ^{T ’01}.

• BCGD-GSeidel is significantly more efficient Meier et al ’06. Every cluster point of the x-sequence is a minimizer ^{T, Yun ’06}. Linear convergence?

Application III: Sparse Inverse Covariance Estimation

X∈Smin₊ⁿ f (X) + ckXk₁

Meinshausen, B ¨uhlmann ’06 Yuan, Lin ’07

Banerjee, El Ghaoui, d’Aspremont ’07 Friedman, Hastie, Tibshirani ’07

c > 0, kXk₁ = P

ij |X_ij|,

f (X) = − log detX + tr(XS) (S ∈ S₊ⁿ is empirical covariance matrix)

• f is strictly convex, cont. diff. on its domain, O(n³) ops to evaluate. k · k₁ is convex, nonsmooth. In applications, n can exceed 6000.

The Fenchel dual problem Rockafellar ’70 is a bound-constrained convex program:

min

W ∈S₊ⁿ,kW −Sk_∞≤c − log det(W ) kY k_∞ = max_ij |Y_ij|.

• IP method requires O(n log(1/)) ops to find -optimal soln. Impractical!

Nesterov’s first-order smoothing method requires O(n^4.5/) ops Banerjee et al ’07.

• Use BCM-GSeidel to solve the dual problem, cycling thru columns

i = 1, ..., n of W. Each iteration reduces (via determinant property & duality) to min

ξ∈<ⁿ⁻¹

1

2ξ^TW_i^¬_i^¬ξ − S_i^T¬_iξ + ckξk₁.

Solve this using IP method (O(n³) ops) Banerjee et al ’07 or BCM-GSeidel Friedman et al ’07.

• Can apply BCGD-GSeidel to either primal or dual problem. More efficient?

Applied to the primal, each iteration entails

u∈<minⁿ



tr((−X⁻¹ + S)D) + 1

2u^THu + ckX + Dk₁



D=u^Te_i+e_iu^T

.

For diagonal H, the minimizing D has closed form! For each trial α in the Armijo LS, det(X + αD) can be evaluated from detX and X⁻¹ in O(n²) ops.

Update X⁻¹ in O(n²) ops. Similar application to the dual. Global convergence, asymptotic linear convergence, complexity analysis... Toh, T, Yun, forthcoming.

Conclusions & Future Work

1. Nonsmooth regularization induces sparsity in the solution, avoids oversmoothing signals, and is useful for variable selection.

2. The regularized problem can be solved effectively by BCM or BCGD, exploiting the problem structure.

3. Extension of BCM, BCGD to handle linear constraints Ax = b is possible, including Support Vector Machine training T, Yun, ’07, ’08. Some open questions on efficient implementation and convergence analysis remain.

4. Many other applications, including stochastic volatility models Neto, Sardy, T, forthcoming.

5. Extension of BCGD to nonconvex nonsmooth regularization is possible (e.g. `_p-regularization, 0 < p < 1) Sardy, T, forthcoming.