Parameter Choices for Sparse Multi-Parameter Regularization with the $\ell_{1}$ Norm

Multi-parameter regularization is a widely used technique for addressing ill-posed problems [4, 6, 8, 9, 10, 18, 19, 28]. Motivated by the challenges posed by big data in practical applications, sparse multi-parameter regularization using the $\ell_{1}$ norm has become a prominent area of research [1, 11, 21, 22, 23, 26]. Empirical studies have shown that, compared to single-parameter regularization with the $\ell_{1}$ norm, the multi-parameter approach offers greater flexibility in promoting sparsity in the regularized solutions while effectively mitigating the ill-posedness inherent in practical problems.

This paper aims to develop a practical parameter selection strategy for multi-parameter regularization with the $\ell_{1}$ norm, enabling the regularized solution to achieve a specified level of sparsity. To achieve this, we first establish a theoretical understanding of how the regularization parameters affect the sparsity of the solution under transform matrices. Unlike single-parameter regularization, the multi-parameter framework allows independent control over the sparsity of the solution corresponding to each transform matrix. Using tools from convex analysis, we characterize the relationship between the choice of multiple regularization parameters and the sparsity of the solution under the transform matrices. In the special case where the transform matrices are degenerate identities and the fidelity term exhibits block separability, this characterization provides a direct multi-parameter selection strategy that ensures each sub-vector of the solution achieves a desired level of block sparsity. In the general case, where such direct adoption is not feasible, we propose an iterative scheme to select the multiple regularization parameters. This scheme ensures the regularized solution attains the prescribed sparsity levels under the transform matrices. Unlike single-parameter regularization models, where only one parameter is optimized, our iterative algorithm considers the interplay between multiple parameters, thereby enhancing the flexibility and effectiveness of the parameter selection process.

In our recent papers [17] and [23], we investigated parameter selection strategies for sparse regularization using the $\ell_{1}$ norm. Compared to these studies, this paper makes two significant contributions: First, we have established strategies for selecting multiple regularization parameters in a model where the objective function comprises a convex fidelity term and multiple $\ell_{1}$-based regularization terms, each composed with a linear transform. Unlike the iterative algorithm in [17], which focused on selecting a single parameter, the algorithm proposed here addresses multiple parameters and explicitly incorporates their interdependencies. While [23] also explored iterative schemes for multi-parameter selection in a highly nonconvex optimization problem with numerous network parameters and extensive training data, the model in this paper facilitates a more sophisticated iterative algorithm. Second, we have developed iterative schemes for selecting multiple regularization parameters in two distinct scenarios: (1) when the fidelity term is differentiable, and the transform matrix has full row rank, and (2) when the fidelity term is nondifferentiable, and the transform matrix lacks full row rank. In contrast, the approaches in [17] and [23] address parameter selection only for the case of a differentiable fidelity term and a full-rank transform matrix.

For the nondifferentiable fidelity term and non-full-rank transform matrix, each iteration requires determining not only the regularized solution but also two auxiliary vectors. To address this, we have characterized these three components using a system of fixed-point equations involving proximity operators. Based on this characterization, we have devised a fixed-point proximity algorithm capable of simultaneously computing all three vectors. The convergence of this algorithm is rigorously analyzed in this paper. By combining the fixed-point proximity algorithm with the sparsity characterization of the regularized solution, we propose a robust multi-parameter selection strategy that ensures the solution achieves a specified level of sparsity.

This paper is organized into seven sections and two appendices. Section 2 introduces the general multi-parameter regularization problem under investigation and highlights several examples of practical significance. In Section 3, we establish a characterization of how each regularization parameter affects the sparsity of the regularized solution under its associated transform matrix. This characterization is further specialized in Section 4 to the case where the transform matrices are degenerate identities. Building on the insights from Sections 3 and 4, Section 5 presents iterative schemes for simultaneously determining multiple parameters and solutions with prescribed sparsity levels under the transform matrices. Two scenarios are addressed: (1) when the fidelity term is differentiable and the transform matrix has full row rank, and (2) when the fidelity term is nondifferentiable and the transform matrix does not have full row rank. Section 6 demonstrates the effectiveness of the proposed parameter selection strategies through a series of numerical examples. Section 7 concludes the paper by summarizing the findings. Appendix A contains the proof of the convergence analysis for the fixed-point proximity algorithm introduced in Section 5, while Appendix B provides closed-form expressions for the proximity operators used in the algorithm.

In this section, we introduce the multi-parameter regularization problem considered in this paper. Additionally, we review several optimization models of practical significance and illustrate how they can be formulated within this general framework.

We begin with describing the multi-parameter regularization problem. For each $s\in\mathbb{N}$, let $\mathbb{N}_{s}:=\{1,2,\ldots,s\}$ and set $\mathbb{N}_{0}:=\emptyset$. Suppose that $\bm{\psi}:\mathbb{R}^{n}\to\mathbb{R}_{+}:=[0,+\infty)$ is a convex function and for each $j\in\mathbb{N}_{d}$, $m_{j}\in\mathbb{N}$ and $\mathbf{B}_{j}$ is an $m_{j}\times n$ real matrix. For each $m\in\mathbb{N}$ and $\mathbf{x}:=[x_{j}:j\in\mathbb{N}_{m}]\in\mathbb{R}^{m}$, we define its $\ell_{1}$ norm by $\|\mathbf{x}\|_{1}:=\sum_{j\in\mathbb{N}_{m}}|x_{j}|$. We consider the multi-parameter regularization problem

where $\lambda_{j}>0,$ $j\in\mathbb{N}_{d},$ is a sequence of positive regularization parameters.

The multi-parameter regularization problem (1) appears in many application areas. Below, we present several examples of it. In imaging reconstruction problems, a combination of multiple regularizers was used to encourage their solution to simultaneously exhibit the characteristics enforced by each of them. For example, a combination of frame-based synthesis and analysis $\ell_{1}$ norm regularizers was proposed in [1] for an imaging deblurring problem. Specifically, for $m,n,s\in\mathbb{N}$, we assume that $\mathbf{A}\in\mathbb{R}^{m\times s}$ represents a periodic convolution, $\mathbf{W}\in\mathbb{R}^{s\times n}$ is a synthesis operator whose columns contain the elements of a frame, and $\mathbf{P}\in\mathbb{R}^{m\times s}$ is an analysis operator of a tight Parseval frame satisfying $\mathbf{P}^{\top}\mathbf{P}=\mathbf{I}_{s}$, where $\mathbf{I}_{s}$ denotes the identity matrix of order $s$. Let $\mathbf{y}\in\mathbb{R}^{m}$ be observed data. The regularization problem combining the synthesis and analysis $\ell_{1}$ norm regularizers has the form

Clearly, problem (2) may be identified as a special case of (1) with

As a generalization of the lasso regularized model [25], the fused lasso regularized model was proposed in [26] for problems with features that can be ordered in some meaningful way. Let $p,n\in\mathbb{N}$. Suppose that a prediction problem with $p$ cases has outcomes $y_{i}$, $i\in\mathbb{N}_{p}$ and features $x_{ij}$, $i\in\mathbb{N}_{p}$, $j\in\mathbb{N}_{n}$. Let $\mathbf{X}:=[x_{ij}:i\in\mathbb{N}_{p},j\in\mathbb{N}_{n}]$ be the $p\times n$ matrix of features and $\mathbf{y}:=[y_{i}:i\in\mathbb{N}_{p}]\in\mathbb{R}^{p}$, $\mathbf{u}\in\mathbb{R}^{n}$ be the vectors of outcomes and coefficients, respectively. The fused lasso regularized model is formulated as

where $\mathbf{D}:=[d_{ij}:i\in\mathbb{N}_{n-1},j\in\mathbb{N}_{n}]$ is the $(n-1)\times n$ first order difference matrix with $d_{ii}=-1,$ $d_{i,i+1}=1$ for $i\in\mathbb{N}_{n-1}$ and $0$ otherwise. By introducing

the fused lasso model (3) can be rewritten in the form of (1). By penalizing the $\ell_{1}$ norm of both the coefficients and their successive differences, the fused lasso regularized model encourages the sparsity of the coefficients and also the sparsity of their differences. As a special case, the fused lasso signal approximation [11] has the form (3) with the feature matrix $\mathbf{X}$ being the identity matrix.

Filtering noisy data was considered in [22] for the case where the underlying signal comprises a low-frequency component and a sparse or sparse-derivative component. Specifically, assume that the noisy data $y(t)$ can be modeled as

where $f$ is a low-pass signal, $u$ is a sparse and sparse-derivative signal and $\eta$ is stationary white Gaussian noise. Given noisy data of the form (4), one seeks the estimate of $f$ and $u$ individually. For this purpose, we first solve the compound sparse denoising problem

to obtain the estimate $\mathbf{u}^{*}$ of $u$. Here, $\mathbf{y}:=[y(t):t\in\mathbb{N}_{n}]$ and $\mathbf{H}$ is the high-pass filter matrix with the form $\mathbf{H}:=\mathbf{A}^{-1}\mathbf{C}$, where $\mathbf{A}$ and $\mathbf{C}$ are banded matrices. We then get the estimate $\mathbf{f}^{*}$ of $f$ as $\mathbf{f}^{*}:=(\tilde{\mathbf{I}}-\mathbf{A}^{-1}\mathbf{C})(\mathbf{y}-\mathbf{u}^{*})$. It is clear that the compound sparse denoising model (5) has the form (1) with

and $\mathbf{B}_{1}:=\mathbf{I}_{n}$, $\mathbf{B}_{2}:=\mathbf{D}$.

Using a technique similar to that used in the fused lasso regularized model, the fused SVM was proposed for classification of array-based comparative genomic hybridization (arrayCGH) data [21, 27]. Given training data $\{(\mathbf{x}_{j},y_{j}):j\in\mathbb{N}_{p}\}$ composed of sample points $\{\mathbf{x}_{j}:j\in\mathbb{N}_{p}\}\subset\mathbb{R}^{n}$ and labels $\{y_{j}:j\in\mathbb{N}_{p}\}\subset\{-1,1\}$. The aim of binary classification is to find a decision function $f(\mathbf{x}):=\mathbf{u}^{\top}\mathbf{x}$, $\mathbf{x}\in\mathbb{R}^{n}$ predicting the class $y_{j}=-1$ or $y_{j}=1$. The class prediction for a profile $\mathbf{x}$ is then $1$ if $f(\mathbf{x})\geq 0$ and $-1$ otherwise. The fused SVM based on the hinge loss function has the form

We rewrite model (7) in the form (1) as follows. We define the matrix $\mathbf{X}:=[\mathbf{x}_{j}:j\in\mathbb{N}_{p}]^{\top}$, the matrix $\mathbf{Y}:=\mathrm{diag}(y_{j}:j\in\mathbb{N}_{p})$ and the function $\bm{\phi}(\mathbf{z}):=\sum_{j\in\mathbb{N}_{p}}\mathrm{max}\{0,1-z_{j}\}$, for all $\mathbf{z}:=[z_{j}:j\in\mathbb{N}_{p}]\in\mathbb{R}^{p}$. Then by introducing $\bm{\psi}(\mathbf{u}):=\bm{\phi}(\mathbf{YXu})$, for $\mathbf{u}\in\mathbb{R}^{n}$, $\mathbf{B}_{1}:=\mathbf{I}_{n}$, and $\mathbf{B}_{2}:=\mathbf{D}$, the fused SVM model (7) can be represented in the form (1).

In this section, we characterize the relation between the multiple regularization parameters and the sparsity of the regularized solution under the transform matrices $\mathbf{B}_{j}$, $j\in\mathbb{N}_{d}$, of the regularization problem (1). Unlike the single-parameter regularization problem, the use of multiple regularization parameters in problem (1) allows us to separately consider the sparsity of the solution under each transform matrix $\mathbf{B}_{j}$.

We begin by recalling the definition of the level of sparsity for a vector in $\mathbb{R}^{s}$. For $s\in\mathbb{N}$, we set $\mathbb{Z}_{s}:=\{0,1,\ldots,s-1\}$. A vector $\mathbf{x}\in\mathbb{R}^{s}$ is said to have sparsity of level $l\in\mathbb{Z}_{s+1}$ if it has exactly $l$ nonzero components.

We now reformulate the regularization problem (1) into an equivalent form to facilitate the characterization of the sparsity of the solution under the transform matrices $\mathbf{B}_{j}$, $j\in\mathbb{N}_{d}$. Let $p_{0}:=0$ and $p_{j}:=\sum_{i\in\mathbb{N}_{j}}m_{i}$, $j\in\mathbb{N}_{d}$. We decompose a vector $\mathbf{z}:=[z_{k}:k\in\mathbb{N}_{p_{d}}]\in\mathbb{R}^{p_{d}}$ into $d$ sub-vectors by setting

By introducing a column block matrix

we write $\mathbf{B}\mathbf{u}$ for the block column vector $[\mathbf{B}_{j}\mathbf{u}:j\in\mathbb{N}_{d}]$. Thus, we may rewrite the regularization problem (1) as

We further convert problem (9) into an equivalent form by utilizing a change of variables. To achieve this, we consider inverting the linear system

Here, $\mathcal{R}(\mathbf{B})$ denotes the range of $\mathbf{B}$. It is known from [5, 12] that the general solution of the linear system (10) can be represented by the pseudoinverse of $\mathbf{B}$. An alternative form of the general solution was provided in [17]. To describe this result, we recall that if $\mathbf{B}$ has the rank $r$ satisfying $0<r\leq\mathrm{min}\{p_{d},n\}$, then $\mathbf{B}$ has the SVD as $\mathbf{B}=\mathbf{U}\mathbf{\Lambda}\mathbf{V}^{\top}$, where $\mathbf{U}$ and $\mathbf{V}$ are $p_{d}\times p_{d}$ and $n\times n$ orthogonal matrices, respectively, and $\mathbf{\Lambda}$ is a $p_{d}\times n$ diagonal matrix with the nonzero diagonal entries $\sigma_{1}\geq\cdots\geq\sigma_{r}>0$, which are the nonzero singular values of $\mathbf{B}$. In order to represent the general solution of linear system (10), we define an $n\times(p_{d}+n-r)$ matrix by employing the SVD of $\mathbf{B}$. Specifically, we denote by $\widetilde{\mathbf{U}}_{r}\in\mathbb{R}^{p_{d}\times r}$ the matrix composed of the first $r$ columns of $\mathbf{U}$ and define an $n\times(p_{d}+n-r)$ block diagonal matrix by setting

We also introduce a diagonal matrix of order $n$ by

Using these matrices, we define an $n\times(p_{d}+n-r)$ matrix by $\mathbf{B}^{\prime}:=\mathbf{V}\mathbf{\Lambda}^{\prime}\mathbf{U}^{\prime}.$ As has been shown in [17], for each solution $\mathbf{u}$ of system (10), there exists a unique vector $\mathbf{v}\in\mathbb{R}^{n-r}$ such that

As a result, the mapping $\mathcal{B}$, defined for each $\mathbf{u}\in\mathbb{R}^{n}$ by $\mathcal{B}\mathbf{u}:=\scriptsize{\begin{bmatrix}\mathbf{z}\\ \mathbf{v}\end{bmatrix}},$ where $\mathbf{z}:=\mathbf{B}\mathbf{u}$ and $\mathbf{v}\in\mathbb{R}^{n-r}$ satisfies equation (11), is bijective from $\mathbb{R}^{n}$ onto $\mathcal{R}(\mathbf{B})\times\mathbb{R}^{n-r}$.

By making use of the change of variables defined by equation (11), we reformulate problem (9) as an equivalent multi-parameter regularization problem with for each $j\in\mathbb{N}_{d}$, $\mathbf{B}_{j}$ being a degenerated identity. We set $\mathbb{M}:=\mathcal{R}(\mathbf{B})\times\mathbb{R}^{n-r}$ and let $\iota_{\mathbb{M}}:\mathbb{R}^{p_{d}+n-r}\rightarrow\overline{\mathbb{R}}:=\mathbb{R}\cup\{+\infty\}$ denote the indicator function of $\mathbb{M}$, that is, $\iota_{\mathbb{M}}(\mathbf{w})=0$ if $\mathbf{w}\in\mathbb{M}$, and $+\infty$ otherwise. For each $j\in\mathbb{N}_{d}$, we introduce a degenerated identity matrix by

where $\mathbf{0}_{s\times t}$ denotes the zero matrix of order $s\times t$. We show in the following lemma that the regularization problem (1) is equivalent to the regularization problem

Below, we consider how the regularization parameters $\lambda_{j}$, $j\in\mathbb{N}_{d}$ influence the sparsity of the solution of problem (13). To characterize the solution of problem (13), we need the notion of the subdifferential of a convex function. Suppose that $f:\mathbb{R}^{s}\to\overline{\mathbb{R}}$ is a proper convex function. The subdifferential of $f$ at $\mathbf{x}\in{\rm dom}(f):=\{\mathbf{y}\in\mathbb{R}^{s}:f(\mathbf{y})<+\infty\}$ is defined by

It is known [30] that for two convex functions $f$ and $g$ on $\mathbb{R}^{s}$, if $g$ is continuous on $\mathbb{R}^{s}$ then $\partial(f+g)(\mathbf{x})=\partial f(\mathbf{x})+\partial g(\mathbf{x})$, for all $\mathbf{x}\in\mathbb{R}^{s}$. We also describe the chain rule of the subdifferential [24]. Suppose that $f:\mathbb{R}^{s}\rightarrow\overline{\mathbb{R}}$ is a convex function and $\mathbf{M}$ is an $s\times t$ matrix. If $f$ is continuous at some point of the range of $\mathbf{M}$, then for all $\mathbf{x}\in\mathbb{R}^{t}$

The Fermat rule [30] states that if a proper convex function $f:\mathbb{R}^{s}\rightarrow\overline{\mathbb{R}}$ has a minimum at $\mathbf{x}\in\mathbb{R}^{s}$ if and only if $\mathbf{0}\in\partial f(\mathbf{x})$.

For the purpose of characterizing sparsity of vectors in $\mathbb{R}^{s}$, we also recall the notion of the sparsity partition of $\mathbb{R}^{s}$ introduced in [29]. Specifically, by using the canonical basis $\mathbf{e}_{s,j},j\in\mathbb{N}_{s},$ for $\mathbb{R}^{s}$, we introduce $s+1$ numbers of subsets of $\mathbb{R}^{s}$ by

It is clear that the sets $\Omega_{s,l},l\in\mathbb{Z}_{s+1}$, form a partition for $\mathbb{R}^{s}$ and for each $l\in\mathbb{Z}_{s+1}$, $\Omega_{s,l}$ coincides with the set of all vectors in $\mathbb{R}^{s}$ having sparsity of level $l$.

We provide in the next lemma a characterization of the sparsity of the solution of problem (13). For an $s\times t$ matrix $\mathbf{M}$, we denote by $\mathcal{N}(\mathbf{M})$ the null space of $\mathbf{M}$ and for each $i\in\mathbb{N}_{t}$, denote by $\mathbf{M}_{(i)}$ the $i$th column of $\mathbf{M}$.

Combining Lemmas 3.1 with 3.2, we establish the relation between the regularization parameters and the sparsity of the solution of problem (1).

Theorem 3.3 provides a characterization of the multiple regularization parameter $\lambda_{j}$, $j\in\mathbb{N}_{d}$, with which problem (1) has a solution with sparsity of a certain level $l_{j}$ under each transform matrix $\mathbf{B}_{j}$. By specifying the sparsity level of the solution under each transform matrix $\mathbf{B}_{j}$ to be $l_{j}^{*}$, our goal is to find regularization parameters $\lambda_{j}^{*}$, $j\in\mathbb{N}_{d}$, satisfying conditions (16), (17) and (18). However, since these conditions depend on the corresponding solution, the characterization stated in Theorem 3.3 can not be used directly as a multi-parameter choice strategy. Motivated by Theorem 3.3, an iterative scheme to be developed in section 5 will enable us to choose multiple regularization parameters with which a minimizer of problem (1) has a prescribed sparsity level under each transform matrix.

The next result concerns the special case that matrix $\mathbf{B}$ defined by (8) has full row rank, that is, $\mathrm{rank}(\mathbf{B})=p_{d}$.

Note that if $\bm{\psi}$ is differentiable, then the subdifferential of $\bm{\psi}$ at $\mathbf{u}^{*}$ is the singleton $\nabla\bm{\psi}(\mathbf{u}^{*})$. In this case, Theorem 3.3 and Corollary 3.4 have the following simple form.

In this section, we consider the special case where $p_{d}=n$ and for each $j\in\mathbb{N}_{d},$ the transform matrix $\mathbf{B}_{j}$ takes the form

In this scenario, the multi-parameter regularization problem (1) assumes the special form

We specialize the characterizations of the sparsity of regularized solutions established in the previous section to this special case. Moreover, particular attention is given to scenarios where the fidelity term $\bm{\psi}$ is block-separable.

We first characterize the sparsity of the solution of problem (25). It follows from equation (24) that matrix $\mathbf{B}$ defined by (8) coincides with the identity matrix of order $n$. Given that matrix $\mathbf{B}$ has full row rank, we specialize Corollary 3.4 to the regularization problem (25). In addition, the transform matrices $\mathbf{B}_{j}$, $j\in\mathbb{N}_{d}$, with the form (24) enable us to consider the sparsity of each sub-vector of the regularized solution separately.