# Enhancement of polarization and magnetization in polycrystalline
magnetoelectric composite

K.P. Jayachandran<EMAIL_ADDRESS>J.M. Guedes  H.C. Rodrigues
IDMEC, Instituto Superior Técnico, University of Lisbon, Av. Rovisco Pais,
1049-001 Lisbon , Portugal

###### Abstract

Electrical control of magnetization or magnetic control of polarization offers
an extra degree of freedom in materials possessing both electric and magnetic
dipole moments _viz.,_ magnetoelectric multiferroics. Microstructure with
polycrystalline configurations that enhances the overall
polarization/magnetization and that outperform single crystalline
configurations are identified. The characterization of local fields
corresponding to the polycrystal configuration underlines nontrivial role
played by randomness in better cross-coupling mediated by anisotropic and
asymmetric strains.

The magnetoelectric (ME) effect manifests in the linear relation between the
magnetic and electric fields in matter and it causes, for instance, either a
magnetization or electrical polarization proportional to applied
electric/magnetic fields [1]. Yet, the ME effect is possible only for certain
magnetic symmetry classes, akin to piezomagnetism. Set of materials classified
as magnetoelectric multiferroics possess both the magnetic and ferroelectric
orders in the same phase[2]. Apart from exhibiting functionalities of both the
orders, a coupling between the ferri-/ferro-magnetic (FM) and ferroelectric
(FE) states enable appearance of novel characteristics not present in either
of the states [3]. Nonetheless, simultaneous occurrence of magnetism and
ferroelectricity in materials is constrained by the conflicting classic
chemical requirements regarding electronic orbital occupancy [4, 5]. In such a
situation, one should envisage only a moderate ME coupling though [6]. For
instance, only few examples exist of materials possessing intrinsic ME
coupling at room temperature [7, 8]. An unambiguous resolution to circumvent
this chemical ’contraindication’ is by juxtaposing an FE and FM material
artificially into strain-mediated multiphase materials such as composites
combining the two phases, that thereby yields substantial ME effect [5, 6, 9].
Moreover, the composites do not need to adhere to the symmetry restrictions of
single phase materials and consequently do possess the freedom of choice from
a wide variety of ferroics existing above room temperature. Crystallites of
the two phases in an ME composite are assumed to be in good mechanical contact
if the two phases are polycrystalline. The change in shape of the FE grains in
response to an applied electric field causes the ferromagnetic grains to
deform, resulting in a change in magnetization [10].

Figure 1:  Schematic hierarchical diagram of the ME multiferroic
microstructure browsing through the succession of modeling tasks in the study.
The microstructure $Y$ is a fibre-reinforced ME composite where an FM fibre
made of CoFe2O4 is embedded in a perovskite FE matrix of BaTiO3. BaTiO3 is a
polycrystal – an aggregate of grains of various orientations (depicted in
different colours). The orientation of the grains quantified by the Euler
angles $(\phi,\theta,\psi)$ is in fact that of the underlying crystallographic
unit cell and is measured against the microscopic ($y_{1},y_{2},y_{3}$)
coordinate system.

Magneto-elastic effect by means of strain transfer across the interface of an
FM/nonmagnetic bilayer could unleash significant changes in the magnetic
properties of the ferromagnet. When an FE material is used as the nonmagnetic
layer, strain transfer can be reinforced through the induction by an electric
field thanks to the inverse piezoelectric effect prevails in the FE material.
This additional strain transfer allows one to manipulate the magnetisation by
an electric field [11]. The concept of magnetic information recording and
memory device technologies as well as electric-field driven magnetoelastic
effects on magnetic anisotropy have attracted much interest recently from both
fundamental and technological aspects owing to the said modification
accomplished on magnetisation orientation [12, 11, 5]. Among many FM/FE hybrid
structures, perovskite manganites/BTO are the well-suited magnetoelastic
heterostructures since transport, magnetic, and electronic properties of
perovskite manganites are very sensitive to lattice strain and can be easily
integrated with BaTiO3 (BTO) due to the similar crystal structure [13]. For
instance when an electric field is applied to the nanostructured BTO/CFO 1–3
composite, the strain created by the piezoelectric (BTO) matrix couples to the
magnetostrictive CFO(CoFe2O4) pillars and hence the strength of the coupling
between the ferroelectricity and magnetism are critically dependent on the
anisotropy and crystalline structure of the constituents[14, 15].

In this letter, we would be studying the impact of an applied
electrical/magnetic fields on a motley ME composite composed of disparate FE
and FM materials. We underscore the significance of heterogeneity of one of
the constituent phases, _viz.,_ the FE phase by introducing randomness in the
microstructure. Prior research have shown that introducing heterogeneity
locally can improve the piezoelectric properties of FE materials and the
ensuing coupling effect of the ME composite of which the former is a component
[16, 17]. Apart from the effect of applied fields, the sway of crystal
anisotropy on the polarization/magnetization as demonstrated in various
studies [18, 19, 20, 21]on magnetoelectric heterostructures is also shown.
Another important challenge is to identify the optimum texture and stress
level and efficient control of the same to achieve targeted magnetoelectric
properties [22]. We demonstrate the effect of local texture in a 1–3
composite, tantamount to the architecture used in Zheng et al., [14] to map
its local field distributions upon applying an external bias field. Moreover,
the average polarization and magnetization post application of ample bias
fields sufficient to saturate the composite would be computed.

Estimation of the equilibrium macroscopic magneto-electro-electric properties
of the ME composite,have been done using the mathematical homogenization
method [23, 24]. A ME multiferroic composite occupying a volume $\Omega$ of
coordinates $\mathbf{x}$ (or $x_{i}$) is considered and the material
properties are assumed to change periodically and that period $\epsilon$ is
characterized by the dimension of an elementary cell $\mathrm{Y}$ of
coordinates $\mathbf{y}$ (or $y_{i}$) of the body (Fig. 1). This lead us to
examine fields for unknown physical quantities in the form of two-scale
asymptotic expansion in powers of $\epsilon$ as

$\zeta^{\epsilon}(\mathbf{x})=\zeta^{0}(\mathbf{x,y})+\epsilon\zeta^{1}(\mathbf{x,y}),\text{with}~{}\mathbf{y}=\mathbf{x}/\epsilon\\\
$ (1)

and the equivalent macroscopic behaviour is estimated to be first order by the
behaviour of $\zeta^{0}$. Here $\zeta$ stands in lieu of the fields _viz.,_
$\mathbf{u}$, the displacement, $\varphi$ and $\psi$, the scalar
electrical/magnetic potentials respectively of a multiferroic composite [25].
Here the functions
$\zeta^{1}(\equiv\mathbf{u}^{1},\varphi^{1}~{}\text{or}~{}\psi^{1})$ are the
local variations (or _field perturbations_) describing the heterogeneous part
of the solutions and are associated with $\mathbf{y}=\mathbf{x}/\epsilon$.
Applying calculus of variations, utilising the asymptotic Eq. (1) and its
derivatives, the homogenization method essentially culminates in the
characterization of effective magneto-electro-elastic moduli when the
characteristic length $\epsilon$ of the period tends to zero. The electrical,
magnetic and mechanical fields (_degrees of freedom_) govern the constitutive
equations of the linear ME multiferroic solid (see Supplementary material).
The model quantifies all the local fields and potentials besides stress
($\mathbf{\sigma}$) and strain ($\mathbf{\varepsilon}$) fields, magnetization
(via magnetic flux density, $\mathbf{B}$), polarization (via electric
displacement, $\mathbf{D}$ as well as von-Mises stress upon the application of
an external field (either electric/magnetic or mechanical). The _macroscopic
averages_ $\langle\sigma_{ij}\rangle$, $\langle D_{i}\rangle$, and $\langle
B_{i}\rangle$ can be computed once the homogenized solution for
$\mathbf{u}^{0}$, $\varphi^{0}$ and $\psi^{0}$ and that of the corresponding
fields _viz.,_ $\frac{\partial u^{0}_{j}(\mathbf{x})}{\partial x_{k}}$,
$\frac{\partial\varphi^{0}(\mathbf{x})}{\partial x_{j}}$ and
$\frac{\partial\psi^{0}(\mathbf{x})}{\partial x_{j}}$ are prescribed (see
Supplementary material). This _postulate_ is equally applicable for the cases
when magnetic field $\textbf{H}=0$, the applied macroscopic electric field
$\langle E_{i}\rangle$ generates average magnetization [1]

$\langle M_{k}\rangle=\widetilde{\alpha}_{ik}\langle E_{i}\rangle$ (2)

and when electric field $\textbf{E}=0$, the applied macroscopic magnetic field
$\langle H_{i}\rangle$ generates average electrical polarization [1]

$\langle P_{i}\rangle=\widetilde{\alpha}_{ik}\langle H_{k}\rangle$ (3)

Figure 2:  Average fields for the multiferroic ME 1–3 composite BTO–CFO. (a)
shows the average dielectric displacement $\langle D_{i}\rangle$. (b) is the
average magnetic flux density, $\langle B_{i}\rangle$ (c) shows the average
magnetization $\langle M_{i}\rangle$ and (d) show the relevant tensor
components of the average stress $\langle\sigma_{IJ}\rangle$ upon applying an
external electric field $\langle E_{3}\rangle$ along the $y_{3}$–axis of the
composite (i.e., along the normal to the composite plane). The data fall in
the shaded area corresponds to single crystal BTO matrix–polycrystalline CFO
fibre constituting the ME composite. The BTO orientations in Euler angles
$(\phi^{\circ},\theta^{\circ},\psi^{\circ})$ are marked along the horizontal
axis in the shaded area. Rest of the data correspond to polycrystalline ME
composite with BTO orientation distribution
$(\mu_{\phi},\sigma_{\phi},\mu_{\theta},\sigma_{\theta},\mu_{\psi},\sigma_{\psi})$
(along the horizontal axis). Here the x–axes is representative line separating
the quantities plotted on the y–axes. i.e., the data points in the plot have
no abscissae. Figure 3:  Average fields for the multiferroic ME 1–3 composite
BTO–CFO. (a) shows the average dielectric displacement $\langle D_{i}\rangle$.
(b) is the average magnetic flux density, $\langle B_{3}\rangle$ (c) shows the
average magnetization $\langle M_{i}\rangle$ and (d) show the relevant tensor
components of the average stress $\langle\sigma_{IJ}\rangle$ upon applying an
external magnetic field $\langle H_{3}\rangle$ along the $y_{3}$–axis of the
composite (i.e., along the normal to the composite plane). Here the x–axes is
representative line separating the quantities plotted on the y–axes. i.e., the
data points in the plot have no abscissae. Figure 4: Pole figures (contour
plots) of the polycrystalline BTO matrix while the ODPs are
$(\mu_{\phi}=0,\sigma_{\phi}=1,\mu_{\theta}=0,\sigma_{\theta}=1,\mu_{\psi}=0,\sigma_{\psi}=1)$.
At this configuration, the ME composite deliver maximum $\langle M_{3}\rangle$
and $\langle P_{3}\rangle$.

The FE matrix of the ME composite is composed of polycrystalline BaTiO3. The
orientation of the BaTiO3 crystallographic grains in the composite is
characterized using Euler angles. It is seen from the Euler angles
$(\phi,\theta,\psi)$ (Fig. 1) that it measure the rotation of the
crystallographic coordinates (x, y, z) with reference to the microstructure
coordinates ($y_{1},y_{2},y_{3}$). The three Euler angles are modelled as
statistically distributed in a normal distribution after poling obeying the
probability distribution function
$f(\alpha\mid\mu,\sigma)=1/(\sigma\sqrt{2\pi})exp-[(\alpha-\mu)^{2}/2\sigma^{2}]$.
Here $\alpha$ is the random variable representing the orientations (i.e., the
Euler angles) $(\phi,\theta,\psi)$ and $\mu$ and $\sigma$ are the mean and the
standard deviation, respectively. The fibres embedded in the matrix is
composed of the FM material CoFe2O4 which is treated as bulk material without
assigning any orientation whatsoever. The numerical model developed is
implemented in finite element method. The convergence of piezoelectric
properties with unit cell size allows us to determine the simulation-space
independent, macroscopic magnetoelectric properties at various distribution of
grains (See Supplementary material).

Table 1: Average values of magnetization $\langle M_{3}\rangle$ (in A/m) at constant external electric field, and electric polarization $\langle P_{3}\rangle$ (C/m2) at external magnetic field of 1–3 magnetoelectric composite BaTiO3 (single and poly-crystalline)–ceramic CoFe2O4.  BTO phase | Orientation 111Euler angles $(\phi,\theta,\psi)$ (rad) of rotation in single crystal BTO phase | Averages
---|---|---
| or ODP222Orientation Distribution Parameters $(\mu_{\phi},\sigma_{\phi})\|(\mu_{\theta},\sigma_{\theta})\|(\mu_{\psi},\sigma_{\psi})$ (rad) of polycrystal BTO phase | $\langle M_{3}\rangle$ | $\langle P_{3}\rangle$
Single crystalline | (0,0,0) | 2.63$\times 10^{4}$ | 0.03
(0,$\pi/4$,0) | -5.48$\times 10^{4}$ | -0.07
(0,$\pi/2$,0) | -1.39 | -2$\times 10^{-6}$
Polycrystalline | (0,$\frac{1}{2}$)$\|$(0,$\frac{1}{2}$)$\|$(0,$\frac{1}{2}$) | 1.38$\times 10^{5}$ | 0.17
(0,1)$\|$(0,0)$\|$(0,0) | 2.47$\times 10^{4}$ | 0.03
(0,0)$\|$(0,1)$\|$(0,0) | 1.3$\times 10^{5}$ | 0.16
(0,0)$\|$(0,0)$\|$(0,1) | 2.47$\times 10^{4}$ | 0.03
(0,1)$\|$(0,1)$\|$(0,1) | 1.49$\times 10^{5}$ | 0.19
(0,$\pi$)$\|$(0,$\pi$)$\|$(0,$\pi$) | 1.46$\times 10^{4}$ | 0.02
Experiment 333CoFe2O4 nanopillars embedded BaTiO3 matrix in a 1–3 composite from Ref. [14]. | – |  $\approx$ 3.5$\times 10^{5}$ |  $\approx$ 0.23
Experiment 444CoFe2O4–BaTiO3 1–3 nanocomposite from Ref. [26]. | – |  $\approx$ 3.15$\times 10^{5}$ |  $\approx$ 0.1

The average magnetic property response along the $y_{3}$–axis of the local
coordinates (which was set to align along the normal to the composite plane)
consequent to an external biasing electrical field is studied first. Here we
apply an electric field ( $E_{3}$ ) which suffices to saturate the
polycrystalline BTO–CFO composite. The results are summarised in Fig. 2. In
general, the average fields acting through out-of-plane (the components along
$y_{3}$–axis) to the composite plane are greater in magnitude than the other
components. Unlike the common perception of great performance of single
crystal piezoelectrics, we see much better magnetization $\langle
M_{3}\rangle$ while the BTO is still polycrystalline (Fig. 2 and Table 1). In
Fig. 2 the polycrystalline data of the averages are differentiated from single
crystal by painting the single crystal region by a shade. Randomness,
introduced by local microstructural heterogeneity could potentially enhance
piezoelectricity in relaxor ferroelectric ceramics [16]. Here, we have
incorporated randomness in the ferroelectric BTO phase through Euler
orientations of the constituent grains kept at a normal distribution but with
a standard deviation of $\mu$ rad. The maximum magnetization $\langle
M_{3}\rangle$ is seen for BTO polycrystal phase with orientation distribution
parameters
$(\mu_{\phi}=0,\sigma_{\phi}=1,\mu_{\theta}=0,\sigma_{\theta}=1,\mu_{\psi}=0,\sigma_{\psi}=1)$.
Each grain orientation would be dissected and kept separately and the grains
are allowed to chose the combination of three $(\phi,\theta,\psi)$ but
dictated by the corresponding mean ($\mu$) and standard deviation ($\mu$). The
value of $\langle M_{3}\rangle=1.49\times 10^{5}(A/m)$ obtained here compares
with the order of $\langle M_{3}\rangle$ measured by Zheng et al., [14] on
nanostructures with vertically aligned CoFe2O4 nanopillars embedded in BaTiO3
matrix. Zheng et al. value is the saturation magnetization ($M_{s}$) response
against an applied magnetic field (hysteresis) of the composite rather than
the cross coupling magnetization value as is obtained in the present work
(Table 1) and hence the deviation is obvious. Epitaxial CFO films on BaTiO3
single crystal substrate shows $\langle M_{3}\rangle\approx 2.5\times
10^{5}~{}(A/m)$ [27] under a magnetic field.

The property variation at various configurations (single crystalline data in
shaded area) consequent to the application of external magnetic field $H_{3}$
are given in Fig. 3. Here the electric field is off and hence the electrical
responses are dictated solely by the piezomagnetic effect of the CFO component
of the ME composite. More precisely, the emergence of electrical displacement
$\langle D_{j}\rangle$ (Fig. 3) which bears the signature of polarization and
the electrical polarization $\langle P_{j}\rangle$ (Fig. 3 and Table 1)
itself. The average polarization out-of-plane of the composite layer $\langle
P_{3}\rangle$ peaks at the polycrystal phase of BTO as seen in magnetization
$\langle M_{3}\rangle$ in the previous simulation experiment with applied
electric field (Fig. 2). The polarization $\langle P_{3}\rangle$ compares well
with the experiments (Table 1)

The causal relationship between local microstructure and better magnetic
response of polycrystalline BTO–CFO composite can be explored from mapping the
local field distribution in response to the external field. The average or
macroscopic field applied to a magnetoelectric composite would permeate the
material microstructure and would spread unequally into different points owing
to the heterogeneity of the material. The anisotropy due to the underlying
crystalline structure of the constituent phases contribute to this phenomenon.
The associated local field distribution exhibit spacial fluctuations critical
to coupling phenomena. The stress/strain mediating the coupling phenomenon is
spread non-uniformly across the microstructure as is seen (in Fig.
SLABEL:fig:S1) in Supplementary material. The equivalent von Mises stress
indicates a stress concentration across BTO matrix compared to the CFO fibre.
This is reflected and underlined in the variation of microscopic displacement
$\mathbf{u}^{\varepsilon}$ (see Fig. SLABEL:fig:S1). The asymmetry and
anisotropy of the local stress/strain distributed across the microstructure
could be seen in histograms (Fig. SLABEL:fig:S2). The in-plane strain
$\epsilon_{12}^{\varepsilon}$ records values orders of magnitude greater than
other components (Fig. SLABEL:fig:S2) conforming the compression of the
composite in $x$ and $y$ directions [26]. The cross-coupling effect of the
applied electric field results in appreciable local magnetic potential
$\psi^{\varepsilon}$ and associated magnetism (Fig. SLABEL:fig:S3). The
magnetic field induced local profile of fields are shown in the Supplementary
materials. It reinforces the cross-coupling between the magnetic and electric
degrees of freedom through mechanical stress. The pole figure in Fig. 4 shows
the distribution of grain orientation about the $y_{3}$–axis (or $z$–axis) of
the composite structure. The $<001>$ axes of the ferroelectric BTO (where the
spontaneous polarization is oriented ) is aligned mostly along the
$y_{3}$–axis (here the $[001]$ psuedo cubic axes of BTO is directed ou-of-
plane).

In summary, strong magnetoelectric coupling is resulted while a
polycrystalline ME composite of 1–3 BaTiO3–CoFe2O4 is subjected external
electric/magnetic fields. In contrast to single crystal BaTiO3–polycrystal
CoFe2O4 composite, the averages of polarization and magnetization of the
polycrystalline BaTiO3–CoFe2O4 exceeds that of the single crystal version. The
depiction of local fields corresponding to the polycrystal configuration
suggests nontrivial role played by randomness in better cross-coupling
mediated by anisotropic and asymmetric strains.

##  I Acknowledgments

This work was supported by FCT-Fundação para a Ciência e a Tecnologia, through
IDMEC, under LAETA, project UIDB/50022/2020.

## References

  * Landau and Lifshitz [1960] L. D. Landau and E. M. Lifshitz,  _Electrodynamics of Continuous Media_ (Pergamon, Oxford, 1960).
  * Schmid [1994] H. Schmid, “Multiferroic magnetoelectrics,” Ferroelectrics 162, 317–338 (1994).
  * Fiebig  _et al._ [2016] M. Fiebig, T. Lottermoser, D. Meier, and M. Trassin, “The evolution of multiferroics,” Nat. Rev. Mater. 1, 16046 (2016).
  * Hill [2000] N. A. Hill, “Why are there so few magnetic ferroelectrics?” J. Phys. Chem. B 104, 6694–6709 (2000).
  * Spaldin and Ramesh [2019] N. A. Spaldin and R. Ramesh, “Advances in magnetoelectric multiferroics,” Nat Mater 18, 203–212 (2019).
  * Cheong and Mostovoy [2007] S.-W. Cheong and M. Mostovoy, “Multiferroics: a magnetic twist for ferroelectricity,” Nat Mater 6, 13–20 (2007).
  * Mandal  _et al._ [2015] P. Mandal, M. J. Pitcher, J. Alaria, H. Niu, P. Borisov, P. Stamenov, J. B. Claridge, and M. J. Rosseinsky, “Designing switchable polarization and magnetization at room temperature in an oxide,” Nature 525, 363–366 (2015).
  * Long  _et al._ [2020] J. Long, M. S. Ivanov, V. A. Khomchenko, E. Mamontova, J.-M. Thibaud, J. Rouquette, M. Beaudhuin, D. Granier, R. A. S. Ferreira, L. D. Carlos, B. Donnadieu, M. S. C. Henriques, J. A. Paixão, Y. Guari, and J. Larionova, “Room temperature magnetoelectric coupling in a molecular ferroelectric ytterbium(iii) complex,” Science 367, 671–676 (2020).
  * Avellaneda and Harshe [1994] M. Avellaneda and G. Harshe, “Magnetoelectric effect in piezoelectric magnetostrictive multilayer (2–2) composites,” J. Intell. Mater. Syst. Struct. 5, 501–513 (1994).
  * van den Boomgaard and Born [1978] J. van den Boomgaard and R. A. J. Born, “A sintered magnetoelectric composite material BaTiO3-Ni(Co, Mn)Fe2O4,” J. Mater. Sci. 13, 1538–1548 (1978).
  * Taniyama [2015] T. Taniyama, “Electric-field control of magnetism via strain transfer across ferromagnetic/ferroelectric interfaces,” J Phys-Condens Mat 27, 504001 (2015).
  * Nan  _et al._ [2008] C. W. Nan, M. I. Bichurin, S. Dong, D. Viehland, and G. Srinivasan, “Multiferroic magnetoelectric composites: Historical perspective, status, and future directions,” J. Appl. Phys 103, 031101 (2008).
  * Panchal  _et al._ [2018] G. Panchal, D. M. Phase, V. R. Reddy, and R. J. Choudhary, “Strain-induced elastically controlled magnetic anisotropy switching in epitaxial ${\mathrm{la}}_{0.7}{\mathrm{sr}}_{0.3}\mathrm{Mn}{\mathrm{o}}_{3}$ thin films on $\mathrm{BaTi}{\mathrm{o}}_{3}$(001),” Phys. Rev. B 98, 045417 (2018).
  * Zheng  _et al._ [2004] H. Zheng, J. Wang, S. E. Lofland, Z. Ma, L. Mohaddes-Ardabili, T. Zhao, L. Salamanca-Riba, S. R. Shinde, S. B. Ogale, F. Bai, D. Viehland, Y. Jia, D. G. Schlom, M. Wuttig, A. Roytburd, and R. Ramesh, “Multiferroic BaTiO3-CoFe2O4 nanostructures,” Science 303, 661–663 (2004).
  * Zheng  _et al._ [2007] H. Zheng, J. Kreisel, Y.-H. Chu, R. Ramesh, and L. Salamanca-Riba, “Heteroepitaxially enhanced magnetic anisotropy in batio3-cofe2o4 nanostructures,” Appl Phys Lett 90, 113113 (2007).
  * Li  _et al._ [2018] F. Li, D. Lin, Z. Chen, Z. Cheng, J. Wang, C. Li, Z. Xu, Q. Huang, X. Liao, L.-Q. Chen, T. R. Shrout, and S. Zhang, “Ultrahigh piezoelectricity in ferroelectric ceramics by design,” Nat Mater 17, 349–354 (2018).
  * Jayachandran, Guedes, and Rodrigues [2008] K. P. Jayachandran, J. M. Guedes, and H. C. Rodrigues, “Piezoelectricity enhancement in ferroelectric ceramics due to orientation,” Appl. Phys. Lett. 92, 232901 (2008).
  * Sahoo  _et al._ [2007] S. Sahoo, S. Polisetty, C.-G. Duan, S. S. Jaswal, E. Y. Tsymbal, and C. Binek, “Ferroelectric control of magnetism in BaTiO3–Fe heterostructures via interface strain coupling,” Phys. Rev. B 76, 092108 (2007).
  * Han  _et al._ [2014] X. Han, L. Xi, Y. Li, X. Guo, D. Li, Z. Wang, Y. Zuo, and D. Xue, “Demonstration of magnetoelectric memory cell in (110) [pb(mg1/2nb2/3)o3]0.68-[pbtio3]0.32/ru/feco heterostructures,” Appl Phys Lett 105, 122402 (2014).
  * Zhang  _et al._ [2014] S. Zhang, Y. Zhao, X. Xiao, Y. Wu, S. Rizwan, L. Yang, P. Li, J. Wang, M. Zhu, H. Zhang, X. Jin, and X. Han, “Giant electrical modulation of magnetization in Co40Fe40B20/Pb(Mg1/3Nb2/3)0.7Ti0.3O3(011) heterostructure,” Sci Rep 4, 3727 (2014).
  * Zhou  _et al._ [2014] W. P. Zhou, Q. Li, Y. Q. Xiong, Q. M. Zhang, D. H. Wang, Q. Q. Cao, L. Y. Lv, and Y. W. Du, “Electric field manipulation of magnetic and transport properties in SrRuO3/Pb(Mg1/3Nb2/3)O3-PbTiO3 heterostructure,” Sci Rep 4, 6991 (2014).
  * Stephen [2019] S. Stephen, “Multiferroics and the path to the market,”  (2019), Nian Sun, professor at Northeastern University, talks to Nature Materials about the potential applications of multiferroic materials.(https://doi.org/10.1038/s41563-019-0295-6 ).
  * Sanchez-Palencia [1980] E. Sanchez-Palencia,  _Non-homogeneous media and vibration theory, Lecture notes in physics 127_ (Springer-Verlag, Berlin, 1980).
  * Jayachandran, Guedes, and Rodrigues [2020] K. P. Jayachandran, J. M. Guedes, and H. C. Rodrigues, “Homogenization method for microscopic characterization of the composite magnetoelectric multiferroics,” Sci Rep 10, 1276 (2020).
  * Jayachandran, Guedes, and Rodrigues [2014] K. P. Jayachandran, J. M. Guedes, and H. C. Rodrigues, “A generic homogenization model for magnetoelectric multiferroics,” J Intel Mat Syst Str 25, 1243–1255 (2014).
  * Schmitz-Antoniak  _et al._ [2013] C. Schmitz-Antoniak, D. Schmitz, P. Borisov, F. M. F. de Groot, S. Stienen, A. Warland, B. Krumme, R. Feyerherm, E. Dudzik, W. Kleemann, and H. Wende, “Electric in-plane polarization in multiferroic cofe2o4/batio3 nanocomposite tuned by magnetic fields,” Nat Commun 4, 2051 (2013).
  * Chopdekar and Suzuki [2006] R. V. Chopdekar and Y. Suzuki, “Magnetoelectric coupling in epitaxial cofe2o4 on batio3,” Appl Phys Lett 89, 182506 (2006), https://doi.org/10.1063/1.2370881 .


# A General Framework of Nonparametric Feature Selection in High-Dimensional
Data

Hang Yu, Yuanjia Wang, and Donglin Zeng

Nonparametric feature selection in high-dimensional data is an important and
challenging problem in statistics and machine learning fields. Most of the
existing methods for feature selection focus on parametric or additive models
which may suffer from model misspecification. In this paper, we propose a new
framework to perform nonparametric feature selection for both regression and
classification problems. In this framework, we learn prediction functions
through empirical risk minimization over a reproducing kernel Hilbert space.
The space is generated by a novel tensor product kernel which depends on a set
of parameters that determine the importance of the features. Computationally,
we minimize the empirical risk with a penalty to estimate the prediction and
kernel parameters at the same time. The solution can be obtained by
iteratively solving convex optimization problems. We study the theoretical
property of the kernel feature space and prove both the oracle selection
property and the Fisher consistency of our proposed method. Finally, we
demonstrate the superior performance of our approach compared to existing
methods via extensive simulation studies and application to a microarray study
of eye disease in animals.

Keywords: Tensor product kernel ; Reproducing kernel Hilbert space; Fisher
consistency; Oracle property; Variable selection

††footnotetext: Hang Yu is PhD candidate at Department of Statistics and
Operation Research, University of North Carolina, Chapel Hill, NC 27599
(hangyu@live.unc.edu), Yuanjia Wang is Professor at Department of
Biostatistics, Columbia University, New York, NY 10032
(yw2016@cumc.columbia.edu), and Donglin Zeng is Professor at Department of
Biostatistics, University of North Carolina, Chapel Hill, NC 27599
(dzeng@email.unc.edu). This research is supported by U.S. NIH grants NS073671,
GM124104, and MH117458.

##  1 Introduction

With biotechnology advances in modern medicine, biomedical studies collecting
complex data with a large number of features are becoming the norm. High-
dimensional feature selection is an essential tool to allow using such data
for disease prediction or precision medicine, for instance, to discover a set
of diagnostic biomarkers from neuroimaging measures for early prediction of
neurodegenerative diseases, or to determine predictive biomarkers for
effective management of type 2 diabetic patients’ healthcare. Accurately
identifying the subset of true important features is even more crucial and
challenging than before in the fields of statistics and machine learning.

High-dimensional feature selection has been extensively studied for linear or
generalized linear models in the past decades, and many methods have been
developed including Lasso (Tibshirani, 1996), SCAD (Fan and Li, 2001), MCP
(Zhang, 2010) and (Wang and Kulasekera, 2012). In these parametric models, the
importance of individual features is characterized by non-null coefficients
associated with them, so proper penalization can identify those non-null
coefficients with probability tending to one when the sample size increases.
However, parametric model assumptions are likely to be incorrect for many
biomedical data due to potential correlations and higher-order interactions
among feature variables. In fact, applying these approaches to any simple
transformation of feature variables may lead to very different feature
selection results.

More recently, increasing efforts have been devoted to high-dimensional
feature selection when parametric assumptions, especially linearity
assumption, do not hold. Various approaches were proposed to select features
based on measuring certain marginal dependency (Guyon and Elisseeff (2003),
Fan and Lv (2008), Fan et al. (2011), Song et al. (2012), Yamada et al.
(2014), Urbanowicz et al. (2018)). For example, nonparametric association
between each feature and outcome was used for screening (Fan and Lv (2008),
Fan et al. (2011), Song et al. (2012)). LI et al. (2012) adopted a a robust
rank correlation screening method based on marginal Kendall correlation
coefficient. Yamada et al. (2014) considered a feature-wise kernelized Lasso,
namely HSICLasso, for capturing nonlinear dependency between features and
outcomes. In this approach, after a Lasso-type regression of an output kernel
matrix on each feature-wise kernel matrix, unimportant features with small
marginal dependence in terms of a Hilbert-Schmidt independence criterion
(HSIC) would be removed. However, all methods based on marginal dependence may
fail to select truly important variables since marginal dependency does not
necessarily imply the significance of a feature when other features are also
included for prediction, which is the case even for a simple linear model.

Alternatively, other approaches were proposed to relax parametric model
assumptions and perform feature selection and prediction simultaneously. Lin
and Zhang (2006) proposed COmponent Selection and Smoothing Operator (COSSO)
to perform penalized variable selection based on smoothing spline ANOVA.
Ravikumar et al. (2009) studied feature selection in a sparse additive model
(SpAM), which assumed an additive model but allowed arbitrary nonparametric
smoothers such as approximation in a reproducing kernel Hilbert space (RKHS)
for each individual component function. Huang et al. (2010) considered spline
approximation in the same model and adopted an adaptive group Lasso method to
perform feature selection. Although both COSSO and SpAM allowed nonlinear
prediction from each feature, they still imposed restrictive additive model
structures, possible with some higher-order interactions. To allow arbitrary
interactions among the features and perform a fully nonparametric prediction,
Allen (2013) propsed a procedure named KerNel Iterative Feature Extraction
(KNIFE), in which the feature input was constructed in a Gaussian RKHS in
order to perform nonparametric prediction. Different weights were used for
different features in the constructed Gaussian kernel function so that a
larger weight implied a higher importance of the corresponding feature
variable. However, due to high nonlinearity in the kernel function, estimating
the weights was numerically unstable even when the dimension of the features
was moderate.

In this paper, we propose a general framework to perform nonparametric high-
dimensional feature selection. We consider a general loss function which
includes both regression models and classification as special cases. To
perform nonparametric prediction, we construct a novel RKHS based on a tensor
product of kernels for individual features. The constructed tensor product
kernel, as discussed in Gao and Wu (2012), can handle any high-order nonlinear
relationship between the features and outcome and any high-order interactions
among the features. More importantly, each feature kernel depends on a non-
negative parameter which determines the feature importance, so for feature
selection, we further introduce a $l_{1}$-penalty of these parameters in the
estimation. Computationally, coordinate descent algorithms are used for
updating parameters and each step involves simple convex optimization
problems. Thus, our algorithm is numerically stable and can handle high-
dimensional features easily. Theoretically, we first derive the approximation
property of the proposed RKHS and characterize the complexity of the unit ball
in this space in terms of bracket covering numbers. We then show that the
estimated prediction function from our approach is consistent and moreover, we
show that under some regularity conditions, the important features can be
selected with probability tending to one.

The rest of the paper is organized as follows. In Section 2, we introduce our
proposed regularized tensor product kernel and lay out a penalized framework
for both estimation and feature selection. We then provide detailed
computational algorithms to solve the optimization problem. In Section 3, we
provide two theorems studying the property of the proposed RKHS. We then give
the main result of this paper including the consistency of the estimated
prediction function and the oracle property of the feature selection. In
Section 4, two simulation studies for regression and classification problems
are conducted and we compare our method to existing methods. Application to a
microarray study is given in Section 5. We conclude the paper with some
discussion in Section 6.

##  2 Method

Suppose data are obtained from $n$ independent subjects and consist of
$(\mathbf{X}_{i},Y_{i}),i=1,...,n$, where we let $\mathbf{X}$ denote
$p_{n}$-dimensional feature variables and $Y$ be the outcome which can be
continuous, binary or ordinal. Our goal is to use the data to learn a
nonparametric prediction function, $f(\mathbf{X})$, for the outcome $Y$.

We learn $f(\mathbf{X})$ through a regularized empirical risk minimization by
assuming $f(\cdot)$ belongs to a RKHS associated with a kernel function,
$\kappa(\mathbf{X},\widetilde{\mathbf{X}})$, which will be described later.
Specifically, if we denote the RKHS generated by
$\kappa(\mathbf{X},\widetilde{\mathbf{X}})$ by $\mathcal{H}_{\kappa}$,
equipped with norm $\|\cdot\|_{\mathcal{H}_{\kappa}}$, then the empirical
regularized risk minimization on RKHS for estimating $f(\mathbf{X})$ solves
the following optimization problem:

$\min_{f}{\bf
P}_{n}l(Y,f(\mathbf{X}))+\gamma_{n}\|f\|_{\mathcal{H}_{\kappa}}^{2},$

where $l(y,f)$ a pre-specified non-negative and convex loss function to
quantify the prediction performance, ${\bf P}_{n}$ denotes the empirical
measure from $n$ observations, i.e., ${\bf
P}_{n}g(Y,\mathbf{X})=n^{-1}\sum_{i=1}^{n}g(Y_{i},\mathbf{X}_{i})$, and
$\gamma_{n}$ is a tuning parameter to control the complexity of $f$. For a
continuous outcome, $l(y,f)$ is often chosen to be a $L_{2}$-loss given as
$(y-f)^{2}$, while for a binary outcome, it can be one of the large-margin
losses such as $\exp\\{-yf\\}$ in Adaboost. There are many choices of kernel
functions for $\kappa(\cdot,\cdot)$ so that the estimated $f(\mathbf{X})$ is
nonlinear. One of the most commonly used kernel functions in machine learning
is the Gaussian kernel function given by
$\kappa(\mathbf{X},\widetilde{\mathbf{X}})=\exp\left\\{-\|\mathbf{X}-\widetilde{\mathbf{X}}\|^{2}/\sigma^{2}\right\\}$
for some bandwidth $\sigma$, where $\|\cdot\|$ is the Euclidean norm. To
handle high-dimensional features, SpAM considered an additive kernel function
by assuming
$\kappa(\mathbf{X},\widetilde{\mathbf{X}})=\sum_{j=1}^{p_{n}}\exp\left\\{-|X_{j}-\widetilde{X}_{j}|^{2}/\sigma^{2}\right\\}$.
In the KNIFE procedure, the kernel function is defined as
$\kappa_{\bm{\omega}}(\mathbf{X},\widetilde{\mathbf{X}})=\exp\left\\{-{\sum_{j=1}^{p_{n}}\omega_{j}(X_{j}-\widetilde{X}_{j})^{2}}/{\sigma^{2}}\right\\}$,
where $\omega_{j},j=1,..,p_{n}$ are the additional weights to determine the
feature importance.

To achieve the goal of both nonparametric prediction and feature selection, we
propose a tensor product kernel as follows. For any given nonnegative vector
$\bm{\lambda}=(\lambda_{1},\lambda_{2},\cdots\lambda_{p_{n}})^{\intercal}$, we
define a $\bm{\lambda}$-regularized kernel function as

$\kappa_{\bm{\lambda},\sigma_{n}}(\mathbf{X},\widetilde{\mathbf{X}})=\prod_{m=1}^{p_{n}}\left\\{1+\lambda_{m}\kappa_{n}(X_{m},\widetilde{X}_{m})\right\\},$
(1)

where $\kappa_{n}(x,y)=\exp\left\\{-(x-y)^{2}/2\sigma_{n}^{2}\right\\}$ with a
pre-defined bandwidth $\sigma_{n}$ in $\mathcal{R}$. There are two important
observations for this new kernel function. First, it is a product of a
univariate kernel function for each feature variable, which is given by
$1+\lambda_{m}\kappa_{n}(X_{m},\widetilde{X}_{m})$. Thus, the RKHS generated
by $\kappa_{\bm{\lambda},\sigma_{n}}$ is equivalent to the tensor product of
the RKHS generated by each feature-specific space. Second, each univariate
kernel function is essentially the same as the Gaussian kernel function when
$\lambda_{m}\neq 0$. Consequently, the resulting tensor product space is the
same as the RKHS generated by the multivariate Gaussian kernel function from
all features whose $\lambda_{m}$’s are non-zero. Therefore, the closure for
the RKHS generated by $\kappa_{\bm{\lambda},\sigma_{n}}$ consists of all
functions that only depend on feature variables for which $\lambda_{m}\neq 0$.
In other words, non-negative parameters, $\lambda_{m}$, completely capture and
regularize the contribution of each feature $X_{m}$. In this way, the feature
selection can be achieved by estimating the regularization parameters,
$\lambda_{m}$’s, in the kernel function.

More specifically, using the proposed kernel function, we let
${\mathcal{H}_{\bm{\lambda},\sigma_{n}}}$ denote the RKHS corresponding to
$\kappa_{\bm{\lambda},\sigma_{n}}$ so we aim to minimize

$\displaystyle L_{n}(\bm{\lambda},f)\equiv$ $\displaystyle{\bf
P}_{n}l(Y,f(\mathbf{X}))+\gamma_{1n}||f||_{\mathcal{H}_{\bm{\lambda},\sigma_{n}}}^{2}+\gamma_{2n}P(\bm{\lambda})$
(2) subject to $\displaystyle
M\geq\lambda_{1},\lambda_{2},\cdots,\lambda_{p_{n}}\geq 0,$

where $M$ is a pre-specified large constant.
$P(\bm{\lambda})=\sum_{m=1}^{p_{n}}P(\lambda_{m})=\sum_{m=1}^{p_{n}}\lambda_{m}I(\lambda_{m}<M/2),$
which is a truncated Lasso, and $\gamma_{1n}$, $\gamma_{2n}$ are tuning
parameters. Here, we include an $l_{1}$ penalization term on the
regularization vector to perform feature selection and restrict $\lambda_{m}$
to be bounded. The latter bound is useful for numerical convergence to avoid
the situation that some $\lambda_{m}$ can diverge. Since our RKHS contains
constant and based on the representation theory for RKHS, solution for (2)
takes form

$f(\mathbf{X})=\sum_{i=1}^{n}\alpha_{i}\kappa_{\bm{\lambda},\sigma_{n}}(\mathbf{X},\mathbf{X}_{i})$

and

$\|f\|_{\mathcal{H}_{\bm{\lambda},\sigma_{n}}}^{2}=\bm{\alpha}^{T}\mathbf{K}_{\bm{\lambda},\sigma_{n}}\bm{\alpha},$

where $\bm{\alpha}=(\alpha_{1},...,\alpha_{n})^{T}$ and
$\mathbf{K}_{\bm{\lambda},\sigma_{n}}$ is an $n\times n$ matrix with entry
$\kappa_{\bm{\lambda},\sigma_{n}}(\mathbf{X}_{i},\mathbf{X}_{j}).$ Then the
optimization becomes solving

$\displaystyle\min_{\alpha_{1},...,\alpha_{n},\bm{\lambda}}$
$\displaystyle{\bf
P}_{n}l(Y,\sum_{i=1}^{n}\alpha_{i}\kappa_{\bm{\lambda},\sigma_{n}}(\mathbf{X},\mathbf{X}_{i}))+\gamma_{1n}\bm{\alpha}^{T}\mathbf{K}_{\bm{\lambda},\sigma_{n}}\bm{\alpha}+\gamma_{2n}\sum_{m=1}^{p_{n}}\lambda_{m}I(\lambda_{m}<M/2)$
subject to $\displaystyle
M\geq\lambda_{1},\lambda_{2},\cdots,\lambda_{p_{n}}\geq 0.$

We iterate between $\bm{\alpha}$ and $\bm{\lambda}$ to solve the above
optimization problem. At the $k$-th iteration,

$\displaystyle\bm{\alpha}^{k+1}=\min\limits_{\bm{\alpha}}n^{-1}\sum_{j=1}^{n}l(Y_{j},\sum_{i=1}^{n}\alpha_{i}\kappa_{\bm{\lambda}^{k},\sigma_{n}}(\mathbf{X}_{j},\mathbf{X}_{i}))+\gamma_{1n}\bm{\alpha}^{\intercal}\mathbf{K}_{\bm{\lambda}^{k},\sigma_{n}}\bm{\alpha}$
(3) $\displaystyle\bm{\lambda}^{k+1}=\min\limits_{0\leq\bm{\lambda}\leq
M}n^{-1}\sum_{j=1}^{n}l(Y_{j},\sum_{i=1}^{n}\alpha_{i}^{k+1}\kappa_{\bm{\lambda},\sigma_{n}}(\mathbf{X}_{j},\mathbf{X}_{i}))$
$\displaystyle\qquad\qquad\qquad+\gamma_{1n}({\bm{\alpha}^{k+1}})^{\intercal}\mathbf{K}_{\bm{\lambda},\sigma_{n}}\bm{\alpha}^{k+1}+\gamma_{2n}\sum_{m=1}^{p_{n}}\lambda_{m}I(\lambda_{m}<M/2).$
(4)

Since the loss function is a convex loss, the optimization in (3) is a convex
minimization problem, so many optimization algorithms can be applied. To solve
(4) for $\bm{\lambda}$, we adopt a coordinate descent algorithm to update each
$\lambda_{q}\ (q=1,2,\cdots,p_{n})$ in turn. Specifically, to obtain
$\lambda_{q}^{k+1}$, we fix
$\lambda_{1}^{k+1},\lambda_{2}^{k+1},\cdots,\lambda_{q+1}^{k},\lambda_{q+2}^{k},\cdots,\lambda_{p_{n}}^{k}$
and then after simple calculation, the objective function takes the following
form,

$\min\limits_{\lambda_{q}\geq
0}\frac{1}{n}\sum_{i=1}^{n}g(a_{iq}+b_{iq}\lambda_{q})+d_{q}\lambda_{q},$ (5)

where $g(\lambda_{q})$ is equal to
$l(Y_{j},\sum_{i=1}^{n}\alpha_{i}^{k+1}\kappa_{\bm{\lambda},\sigma_{n}}(\mathbf{X}_{j},\mathbf{X}_{i}))$
as a function of $\lambda_{q}$, and $a_{iq},b_{iq},d_{q}$’s are constants. By
the construction of $\kappa_{\bm{\lambda},\sigma_{n}}$, $g(\lambda_{q})$ is a
convex function so each step in the coordinating descent algorithm is a
constrained convex minimization problem in a bounded inteval, which is easy to
solve. Thus, our algorithm guarantees that the objective function decreases
over iterations and converges to a local minimum. We summarize the algorithm
in the following table. At the convergence after $k$ iterations, the final
prediction function is given as

$\widehat{f}_{\bm{\lambda}^{k+1}}(\mathbf{X})=\sum_{i=1}^{n}{\alpha}_{i}^{k+1}\kappa_{\bm{\lambda}^{k+1},\sigma_{n}}(\mathbf{X},\mathbf{X}_{i}).$

For classification problem, the classification rule is
$\textrm{sign}(\widehat{f}_{\bm{\lambda}^{k+1}}(\mathbf{X}))=\textrm{sign}(\sum_{i=1}^{n}{\alpha}_{i}^{k+1}\kappa_{\bm{\lambda}^{k+1},\sigma_{n}}(\mathbf{X},\mathbf{X}_{i})).$
We give details of our algorithm below (Algorithm 1).

Algorithm 1 Algorithm for learning $f(\mathbf{X})$

Data $(\mathbf{X},\mathbf{Y})$; Regularization parameter $\gamma_{1n}$ and
$\gamma_{2n}$; Former updating results,
$\widehat{\bm{\alpha}}^{k},\widehat{\bm{\lambda}}^{k},\widehat{f}_{\widehat{\bm{\lambda}}^{k}}$;
Initialize For regression, $\widehat{\bm{\lambda}}_{0}={\mathbf{0}}$; For
classification, $\widehat{\bm{\lambda}}_{0}=(0,\cdots,1,\cdots,0)$, where all
elements equal to 0, expect the one having largest margin correlation with
outcome. Iterate until convergence
($\delta=|L_{n}(\widehat{\bm{\lambda}}^{k+1},\widehat{f}_{\widehat{\bm{\lambda}}^{k+1}})-L_{n}(\widehat{\bm{\lambda}}^{k},\widehat{f}_{\widehat{\bm{\lambda}}^{k}})|\leq
c_{1}$,
$e={\|\widehat{\bm{\lambda}}^{k+1}-\widehat{\bm{\lambda}}^{k}\|}_{1}\leq
c_{2}$, where $c_{1}$ and $c_{2}$ are given cut points):

  * (i)

Update $\widehat{\bm{\alpha}}^{k+1}$ for fix $\widehat{\bm{\lambda}}^{k}$,
which can be solved explicitly for regression and via fminsearch function for
classification.

  * (ii)

Update $\widehat{\bm{\lambda}}^{k+1}$ for fixed $\widehat{\bm{\alpha}}^{k+1}$
via coordinate descent algorithm.

  * (iii)

$\delta=|L_{n}(\widehat{\bm{\lambda}}^{k+1},\widehat{f}_{\widehat{\bm{\lambda}}^{k+1}})-L_{n}(\widehat{\bm{\lambda}}^{k},\widehat{f}_{\widehat{\bm{\lambda}}^{k}})|$
and $e={\|\widehat{\bm{\lambda}}^{k+1}-\widehat{\bm{\lambda}}^{k}\|}_{1}$ .

$\widehat{\bm{\alpha}}^{k+1},\widehat{\bm{\lambda}}^{k+1},\widehat{f}_{\widehat{\bm{\lambda}}^{k+1}}$.

Remark 1. When updating $\bm{\alpha}$ interatively, for regression, it can be
solved in a closed form as
$\widehat{\bm{\alpha}}^{k+1}=(\mathbf{K}_{\widehat{\bm{\lambda}}^{k},\sigma_{n}}^{\intercal}\mathbf{K}_{\widehat{\bm{\lambda}}^{k},\sigma_{n}}+n\gamma_{1n}\mathbf{K}_{\widehat{\bm{\lambda}}^{k},\sigma_{n}})^{-1}\mathbf{K}_{\widehat{\bm{\lambda}}^{k},\sigma_{n}}^{\intercal}Y$.
For classification, we apply one-step Newton method for updating. Tuning
parameters in the algorithm are chosen via cross-validation over a grid of
$2^{-15},2^{-13},\cdots,2^{-13},2^{15}$. Although the kernel bandwidth,
$\sigma_{n}$, can also be tuned, to save computation cost, we follow Jaakkola
et al. (1999) to set it to be the median value of the paired distances.

##  3 Theoretical Properties

In this section, we present some theoretical properties of our proposed
method. Since our proposed kernel function is new, we first provide two
theorems that describe the properties for the RKHS generated by this kernel
function. In the first theorem, we show that this space is dense in $L_{2}(P)$
subspace consisting of all measurable functions that only depend on the
feature variables for which $\lambda_{m}\neq 0$ in the kernel function. In the
second theorem, we obtain the entropy number for the unit ball in this space.
Both theorems are necessary to establish the asymptotic properties of the
proposed estimator for $f(\mathbf{X})$ as given in the previous section.

To state our results, we define $f_{0}({\bf X})$ as the Bayesian prediction
function, which is assumed to be unique. That is, $E[l(Y,f)]$ attains its
minimum when $f=f_{0}$. We assume that feature variables
$X_{1},X_{2},\cdots,X_{q}$ are important in terms that $f_{0}({\bf X})$ is
only a function of $X_{1},X_{2},...,X_{q}$ and for any $1\leq s\leq q$,

$E\left\\{\Big{(}f_{0}({\bf X})-E\left[f_{0}({\bf
X})\Big{|}X_{1},X_{2},X_{s-1},X_{s+1}\cdots,X_{q}\right]\Big{)}^{2}\right\\}>0.$

Finally, we let $d_{2}(f_{0},\mathcal{H}_{\bm{\lambda},\sigma_{n}})$ denote
the $L_{2}(P)$-distance between $f_{0}$ and the RKHS generated by
$\kappa_{\bm{\lambda},\sigma_{n}}$.

######  Theorem 1.

For a vector $\bm{\lambda}_{n}=(\lambda_{n1},...,\lambda_{np_{n}})$ with
$\lambda_{nm}\geq 0$ for $m=1,...,p_{n}$, the following results hold:

  1. (i)

If $\lambda_{nm}>0$ for $m=1,...,q$, i.e., $\lambda_{n}$’s that are associated
with the important features are strictly positive, then
$d_{2}(f_{0},\mathcal{H}_{\bm{\lambda}_{n},{\sigma}_{n}})\to 0.$

  2. (ii)

If for some $m\leq q$, $\lambda_{nm}=0$, then $\lim\inf
d_{2}(f_{0},\mathcal{H}_{\bm{\lambda}_{n},{\sigma}_{n}})>0.$

Note: The Theorem holds for $\bm{\lambda}$ whose value depends on $n$ and
denoted as $\bm{\lambda}_{n}$.

###### Proof.

To prove (i), we first note that after expansion,
$\kappa_{\bm{\lambda}_{n},\sigma_{n}}(\mathbf{X},\widetilde{\mathbf{X}})$ is
the summation of a number of Gaussian kernels. In particular, one term of this
summation is

$\left\\{\lambda_{n1}\lambda_{n2}\cdots\lambda_{nq}{\kappa_{\sigma_{n}}(X_{1},\widetilde{X}_{1})\kappa_{\sigma_{n}}(X_{2},\widetilde{X}_{2})\cdots\kappa_{\sigma_{n}}(X_{q},\widetilde{X}_{q})}\right\\},$

where $\kappa_{\sigma}(x,y)=\exp\\{-(x-y)^{2}/\sigma^{2}\\}.$ Since
$\lambda_{n1},...,\lambda_{nq}>0$, the kernel function associated with this
term is proportional to the Gaussian kernel in the space of
$(X_{1},\cdots,X_{q})$ with bandwidth $\sigma_{n}$ for each domain $k$.
Therefore, the closure of the RKHS generated by
$\kappa_{\bm{\lambda}_{n},\sigma_{n}}$ includes the RKHS generated by the
Gaussian kernel in the space of $(X_{1},\cdots,X_{q})$. The result in (i)
holds since the latter is asymptotically dense in the subspace of $L_{2}(P)$
consisting of any functions depending on $(x_{1},...,x_{q})$.

To prove (ii), if $\lambda_{m}=0$, then it is clear that any function in
$\mathcal{H}_{\bm{\lambda}_{n},{\sigma}_{n}}$ only depends on the feature
variables except $X_{m}$. Therefore,

$\mathcal{H}_{\bm{\lambda}_{n},{\sigma}_{n}}\subset\left\\{g({\bf
X}_{-m}):g\in L_{2}(P)\right\\},$

where ${\bf X}_{-m}$ denotes all the feature variables excluding $X_{m}$. On
the other hand, the projection of $f_{0}$ on the latter space is $E[f_{0}|{\bf
X}_{-m}]$. Therefore,

$\lim\inf d(f_{0},\mathcal{H}_{\bm{\lambda}_{n},{\sigma}_{n}})\geq
d(f_{0},E[f_{0}|{\bf X}_{-m}])>0$

since $X_{m}$ is one important variable for $f_{0}$. We obtain the result. ∎

Our next theorem studies the bracket covering number for a unit ball in ${\cal
H}_{\bm{\lambda}_{n},\sigma_{n}}$. We consider ${\cal B}_{n}$ as the unit ball
in ${\cal H}_{\bm{\lambda}_{n},\sigma_{n}}$, i.e., ${\cal
B}_{n}\equiv\left\\{f({\bf x}):\|f\|_{{\cal
H}_{\bm{\lambda}_{n},\sigma_{n}}}\leq 1\right\\},$ Then the $\epsilon$-bracket
covering number for ${\cal B}_{n}$, denoted as
$N_{[]}(\epsilon,\mathcal{B}_{n},\|\cdot\|_{L_{2}(P)})$, is defined as the
minimal number of pairs $[l({\bf x}),u({\bf x})]$ such that any function
$\|u(\mathbf{X})-l(\mathbf{X})\|_{L_{2}(P)}\leq\epsilon$ and any function $f$
in ${\cal B}_{n}$ is between one pair, i.e., $l({\bf x})\leq f({\bf x})\leq
u({\bf x})$.

######  Theorem 2.

For a vector $\bm{\lambda}_{n}=(\lambda_{n1},...,\lambda_{np_{n}})$ such that
$\lambda_{nm}$ is uniformly bounded by a constant $M$ for $m=1,...,q$ and
$\lambda_{n(q+1)}=...=\lambda_{np_{n}}=0$, it holds

$\log\mathcal{N}_{[]}(\epsilon,\mathcal{B}_{n},\|\cdot\|_{L_{2}(P)})\leq
C\sigma_{n}^{-(1-v/4)q}\epsilon^{-v},$

where $v$ is any constant within $(0,2)$ and $C$ only depends on $M$ and $q$.

###### Proof.

For any $f\in{\cal B}_{n}$ with form

$f({\bf
x})=\sum_{i=1}^{\infty}\alpha_{i}\kappa_{\bm{\lambda}_{n},\sigma_{n}}({\bf
x},{\bf x}_{i}),$

where ${\bf x}_{1},{\bf x}_{2},...$ are a sequence of given points. Using the
expansion of $\kappa_{\bm{\lambda}_{n},\sigma_{n}}$, we have

$\displaystyle f({\bf x})$ $\displaystyle=$
$\displaystyle\sum_{\\{k_{1},...,k_{s}\\}\subset\\{1,...,q\\}\cup\phi}\lambda_{nk_{1}}\cdots\lambda_{nk_{s}}\sum_{i=1}^{\infty}\alpha_{i}\exp\left\\{{-\frac{(x_{ik_{1}}-x_{k_{1}})^{2}+\cdots+(x_{ik_{s}}-x_{k_{s}})^{2}}{\sigma_{n}^{2}}}\right\\}$
$\displaystyle=$
$\displaystyle\sum_{\\{k_{1},...,k_{s}\\}\subset\\{1,...,q\\}\cup\phi}\sqrt{\lambda_{nk_{1}}\cdots\lambda_{nk_{s}}}f_{k_{1}...k_{s}}(\mathbf{x}),$

where $x_{ik}$ and $x_{k}$ are respectively the $k$th component of ${\bf
x}_{i}$ and ${\bf x}$, and

$f_{k_{1}...k_{s}}(\mathbf{x})=\sum_{i=1}^{\infty}\alpha_{i}\sqrt{{\lambda}_{nk_{1}}\cdots{\lambda}_{nk_{s}}}\exp\left\\{-\frac{(x_{ik_{1}}-x_{k_{1}})^{2}+\cdots+(x_{ik_{s}}-x_{k_{s}})^{2}}{\sigma_{n}^{2}}\right\\}.$

Here, if the index set if empty, then the exponential part in the summation is
replaced by 1.

Clearly, if we denote ${\cal H}_{k_{1}...k_{s}}$ as the reproducing kernel
Hilbert space generated by the Gaussian kernel
$\exp\left\\{-{[(\widetilde{x}_{k_{1}}-x_{k_{1}})^{2}+\cdots+(\widetilde{x}_{k_{s}}-x_{k_{s}})^{2}]}/{\sigma_{n}^{2}}\right\\},$
then $f_{k_{1}...k_{s}}(\mathbf{x})\in{\cal H}_{k_{1}...k_{s}}$ and moreover,

$\displaystyle\|f\|_{\mathcal{H}_{{\bm{\lambda}_{n}},\sigma_{n}}}^{2}$
$\displaystyle=\sum_{i=1}^{\infty}\sum_{j=1}^{\infty}\alpha_{i}\alpha_{j}\kappa_{{\bm{\lambda}_{n},\sigma_{n}}}(\mathbf{x}_{i},\mathbf{x}_{j})$
$\displaystyle=\sum_{i=1}^{n}\sum_{j=1}^{n}\alpha_{i}\alpha_{j}\sum_{\\{k_{1},...,k_{s}\\}\subset\\{1,...,q\\}\cup\phi}{\lambda}_{nk_{1}}\cdots{\lambda}_{nk_{s}}\exp\left\\{-\frac{(x_{ik_{1}}-x_{jk_{1}})^{2}+\cdots+(x_{ik_{s}}-x_{jk_{s}})^{2}}{\sigma_{n}^{2}}\right\\}$
$\displaystyle=\sum_{\\{k_{1},...,k_{s}\\}\subset\\{1,...,q\\}\cup\phi}\sum_{i=1}^{\infty}\sum_{j=1}^{\infty}\alpha_{i}\alpha_{j}{\lambda}_{nk_{1}}\cdots{\lambda}_{nk_{s}}\exp\left\\{-\frac{(x_{ik_{1}}-x_{jk_{1}})^{2}+\cdots+(x_{ik_{s}}-x_{jk_{s}})^{2}}{\sigma_{n}^{2}}\right\\}$
$\displaystyle=\sum_{\\{k_{1},...,k_{s}\\}\subset\\{1,...,q\\}\cup\phi}\|f_{k_{1}...k_{s}}\|_{\mathcal{H}_{k_{1}...k_{s}}}^{2}.$

Thus, $\|f\|_{\mathcal{H}_{\bm{\lambda}_{n},\sigma_{n}}}\leq 1$ implies
$\|f_{k_{1}...k_{s}}\|_{\mathcal{H}_{k_{1}...k_{s}}}\leq 1$ for any
$k_{1},...,k_{s}$.

Consequently, since such $f$ is dense in $\mathcal{B}_{n}$, we conclude

$\mathcal{B}_{n}\subseteq\overline{\left\\{\sum_{\\{k_{1},...,k_{s}\\}\subset\\{1,...,q\\}\cup\phi}f_{k_{1}...k_{s}}({\bf
x})\sqrt{{\lambda}_{nk_{1}}\cdots{\lambda}_{nk_{s}}}:\|f_{k_{1}...k_{s}}\|_{\mathcal{H}_{k_{1}...k_{s}}}^{2}\leq
1\right\\}}.$

Thus, there exists a constant $C$ only depending on $M$ and $q$ such that

$\log\mathcal{N}_{[]}(2^{q}M^{q/2}\epsilon,\mathcal{B}_{n},\|\cdot\|_{L_{2}(P)})\leq{\sum_{\\{k_{1},...,k_{s}\\}\subset\\{1,...,q\\}\cup\phi}}\log\mathcal{N}_{[]}(\epsilon,\\{f_{k_{1}...k_{s}}({\bf
x}),\|f_{k_{1}...k_{s}}\|_{\mathcal{H}_{k_{1}...k_{s}}}\leq
1\\},\|\cdot\|_{L_{2}(P)})$

According to (Steinwart and Scovel (2007)), we know

$\log\mathcal{N}_{[]}(\epsilon,\\{f_{k_{1}...k_{s}}({\bf
x}),\|f_{k_{1}...k_{s}}\|_{\mathcal{H}_{k_{1}...k_{s}}}^{2}\leq
1\\},\|\cdot\|_{L_{2}(P)})\leq C\sigma_{n}^{-(1-v/4)s}\epsilon^{-v},$

for any constant $v\in(0,2)$ and a constant $C$ only depending on $s$.
Therefore,

$\log\mathcal{N}(\epsilon,\mathcal{B}_{n},\|\cdot\|_{L_{2}(P)})\leq
C(M,q)\sum_{\\{k_{1},...,k_{s}\\}\subset\\{1,...,q\\}\cup\phi}\sigma_{n}^{-(1-v/4)s}\epsilon^{-v}\leq
C(M,q)\sigma_{n}^{-(1-v/4)q}\epsilon^{-v}$

for a constant $C(M,q)$. We have proved Theorem 2. ∎

Our next theorem gives the main properties of the estimated prediction
function. We show that the resulting prediction function from our method leads
to Bayesian risk asymptotically. Moreover, with probability tending to one,
the variable selection based on non-zero $\lambda_{n}$’s is oracle as if we
knew which variables were important. Recall that
$(\widehat{\bm{\lambda}_{n}},\widehat{f})$ is the optimal solution of the
objective function

$L_{n}(\bm{\lambda}_{n},f)={\bf
P}_{n}\mathnormal{l}(Y,f(\mathbf{X}))+\gamma_{1n}\|f\|_{\mathcal{H}_{\bm{\lambda}_{n},\bm{\sigma}_{n}}}^{2}+\gamma_{2n}P(\bm{\lambda}_{n}),$
(6)

where $P(\bm{\lambda}_{n})$ is the truncated Lasso penalty for
$\bm{\lambda}_{n}$. Equivalently, if we define for any $\bm{\lambda}_{n}$,

$\widehat{f}_{\bm{\lambda}_{n}}=\textrm{arg}\min_{f}L_{n}(\bm{\lambda}_{n},f),$

which exists due to the convexity of $L_{n}(\bm{\lambda}_{n},f)$ in $f$, then
$\widehat{\bm{\lambda}}$ minimizes
$L_{n}(\widehat{\bm{\lambda}}_{n},\widehat{f}_{\bm{\lambda}_{n}})$ and
$\widehat{f}=\widehat{f}_{\widehat{\bm{\lambda}}_{n}}$.

For the main theorem, we assume $(Y,\mathbf{X})$ to have a bounded support and
need the following conditions.
(C1). The loss function $l(y,f)$ is convex and is Lipschtisz continuous with
respect to $f$ in any bounded set.
(C2). There exit $\delta>0$ and a constant $c_{1}>0$ such that

$E[l(Y,f({\bf X}))-l(Y,f_{0}({\bf X}))]\geq c_{1}\|f({\bf X})-f_{0}({\bf
X})\|_{L_{2}(P)}^{2}$

whenever $E[l(Y,f({\bf X}))-l(Y,f_{0}({\bf X}))]$ is smaller than $\delta$.
(C3). Assume
$\|l_{2}(Y,f(\mathbf{X}))-l_{2}(Y,f_{0}(\mathbf{X}))\|_{L_{2}(P)}\leq
c_{2}\|f({\bf X})-f_{0}({\bf X})\|_{L_{2}(P)}$ for a constant $c_{2}$, where
$l_{2}(y,x)=\partial l(y,x)/\partial x$.
(C4). For any
$\widetilde{\bm{\lambda}_{n}}=(\lambda_{n1},...,\lambda_{np_{n}})$ such that
$\lambda_{nk}=0$ for $k>q$, let $\Lambda_{\max}(\mathbf{X}_{-q})$ and
$\Lambda_{\min}(\mathbf{X}_{-q})$ be the largest and smallest eigenvalues of
the matrix
$\left(E[K_{\widetilde{\bm{\lambda}}_{n}}(\mathbf{X}_{j},\mathbf{X})K_{\widetilde{\bm{\lambda}}_{n}}(\mathbf{X}_{l},\mathbf{X})|\mathbf{X}_{-q}]\right)$
where $\mathbf{X}_{-q}$ denotes all unimportant variables. We assume that with
probability one, there exists one constant $c$ such that
$\Lambda_{\max}(\mathbf{X}_{-q})/\Lambda_{\min}(\mathbf{X}_{-q})\leq
c\sigma_{n}^{-1/2}$ and
$E[\Lambda_{\min}(\mathbf{X}_{-q})\kappa_{n}(x,X_{m})^{2}]\leq
c\sigma_{n}^{1/2}$ for any $m>q$.
(C5). Assume $\log p_{n}=o(n^{1-(2+q)\alpha_{1}-\alpha_{2}-\alpha_{3}})$.
Moreover, we assume $\sigma_{n}=n^{-\alpha_{1}}$,
$\gamma_{1n}=n^{-\alpha_{2}},\gamma_{2n}=n^{-\alpha_{3}}$, where
$\alpha_{k}>0$ for $k=1,2,3$ and they satisfy

(i) $1-(2+q)\alpha_{1}-\alpha_{2}>0$

(ii)
$0<\alpha_{3}<\min{\Big{(}\frac{1}{4}(1+\frac{\alpha_{1}q}{2}+\alpha_{2}),1-(2+q)\alpha_{1}-\alpha_{2},\frac{\alpha_{1}}{2},\frac{\alpha_{2}}{2}\Big{)}}$.

Conditions (C1)-(C3) give the assumptions for the loss functions. It can be
verified that they hold for $l(y,f)=(y-f)^{2}$ for a continuous $y$ and for
$l(y,f)=\exp(-yf)$ for a binary $y$. Condition (C4) implies the equivalence
between the Euclidean norm of the coefficients and the reproducing kernel
Hilbert space norm, up to a scale proportional to $\sigma_{n}^{-1/2}$. The
second half of the condition in (C4) holds automatically if the important
variables are independent of the unimportant variable when
$\Lambda_{\min}(\mathbf{X}_{-q})$ does not depend on $\mathbf{X}_{-q}$. We
note that such a condition is analogue to the design matrix condition assumed
in high dimensional linear model literature. Finally, condition (C5) allows
the dimensionality of the feature variable to be ultra-high and imposes
additional constraints for the choices of the bandwidth and two tuning
parameters.

######  Theorem 3.

Under Conditions (C1)-(C5), there exists a local minimizer
$\widehat{\bm{\lambda}}_{n}$ for
$L_{n}(\bm{\lambda}_{n},\widehat{f}_{\bm{\lambda}_{n}})$ such that with
probability tending to one,
(a) $E[l\big{(}Y,\widehat{f}_{\widehat{\bm{\lambda}}_{n}}\big{)}]$ converges
to $E[l\big{(}Y,f_{0}\big{)}]$.
(b) For $m=1,...,q,\widehat{\lambda}_{nm}>0$.
(c) For $m=q+1,q+2,\cdots,p_{n},\ \widehat{\lambda}_{nm}=0$.

The first part of Theorem 3 implies that the loss of the estimated prediction
function converges to the Bayes risk. The last two conclusions in Theorem 3
show that the $\widehat{\lambda}_{nm}$’s associated with important feature
variables should be non-zero, i.e., the estimated function does depend on
important variables. More importantly, the proposed method can estimate the
predicted function as if we knew which variables are important in the truth.
The proof for Theorem 3 is given in the supplementary file. The proof of
Theorem 3(a) entails careful examination of the stochastic variability of
$L_{n}(\bm{\lambda}_{n},\widehat{f}_{\bm{\lambda}_{n}})$, for which we first
establish a preliminary bound for $\widehat{f}_{\bm{\lambda}_{n}}$ and then
appeal to some concentration inequalities for empirical processes with metric
entropy as derived from Theorem 2. To prove Theorem 3(b) and (c) in the
theorem, we examine the KKT conditions to show that the oracle estimators,
i.e., $\lambda_{nm}$ is known to be zero for $m>q$, satisfies the KKT
conditions with probability tending to one. Again, concentration inequalities
for empirical processes are needed in technical arguments in the proof.

##  4 Simulation Study

We conducted two simulation studies, one for a regression problem with
continuous $Y$ and the other for classification with binary $Y$. In the first
simulation study, we considered a continuous outcome model with total number
of $p$ correlated feature variables, which were generated from a multivariate
normal distribution, each with mean zero and variance one. Furthermore,
$X_{1},X_{2},X_{3},X_{4}$ were correlated with $corr(X_{1},X_{2})=0.4,\
corr(X_{1},X_{3})=-0.3,\ corr(X_{2},X_{3})=0.5$ and $corr(X_{3},X_{4})=0.2$,
while the others were all independent. The outcome variable, $Y$, was
simulated from a linear model

$Y=0.9X_{5}^{3}+4X_{1}X_{2}X_{3}+2.3\exp(-X_{3})+4X_{4}+\epsilon,$

where $\epsilon\sim N(0,1)$. Thus, $X_{1}$ to $X_{5}$ were important variables
but not any others. In the second simulation study, $X$’s were generated
similarly but with some different correlations: $corr(X_{1},X_{2})=-0.2,\
corr(X_{1},X_{4})=0.2,\ corr(X_{2},X_{3})=0.5,\ corr(X_{3},X_{4})=0.3$ and
$corr(X_{3},X_{4})=-0.4$. The binary outcome, $Y$, with values $-1$ and $1$,
were generated from a Bernoulli distribution with the probability of being one
given by

$\Big{\\{}1+e^{-0.25+(X_{2}-1.1X_{3}+0.3X_{4})^{3}}\Big{\\}}^{-1},$

so only $X_{2}$ to $X_{4}$ were important variables. Since many biomedical
applications (as well as our application in this work) have small to moderate
sample sizes, in both simulation studies, we considered sample size
$n=100,200$ and $400$ and varied the feature dimension from $p=200,400$ to
$1000$. Each simulation setting was repeated 500 times.

For each simulated data, we used the proposed method to learn the prediction
function. Initial values, tuning parameters and optimization package used for
binary case are chosen as in Remark 1 of Section 2, where 3-fold cross-
validation was used for selecting the tuning parameters. The bound of
regularized parameter $M$ was chosen to be $10^{5}$. We also centerized
continuous outcome and re-weighted class label controlled to be balanced
before iteration to make numerical stable. We reported the true positive
rates, true negative rates and the average number of the selected variables
for feature selection. We also reported the prediction errors or
misclassification rates using a large and independent validation data. For
comparison, we compared our proposed method with HSICLasso and SpAM since both
methods were able to estimate nonlinear functions in high dimensional
settings. In addition, we also compared the performance with LASSO in the
first simulation study and $l_{1}$-SVM in the second simulation study, in
order to study the impact due to model misspecification.

The results based on 500 replicates are summarized in Tables 1 and 2. From
these tables, we observe that for fixed dimension, the performance of our
method improves as sample size $n$ becomes large in terms of the improved true
positive and true negative rates for feature selection as well as decreasing
prediction errors. In almost all cases, our true negative rate is close to
100$\%$, which shows that noise variables can be identified with a very high
chance. As expected, the performance deteriorates as the dimensionality
increases. Interestingly, our method continues to select only a small number
of feature variables. Comparatively, HSICLasso selected many more noise
variables and had larger prediction errors, while SpAM also tended to select
more features than our method. The performance of these methods become much
worse when the feature dimension is 1000. Clearly, LASSO and $l_{1}$-SVM did
not yield reasonable variable selection results and their prediction errors
are much higher due to model misspecification. We also give boxplots to
visualize prediction performance of 500 replications in Figures 1 and 2. Since
Lasso cannot provide stable prediction errors, its prediction errors from many
replicates are out of the bound as shown in Figure 1. Figure 1 and 2 further
confirm that our method is superior to all other methods, even when the
dimension is as large as 1000 and the sample size is as small as $n=100$,
which is of similar size as our real data analysis example in Section 5.

Table 1: Results from The Simulation Study with Continuous Outcome

| (a) Summary of Feature Selection Performance
---|---
|  | Proposed Method |  | HSICLasso |  | SPAM |  | LASSO
$p$ | $n$ | TPR | TNR | Avg$\\#$ |  | TPR | TNR | Avg$\\#$ |  | TPR | TNR | Avg$\\#$ |  | TPR | TNR | Avg$\\#$
100 | 100 | 60.9$\%$ | 97.3$\%$ | 5.6 |  | 81.5$\%$ | 78.5$\%$ | 24.5 |  | 99.6$\%$ | 34.6$\%$ | 67.1 |  | 98.8$\%$ | 1.3$\%$ | 98.8
| 200 | 71.2$\%$ | 99.0$\%$ | 4.5 |  | 98.0$\%$ | 60.4 $\%$ | 42.5 |  | 100.0$\%$ | 4.4$\%$ | 95.8 |  | 100.0$\%$ | 0.1$\%$ | 99.9
| 400 | 82.7$\%$ | 98.4$\%$ | 5.7 |  | 99.6 $\%$ | 78.0 $\%$ | 25.8 |  | 100.0$\%$ | 0.3$\%$ | 99.7 |  | 100.0$\%$ | 0.1$\%$ | 99.9
200 | 100 | 57.2$\%$ | 98.7$\%$ | 5.5 |  | 75.6$\%$ | 88.8$\%$ | 25.6 |  | 99.1$\%$ | 63.0$\%$ | 77.1 |  | 84.0$\%$ | 52.2$\%$ | 97.5
| 200 | 66.6$\%$ | 99.5$\%$ | 4.2 |  | 94.0$\%$ | 75.2$\%$ | 53.1 |  | 100.0$\%$ | 33.7$\%$ | 134.1 |  | 99.1$\%$ | 0.0$\%$ | 198.6
| 400 | 78.1$\%$ | 99.4$\%$ | 5.0 |  | 99.8 $\%$ | 84.2$\%$ | 35.8 |  | 100.0$\%$ | 5.4$\%$ | 189.5 |  | 100.0$\%$ | 0.12$\%$ | 199.8
400 | 100 | 47.3$\%$ | 99.3$\%$ | 5.2 |  | 68.5$\%$ | 90.4$\%$ | 41.5 |  | 98.2$\%$ | 80.8$\%$ | 80.8 |  | 79.4$\%$ | 76.4$\%$ | 97.1
| 200 | 65.0$\%$ | 99.7$\%$ | 4.5 |  | 86.3$\%$ | 89.0$\%$ | 47.6 |  | 100.0$\%$ | 62.1$\%$ | 154.6 |  | 90.7$\%$ | 51.4$\%$ | 196.6
| 400 | 73.1$\%$ | 99.8$\%$ | 4.4 |  | 99.7$\%$ | 87.6$\%$ | 54.0 |  | 100.0$\%$ | 34.0$\%$ | 265.8 |  | 99.1$\%$ | 0.7$\%$ | 397.3
1000 | 100 | 40.7$\%$ | 99.7$\%$ | 5.0 |  | 56.0$\%$ | 91.8$\%$ | 84.5 |  | 93.7$\%$ | 92.2$\%$ | 82.6 |  | 73.6$\%$ | 90.6$\%$ | 97.2
| 200 | 61.2$\%$ | 99.9$\%$ | 4.5 |  | 78.2$\%$ | 98.6$\%$ | 18.4 |  | 99.9$\%$ | 84.2$\%$ | 162.0 |  | 85.5$\%$ | 80.7$\%$ | 196.1
| 400 | 70.7$\%$ | 99.9$\%$ | 4.0 |  | 99.4$\%$ | 91.0$\%$ | 94.9 |  | 100.0$\%$ | 68.7$\%$ | 316.4 |  | 94.5$\%$ | 60.7$\%$ | 395.6

|  |  | (b) Summary of Prediction Errors
---|---|---|---
$p$ | $n$ |  | Proposed Method | HSICLasso | SPAM | LASSO
100 | 100 |  | 7.405 (0.527) | 7.695 (0.291) | 6.985 (0.291) | 41.663 (11.325)
| 200 |  | 5.929 (0.950) | 7.323 (0.098) | 7.299 (0.389) | 9.508 (0.575)
| 400 |  | 4.424 (0.777) | 7.115 (0.053) | 6.868 (0.292) | 7.840 (0.183)
200 | 100 |  | 7.567 (0.493) | 7.603 (0.286) | 6.672 (0.336) | 10.176 (0.794)
| 200 |  | 6.623 (0.412) | 7.313 (0.130) | 6.404 (0.279) | 44.464 (9.125)
| 400 |  | 5.661(0.580) | 6.946 (0.054) | 6.767 (0.305) | 9.370 (0.433)
400 | 100 |  | 7.920 (0.670) | 8.001 (0.284) | 6.815 (0.399) | 9.091 (0.532)
| 200 |  | 7.008 (0.346) | 7.563 (0.233) | 6.222 (0.218) | 10.151 (0.722)
| 200 |  | 6.444 (0.199) | 7.061 (0.049) | 6.079 (0.192) | 40.190 (6.402)
1000 | 100 |  | 8.215 (0.764) | 8.638 (0.263) | 7.067 (0.372) | 8.851 (0.406)
| 200 |  | 7.324 (0.368) | 7.539 (0.252) | 6.214 (0.242) | 8.871 (0.379)
| 400 |  | 6.818 (0.250) | 7.376 (0.068) | 5.870 (0.161) | 9.652 (0.429)

Note. In (a), “TPR" is the true positive rate, “TNR" is the true negative
rate, and “Avg$\\#$" is the average number of the selected variables from 500
replicates. In (b), the numbers are the mean squared errors from prediction,
and the numbers within parentheses are the median absolute deviations from 500
replicates.

Table 2: Results from The Simulation Study with Binary Outcome

| (a) Summary of Feature Selection Performance
---|---
|  | Proposed Method |  | HSICLasso |  | SPAM |  |  $l_{1}$-SVM
$p$ | $n$ | TPR | TNR | Avg$\\#$ |  | TPR | TNR | Avg$\\#$ |  | TPR | TNR | Avg$\\#$ |  | TPR | TNR | Avg$\\#$
100 | 100 | 74.7$\%$ | 99.0$\%$ | 3.3 |  | 71.1$\%$ | 79.4$\%$ | 22.1 |  | 64.5$\%$ | 89.6$\%$ | 12.1 |  | 76.2$\%$ | 75.1$\%$ | 26.5
| 200 | 83.9$\%$ | 99.9$\%$ | 2.6 |  | 80.7$\%$ | 89.7 $\%$ | 12.4 |  | 53.4$\%$ | 98.9$\%$ | 2.6 |  | 92.5$\%$ | 80.4$\%$ | 21.8
| 400 | 86.0$\%$ | 99.9$\%$ | 2.6 |  | 87.8$\%$ | 90.3$\%$ | 12.1 |  | 50.6$\%$ | 99.9$\%$ | 1.5 |  | 98.8$\%$ | 71.3$\%$ | 30.8
200 | 100 | 70.4$\%$ | 99.3$\%$ | 3.5 |  | 71.3$\%$ | 80.1$\%$ | 41.3 |  | 63.6$\%$ | 91.2$\%$ | 19.2 |  | 71.3$\%$ | 85.4$\%$ | 31.0
| 200 | 84.1$\%$ | 99.8$\%$ | 2.9 |  | 78.3$\%$ | 95.0 $\%$ | 12.2 |  | 54.5$\%$ | 98.7$\%$ | 4.1 |  | 90.7$\%$ | 80.2$\%$ | 41.8
| 400 | 87.0$\%$ | 100.0$\%$ | 2.7 |  | 83.1 $\%$ | 96.5$\%$ | 9.3 |  | 50.6$\%$ | 99.9$\%$ | 1.6 |  | 89.3$\%$ | 73.4$\%$ | 55.0
400 | 100 | 68.5$\%$ | 99.5$\%$ | 3.9 |  | 70.9$\%$ | 79.4$\%$ | 84.0 |  | 63.7$\%$ | 92.8$\%$ | 30.6 |  | 65.9$\%$ | 86.7$\%$ | 54.7
| 200 | 84.5$\%$ | 99.9$\%$ | 3.0 |  | 76.9$\%$ | 95.5 $\%$ | 20.2 |  | 57.7$\%$ | 98.3$\%$ | 8.3 |  | 87.0$\%$ | 91.0$\%$ | 38.1
| 400 | 87.0$\%$ | 100.0$\%$ | 2.6 |  | 79.6 $\%$ | 98.9$\%$ | 6.8 |  | 51.9$\%$ | 100.0$\%$ | 1.8 |  | 99.1$\%$ | 82.3$\%$ | 73.0
1000 | 100 | 61.3$\%$ | 99.8$\%$ | 4.1 |  | 72.2$\%$ | 77.4$\%$ | 227.4 |  | 61.0$\%$ | 95.7 $\%$ | 45.0 |  | 58.4$\%$ | 90.3 $\%$ | 98.9
| 200 | 86.3$\%$ | 99.9$\%$ | 3.3 |  | 75.5 $\%$ | 95.9 $\%$ | 43.6 |  | 54.3$\%$ | 98.9 $\%$ | 12.8 |  | 79.4$\%$ | 91.4$\%$ | 87.4
| 400 | 87.7$\%$ | 100.0$\%$ | 2.8 |  | 73.9 $\%$ | 99.6 $\%$ | 6.6 |  | 50.0$\%$ | 100.0$\%$ | 1.9 |  | 96.8$\%$ | 90.1 $\%$ | 101.6

|  |  | (b) Summary of Misclassification Errors
---|---|---|---
$p$ | $n$ |  | Proposed Method | HSICLasso | SPAM |  $l_{1}$-SVM
100 | 100 |  | 0.314 (0.017) | 0.345 (0.028) | 0.343 (0.018) | 0.359 (0.032)
| 200 |  | 0.290 (0.009) | 0.307 (0.012) | 0.312 (0.002) | 0.305 (0.011)
| 400 |  | 0.283 (0.004) | 0.292 (0.012) | 0.297 (0.002) | 0.292 (0.007)
200 | 100 |  | 0.316 (0.019) | 0.351 (0.042) | 0.344 (0.034) | 0.352 (0.031)
| 200 |  | 0.280 (0.008) | 0.302 (0.015) | 0.302 (0.003) | 0.321 (0.028)
| 400 |  | 0.270 (0.004) | 0.282 (0.014) | 0.297 (0.002) | 0.326 (0.025)
400 | 100 |  | 0.331 (0.024) | 0.372 (0.047) | 0.369 (0.046) | 0.390 (0.031)
| 200 |  | 0.286(0.010) | 0.319 (0.018) | 0.311 (0.003) | 0.327 (0.026)
| 200 |  | 0.277 (0.004) | 0.288 (0.014) | 0.305 (0.001) | 0.295 (0.010)
1000 | 100 |  | 0.352 (0.027) | 0.397 (0.037) | 0.390 (0.036) | 0.416 (0.027)
| 200 |  | 0.287 (0.008) | 0.335 (0.024) | 0.315 (0.003) | 0.381 (0.020)
| 400 |  | 0.277 (0.004) | 0.294 (0.008) | 0.305 (0.001) | 0.353 (0.016)

Note. See Table 1.

Figure 1: Boxplots of Prediction Errors for Continuous Outcome

Note. The plots give the distribution of prediction errors among four
competing methods. The comparing methods from left to right in each plot are
our proposed method, HSICLasso, SpAM and Lasso.

Figure 2: Boxplots of Misclassification Errors for Binary Outcome

Note. The plots give the distribution of misclassification rates among four
competing methods. The comparing methods from left to right in each plot are
our proposed method, HSICLasso, SpAM and $l_{1}$-SVM.

##  5 Application

We applied our proposed method to analyze a gene expression study in Scheetz
et al. (2006). This study analyzed microarrays RNAs of eye disease from 120
male rats, containing the expression levels from about $31,000$ gene probes.
One interesting question was to determine which probes might be associated
with the expression of gene TRIM32, which had been implicated in a number of
diverse biological pathways and also known to be one of 14 genes linked to
Bardet-Biedl syndrome (Locke et al., 2009). For this purpose, we dichotomized
TRIM32 based on whether it was over expressed as compared to a reference
sample in the dataset. We further restricted our feature variables to the top
1000 probe sets that were most correlated with TRIM32. All feature variables
were on a log-scale and standardized in the analysis. To examine the
performance of our method, we randomly divided the whole sample so that $70\%$
were used for training and the rest were used for testing. This random
splitting was then repeated 500 times to obtain reliable results. For each
training data, we used 3-fold cross validation to choose tuning parameters. We
also applied HSICLasso, SpAM and $l_{1}$-SVM for comparison.

The analysis results are shown in Table 3. We notice that our method gives
almost the same classification error as $l_{1}$-SVM, which is the smallest on
average. However, our method selects a much smaller set of feature variables
with an average of 5 variables. SpAM selects 13 variables on average but its
classification error is higher. In Table 4, we report the top 10 most-frequent
selected features among all 500 replications for each method. We notice that
some features such as Fbxo7 and LOC102555217 were selected by at least three
methods. In addition, Gene Sirt 3 was identified by all three nonlinear
feature selection methods, but not $l_{1}$-SVM, indicating some possible
nonlinear relationship between Sirt 3 and TRIM32. In fact, Figure 3 reveals
some nonlinear relationship between Sirt 3 and Fbxo7 using 5-Nearest-Neighbors
model. Our method also selected some genes that were not identified by any
other method. We applied our method to analyze the whole sample and obtained a
training error of $21.9\%$ along five 5 genes identified (Fbxo7, Plekha6,
Nfatc4, 1375872 and 1388656), which were all selected as the top 10 genes in
the previous random splitting experiment.

Table 3: Summary of Feature Selection Results in The Real Data Application

| min # | max # | avg # | classification error
---|---|---|---|---
Proposed Method | 2 | 13 | 5.1 | 0.286 (0.057)
HSICLasso | 1 | 1000 | 250.3 | 0.293 (0.046)
SpAM | 1 | 26 | 12.3 | 0.316 (0.057)
$l_{1}$-SVM | 7 | 990 | 448.7 | 0.283 (0.058)

Note. The numbers are the mean of misclassification rates from 500 replicates.
The numbers within parentheses are the median absolute deviations from 500
replicates. “min$\\#$" is the minimum number of the selected features,
“max$\\#$" is the max number of the selected features, and “avg.$\\#$" is the
average number of the selected features.

Table 4: Top 10 Most Selected Genes for Each Method Based on 500 Random
Splittings

Proposed Method | HSICLasso | SpAM |  $l_{1}$-SVM
---|---|---|---
Fbxo7 (67.5$\%$) |  Ska1 (76.6$\%$) |  1388491 (46.2$\%$) |  1376747 (99.1$\%$)
Plekha6 (47.3$\%$) |  Sirt3 (76.2$\%$) |  Fbxo7 (37.8$\%$) |  1390538 (98.9$\%$)
LOC102555217 (24.5$\%$) |  Ddx58 (76.2$\%$) |  Slco1c1 (36.6$\%$) |  RragB (98.6$\%$)
Nfatc4 (22.7$\%$) |  1371610 (76.0$\%$) |  Stmn1 (35.4$\%$) |  Atl1 (97.9$\%$)
1390538 (20$\%$) |  LOC100912578 (73.2$\%$) |  1373944 (32.4$\%$) |  Fbxo7 (97.3$\%$)
1375872 (20$\%$) |  Ttll7 (70.4$\%$) |  Ufl1 (32.2$\%$) |  Plekha6 (95.1$\%$)
RGD1306148 (13.4$\%$) |  Decr1 (70.4$\%$) |  LOC100912578 (31.0$\%$) |  1375872 (94.8$\%$)
Sirt3 (11.6$\%$) |  Mff (68.0$\%$) |  LOC100911357 (28.6$\%$) |  RGD1306148 (94.1$\%$)
Prpsap2 (11.4$\%$) |  Pkn2 (67.0$\%$) |  LOC102555217 (26.8$\%$) |  Ska1 (93.6$\%$)
1388656 (10.2$\%$) |  Taf11 (65.0$\%$) |  Sirt3 (22.4$\%$) |  LOC102555217 (93.2$\%$)

Note. The numbers within parentheses are the frequencies to be selected in 500
random splittings. The genes also selected by the proposed method are
highlighted in boldface.

Figure 3: 5-Nearest-Neighbor Plot of Sirt3 versus Fbxo7 in Real Data Study

##  6 Discussion

In this work, we have proposed a general framework for nonparametric feature
selection for both regression and classification in high dimensional settings.
We introduced a novel tensor product kernel for empirical risk minimization.
This kernel led to fully nonparametric estimation for the prediction function
but allowed the importance of each feature to be captured by a non-negative
parameter in the kernel function. Our approach is computationally efficient
because it iteratively solves a convex optimization problem in a coordinate
descent manner. We have shown that the proposed method has theoretical oracle
property for variable selection. The superior performance of the proposed
method was demonstrated via simulation studies and a real data application
with a large number of feature variables.

We considered $l_{2}$ loss function for regression and exponential loss
function for classification as examples. Clearly, the proposed framework
applies to feature selection under many different loss functions in machine
learning field. Another extension is to incorporate structures of feature
variables in constructing the kernel function. For example, in integrative
data analysis, feature variables arise from many different domains such as
clinical domain, DNA, RNA, imaging and nutrition. It will be interesting to
construct a hieachical kernel function which can not only identify feature
variables within each domain but also identify important domains at the same
time.

Our framework of nonparametric feature selection can be generalized to
precision medicine where one of the main goals is to identify predictive
biomarkers for treatment response. We can adopt loss functions used for
precision medicine in our proposed method to simultaneously accomplish
variable selection and discovering optimal individual treatment rules.
Extensions to categorical outcomes and multi-stage treatment rule estimation
are also possible under our general framework, which can be pursued in future
work.

## References

  * Allen (2013) Allen, G. I. (2013).  Automatic feature selection via weighted kernels and regularization.  Journal of Computational and Graphical Statistics 22, 284–299.
  * Fan et al. (2011) Fan, J., Feng, Y., and Song, R. (2011).  Nonparametric independence screening in sparse ultra-high dimensional additive models.  Journal of the American Statistical Association 106, 544–557.
  * Fan and Li (2001) Fan, J. and Li, R. (2001).  Variable selection via nonconcave penalized likelihood and its oracle properties.  Journal of the American Statistical Association 96, 1348–1360.
  * Fan and Lv (2008) Fan, J. and Lv, J. (2008).  Sure independence screening for ultra-high dimensional feature space.  Journal of the Royal Statistical Society: Series B 70, 849–911.
  * Gao and Wu (2012) Gao, C. and Wu, X. (2012).  Kernel support tensor regression.  2012 International Workshop on Information and Electronics Engineering (IWIEE) 29, 3986–3990.
  * Guyon and Elisseeff (2003) Guyon, I. and Elisseeff, A. (2003).  An introduction to variable and feature selection.  Journal of Machine Learning Research 3, 1157–1182.
  * Huang et al. (2010) Huang, J., Horowitz, J. L., and Wei, F. (2010).  Variable selection in nonparametric additive model.  The Annals of Statistics 38, 2282–2313.
  * Jaakkola et al. (1999) Jaakkola, T., Diekhans, M., and Haussler, D. (1999).  Using the fisher kernel method to detect remote protein.  ISMB 99, 149–158.
  * LI et al. (2012) LI, G., PENG, H., ZHANG, J., and ZHU, L. (2012).  Robust rank correlation based screening.  Annals of Statistics 40, 1846–1877.
  * Lin and Zhang (2006) Lin, Y. and Zhang, H. H. (2006).  Component selection and smoothing in multivariate nonparametric regression.  The Annals of Statistics 34, 2272–2297.
  * Locke et al. (2009) Locke, M., Tinsley, C. L., Benso, M. A., and Blake, D. J. (2009).  Trim32 is an e3 ubiquitin ligase for dysbindin.  Human Molecular Genetics 18, 2344–2358.
  * Ravikumar et al. (2009) Ravikumar, P., Lafferty, J., Liu, H., and Wasserman, L. (2009).  Sparse additive models.  Journal of the Royal Statistical Society: Series B 101,.
  * Scheetz et al. (2006) Scheetz, T. E., Kim, K.-Y. A., Swiderski, R. E., and Philp, A. R. (2006).  Regulation of gene expression in the mammalian eye and its relevance to eye disease.  Proc Natl Acad Sci U S A. 103, 14429–14434.
  * Song et al. (2012) Song, L., Smola, A., Gretton, A., Bedo, J., and Borgwardt, K. (2012).  Feature selection via dependence maximization.  Journal of Machine Learning Research 13, 1393–1434.
  * Steinwart and Scovel (2007) Steinwart, I. and Scovel, C. (2007).  Fast rates for support vector machines using gaussian kernels.  Annals of Statistics 35, 575–607.
  * Tibshirani (1996) Tibshirani, R. (1996).  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B 58, 267–288.
  * Urbanowicz et al. (2018) Urbanowicz, R. J., Meeker, M., Cava, W. L., and Olson, R. S. (2018).  Relief-based feature selection: Introduction and review.  Journal of Biomedical Informatics 85, 189–203.
  * Wang and Kulasekera (2012) Wang, D. and Kulasekera, K. (2012).  Parametric component detection and variable selection in varying-coefficient partially linear models.  Journal of Multivariate Analysis 112, 118–129.
  * Yamada et al. (2014) Yamada, M., Jitkrittum, W., Sigal, L., Xing, E. P., and Sugiyama, M. (2014).  High-dimensional feature selection by feature-wise non-linear lasso.  Neural Computation 26, 185–207.
  * Zhang (2010) Zhang, C. (2010).  Nearly unbiased variable selection under minimax concave penalty.  The Annals of Statistics 38, 894–942.