Comparative study of computational algorithms for the Lasso with high-dimensional, highly correlated data

Abstract

Variable selection is important in high-dimensional data analysis. The Lasso regression is useful since it possesses sparsity, soft-decision rule, and computational efficiency. However, since the Lasso penalized likelihood contains a nondifferentiable term, standard optimization tools cannot be applied. Many computation algorithms to optimize this Lasso penalized likelihood function in high-dimensional settings have been proposed. To name a few, coordinate descent (CD) algorithm, majorization-minimization using local quadratic approximation, fast iterative shrinkage thresholding algorithm (FISTA) and alternating direction method of multipliers (ADMM). In this paper, we undertake a comparative study that analyzes relative merits of these algorithms. We are especially concerned with numerical sensitivity to the correlation between the covariates. We conduct a simulation study considering factors that affect the condition number of covariance matrix of the covariates, as well as the level of penalization. We apply the algorithms to cancer biomarker discovery, and compare convergence speed and stability.

Appendix:

1.1 A Preconditioned conjugate gradient (PCG) method

Conjugate gradient (CG) is a method to resolve positive definite linear equations, A x = b, applied to sparse system that has too large data to solve with Cholesky Decomposition. Instead of directly solving linear equations, this method is to minimize the following function f(x),

$$f(x) = \frac{1}{2}x^{T}Ax -b^{T}x. $$

For a positive definite A, two nonzero vectors u, v are said to be conjugate with respect to A, if they satisfy

$$\langle u,v\rangle_{A} \triangleq u^{T}Av=0. $$

If P is defined as

$$P=\left\lbrace p_{k} : \forall i \neq k, \langle p_{i},p_{k}\rangle_{A}=0 \right\rbrace, $$

it means the set of n number of mutually conjugate directional vectors. Thus, the set P becomes a basis of $\mathbb {R}^{n}$ and x is represented in the form of

$$x=\sum\limits_{i=1}^{n} \alpha_{i} p_{i}. $$

By multiplying both sides by matrix A, b is decomposed by

$$b=Ax=\sum\limits_{i=1}^{n} \alpha_{i} Ap_{i}. $$

Multiplying an arbitrary directional vector p _k ∈P,

$${p_{k}^{T}}b={p_{k}^{T}}Ax=\sum\limits_{i=1}^{n} \alpha_{i} {p_{k}^{T}} A p_{i} = \alpha_{k} {p_{k}^{T}} A p_{k}. $$

Accordingly, the explicit form of α _k can be derived as followed,

$$\alpha_{k}=\frac{{p_{k}^{T}}b}{{p_{k}^{T}}Ap_{k}}=\frac{\langle p_{k},b\rangle} {\|p_{k}\|_{A}^{2}}. $$

If mutually conjugate directional vectors are not given, conjugate gradient (CG) solves the problem iteratively. Set x ₀ as an initial value of x, and a linear equation given by

$$Az=b-Ax_{0} $$

becomes a target function to solve. If we regarding r _k = b−A x _k as k-th residual, r _k becomes a negative gradient of convex function x = x _k , ∇f(x) given by,

$$\nabla f(x_{k}) = Ax_{k} -b, $$

which means that conjugate gradient method moves toward the direction of r _k . Since all directional vectors should satisfy the condition that all vectors are conjugate with respect to A, then k-th direction p _k is given by,

$$p_{k}=r_{k}-\sum\limits_{i>k}\frac{{p_{i}^{T}}Ar_{k}}{{p_{i}^{T}}Ap_{i}}p_{i}. $$

Following this direction, next value of x is updated as followed,

$$x_{k+1}=x_{k} +\alpha_{k} p_{k}, $$

where

$$\alpha_{k}=\frac{{p_{k}^{T}}b}{{p_{k}^{T}}Ap_{k}}=\frac{{p_{k}^{T}}r_{k-1}}{{p_{k}^{T}}Ap_{k}}. $$

Convergence rate of conjugate gradient method depends on condition number of A and especially eigenvalues of A [21]. Accordingly, A x = b problem can be regarded same as linear equation that multiply by inverse matrix of preconditioner given by

$$M^{-1}Ax=M^{-1}b. $$

In choosing an appropriate preconditioner, it should satisfy some necessary conditions.

• M is both symmetric and positive definite matrix.

• M ⁻¹ A is well conditioned and hardly has extreme eigenvalues.

• M x = b is easy to solve.

Widely used preconditioners that satisfy these conditions are the followings;

• 1) Diagonal: M=diag(1/A ₁₁,...,1/A _{n n} ),

• 2) Incomplete(approximate) Cholesky factorization: $M=\hat {A}^{-1}$, where $\hat {A}=\hat {L}\hat {L}^{T}$.