# Linearity Identification for General Partial Linear Single-Index Models

Authors: Shaogao Lv; Luhong Wang
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Partial linear models, a family of popular semiparametric models, provide us with an interpretable and flexible assumption for modelling complex data. One challenging question in partial linear models is the structure identification for the linear components and the nonlinear components, especially for high dimensional data. This paper considers the structure identification problem in the general partial linear single-index models, where the link function is unknown. We propose two penalized methods based on a modern dimension reduction technique. Under certain regularity conditions, we show that the second estimator is able to identify the underlying true model structure correctly. The convergence rate of the new estimator is established as well.


---

## Body

## 1. Introduction

Partially linear models (PLMs), containing both linear and nonlinear additive components, provide a useful class of tools for modelling complex data. PLMs have wide applications in practice due to their flexibility and interpretability.Given the data set{
y
i
,
x
i
(
1
)
,
x
i
(
2
)
}
i
=
1
n, where y
i is the response and x
i
(
1
)
=
(
x
i
1
,
…
,
x
i
p
) and x
i
(
2
)
=
(
x
i
(
p
+
1
)
,
…
,
x
i
(
p
+
q
)
) are vectors of covariates, the basic PLM has the following form: (1)
y
i
=
b
+
β
o
′
x
i
1
+
∑
j
=
1
q
f
o
,
j
x
i
p
+
j
+
ε
i
,where b is the intercept, β
o is a vector of unknown parameters for linear terms, each f
o
,
j is an unknown nonlinear function from R
q to R, and ε
i’s are i.i.d. random errors with zero conditional mean and a finite variance.Given that the linear term and nonlinear term are known in advance, estimation and inference for various PLMs have been well studied in literature, such as [1–3]. Under the same settings, motivated by sparse models and the Lasso idea developed rapidly in recent years, there exists a substantial work on variable selection and estimation for sparse PLMs, including [4–6] and others. However, faced with complicated data-generating processes, a key and challenging question to use PLM is how to decide which features should be assigned to be linear and which are nonlinear, especially for high dimensional problems such as text categorization and biology. Furthermore, the difficulty level of model search increases dramatically as the data dimension grows due to the curse of dimensionality. In the last several years, a number of methods have been proposed to address various aspects of this problem. Two commonly used methods are the screening and hypothesis testing procedures. The screening method is sometimes useful in practice but lacks theoretical justifications, while, for hypothesis testing, it is often hard to construct proper test statistics and the tests may have a poor performance when the number of covariates is large. To handle estimation and component identification problems efficiently for PLM, Zhang et al. [7] originally proposed a penalized method called LAND for automatically identifying linear and nonlinear terms for PLM. This method is based on the orthogonal decomposition of Sobolev space of order 2, consisting of linear subspace and nonlinear subspaces. Another main stream of research is based on spline approximation for nonparametric functions and penalized likelihood. For example, Huang et al. [8] start with a nonparametric additive model but use spline series expansions to approximate the component functions. This spline series approximation approach allows them to prove selection consistency for more general functions. Lian et al. [9] proposed a doubly SCAD-penalized approach for partial linear structure identification problems of nonpolynomial (NP) dimensionality. These above methods can be easily extended to be generalized partial linear models; that is, the added link function is known in advance. If the link function is unknown, a standard technique in the literature is the two-step estimation. Firstly given that the component functions are fixed, a rough estimator of the link function is obtained; then the derived link function is fixed, estimating the component functions again. Repeat the process until convergence. However, this method is unstable and quite sensitive to the choice of penalty parameters. Moreover, this method involves a nonconvex optimization, which can not guarantee satisfactory numerical results.In this paper we consider a general single-index models with the partial linear form(2)
y
i
=
G
β
o
′
x
i
1
+
∑
j
=
1
q
f
o
,
j
x
i
p
+
j
,
ε
i
,where the error term ε
i is assumed to be independent of covariates, besides satisfying the standard moment conditions mentioned above. The link function G
(
·
) is typically unknown. Note that additive structure of the unknown link function and the error is not assumed in model (2). This model is quite general; it includes the classical generalized partial linear model with the form y
i
=
G
1
(
β
o
′
x
i
(
1
)
+
∑
j
=
1
q
f
o
,
j
(
x
i
(
p
+
j
)
)
)
+
ε
i, as well as the partial linear transformation model G
2
(
y
i
)
=
β
o
′
x
i
(
1
)
+
∑
j
=
1
q
f
o
,
j
(
t
i
(
p
+
j
)
)
+
ε
i. Clearly, when the link function G
(
·
) is not specified, the partial linear term is identifiable only up to a multiplicative scalar because any location-scale change in {
β
o
′
x
i
(
1
)
+
∑
j
=
1
q
f
o
,
j
(
x
i
(
p
+
j
)
)
} can be absorbed into the link function.By making use of the cubic spline approximation to smooth functions and the regularization technique, we develop a new penalized method to identify model structure for the general partial linear models (2). The proposed method is formulated by the theory of sufficient dimension reduction and two different variant forms of eigenvalues problems. We specially use a modern sufficient dimension technique for estimating the linear and nonlinear terms, called gradient-based kernel dimension reduction (gKDR), proposed by Fukumizu and Leng [10]. This method shows some specific advantages compared with the above-mentioned techniques; for example, the gKDR method can handle any type of variable for y including multivariate or nonvectorial one in the same way. Besides, the linear condition and constant variance conditions are avoided when the gKDR is adopted.Based on the gKDR and two variant forms for the eigenvalues problem, we design two different penalized approaches by adding twofold Lasso-type penalties. The first algorithm is a nonconvex optimization problem with equality constraint, which can be solved efficiently by analternating direction method of multipliers (ADMM). To provide better statistical properties, a convex relaxation of the eigenvalues problem is used to formulate our other method, a semidefinite programme. Using this convex approach, we can establish the estimation consistency and support recovery for the general partial linear models. The convex optimization problem can also be solved in a polynomial time algorithm by the ADMM.The rest of the article is organized as follows. In Section2 we transform model (2) into a surrogate model by the cubic spline approximation. Next we introduce the gKDR in the literature of dimension reduction. Then we propose two different penalized approaches based on two variants of eigenvalues problems. Statistical properties of the new estimators, including its convergence rate and selection consistency, are established in Section 3. All the proofs are relegated to the Appendix.


## 2. Problem Formulation

Notation. We collect here some standard notation used throughout the paper. For matrices A
,
B with the same dimension A
,
B
≔
tr
⁡
(
A
T
B
) is the Frobenius inner product, and A
F
≔
〈
A
,
A
〉 is the square Frobenius norm. x
q is the usual norm in the Euclidean spaces. A
a
,
b is (
a
,
b
)-norm defined to be l
b norm of the vector of rowwise l
a norms of A. We denote A
⪯
B, meaning that B
-
A is semidefinite positive.We mainly consider the problem of linear identification of the general single-index model (2) and estimate relevant variables to the response. For nonparametric components in model (2), one often uses spline or wavelet approximation as finite representations, because of their good approximation capabilities to functions in various nonparametric classes and computational advantages. In our specific setting, we adopt the cubic spline approximation, because two important properties of cubic spline are critical for our analysis. One property is that all the its bases are linearly independent of each other. Another property is that the function x is contained in the cubic spline bases. We refer the reader to [11] for a detailed description for spline functions. Thus, benefited from the two properties of cubic spline, each function of this spline space has a unique finite representation, and the coefficient of the base function x corresponds to the linear part.Letx
=
(
x
(
1
)
,
x
(
2
)
) be the total covariate defined on X
⊂
R
p
+
q, and Y
∈
R is the response. For each covariate x
j
,
j
=
1
,
…
,
p
+
q, let ϕ
(
x
j
)
=
[
ϕ
1
(
x
j
)
,
…
,
ϕ
M
(
x
j
)
] be an M-dimensional basis expansion, where ϕ
1
(
x
j
)
=
x
j. For any vector β
n
j
=
(
β
j
1
,
…
,
β
j
M
), we define f
j
(
x
)
=
β
n
j
′
ϕ
(
x
) as cubic spline-based function we consider. With this spline finite representation, each f
o
,
j in our model (2) can be approximated well by this form of f
j above. Let Φ
n
(
x
)
=
(
ϕ
(
x
1
)
,
…
,
ϕ
(
x
p
+
q
)
) be M
×
(
p
+
q
) vector-vauled base function for any x
=
(
x
1
,
…
,
x
p
+
q
)
∈
X. To represent our model by a unified form, we denote β
s
j
=
(
β
o
,
j
,
0
,
…
,
0
) as the coefficient representation of the jth linear-form variable. Let β
p
=
β
s
1
′
,
…
,
β
s
p
′
′ and β
q
=
β
o
,
1
′
,
β
o
,
q
′
′, where β
o
,
j
(
j
=
1
,
…
,
q
) is the coefficient corresponding to the jth true nonlinear component. Thus, let θ
o
=
β
p
′
,
β
q
′
′ be M
×
(
p
+
q
) approximating coefficient, our focus on the original model (2) is now transformed to the following one: (3)
y
i
=
G
θ
o
′
Φ
n
x
i
,
ε
i
.This semiparametric model (3) indicates essentially that the response y depends solely on the generalized predictor vector Φ
n through a linear combination θ
o
′
Φ
n, or, equivalently, y is independent of the original covariate x when θ
o
′
Φ
n is given. In the literature, the theory of sufficient dimension reduction provides an effective starting point to estimate θ
o
′ without loss of regression information of y on x and without assuming the specific link function. The popular methods for sufficient dimension reduction include sliced inversion regression [12], sliced average variance estimation [13], contour regression, directional regression [14], and references therein.In this article, we use the gKDR formulated by the reproducing kernel techniques. This method shows some specific advantages compared with the above-mentioned techniques; for example, the gKDR method can handle any type of variable fory including multivariate or nonvectorial one in the same way. Besides, the nonparametric nature of the kernel method avoids making strong assumptions on the distribution of covariates, the response, or the conditional probability, which are required usually in many classical dimension reduction methods such as SIR, pHd, and contour regression.


### 2.1. gKDR Method

To give the gKDR method, we need some technical notations. For a compact setΩ in some metric space, we denote apositive-definite kernel on Ω
×
Ω, that is, a symmetric function k
:
Ω
×
Ω
→
R satisfying the following semidefinite positive condition: ∑
i
,
j
=
1
n
c
i
c
j
k
(
x
i
,
x
j
)
≥
0 for any x
1
,
…
,
x
n in Ω and c
1
,
…
,
c
n
∈
R. It is known [15] that a positive-definite kernel on Ω is uniquely associated with a Hilbert space H consisting of functions on Ω, and H is called areproducing kernel Hilbert space (RKHS) associated with the kernel k.For short, we denoteZ
≔
{
Φ
n
x
,
x
∈
X
} as a subset of R
M
(
p
+
q
). Let (
Z
,
B
Z
,
μ
Z
) and (
Y
,
B
Y
,
μ
Y
) be measurable spaces and (
Z
,
Y
) be a random variable on Z
×
Y with probability P. Let k
Z and k
Y be positive-definite kernels on Z and Y, respectively, and the corresponding RKHSs denoted by H
Z and H
Y. It is assumed that E
(
k
Z
(
Z
,
Z
)
) and E
(
k
Y
(
Y
,
Y
)
) are finite. Now we introduce two integral operators, so that an appropriate variant of the problem of estimating θ
o in (3) can be established. Thecross-covariance operator
C
Y
Z
:
H
Z
→
H
Y is defined as the operator as follows: (4)
C
Y
Z
f
y
=
∫
k
y
,
y
~
f
z
~
d
P
z
~
,
y
~
,
∀
f
∈
H
X
.Similarly, C
Z
Z denotes theself-covariance operator on H
Z; that is, (5)
C
Z
Z
f
z
=
∫
k
z
,
z
~
f
z
~
d
P
Z
z
~
,
∀
f
∈
H
Z
.Remark that these definitions are natural extensions of the ordinary covariance metrics on Euclidean spaces. Although C
Y
Z and C
Z
Z depend on the kernels, we omit the dependence for notational simplicity.With these definitions, we now define the populationM
(
p
+
q
)
×
M
(
p
+
q
) matrix M
(
z
)
=
(
M
i
j
(
z
)
) by (6)
M
i
j
z
=
C
Y
Z
C
X
X
-
1
∂
k
Z
·
,
z
∂
z
i
,
C
Y
Z
C
Z
Z
-
1
∂
k
Z
·
,
z
∂
z
j
H
Y
,
i
,
j
=
1
,
…
,
M
p
+
q
.Reference [10] has shown that the eigenvector corresponding to the largest eigenvalue of the M
(
z
) is equal to the true vector θ
o for any z, provided that the maximal eigenvalue is simple. To be precise, the following equality holds: (7)
θ
o
=
v
→
m
a
x
E
M
Z
,where v
→
m
a
x
(
A
) is denoted to be the maximum eigenvector of the matrix A. Thus, the problem of estimating θ
o of (3) is equivalent to the eigenvalue problem of solving M
(
z
).We now turn to the empirical version ofM
(
z
), so as to estimate θ
o based on available sample. Denote n
×
n Gram matrices (
k
Z
(
z
i
,
z
j
)
) and (
k
Y
(
y
i
,
y
j
)
) by G
Z and G
Y, respectively. Let ∇
k
Z
(
z
)
=
(
∂
k
Z
(
z
1
,
z
)
/
∂
z
,
…
,
∂
k
Z
(
z
n
,
z
)
/
∂
z
)
∈
R
n
×
M
(
q
+
p
). Then, the gKDR proposes the eigenvectors of M
(
p
+
q
)
×
M
(
p
+
q
) symmetric matrix: (8)
M
n
z
≔
∇
k
Z
z
′
G
Z
+
n
ϵ
n
I
n
-
1
G
Y
G
Z
+
n
ϵ
n
I
n
-
1
·
∇
k
Z
z
,where ϵ
n is a regularization parameter in Thikonov-type regularization. Then, we define (9)
M
^
n
≔
1
n
∑
i
=
1
n
M
n
z
ias an empirical estimator of θ
o. In the next two subsections, we introduce two variational forms with respect to the eigenvalue problems, so that our proposed methods based on Lasso idea can be formulated naturally.


### 2.2. A Nonconvex Variational Algorithm

In order to employ the popular Lasso idea, we need to transform the eigenvalue problem (7) to an optimization scheme. Recall the followingCourant-Fischer variational representation of the maximal eigenvalue and eigenvector: (10)
m
a
x
θ
∈
R
M
p
+
q
θ
′
E
M
Z
θ
,
s
.
t
.
θ
2
2
=
1
.To identify the linear and nonlinear terms in model (3), we define an estimator by adding twofold convex penalties to (10) such that some sparse estimators might be generated. To be precise, let θ
=
(
θ
1
,
…
,
θ
p
+
q
), where each θ
j
∈
R
M. The penalized estimation method is defined as follows: (11)
θ
^
n
=
a
r
g
m
i
n
θ
∈
R
M
p
+
q
-
θ
′
M
^
n
θ
+
λ
∑
j
=
1
p
+
q
α
θ
j
,
1
2
+
1
-
α
θ
j
,
-
1
2
,
s
.
t
.
θ
2
2
=
1
,where θ
j
,
1 is the first element of θ
j and θ
j
,
-
1 is the last M
-
1 component of θ
j
,
j
=
1
,
…
,
p
+
q. Note that λ is the penalty parameter and α is a tuning parameter, typically, α
=
M
/
(
1
+
M
) for normalization. Our penalized method is motivated by the following fact. For any given function f
(
x
)
=
∑
k
=
1
M
β
k
ϕ
k
(
x
), by the cubic spline representation, as mentioned earlier, the first spline base is the single function x and all the bases are linearly independent of each other. Consequently, β
1
=
0 and β
k
≠
0 with k
>
1 if and only if f is a nonlinear function. Based on this and the lasso constraint in (11), the proposed method can identify both the linear and nonlinear terms, of course, also finding the relevant feature simultaneously.We observe that optimization of  (11) is complicated, since it belongs to a class of constrained nonconvex programmes. In this article, we consider an alternative formulation of (11); precisely, we use the ADMM [16]. ADMM is an algorithm that is intended to blend the decomposability of dual ascent with the superior convergence properties of the method of multipliers. Let ρ
(
θ
)
=
∑
j
=
1
p
+
q
(
α
θ
j
,
1
2
+
1
-
α
θ
j
,
-
1
2
), and by introducing an additional constraint, we can reformulate algorithm (11) as (12)
m
i
n
θ
,
γ
∈
R
M
p
+
q
∞
×
Π
2
γ
-
θ
′
M
^
n
θ
+
λ
ρ
θ
,
s
.
t
.
θ
=
γ
,where Π
2 is the 0
-
1 indicator function for the unit l
2-ball surface in R
M
(
p
+
q
) and we adopt the convention ∞
×
0
=
0. As in the method of multipliers, we form the augmented Lagrangian with a parameter τ: (13)
L
γ
,
θ
,
ν
,
τ
=
∞
×
Π
1
γ
-
θ
′
M
^
n
θ
+
λ
ρ
θ
+
ν
′
θ
-
γ
+
τ
2
θ
-
γ
2
2
,where (
ν
,
τ
) is the Lagrangian multipliers. Defining the residual r
=
θ
-
γ, we have (14)
L
γ
,
θ
,
ν
,
τ
=
∞
×
Π
1
γ
-
θ
′
M
^
n
θ
+
λ
ρ
θ
+
τ
2
r
+
u
2
2
-
τ
2
u
2
2
,where u
=
ν
/
τ is called the scaled dual variable. Using the scaled dual variable, we can express ADMM as (15)
θ
k
+
1
=
arg
⁡
min
θ
⁡
λ
ρ
θ
-
θ
′
M
^
n
θ
+
τ
2
θ
-
γ
k
+
u
k
2
2
,
(16)
γ
k
+
1
=
Π
B
2
θ
k
+
1
+
u
k
,
(17)
u
k
+
1
=
u
k
+
θ
k
+
1
-
γ
k
+
1
,where Π
B
2 is Euclidean projection onto the unit l
2-ball surface.There are many convergence results for ADMM discussed in the literature. ADMM consists of iteratively minimizingL
(
γ
,
θ
,
ν
,
τ
) with respect to γ, minimizing L
(
γ
,
θ
,
ν
,
τ
) with respect to θ, and then updating the dual variable u. The θ-update involves iterative evaluation of the proximal operator associated with (18)
F
θ
≔
-
θ
′
M
^
n
θ
+
τ
2
θ
-
γ
k
+
u
k
2
2
.To this end, we use theblock coordinate descent (BCD) [17] for l
1
/
l
2-regularization to compute θ
k
+
1 of (15). Specially at tth iteration, the BCD solves a problem of the form (19)
m
i
n
θ
g
∇
g
F
θ
t
k
′
θ
g
-
θ
g
,
t
k
+
1
2
θ
g
-
θ
g
,
t
k
′
H
g
g
θ
g
-
θ
g
,
t
k
+
μ
θ
g
2
,where θ
g is the subelement of θ with group g
∈
G forming a partition of {
1
,
…
,
M
(
p
+
q
)
} and H
g
g equals or approximates ∇
g
g
2
F
(
θ
t
k
). Note that μ is taken as λ
α or λ
(
1
-
α
) according to the specified group. Obviously, the BCD can be viewed as a generalization of the classical coordinate descent for l
1-regularization. In our specific setting, a direct computation on F
(
θ
) yields that (20)
∇
g
F
θ
=
-
2
θ
′
M
^
n
,
g
1
,
…
,
θ
′
M
^
n
,
g
k
+
τ
θ
+
u
k
-
γ
k
g
,
(21)
∇
2
F
θ
g
g
=
-
2
M
^
n
g
g
+
τ
I
g
,where we write g
=
(
g
1
,
…
,
g
k
) and M
^
n
,
g
l is the lth column of M
^
n
,
l
=
1
,
…
,
k. I
g
g is the identity and (
M
^
n
)
g
g is the submatrix of M
^
n with the index g. Since (
M
^
n
)
g
g in (21) is usually not a multiple of the identity, the solution minimizing (19) has no closed form. We here use a simple but efficient replacement; that is, H
g
g
=
h
g
g
I
|
g
|, where h
g
g is the largest eigenvalue of -
2
(
M
^
n
)
g
g
+
τ
I
|
g
|, which can be computed efficiently, since -
2
(
M
^
n
)
g
g
+
τ
I
|
g
| is only a submatrix indexed by g and independent of θ. In this case, the iterative solution θ
t
+
1
k of (19) is obtained bygroup soft thresholding of the Newton step: (22)
θ
g
,
t
+
1
k
=
P
r
o
x
μ
·
2
θ
g
,
t
k
-
h
g
g
-
1
∇
g
F
θ
t
k
,where P
r
o
x
μ
·
2 is called the soft thresholding operator associated with ·
2, defined as P
r
o
x
μ
·
2
(
θ
g
)
=
1
-
μ
/
θ
g
2
+
θ
g for any θ
g. Under mild conditions, θ
t
+
1
k converges to some vector as t goes to infinity, denoted by θ
k
+
1. It is worthwhile to note that, by (20), we find that the iteration in (22) can be computed parallel between groups.With regard to projection computation of (16), it is known that γ
k
+
1
=
θ
k
+
1
+
u
k
/
θ
k
+
1
+
u
k
2 if θ
k
+
1
+
u
k
2
≠
0; otherwise, γ
k
+
1 is set to any given value on l
2-ball surface. Thus, An efficient implementation for ADMM is obtained; to be precise, the overall procedure can be stated as in Algorithm 1.Algorithm 1:ADMM algorithm.
given: parameters λ, α, ϵ
n,  τ, the group indexes g and the kernels k
Z, k
Y.
initialize: θ
0
=
0, γ
0
=
0, u
0
=
0,
for  k
≥
0
repeat
for  t
≥
0  repeat
θ
g
,
t
+
1
k
=
P
r
o
x
μ
⋅
2
(
θ
g
,
t
k
-
h
g
g
-
1
∇
g
F
(
θ
t
k
)
), for each group g,
t
←
t
+
1
until  θ
g
,
t
k converges to some θ
g
∗
θ
g
k
+
1
←
θ
g
∗, then compute γ
k
+
1
=
Π
B
2
(
θ
k
+
1
+
u
k
),
andu
k
+
1
=
u
k
+
θ
k
+
1
-
γ
k
+
1
thenk
←
k
+
1
until
θ
k
+
1 converges to θ
^
n.
then return
(
θ
^
n
)


### 2.3. A Semidefinite Programming Relaxation

This section intrduces a convex relaxation of (10). Since estimate of individual eigenvectors of θ
o in (7) is unstable if the gap between their eigenvalues is small, it seems reasonable to instead focus on their span, that is, the principal subspace of variation. Let Λ
o
=
θ
o
θ
o
′ as the projection onto the subspace spanned by θ
o; a lesser known but equivalent variational representation of (10) is of the convex program: (23)
m
a
x
Θ
∈
S
+
M
p
+
q
,
tr
⁡
Θ
=
1
tr
⁡
E
M
Z
Θ
,where S
+
M
(
p
+
q
)
=
{
Θ
∈
R
M
(
p
+
q
)
×
M
(
p
+
q
)
∣
0
⪯
Θ
⪯
I
}. The optimization problem (23) is a semidefinite program (SDP), which can be solved exactly in polynomial time. If the maximal eigenvalue is simple, the optimum always achieved at a rank-one matrix, denoted by Θ
-, corresponds to the maximal eigenvector; that is, Θ
-
=
θ
o
θ
o
′. Motivated by this nice conclusion, we can estimate the sparse vector θ
o by adding l
1-type regularization, (24)
Θ
^
≔
arg
⁡
m
a
x
Θ
∈
S
+
M
p
+
q
tr
⁡
M
^
n
Θ
-
λ
Θ
1,1
,
s
.
t
.
tr
⁡
Θ
=
1
,and computing the maximal eigenvector θ
^
=
v
→
m
a
x
(
Θ
^
).Using the ADMM again, we rewrite (11) as the equivalent problem with equality constrained (25)
m
i
n
∞
×
1
F
1
Θ
-
M
^
n
,
Θ
+
λ
Δ
1,1
s
.
t
.
Θ
-
Δ
=
0
,where F
1
=
{
Θ
∈
S
M
p
+
q
,
tr
⁡
Θ
=
1
} and 1
F
1 is the 0
-
1 indicator function for F
1. The augmented Lagrangian associated with (25) has the form (26)
L
τ
Θ
,
Δ
,
U
≔
∞
×
1
F
1
Θ
-
M
^
n
,
Θ
+
λ
Δ
1,1
+
τ
2
Θ
-
Δ
+
U
F
2
-
U
F
2
.
Θ and Δ updates mainly involve computing the projection operator and the proximal operator, respectively; that is, (27)
Π
F
1
Δ
-
U
+
M
^
n
τ
≔
arg
⁡
min
Θ
∈
F
1
⁡
1
2
Θ
-
Δ
-
U
+
M
^
n
τ
F
2
,
S
λ
/
τ
Θ
+
U
≔
arg
⁡
min
Δ
⁡
λ
τ
Δ
1,1
+
1
2
Θ
+
U
-
Δ
F
2
,where S
λ
/
τ is the elementwise soft thresholding operator defined as (28)
S
λ
/
τ
x
=
sign
⁡
x
max
⁡
x
-
λ
τ
,
0
.On the other hand, Π
F
1 is the Euclidean projection onto F
1 and has a closed form as follows.Lemma 1 (Fantope projection).
IfU
=
∑
j
γ
j
u
j
u
j
′ is a spectral decomposition of U, then Π
F
1
(
U
)
=
∑
j
γ
j
+
θ
u
j
u
j
′, where γ
j
+
θ
=
min
⁡
max
⁡
γ
j
-
θ
,
0
,
1 and θ satisfies the equation ∑
j
γ
j
+
(
θ
)
=
1.We notice that computing ofΠ
F
1
(
U
) involves an eigendecomposition of U and then modifying the eigenvalues by solving a monotone, piecewise linear equation.


## 2.1. gKDR Method

To give the gKDR method, we need some technical notations. For a compact setΩ in some metric space, we denote apositive-definite kernel on Ω
×
Ω, that is, a symmetric function k
:
Ω
×
Ω
→
R satisfying the following semidefinite positive condition: ∑
i
,
j
=
1
n
c
i
c
j
k
(
x
i
,
x
j
)
≥
0 for any x
1
,
…
,
x
n in Ω and c
1
,
…
,
c
n
∈
R. It is known [15] that a positive-definite kernel on Ω is uniquely associated with a Hilbert space H consisting of functions on Ω, and H is called areproducing kernel Hilbert space (RKHS) associated with the kernel k.For short, we denoteZ
≔
{
Φ
n
x
,
x
∈
X
} as a subset of R
M
(
p
+
q
). Let (
Z
,
B
Z
,
μ
Z
) and (
Y
,
B
Y
,
μ
Y
) be measurable spaces and (
Z
,
Y
) be a random variable on Z
×
Y with probability P. Let k
Z and k
Y be positive-definite kernels on Z and Y, respectively, and the corresponding RKHSs denoted by H
Z and H
Y. It is assumed that E
(
k
Z
(
Z
,
Z
)
) and E
(
k
Y
(
Y
,
Y
)
) are finite. Now we introduce two integral operators, so that an appropriate variant of the problem of estimating θ
o in (3) can be established. Thecross-covariance operator
C
Y
Z
:
H
Z
→
H
Y is defined as the operator as follows: (4)
C
Y
Z
f
y
=
∫
k
y
,
y
~
f
z
~
d
P
z
~
,
y
~
,
∀
f
∈
H
X
.Similarly, C
Z
Z denotes theself-covariance operator on H
Z; that is, (5)
C
Z
Z
f
z
=
∫
k
z
,
z
~
f
z
~
d
P
Z
z
~
,
∀
f
∈
H
Z
.Remark that these definitions are natural extensions of the ordinary covariance metrics on Euclidean spaces. Although C
Y
Z and C
Z
Z depend on the kernels, we omit the dependence for notational simplicity.With these definitions, we now define the populationM
(
p
+
q
)
×
M
(
p
+
q
) matrix M
(
z
)
=
(
M
i
j
(
z
)
) by (6)
M
i
j
z
=
C
Y
Z
C
X
X
-
1
∂
k
Z
·
,
z
∂
z
i
,
C
Y
Z
C
Z
Z
-
1
∂
k
Z
·
,
z
∂
z
j
H
Y
,
i
,
j
=
1
,
…
,
M
p
+
q
.Reference [10] has shown that the eigenvector corresponding to the largest eigenvalue of the M
(
z
) is equal to the true vector θ
o for any z, provided that the maximal eigenvalue is simple. To be precise, the following equality holds: (7)
θ
o
=
v
→
m
a
x
E
M
Z
,where v
→
m
a
x
(
A
) is denoted to be the maximum eigenvector of the matrix A. Thus, the problem of estimating θ
o of (3) is equivalent to the eigenvalue problem of solving M
(
z
).We now turn to the empirical version ofM
(
z
), so as to estimate θ
o based on available sample. Denote n
×
n Gram matrices (
k
Z
(
z
i
,
z
j
)
) and (
k
Y
(
y
i
,
y
j
)
) by G
Z and G
Y, respectively. Let ∇
k
Z
(
z
)
=
(
∂
k
Z
(
z
1
,
z
)
/
∂
z
,
…
,
∂
k
Z
(
z
n
,
z
)
/
∂
z
)
∈
R
n
×
M
(
q
+
p
). Then, the gKDR proposes the eigenvectors of M
(
p
+
q
)
×
M
(
p
+
q
) symmetric matrix: (8)
M
n
z
≔
∇
k
Z
z
′
G
Z
+
n
ϵ
n
I
n
-
1
G
Y
G
Z
+
n
ϵ
n
I
n
-
1
·
∇
k
Z
z
,where ϵ
n is a regularization parameter in Thikonov-type regularization. Then, we define (9)
M
^
n
≔
1
n
∑
i
=
1
n
M
n
z
ias an empirical estimator of θ
o. In the next two subsections, we introduce two variational forms with respect to the eigenvalue problems, so that our proposed methods based on Lasso idea can be formulated naturally.


## 2.2. A Nonconvex Variational Algorithm

In order to employ the popular Lasso idea, we need to transform the eigenvalue problem (7) to an optimization scheme. Recall the followingCourant-Fischer variational representation of the maximal eigenvalue and eigenvector: (10)
m
a
x
θ
∈
R
M
p
+
q
θ
′
E
M
Z
θ
,
s
.
t
.
θ
2
2
=
1
.To identify the linear and nonlinear terms in model (3), we define an estimator by adding twofold convex penalties to (10) such that some sparse estimators might be generated. To be precise, let θ
=
(
θ
1
,
…
,
θ
p
+
q
), where each θ
j
∈
R
M. The penalized estimation method is defined as follows: (11)
θ
^
n
=
a
r
g
m
i
n
θ
∈
R
M
p
+
q
-
θ
′
M
^
n
θ
+
λ
∑
j
=
1
p
+
q
α
θ
j
,
1
2
+
1
-
α
θ
j
,
-
1
2
,
s
.
t
.
θ
2
2
=
1
,where θ
j
,
1 is the first element of θ
j and θ
j
,
-
1 is the last M
-
1 component of θ
j
,
j
=
1
,
…
,
p
+
q. Note that λ is the penalty parameter and α is a tuning parameter, typically, α
=
M
/
(
1
+
M
) for normalization. Our penalized method is motivated by the following fact. For any given function f
(
x
)
=
∑
k
=
1
M
β
k
ϕ
k
(
x
), by the cubic spline representation, as mentioned earlier, the first spline base is the single function x and all the bases are linearly independent of each other. Consequently, β
1
=
0 and β
k
≠
0 with k
>
1 if and only if f is a nonlinear function. Based on this and the lasso constraint in (11), the proposed method can identify both the linear and nonlinear terms, of course, also finding the relevant feature simultaneously.We observe that optimization of  (11) is complicated, since it belongs to a class of constrained nonconvex programmes. In this article, we consider an alternative formulation of (11); precisely, we use the ADMM [16]. ADMM is an algorithm that is intended to blend the decomposability of dual ascent with the superior convergence properties of the method of multipliers. Let ρ
(
θ
)
=
∑
j
=
1
p
+
q
(
α
θ
j
,
1
2
+
1
-
α
θ
j
,
-
1
2
), and by introducing an additional constraint, we can reformulate algorithm (11) as (12)
m
i
n
θ
,
γ
∈
R
M
p
+
q
∞
×
Π
2
γ
-
θ
′
M
^
n
θ
+
λ
ρ
θ
,
s
.
t
.
θ
=
γ
,where Π
2 is the 0
-
1 indicator function for the unit l
2-ball surface in R
M
(
p
+
q
) and we adopt the convention ∞
×
0
=
0. As in the method of multipliers, we form the augmented Lagrangian with a parameter τ: (13)
L
γ
,
θ
,
ν
,
τ
=
∞
×
Π
1
γ
-
θ
′
M
^
n
θ
+
λ
ρ
θ
+
ν
′
θ
-
γ
+
τ
2
θ
-
γ
2
2
,where (
ν
,
τ
) is the Lagrangian multipliers. Defining the residual r
=
θ
-
γ, we have (14)
L
γ
,
θ
,
ν
,
τ
=
∞
×
Π
1
γ
-
θ
′
M
^
n
θ
+
λ
ρ
θ
+
τ
2
r
+
u
2
2
-
τ
2
u
2
2
,where u
=
ν
/
τ is called the scaled dual variable. Using the scaled dual variable, we can express ADMM as (15)
θ
k
+
1
=
arg
⁡
min
θ
⁡
λ
ρ
θ
-
θ
′
M
^
n
θ
+
τ
2
θ
-
γ
k
+
u
k
2
2
,
(16)
γ
k
+
1
=
Π
B
2
θ
k
+
1
+
u
k
,
(17)
u
k
+
1
=
u
k
+
θ
k
+
1
-
γ
k
+
1
,where Π
B
2 is Euclidean projection onto the unit l
2-ball surface.There are many convergence results for ADMM discussed in the literature. ADMM consists of iteratively minimizingL
(
γ
,
θ
,
ν
,
τ
) with respect to γ, minimizing L
(
γ
,
θ
,
ν
,
τ
) with respect to θ, and then updating the dual variable u. The θ-update involves iterative evaluation of the proximal operator associated with (18)
F
θ
≔
-
θ
′
M
^
n
θ
+
τ
2
θ
-
γ
k
+
u
k
2
2
.To this end, we use theblock coordinate descent (BCD) [17] for l
1
/
l
2-regularization to compute θ
k
+
1 of (15). Specially at tth iteration, the BCD solves a problem of the form (19)
m
i
n
θ
g
∇
g
F
θ
t
k
′
θ
g
-
θ
g
,
t
k
+
1
2
θ
g
-
θ
g
,
t
k
′
H
g
g
θ
g
-
θ
g
,
t
k
+
μ
θ
g
2
,where θ
g is the subelement of θ with group g
∈
G forming a partition of {
1
,
…
,
M
(
p
+
q
)
} and H
g
g equals or approximates ∇
g
g
2
F
(
θ
t
k
). Note that μ is taken as λ
α or λ
(
1
-
α
) according to the specified group. Obviously, the BCD can be viewed as a generalization of the classical coordinate descent for l
1-regularization. In our specific setting, a direct computation on F
(
θ
) yields that (20)
∇
g
F
θ
=
-
2
θ
′
M
^
n
,
g
1
,
…
,
θ
′
M
^
n
,
g
k
+
τ
θ
+
u
k
-
γ
k
g
,
(21)
∇
2
F
θ
g
g
=
-
2
M
^
n
g
g
+
τ
I
g
,where we write g
=
(
g
1
,
…
,
g
k
) and M
^
n
,
g
l is the lth column of M
^
n
,
l
=
1
,
…
,
k. I
g
g is the identity and (
M
^
n
)
g
g is the submatrix of M
^
n with the index g. Since (
M
^
n
)
g
g in (21) is usually not a multiple of the identity, the solution minimizing (19) has no closed form. We here use a simple but efficient replacement; that is, H
g
g
=
h
g
g
I
|
g
|, where h
g
g is the largest eigenvalue of -
2
(
M
^
n
)
g
g
+
τ
I
|
g
|, which can be computed efficiently, since -
2
(
M
^
n
)
g
g
+
τ
I
|
g
| is only a submatrix indexed by g and independent of θ. In this case, the iterative solution θ
t
+
1
k of (19) is obtained bygroup soft thresholding of the Newton step: (22)
θ
g
,
t
+
1
k
=
P
r
o
x
μ
·
2
θ
g
,
t
k
-
h
g
g
-
1
∇
g
F
θ
t
k
,where P
r
o
x
μ
·
2 is called the soft thresholding operator associated with ·
2, defined as P
r
o
x
μ
·
2
(
θ
g
)
=
1
-
μ
/
θ
g
2
+
θ
g for any θ
g. Under mild conditions, θ
t
+
1
k converges to some vector as t goes to infinity, denoted by θ
k
+
1. It is worthwhile to note that, by (20), we find that the iteration in (22) can be computed parallel between groups.With regard to projection computation of (16), it is known that γ
k
+
1
=
θ
k
+
1
+
u
k
/
θ
k
+
1
+
u
k
2 if θ
k
+
1
+
u
k
2
≠
0; otherwise, γ
k
+
1 is set to any given value on l
2-ball surface. Thus, An efficient implementation for ADMM is obtained; to be precise, the overall procedure can be stated as in Algorithm 1.Algorithm 1:ADMM algorithm.
given: parameters λ, α, ϵ
n,  τ, the group indexes g and the kernels k
Z, k
Y.
initialize: θ
0
=
0, γ
0
=
0, u
0
=
0,
for  k
≥
0
repeat
for  t
≥
0  repeat
θ
g
,
t
+
1
k
=
P
r
o
x
μ
⋅
2
(
θ
g
,
t
k
-
h
g
g
-
1
∇
g
F
(
θ
t
k
)
), for each group g,
t
←
t
+
1
until  θ
g
,
t
k converges to some θ
g
∗
θ
g
k
+
1
←
θ
g
∗, then compute γ
k
+
1
=
Π
B
2
(
θ
k
+
1
+
u
k
),
andu
k
+
1
=
u
k
+
θ
k
+
1
-
γ
k
+
1
thenk
←
k
+
1
until
θ
k
+
1 converges to θ
^
n.
then return
(
θ
^
n
)


## 2.3. A Semidefinite Programming Relaxation

This section intrduces a convex relaxation of (10). Since estimate of individual eigenvectors of θ
o in (7) is unstable if the gap between their eigenvalues is small, it seems reasonable to instead focus on their span, that is, the principal subspace of variation. Let Λ
o
=
θ
o
θ
o
′ as the projection onto the subspace spanned by θ
o; a lesser known but equivalent variational representation of (10) is of the convex program: (23)
m
a
x
Θ
∈
S
+
M
p
+
q
,
tr
⁡
Θ
=
1
tr
⁡
E
M
Z
Θ
,where S
+
M
(
p
+
q
)
=
{
Θ
∈
R
M
(
p
+
q
)
×
M
(
p
+
q
)
∣
0
⪯
Θ
⪯
I
}. The optimization problem (23) is a semidefinite program (SDP), which can be solved exactly in polynomial time. If the maximal eigenvalue is simple, the optimum always achieved at a rank-one matrix, denoted by Θ
-, corresponds to the maximal eigenvector; that is, Θ
-
=
θ
o
θ
o
′. Motivated by this nice conclusion, we can estimate the sparse vector θ
o by adding l
1-type regularization, (24)
Θ
^
≔
arg
⁡
m
a
x
Θ
∈
S
+
M
p
+
q
tr
⁡
M
^
n
Θ
-
λ
Θ
1,1
,
s
.
t
.
tr
⁡
Θ
=
1
,and computing the maximal eigenvector θ
^
=
v
→
m
a
x
(
Θ
^
).Using the ADMM again, we rewrite (11) as the equivalent problem with equality constrained (25)
m
i
n
∞
×
1
F
1
Θ
-
M
^
n
,
Θ
+
λ
Δ
1,1
s
.
t
.
Θ
-
Δ
=
0
,where F
1
=
{
Θ
∈
S
M
p
+
q
,
tr
⁡
Θ
=
1
} and 1
F
1 is the 0
-
1 indicator function for F
1. The augmented Lagrangian associated with (25) has the form (26)
L
τ
Θ
,
Δ
,
U
≔
∞
×
1
F
1
Θ
-
M
^
n
,
Θ
+
λ
Δ
1,1
+
τ
2
Θ
-
Δ
+
U
F
2
-
U
F
2
.
Θ and Δ updates mainly involve computing the projection operator and the proximal operator, respectively; that is, (27)
Π
F
1
Δ
-
U
+
M
^
n
τ
≔
arg
⁡
min
Θ
∈
F
1
⁡
1
2
Θ
-
Δ
-
U
+
M
^
n
τ
F
2
,
S
λ
/
τ
Θ
+
U
≔
arg
⁡
min
Δ
⁡
λ
τ
Δ
1,1
+
1
2
Θ
+
U
-
Δ
F
2
,where S
λ
/
τ is the elementwise soft thresholding operator defined as (28)
S
λ
/
τ
x
=
sign
⁡
x
max
⁡
x
-
λ
τ
,
0
.On the other hand, Π
F
1 is the Euclidean projection onto F
1 and has a closed form as follows.Lemma 1 (Fantope projection).
IfU
=
∑
j
γ
j
u
j
u
j
′ is a spectral decomposition of U, then Π
F
1
(
U
)
=
∑
j
γ
j
+
θ
u
j
u
j
′, where γ
j
+
θ
=
min
⁡
max
⁡
γ
j
-
θ
,
0
,
1 and θ satisfies the equation ∑
j
γ
j
+
(
θ
)
=
1.We notice that computing ofΠ
F
1
(
U
) involves an eigendecomposition of U and then modifying the eigenvalues by solving a monotone, piecewise linear equation.


## 3. Statistical Properties

Before giving statistical results of the estimator (24) in terms of consistency and model selection properties, we recall the existing error bound [10] for M
n
(
z
)
-
M
(
z
) in Frobenius norm for any z.Recall that the following technique conditions are required to derive the gradient-based estimator. For shorthand, letm
=
M
(
p
+
q
).(i)
H
Z and H
Y are separable.
(ii)
k
Z and k
Y are measurable, and E
[
k
Z
(
Z
,
Z
)
]
<
∞, E
[
k
Y
(
Y
,
Y
)
]
<
∞.
(iii)
k
Z
(
z
,
z
~
) is continuously differentiable and ∂
k
Z
(
·
,
z
)
/
∂
z
i
∈
R
(
C
Z
Z
) for i
=
1
,
…
,
m.
(iv)
E
[
k
Y
(
y
,
Y
)
∣
Z
=
·
]
∈
H
Z for any y
∈
Y.
(v)
DefineE
[
h
(
Y
)
∣
Z
]
=
φ
h
(
β
′
Z
) for any given function h. Suppose that φ
h is differentiable with respect to z, and the linear function h
→
∂
φ
h
(
z
)
/
∂
z
a is continuous for any z
∈
R
m and a
=
1
,
…
,
m.
(vi)
There isϑ
m
≥
0 and L
m
≥
0 such that some h
a
,
z
m
∈
H
Z
m satisfies (29)
∂
k
Z
·
,
z
∂
z
a
=
C
Z
Z
ϑ
m
+
1
h
a
,
z
m
,
a
=
1
,
…
,
m
,
andE
[
h
Z
a
H
Z
]
≤
L
m
.
(vii)
Letα
m
≔
(
E
[
k
m
(
Z
,
Z
)
]
-
E
[
k
m
(
Z
,
Z
~
)
]
)
1
/
2, where Z
~ is an i.i.d. copy of Z. Then α
m
/
n
→
0
(
n
→
∞
).Theorem 2.
Under Assumptions (i)–(vii), for the choiceε
n
=
n
m
a
x
{
1
/
3
,
1
/
2
β
+
2
}, we have (30)
M
n
z
-
M
z
F
=
O
p
m
L
m
2
α
m
2
n
m
i
n
1
/
3
,
2
ϑ
m
+
1
/
4
ϑ
m
+
4
,for every z
∈
Z as n
→
∞. If additionally, m
L
m
2
/
n
→
0  (
n
→
∞
), then M
^
n converges in probability to E
(
M
(
Z
)
) of the order O
p
m
L
m
2
/
n
+
m
L
m
α
m
2
/
n
m
i
n
{
1
/
3
,
2
ϑ
m
+
1
/
4
ϑ
m
+
4
} in Frobenius norm.For simplicity, we writeδ
o
=
v
→
1
(
E
[
M
(
Z
)
]
)
-
v
→
2
(
E
[
M
(
Z
)
]
). We write the second error bound of Theorem 2 by Ω
(
n
,
m
)
≔
O
p
m
L
m
2
/
n
+
m
L
m
α
m
2
/
n
m
i
n
{
1
/
3
,
2
ϑ
m
+
1
/
4
ϑ
m
+
4
}.Theorem 3.
Under the assumptions of Theorem2, let s
=
Λ
o
2,0 and λ is chosen to satisfy λ
=
Ω
(
n
,
m
); then (31)
Θ
^
-
Λ
o
F
≤
4
s
δ
o
Ω
n
,
m
.The proof is shown in the Appendix. The following theorem provides a sufficient condition for support recovery by diagonal thresholdingΘ
^.Theorem 4.
Under the assumptions of Theorem3, for any t
>
0, (32)
j
:
Λ
o
j
j
=
0
,
Θ
^
j
j
≥
t
+
j
:
Λ
o
j
j
≥
2
t
,
Θ
^
j
j
≤
t
≤
4
s
δ
o
t
Ω
n
,
m
.Therefore, the feature selection procedure {
j
:
(
Θ
^
)
j
j
≥
t
} succeeds with high probability if m
i
n
j
:
(
Λ
o
)
j
j
≠
0
(
Λ
o
)
j
j
≥
2
t
>
8
s
/
δ
o
Ω
(
n
,
m
).The above theorem follows immediately from the conclusion of Theorem3.



---

*Source: 3537564-2016-09-27.xml*