# Scalable variable selection for two-view learning tasks with projection
operators

Sandor Szedmak
Department of Computer Science
Aalto University
Espoo, Finland
<EMAIL_ADDRESS>
&Riikka Huusari
Department of Computer Science
Aalto University
Espoo, Finland
<EMAIL_ADDRESS>
&Tat Hong Duong Le
Department of Computer Science
Aalto University
Espoo, Finland
<EMAIL_ADDRESS>
&Juho Rousu
Department of Computer Science
Aalto University
Espoo, Finland
<EMAIL_ADDRESS>

###### Abstract

In this paper we propose a novel variable selection method for two-view
settings, or for vector-valued supervised learning problems. Our framework is
able to handle extremely large scale selection tasks, where number of data
samples could be even millions. In a nutshell, our method performs variable
selection by iteratively selecting variables that are highly correlated with
the output variables, but which are not correlated with the previously chosen
variables. To measure the correlation, our method uses the concept of
projection operators and their algebra. With the projection operators the
relationship, correlation, between sets of input and output variables can also
be expressed by kernel functions, thus nonlinear correlation models can be
exploited as well. We experimentally validate our approach, showing on both
synthetic and real data its scalability and the relevance of the selected
features.
Keywords: Supervised variable selection, vector-valued learning, projection-
valued measure, reproducing kernel Hilbert space

##  1 Introduction

Vector-valued, or more generally structured output learning tasks arising from
various domains have attracted much research attention in recent years
[Micchelli and Pontil, 2005, Deshwal et al., 2019, Brogat-Motte et al., 2022].
For both supervised but also unsupervised learning approaches, multi-view data
has been of interest [Hotelling, 1936, Xu et al., 2013, Minh et al., 2016a].
Despite many successful approaches for various multi-view and vector-valued
learning settings, including interpretability to these models has received
less attention. While there are various feature selection and dimensionality
reduction methods either for scalar-valued learning tasks, or unsupervised
methods for data represented in a single view [Zebari et al., 2020, Li et al.,
2017, Anette and Nokto, 2018, Bommert et al., 2020], there is scarcity of
methods suitable for when data is represented in two views, or arises from a
vector-valued learning task. From the point of view of interpretability,
especially feature selection methods are advantageous over dimensionality
reduction since the relevant features are directly obtained as a result and
not given only in (linear) combinations.

Recently, some feature selection methods have been proposed for structured
output learning tasks. [Brouard et al., 2022] proposed kernel-based non-linear
feature selection model relying on sparsity regularization. Another supervised
feature selection approach based on kernel methods is introduced in [Song et
al., 2012], this one relying instead on forward- and backward-selection
ideology. In addition, [Jordan et al., 2021] discusses feature selection in
conjunction with kernel-based models, obtaining sparsity implicitly via loss
function without explicit regularization term. An alternative, spline based,
approach to the non-linear feature selection is proposed by [Boyd et al.,
2018]. These methods, relying on the kernel evaluations between data samples
for both inputs and outputs, tend not to scale very well to large sample
sizes.

In this paper, we introduce a novel variable selection approach for vector-
valued, or two-view learning tasks, including CCA. Our method is based on
efficient iterative computation of projections of input variables to the
vector space intersection between the space spanned by the output variables
and the one of the previously selected input variables. In this space, the
input variables are then selected by a correlation-based criterion. Going one
step further, we also exploit a kernel-based representation of the variables,
allowing us to capture complex, nonlinear relationships. Here, we consider the
kernelised representation of the variables instead of data samples – in
essence, we model the co-variance on the features in the Hilbert space induced
by the kernel. Notably, both input and output features are captured with the
same kernelisation. This is in stark contrast to other proposed kernel-based
feature selection approaches in literature, where separate kernels are used
for data samples in input and output spaces [Brouard et al., 2022, Song et
al., 2012, Jordan et al., 2021]. We can more readily draw comparisons to
canonical correlation analysis (CCA) and its kernelized version, where the
correlations are computed between two sets of variables instead of pairs of
individual ones [Bie et al., 2005].

For many approaches, scalability in feature selection can be a challenge for
when the data dimensionality is extremely large. Some supervised linear
feature selection models adapted to this setting are proposed in [Fan et al.,
2009, Aghazadeh et al., 2018, Valcárcel et al., 2022]. We note, that all these
methods are for the supervised setting, but with scalar-valued output
variables. While scalability w.r.t the feature dimensionality is often
considered due to motivations arising from fat data, the scalability to large
sample sizes is less focused on. Traditionally, kernelized algorithms, while
powerful, are very poorly scalable due to the dependence to the kernel matrix,
especially if its inverse is required. Contrary to the usual, by leveraging
the recursive formulation of our algorithm and a trick with singular value
decomposition on the variable representation, our approach is extremely
scalable to large sample sizes - which we also demonstrate experimentally in
Section 5.

To summarize, our main contributions in this paper are as follows:

  * •

we propose projective selection (ProjSe) algorithm, a novel approach for
variable selection for vector-valued or two-view learning problems that is
based on projection operators. In ProjSe the result of the feature selection
only depends on the subspace spanned by the outputs, not on the specific
values (invariance).

  * •

our proposed iterative method offers high scalability even for the kernelised
formulation capturing non-linearities in the data, due to a trick with
singular value decomposition applied to the feature representation.

  * •

we experimentally validate the proposed approach, showing both relevance of
the selected features and the efficiency of the algorithm.

Table 1: Some of the frequently used notation in this paper. $\mathcal{H}$ |  is a Hilbert space - unless otherwise noted, it has finite dimension $d$, in which case $\mathcal{H}=\mathbb{R}^{d}$.
---|---
$\oplus$ |  denotes the direct sum of subspaces
$\mathbf{I}$ |  is the identity operator acting on $\mathcal{H}$.
$\mathcal{L}$ |  is an arbitrary subspace of $\mathcal{H}$.
$\mathcal{L}_{\mathbf{X}}$ |  is a subspace of $\mathcal{H}$, spanned by the columns of matrix $\mathbf{X}$.
$\mathcal{L}^{\perp}$ |  is a subspace of $\mathcal{H}$, the orthogonal complement of $\mathcal{L}$.
$\mathbf{P}_{\mathcal{L}}$ |  is an orthogonal projection operator into subspace $\mathcal{L}$, $\mathbf{P}_{\mathcal{L}}:\mathcal{H}\rightarrow\mathcal{L}$
$\mathcal{L}_{\mathbf{P}}$ |  is the subspace corresponding to the projection operator $\mathbf{P}$.
$\mathbf{P}_{\mathcal{L}^{\perp}}$ |  is an orthogonal projection operator into the orthogonal complement of subspace $\mathcal{L}$.
$\mathbf{P}_{\mathbf{X}}$ |  is an orthogonal projection operator into the subspace of $\mathcal{H}$, spanned by the columns of matrix $\mathbf{X}$. It is the same as $\mathbf{P}_{\mathcal{L}_{X}}$.
$\mathbf{P}_{\mathbf{X}^{\perp}}$ |  is an orthogonal projection operator into the subspace of $\mathcal{H}$ orthogonal to the subspace spanned by the columns of matrix $\mathbf{X}$. It is the same as $\mathbf{P}_{\mathcal{L}_{X}^{\perp}}$.
$\mathbf{A}^{+}$ |  denotes the Moore-Penrose inverse of matrix $\mathbf{A}$.
$[n]$ |  is a short hand notation for the set $\\{1,\dots,n\\}.$
$\mathbf{A}\circ\mathbf{B}$ |  denotes pointwise(Hadamard) product of matrices $\mathbf{A}$ and $\mathbf{B}$.
$\mathbf{A}^{\circ n}$ |  is the pointwise power of matrix $\mathbf{A}$.
$\mathbf{A}[:,\mathcal{I}]$ |  selects the subset of columns of matrix $\mathbf{A}$ with indices in set $\mathcal{I}$.

The paper is organised as follows. In the next section we give overview of the
method, before moving to more rigorous treatment in Section 3. There we give a
brief introduction to projection operators and their matrix representation,
and discuss the key detail of our approach, expressing the projector into
intersection. We then move on to describing our large-scale kernelized
adaptation of the algorithm in Section 4. We validate our approach
experimentally in Section 5 before concluding.

##  2 Method overview

Figure 1: Illustration of the main steps of the algorithm

Our algorithm is designed to perform variable selection when there are
multiple dependent variables of interest.We denote the matrix containing the
data from which the variables are selected as
$\mathbf{X}\in\mathbb{R}^{m\times n_{x}}$, and the reference data as
$\mathbf{Y}\in\mathbb{R}^{m\times n_{y}}$ – the sample size is $m$, and the
number of features/variables are $n_{x}$ and $n_{y}$ (see other frequently
used notation in Table 1). Here $\mathbf{X}$ and $\mathbf{Y}$ could also
correspond to vector-valued inputs and outputs of some supervised learning
task. Our method is based on defining correlation via projection operators: we
define the correlation between a variable vector $\mathbf{x}\in\mathbb{R}^{m}$
(a column vector from $\mathbf{X}$ containing the values of a single input
variable for all data points) and a set of variables in columns of matrix
$\mathbf{Y}$, as

$\text{corr}(\mathbf{x},\mathbf{Y})=\left\|\mathbf{P}_{\mathcal{L}_{Y}}\frac{\mathbf{x}}{||\mathbf{x}||}\right\|=\left\langle\mathbf{P}_{\mathcal{L}_{Y}}\frac{\mathbf{x}}{||\mathbf{x}||},\mathbf{P}_{\mathcal{L}_{Y}}\frac{\mathbf{x}}{||\mathbf{x}||}\right\rangle^{\frac{1}{2}}=\left\langle\frac{\mathbf{x}}{||\mathbf{x}||},\mathbf{P}_{\mathcal{L}_{Y}}\frac{\mathbf{x}}{||\mathbf{x}||}\right\rangle^{\frac{1}{2}}$
(1)

where $\mathbf{P}_{\mathcal{L}_{Y}}$ (or $\mathbf{P}_{Y}$ in shorthand) is the
orthogonal projection operator into a subspace $\mathcal{L}_{Y}$ spanned by
the columns of $\mathbf{Y}$. This definition is motivated by the concept of
Projection-Valued Measure which plays a significant role in quantum mechanics
theory (see for example [Nielsen and Chuang, 2000]). Our approach selects
variables from input data $\mathbf{X}$ iteratively, such that correlation
between the selected variable and the outputs is high, while correlation to
the previously selected variables is low.

######  Remark 1.

For sake of simplicity, we assume that for all $\mathbf{x}\in\mathbb{R}^{m}$,
$\|\mathbf{x}\|=1$.

Our variable selection algorithm, ProjSe, is illustrated in Figure 1. The
first set of variables is is chosen simply to maximize the projection onto the
subspace spanned by columns of $\mathbf{Y}$, $\mathcal{L}_{Y}$. This is
illustrated with $\mathbf{x}_{1}$, which is projected with $\mathbf{P}_{Y}$ as
$\mathbf{P}_{Y}\mathbf{x}_{1}$. The second set of features chosen,
$\mathbf{x}_{2}$ in the figure, is projected into the intersection of
$\mathcal{L}_{Y}$, and the orthogonal complement of the chosen feature
$\mathbf{x}_{1}$, $\mathcal{L}_{\mathbf{x_{1}}^{\perp}}$. At this step, the
correlation is measured with the projection operator
$\mathbf{P}_{\mathcal{L}_{Y}\cap\mathcal{L}_{\mathbf{x}_{1}^{\perp}}}$.
Interestingly, it turns out that this projected feature,
$\mathbf{P}_{\mathbf{Y}\cap\mathbf{x}_{1}^{\perp}}\mathbf{x}_{2}$, lies also
in the intersection of $\mathcal{L}_{Y}$ and
$\mathcal{L}_{(\mathbf{P}_{Y}\mathbf{x}_{1})^{\perp}}$. This observation paves
the way for building our efficient, recursive algorithm for the feature
selection with projection operators.

The pseudo-code of the basic form of our proposed variable selection by
projection, ProjSe, algorithm is displayed in Figure 2. The approach is fully
deterministic without randomness, and thus practical to apply. Similarly to
CCA, our variable selection algorithm in a sense joins the variable spaces of
the inputs and outputs – both of them are considered in the same space. At the
same time, in order to our selection approach to work,
$\mathcal{L}_{\mathbf{X}}$ should not be fully orthogonal
$\mathcal{L}_{\mathbf{Y}}$. Additionally, due to the properties of the
projection operators, our approach promotes invariance: the selected
explanatory variables (input features) depend only on the subspace spanned by
the response variables (output features), and are independent on any
transformation on the response variables that would span the same subspace.
These transformations can be singular or even nonlinear, as long as they are
automorphisms of the output space.

In this basic form the algorithm is is scalable to medium-scale data, as it is
limited memory required to store the projection matrix.In the following
sections we present techniques that allow scaling to very large datasets, e.g.
$m>1000000$ and $m\gg n_{x},n_{y}$. A recursive representation of the
projection operators (see Section 3.2), and especially the singular vector
based form, (eq. (9)), significantly reduces the demand for resources, both
for memory and for computation time.

Generic algorithm of the variable selection by projection 1. Input: (a) A set
of output variables $\\{\mathbf{y}_{1},\dots,\mathbf{y}_{n_{y}}\\}$ of
$\mathbb{R}^{m}$, in $\mathbf{Y}\in\mathbb{R}^{m\times n_{y}}$. (b) A set of
input variables $\\{\mathbf{x}_{1},\dots,\mathbf{x}_{n_{x}}\\}$ of
$\mathbb{R}^{m}$, collected to $\mathbf{X}\in\mathbb{R}^{m\times n_{x}}$. (c)
$D\leq n_{y}$: number of variables to be chosen from $\mathbf{X}$. 2. Output:
Set $\mathcal{I}_{D}$ of indices of selected variables from $\mathbf{X}$ in
the selection order. The algorithm: 1. Let $t=0$ and
$\mathcal{I}_{t}=\emptyset$.
$\tilde{\mathbf{X}}_{t}=\mathbf{X}[:,\mathcal{I}_{t}]$; since
$\mathcal{I}_{t}=\emptyset$, $\tilde{\mathbf{X}}_{t}^{\perp}=\mathbb{R}^{m}$.
2. For $t\in\\{1,\dots,D\\}$ Do (a) Let
$\mathbf{P}_{\mathcal{L}_{Y}\cap\mathcal{L}_{\tilde{X}_{t}^{\perp}}}$ be the
projection into the intersection of $\mathcal{L}_{Y}$ and
$\mathcal{L}_{\tilde{X}_{t}^{\perp}}$. (b) Choose an index by
$k_{*}=\operatorname*{arg\,max}_{k\in\\{1,\dots,n_{x}\\}\setminus\mathcal{I}_{t}}\left\|\mathbf{P}_{\mathcal{L}_{Y}\cap\mathcal{L}_{\tilde{X}_{t}^{\perp}}}\tfrac{\mathbf{x}_{k}}{\|\mathbf{x}_{k}\|}\right\|^{2}$
(c) $\mathcal{I}_{t+1}=\mathcal{I}_{t}\cup\\{k_{*}\\}$, and
$\tilde{\mathbf{X}}_{t+1}=\mathbf{X}[:,\mathcal{I}_{t+1}]$.

Figure 2: The generic algorithm of supervised variable selection by projection

##  3 Projection operators

This section first introduces relevant background on projection operators and
their algebra. Then, two key points for our algorithm are discussed: the
matrix representation of the projectors, and how the projection into the
intersection can be expressed.

###  3.1 Projection operators, projectors

We now briefly introduce the mathematical framework describing the projection
operators of a Hilbert space. The proofs of the statements mentioned, as well
as further details, are presented for example by [Kreyszig, 1989].

Let $\mathbf{T}$ be a linear operator
$\mathbf{T}:\mathcal{H}\rightarrow\mathcal{H}$. Its adjoint
$\mathbf{T}^{*}:\mathcal{H}\rightarrow\mathcal{H}$ is defined by
$\braket{\mathbf{y},\mathbf{T}^{*}\mathbf{x}}=\braket{\mathbf{T}\mathbf{y},\mathbf{x}}$
for all $\mathbf{x},\mathbf{y}\in\mathcal{H}$. A linear operator $\mathbf{T}$
is self-adjoint, or Hermitian if $\mathbf{T}=\mathbf{T}^{*}$, unitary if
$\mathbf{T}^{*}=\mathbf{T}^{-1}$ and normal if
$\mathbf{T}\mathbf{T}^{*}=\mathbf{T}^{*}\mathbf{T}$. On the set of self-
adjoint operators of $\mathcal{H}$ one can define a partial order $\preceq$ by

$\mathbf{T}_{1}\preceq\mathbf{T}_{2}\Leftrightarrow\braket{\mathbf{T}_{1}\mathbf{x},\mathbf{x}}\leq\braket{\mathbf{T}_{2}\mathbf{x},\mathbf{x}},$
(2)

for all $\mathbf{x}\in\mathcal{H}$. An operator $\mathbf{T}$ is positive if

$\mathbf{0}\preceq\mathbf{T}\Leftrightarrow
0\leq\braket{\mathbf{T},\mathbf{x}},$ (3)

for all $\mathbf{x}\in\mathcal{H}$. As a consequence we have

$\mathbf{T}_{1}\preceq\mathbf{T}_{2}\Leftrightarrow\mathbf{0}\preceq\mathbf{T}_{2}-\mathbf{T}_{1}.$
(4)

Let $\mathcal{L}$ be a subspace of $\mathcal{H}$, the orthogonal complement of
$\mathcal{L}$ is given by
$\mathcal{L}^{\perp}=\\{\mathbf{x}|\mathbf{x}\perp\mathbf{z},\forall\mathbf{z}\in\mathcal{L},\mathbf{x}\in\mathcal{H}\\}$.

######  Theorem 2.

For any subspace $\mathcal{L}\subseteq\mathcal{H}$,
$\mathcal{H}=\mathcal{L}\oplus\mathcal{L}^{\perp}$.

A linear operator $\mathbf{P}$ is a projection operator if
$\mathbf{P}:\mathcal{H}\rightarrow\mathcal{L}$ for a subspace $\mathcal{L}$ of
$\mathcal{H}$. To highlight the connection between the subspace and
projection, they can be also denoted as $\mathcal{L}_{P}$ and
$\mathbf{P}_{L}$.

An operator $\mathbf{P}$ is idempotent if
$\mathbf{P}\mathbf{P}\mathbf{x}=\mathbf{P}\mathbf{x},\ \text{or}\
\mathbf{P}\mathbf{P}=\mathbf{P}$ holds for any $\mathbf{x}\in\mathcal{H}$. The
projection operators can be characterized by the following statements.

######  Theorem 3.

A linear operator $\mathbf{P}:\mathcal{H}\rightarrow\mathcal{H}$ is a
projection if it is self adjoint, $\mathbf{P}=\mathbf{P}^{*}$, and idempotent
$\mathbf{P}\mathbf{P}=\mathbf{P}$.

######  Proposition 4.

The map connecting the set of closed subspaces111In a finite dimensional
Hilbert space all subspaces are closed. of $\mathcal{H}$ and the set of the
corresponding orthogonal projections is bijective.

As a consequence of the idempotent and self-adjoint properties we have that
the range $\mathcal{R}(\mathbf{P})$ and the null space
$\mathcal{N}(\mathbf{P})$ of $\mathbf{P}$ are orthogonal, namely for any
$x,y\in\mathcal{H}$

$\braket{\mathbf{P}x,y-\mathbf{P}y}=\braket{\mathbf{P}^{2}x,y-\mathbf{P}y}=\braket{\mathbf{P}x,(\mathbf{P}-\mathbf{P}^{2})y}=0.$
(5)

The following theorems describe some algebraic properties of projection
operators we are going to exploit.

######  Theorem 5.

(Product of projections) Let $\mathbf{P}_{1}$ and $\mathbf{P}_{2}$ be
projections on $\mathcal{H}$. $\mathbf{P}=\mathbf{P}_{1}\mathbf{P}_{2}$ is
projection if and only if
$\mathbf{P}_{1}\mathbf{P}_{2}=\mathbf{P}_{2}\mathbf{P}_{1}$. Then
$\mathbf{P}:\mathcal{H}\rightarrow\mathcal{L}_{P_{1}}\cap\mathcal{L}_{P_{2}}$.

######  Theorem 6.

(Sum of projections) Let $\mathbf{P}_{1}$ and $\mathbf{P}_{2}$ be projections
on $\mathcal{H}$. $\mathbf{P}=\mathbf{P}_{1}+\mathbf{P}_{2}$ is projection if
and only if $\mathcal{L}_{P_{1}}\perp\mathcal{L}_{P_{2}}$. Then
$\mathbf{P}:\mathcal{H}\rightarrow\mathcal{L}_{P_{1}}\oplus\mathcal{L}_{P_{2}}$.

######  Theorem 7.

(Partial order) Let $\mathbf{P}_{1}$ and $\mathbf{P}_{2}$ be projections on
$\mathcal{H}$, and $\mathcal{N}(\mathbf{P}_{1})$ and
$\mathcal{N}(\mathbf{P}_{2})$ the corresponding null spaces. Then the
following statements are equivalent.

$\begin{array}[]{ll}\mathbf{P}_{1}\mathbf{P}_{2}&=\mathbf{P}_{2}\mathbf{P}_{1}=\mathbf{P}_{1},\\\
\mathcal{L}_{P_{1}}&\subseteq\mathcal{L}_{P_{2}},\\\
\mathcal{N}(\mathbf{P}_{1})&\supseteq\mathcal{N}(\mathbf{P}_{2}),\\\
||\mathbf{P}_{1}\mathbf{x}||&\leq||\mathbf{P}_{2}\mathbf{x}||,\\\
\mathbf{P}_{1}&\preceq\mathbf{P}_{2}.\\\ \end{array}$ (6)

######  Theorem 8.

(Difference of projections) Let $\mathbf{P}_{1}$ and $\mathbf{P}_{2}$ be
projections on $\mathcal{H}$. $\mathbf{P}=\mathbf{P}_{2}-\mathbf{P}_{1}$ is
projection if and only $\mathcal{L}_{P_{1}}\subseteq\mathcal{L}_{P_{2}}$. Then
$\mathbf{P}:\mathcal{H}\rightarrow\mathcal{L}_{P}$, where
$\mathcal{L}_{P_{2}}=\mathcal{L}_{P_{1}}\oplus\mathcal{L}_{P}$, namely
$\mathcal{L}_{P}$ is the orthogonal complement of $\mathcal{L}_{P_{1}}$ in
$\mathcal{L}_{P_{2}}$.

From the theorems above we can derive a simple corollary: if $\mathcal{L}$ is
a subspace, then the projection into its complement is equal to
$\mathbf{P}_{\mathcal{L}^{\perp}}=\mathbf{I}-\mathbf{P}_{\mathcal{L}}$.

######  Theorem 9.

(Monotone increasing sequence) Let $(\mathbf{P}_{n})$ be monotone increasing
sequence of projections defined on $\mathcal{H}$. Then:

  1. 1.

$(\mathbf{P}_{n})$ is strongly operator convergent, and the limit
$\mathbf{P}$, $\mathbf{P}_{n}\rightarrow\mathbf{P}$, is a projection.

  2. 2.

$\mathbf{P}:\mathcal{H}\rightarrow\cup_{n=1}^{\infty}\mathcal{L}_{\mathbf{P}_{n}}$.

  3. 3.

$\mathcal{N}(\mathbf{P})=\cap_{n=1}^{\infty}\mathcal{N}(\mathbf{P}_{n})$.

If $\mathbf{S}$ is a self-adjoint operator and $\mathbf{P}$ is a projection
into the range of $\mathbf{S}$ then
$\mathbf{S}\mathbf{P}=\mathbf{P}\mathbf{S}$, see [Conway, 1997] for further
details.

Let $\mathbf{I}$ be projection into the entire space, and $\mathbf{0}$ its
complement. If $\mathbf{0}\leq\mathbf{S}\leq\mathbf{I}$, and
$\mathbf{T}\geq\mathbf{0}$ operators. If $\mathbf{P}$ is a projection into
$\mathbf{S}+\mathbf{T}$ then $\mathbf{P}$ commutes both $\mathbf{S}$ and
$\mathbf{T}$. See in [Hayashi and Nagaoka, 2002].

###  3.2 Matrix representation of projectors

If a basis of the Hilbert space $\mathcal{H}$ is fixed, then every linear
operator acting on $\mathcal{H}$ can be represented by a matrix. Let the
subspace $\mathcal{L}$ of $\mathcal{H}$ be spanned by the vectors
$\mathbf{a}_{1},\dots,\mathbf{a}_{k}$ of $\mathcal{H}$. Let us construct a
matrix $\mathbf{A}$ whose columns are equal to the vectors
$\mathbf{a}_{1},\dots,\mathbf{a}_{k}$. Here the linear independence of those
vectors is not assumed. The corresponding subspace is denoted by
$\mathcal{L}_{A}$

The matrix representing the orthogonal projection operator into to subspace
$\mathcal{L}_{A}$ can be expressed by a well-known minimization problem [Golub
and Loan, 2013],

$\arg\min_{\mathbf{w}}\|\mathbf{x}-\mathbf{A}\mathbf{w}\|^{2}=\mathbf{w}^{*}=(\mathbf{A}^{T}\mathbf{A})^{+}\mathbf{A}^{T}\mathbf{x},$
(7)

where + denotes the Moore-Penrose pseudo-inverse. Based on eq. (5) the vector
$\mathbf{A}^{T}\mathbf{w}^{*}$ is the orthogonal projection of $\mathbf{x}$
into $\mathcal{L}$. The orthogonal projection of $\mathbf{x}$ is equal to
$\mathbf{P}_{A}\mathbf{x}=\mathbf{A}(\mathbf{A}^{T}\mathbf{A})^{+}\mathbf{A}^{T}\mathbf{x}$.
Since this is true for any $\mathbf{x}\in\mathcal{H}$, the matrix
representation of the orthogonal projection operator $\mathbf{P}_{A}$ is given
by

$\mathbf{P}_{A}=\mathbf{A}(\mathbf{A}^{T}\mathbf{A})^{+}\mathbf{A}^{T}.$ (8)

This formula can be simplified by exploiting the properties of the Moore-
Penrose pseudo-inverse, see for example [Ben-Israel and Greville, 2003], via
the singular value decomposition
$\mathbf{U}_{A}\mathbf{S}_{A}\mathbf{V}_{A}^{T}$ of the matrix $\mathbf{A}$.
Here we assume that the matrix $\mathbf{A}\in\mathbb{R}^{m\times n_{A}}$,
$m>n_{A}$, and $\mathbf{V}_{A}$ is a square matrix, but $\mathbf{U}_{A}$
contains only those left singular vectors where the corresponding singular
values are not equal to zero. We have

$\framebox{$\mathbf{P}_{A}$}=\mathbf{A}(\mathbf{A}^{T}\mathbf{A})^{+}\mathbf{A}^{T}=\mathbf{A}\mathbf{A}^{+}=\mathbf{U}_{A}\mathbf{S}_{A}\mathbf{V}_{A}^{T}\mathbf{V}_{A}\mathbf{S}_{A}^{+}\mathbf{U}_{A}^{T}=\framebox{$\mathbf{U}_{A}\mathbf{U}_{A}^{T}$}.$
(9)

This representation of the projection operator plays a central role in our
variable selection algorithm. The following proposition ensures that the
projection operator does not depend on its representation.

######  Proposition 10.

Assume that two different matrices $\mathbf{A}$ and $\mathbf{B}$ span the same
subspace $\mathcal{L}$ of dimension $k$. Then the two representations
$\mathbf{P}_{A}=\mathbf{U}_{A}\mathbf{U}_{A}^{T}$ and
$\mathbf{P}_{B}=\mathbf{U}_{B}\mathbf{U}_{B}^{T}$ yield the same projection
operator.

###### Proof.

Since the columns of $\mathbf{U}_{B}$ as linear combinations of $\mathbf{B}$
are in the $\mathcal{L}$, thus $\mathbf{P}_{A}\mathbf{U}_{B}=\mathbf{U}_{B}$.
Multiplying both sides with $\mathbf{U}_{B}^{T}$ we obtain that
$\mathbf{P}_{A}\mathbf{U}_{B}\mathbf{U}_{B}^{T}=\mathbf{U}_{B}\mathbf{U}_{B}^{T}$
which is $\mathbf{P}_{A}\mathbf{P}_{B}=\mathbf{P}_{B}$. Because the right hand
side, $\mathbf{P}_{B}$, is a projection, the left hand side
$\mathbf{P}_{A}\mathbf{P}_{B}$ is also one. Following the same line we have
$\mathbf{P}_{B}\mathbf{P}_{A}=\mathbf{P}_{A}$ as well. From Theorem 5 we know
that if the product of projections is a projection, then the product of
projections is commutative,
$\mathbf{P}_{B}\mathbf{P}_{A}=\mathbf{P}_{A}\mathbf{P}_{B}$. Finally we can
conclude that

$\mathbf{P}_{A}=\mathbf{P}_{B}\mathbf{P}_{A}=\mathbf{P}_{A}\mathbf{P}_{B}=\mathbf{P}_{B}.$

∎

We also exploited that if $\mathcal{H}$ is finite dimensional and the
corresponding field is $\mathbb{R}$ then the adjoint of $\mathbf{P}^{*}$ is
represented by the transpose $\mathbf{P}^{T}$ of the matrix $\mathbf{P}$.

####  3.2.1 Projection onto the intersection of subspaces - General view

Our algorithm hinges on the orthogonal projector of the intersection of a set
of subspaces
$\\{\mathcal{L}_{1},\mathcal{L}_{2},\dots,\mathcal{L}_{n_{L}}\\}$. To
introduce this concept, here we mainly follow the line presented by [Ben-
Israel, 2015]. We can start with some classical result, first we can recall
[von Neumann, 1950], who derived a solution in case of two subspaces as a
limit:

$\mathbf{P}_{\mathcal{L}_{1},\cap\mathcal{L}_{2}}=\lim_{n\rightarrow\infty}(\mathbf{P}_{\mathcal{L}_{1}}\mathbf{P}_{\mathcal{L}_{2}})^{n}.$
(10)

That result has been extended to arbitrary finite sets of subspaces by
[Halperin, 1962],:

$\mathbf{P}_{\mathcal{L}_{1},\cap\dots\cap\mathcal{L}_{n_{L}}}=\lim_{n\rightarrow\infty}(\mathbf{P}_{\mathcal{L}_{1}}\dots\mathbf{P}_{L_{n}})^{n}.$
(11)

[Anderson. and Duffin, 1969], gave an explicit formula for the case of two
subspaces by

$\mathbf{P}_{\mathcal{L}_{1},\cap\mathcal{L}_{2}}=2\mathbf{P}_{\mathcal{L}_{1}}(\mathbf{P}_{\mathcal{L}_{1}}+\mathbf{P}_{\mathcal{L}_{2}})^{\dagger}\mathbf{P}_{\mathcal{L}_{2}}.$
(12)

[Ben-Israel, 2015], provides an alternative to compute
$\mathbf{P}_{\mathcal{L}_{1},\cap\dots\cap\mathcal{L}_{n_{L}}}$ Here we rely
on the Lemma 4.1 and Corollary 4.2 of his work:

######  Proposition 11.

For $i=1,\dots,n_{L}$, let $\mathcal{L}_{i}$ be subspaces of $\mathcal{H}$,
$\mathbf{P}_{i}$ be the corresponding projectors,
$\mathbf{P}^{\perp}_{i}=\mathbf{I}-\mathbf{P}_{i}$, and $\lambda_{i}>0$.
Define $\mathbf{Q}:=\sum_{i=1}^{n_{L}}\lambda_{i}\mathbf{P}_{i}^{\perp}.$ then
we have
$\mathbf{P}_{\mathcal{L}_{1},\cap\dots\cap\mathcal{L}_{n_{L}}}=\mathbf{I}-\mathbf{Q}^{\dagger}\mathbf{Q}.$
With the particular choice $\sum_{i=1}^{n_{L}}\lambda_{i}=1$, $\mathbf{Q}$
might be written as
$\mathbf{Q}:=\mathbf{I}-\sum_{i=1}^{n_{L}}\lambda_{i}\mathbf{P}_{i},$
eliminating all the complements of the projectors.

By exploting that for any projector $\mathbf{P}$
$\mathbf{P}^{\perp}=\mathbf{I}-\mathbf{P}$, the $\mathbf{Q}_{t}$ corrsponding
to $\mathbf{P}_{\mathcal{L}_{V}\cap\mathcal{L}_{\tilde{X}_{t}^{\perp}}}$ can
be written as

$\mathbf{Q}_{t}=\lambda_{V}(\mathbf{I}-\mathbf{P}_{\mathcal{L}_{V}})+\sum_{\mathbf{x}\in\tilde{X}_{t}}\lambda_{x}\mathbf{P}_{L_{\\{\mathbf{x}\\}}}.$
(13)

The critical point is the computation of the Moore-Penrose inverse of
$\mathbf{Q}$.

###  3.3 Expressing the projector into intersection

To implement the proposed variable selection algorithm (Figure 2) the
projection into the intersection of an arbitrary subspace
$\mathcal{L}_{\mathbf{P}}$ and the complement of an arbitrary vector
$\mathbf{x}$, $\mathbf{P}_{\mathcal{L}\cap\mathbf{x}^{\perp}}$, has to be
computed. The projector $\mathbf{P}_{\mathcal{L}^{\perp}}$ to the complement
of a subspace $\mathcal{L}$ can be expressed as
$\mathbf{I}-\mathbf{P}_{\mathcal{L}}$, hence the projector into
$\mathbf{P}_{\mathbf{x}^{\perp}}$ is given by
$\mathbf{I}-\dfrac{\mathbf{x}\mathbf{x}^{T}}{||\mathbf{x}||^{2}}$. Since
$\mathcal{L}$ is arbitrary we use $\mathbf{P}$ instead of
$\mathbf{P}_{\mathcal{L}}$ for sake of simplicity.

While we have these two projectors, their product, according to Theorem 5, is
not a projection as it does not commute:

$\mathbf{P}\left(\mathbf{I}-\dfrac{\mathbf{x}\mathbf{x}^{T}}{||\mathbf{x}||^{2}}\right)=\mathbf{P}-\dfrac{\mathbf{P}\mathbf{x}\mathbf{x}^{T}}{||\mathbf{x}||^{2}}\neq\left(\mathbf{I}-\dfrac{\mathbf{x}\mathbf{x}^{T}}{||\mathbf{x}||^{2}}\right)\mathbf{P}=\mathbf{P}-\dfrac{\mathbf{x}\mathbf{x}^{T}\mathbf{P}}{||\mathbf{x}||^{2}},$
(14)

because in the general case
$\mathbf{P}\mathbf{x}\mathbf{x}^{T}\neq\mathbf{x}\mathbf{x}^{T}\mathbf{P}$. To
overcome this problem we can recognize that the intersection
$\mathcal{L}_{\mathbf{P}}\cap\mathcal{L}_{\mathbf{x}^{\perp}}$ can be
expressed after a simple transformation.

######  Lemma 12.

Let $\mathbf{P}$ be a projector and $\mathbf{x}$ be any vector, then the
intersections $\mathcal{L}_{\mathbf{P}}\cap\mathcal{L}_{\mathbf{x}^{\perp}}$
and $\mathcal{L}_{\mathbf{P}}\cap\mathcal{L}_{(\mathbf{P}\mathbf{x})^{\perp}}$
are the same subspaces of $\mathcal{L}_{\mathbf{P}}$.

###### Proof.

Any vector $\mathbf{u}$ is in $\mathcal{L}_{\mathbf{P}}$ if
$\mathbf{P}\mathbf{u}=\mathbf{u}$, $\mathbf{u}$ is in
$\mathcal{L}_{\mathbf{x}^{\perp}}$ if $\braket{\mathbf{x},\mathbf{u}}=0$, and
$\mathbf{u}$ is in $\mathcal{L}_{(\mathbf{P}\mathbf{x})^{\perp}}$ if
$\braket{\mathbf{P}\mathbf{x},\mathbf{u}}=0$. Since
$\mathbf{P}\mathbf{u}=\mathbf{u}$, therefore
$\braket{\mathbf{x},\mathbf{u}}=\braket{\mathbf{P}\mathbf{x},\mathbf{u}}=0$. ∎

By projecting $\mathbf{x}$ into $\mathcal{L}$ first, and then computing the
corresponding intersection, we can compute the projector into
$\mathcal{L}_{\mathbf{P}}\cap\mathcal{L}_{\mathbf{x}^{\perp}}$ in a simple
way.

######  Proposition 13.

Let $\|\mathbf{x}\|=1$ and
$\alpha=\|\mathbf{P}\mathbf{x}\|^{2}=\mathbf{x}^{T}\mathbf{P}\mathbf{P}\mathbf{x}=\mathbf{x}^{T}\mathbf{P}\mathbf{x}$.
If $\alpha>0$ then

$\framebox{$\mathbf{P}_{\mathcal{L}_{\mathbf{P}}\cap\mathcal{L}_{\mathbf{x}^{\perp}}}$}=\mathbf{P}_{\mathcal{L}_{\mathbf{P}}\cap\mathcal{L}_{(\mathbf{P}{\mathbf{x})^{\perp}}}}=\mathbf{P}\left(\mathbf{I}-\tfrac{1}{\alpha}\mathbf{P}\mathbf{x}\mathbf{x}^{T}\mathbf{P}\right)=\framebox{$\mathbf{P}-\tfrac{1}{\alpha}\mathbf{P}\mathbf{x}\mathbf{x}^{T}\mathbf{P}$}.$
(15)

When $\alpha=0$, which means $\mathbf{x}$ orthogonal to
$\mathcal{L}_{\mathbf{P}}$, then we have
$\mathbf{P}_{\mathcal{L}_{\mathbf{P}}\cap\mathcal{L}_{\mathbf{x}^{\perp}}}=\mathbf{P}.$

We can check that
$\mathbf{P}-\tfrac{1}{\alpha}\mathbf{P}\mathbf{x}\mathbf{x}^{T}\mathbf{P}$ is
a real projector. It is idempotent, since

$\begin{array}[]{@{}l@{}}\left(\mathbf{P}-\tfrac{1}{\alpha}\mathbf{P}\mathbf{x}\mathbf{x}^{T}\mathbf{P}\right)^{2}=\mathbf{P}-\tfrac{1}{\alpha}\mathbf{P}\mathbf{x}\mathbf{x}^{T}\mathbf{P}-\tfrac{1}{\alpha}\mathbf{P}\mathbf{x}\mathbf{x}^{T}\mathbf{P}+\tfrac{\alpha}{\alpha^{2}}\mathbf{P}\mathbf{x}\mathbf{x}^{T}\mathbf{P}=\
\mathbf{P}-\tfrac{1}{\alpha}\mathbf{P}\mathbf{x}\mathbf{x}^{T}\mathbf{P}.\end{array}$

This agrees with Theorem 5 which states that the product of projections is a
projection, idempotent and adjoint, and it is map into the intersection of the
corresponding subspaces. Furthermore the orthogonality of $\mathbf{x}$ and the
projection of any $\mathbf{u}\in\mathcal{H}$ by
$\mathbf{P}-\tfrac{1}{\alpha}\mathbf{P}\mathbf{x}\mathbf{x}^{T}\mathbf{P}$ can
also be verified, as

$\begin{array}[]{ll}\left\langle\mathbf{x},\left(\mathbf{P}-\tfrac{1}{\alpha}\mathbf{P}\mathbf{x}\mathbf{x}^{T}\mathbf{P}\right)\mathbf{u}\right\rangle=\mathbf{x}^{T}\mathbf{P}\mathbf{u}-\tfrac{1}{\alpha}\mathbf{x}^{T}\mathbf{P}\mathbf{x}\mathbf{x}^{T}\mathbf{P}\mathbf{u}=\mathbf{x}^{T}\mathbf{P}\mathbf{u}-\tfrac{1}{\alpha}\alpha\mathbf{x}^{T}\mathbf{P}\mathbf{u}=0.\end{array}$

######  Definition 14.

([Tjur, 1984]) Two subspaces $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ are
orthogonal if for any vectors, $\mathbf{x_{1}}\in\mathcal{L}_{1}$ and
$\mathbf{x_{2}}\in\mathcal{L}_{2}$ $\braket{\mathbf{x}_{1},\mathbf{x}_{2}}=0$
holds. Two subspaces $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ are geometrically
orthogonal if
$\mathcal{L}_{1}=(\mathcal{L}_{1}\cap\mathcal{L}_{2})\oplus\mathcal{C}_{1}$
and
$\mathcal{L}_{2}=(\mathcal{L}_{1}\cap\mathcal{L}_{2})\oplus\mathcal{C}_{2}$
and $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ are orthogonal.

######  Lemma 15.

([Tjur, 1984]) Two subspaces $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ are
geometrically orthogonal if and only if
$\mathbf{P}_{\mathcal{L}_{1}}\mathbf{P}_{\mathcal{L}_{2}}=\mathbf{P}_{\mathcal{L}_{2}}\mathbf{P}_{\mathcal{L}_{1}}$
and
$\mathcal{L}_{\mathbf{P}_{\mathcal{L}_{1}}\mathbf{P}_{\mathcal{L}_{2}}}=\mathcal{L}_{1}\cap\mathcal{L}_{2}$.

######  Proposition 16.

The subspaces $\mathcal{L}_{\mathbf{P}}$ and
$\mathcal{L}_{(\mathbf{P}\mathbf{x})^{\perp}}$ are geometrically orthogonal.

###### Proof.

It is a simple Corollary of Proposition 13. ∎

##  4 Selecting variables in RKHS

In order to take into account non-linear correlations in the data, we propose
a kernelized adaptation of the problem. Kernel methods are a group of varied
machine learning models, taking advantage of a symmetric and positive semi-
definite kernel function comparing data samples (sets of features)
$k:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$. The usage of a kernel
function allows including non-linearity to the models implicitly via a feature
map $\varphi:\mathcal{X}\rightarrow\mathcal{F}_{k}$: a kernel evaluated with
two samples corresponds to an inner product in this so-called feature space
(more specifically reproducing kernel Hilbert space, RKHS):
$k(x,z)=\langle\varphi(x),\varphi(z)\rangle_{\mathcal{F}_{k}}$. For more
thorough introduction to traditional kernel methods, we refer the reader e.g.
to [Hofmann et al., 2008].

We here propose to kernelize the variable representation. We consider
$\phi:\mathbb{R}^{m}\rightarrow\mathcal{H}$, where $\mathbb{R}^{m}$ is the
vector space containing all columns of $\mathbf{Y}\in\mathbb{R}^{m\times
n_{y}}$ and $\mathbf{X}\in\mathbb{R}^{m\times n_{x}}$, and $\mathcal{H}$ is a
RKHS. In essence, this corresponds to defining a kernel on the variable
vectors, $\kappa:\mathbb{R}^{m}\times\mathbb{R}^{m}\rightarrow\mathbb{R}$ – in
fact, we assume that the $\phi$ is only given implicitly via $\kappa$. In
mathematical sense, this matrix can equally well be considered to be a kernel
matrix, since distinction between the rows and columns is by convention only.
Usually however the matrix built from inner products between the variables is
referred to as covariance operator. The covariance operators are also extended
to RKHS with various applications in machine learning tasks [Muandet et al.,
2017, Minh et al., 2016b]. Contrary to our approach, there the feature map and
kernel are defined on data space $\mathcal{X}$ instead of variable space
$\mathbb{R}^{m}$. We need to mention here also the Gaussian Process
Regression, [Rasmussen and Williams, 2005] where kernels are also used to
cover the covariance matrix, thus connecting the variables via inner product.

We highlight that as the kernel is defined on variables, we can easily
evaluate $\kappa(\mathbf{x}_{i},\mathbf{y}_{j})$. We use the following
shorthands for feature and kernel evaluations on the available training data:
$\phi(\mathbf{Y})=[\phi(\mathbf{y}_{1}),\dots,\phi(\mathbf{y}_{n_{y}})]$ with
$\phi(\mathbf{y}_{i}),i\in[n_{y}]$ a column vector, and
$\kappa(\mathbf{Y},\mathbf{x})=[\kappa(\mathbf{y}_{1},\mathbf{x}),\dots,\kappa(\mathbf{y}_{n_{y}},\mathbf{x})]^{\top}$
a column vector of kernel evaluations (similarly for $\phi(\mathbf{X})$). Note
that $\kappa(\mathbf{Y},\mathbf{Y})=\phi(\mathbf{Y})^{\top}\phi(\mathbf{Y})$
with this notation. We further denote
$\mathbf{K_{y}}=\kappa(\mathbf{Y},\mathbf{Y})$. We assume that
$\|\phi(\mathbf{x})\|=1$.

###  4.1 Expressing the Projection operator in RKHS

Based on Section 3.2, Equation (9), the projection $\mathbf{P}_{Y}$ is
represented with the left singular vectors of $\mathbf{P}_{Y}$,
$\mathbf{U}_{Y}$. This representation is also needed for the kernelized
algorithm. However calculating directly the singular value decomposition on
$\phi(\mathbf{Y})$,
$\phi(\mathbf{Y})=\mathbf{U}_{Y}\mathbf{S}_{Y}\mathbf{V}_{Y}^{T}$, might not
be feasible if the dimensionality of the feature space is large. Assuming that
$\mathcal{H}$ is finite dimensional222For clarity, we restrict the discussion
to finite dimensions and $\mathcal{H}=\mathbb{R}^{d}$ with $d<\infty$. We note
that the approach is equally valid also with infinite dimensions. with
dimension $d$, we have $\phi(\mathbf{Y}),\mathbf{U}_{Y}\in\mathbb{R}^{d\times
n_{y}}$, and $\mathbf{S}_{Y},\mathbf{V}_{Y}\in\mathbb{R}^{n_{y}\times n_{y}}$.
Therefore we can write

$\mathbf{S}_{Y}=\left[\begin{array}[]{c}\mathbf{D}_{Y}\\\
\boldsymbol{\varnothing}\end{array}\right],\mathbf{D}_{Y}\in\mathbb{R}^{n_{y}\times
n_{y}},\boldsymbol{\varnothing}\in[0]^{m-n_{y},n_{y}},$ (16)

and $\mathbf{D}_{Y}$ diagonal with nonnegative elements of singular values,
thus $\phi(\mathbf{Y})=\mathbf{U}_{Y}\mathbf{D}_{Y}\mathbf{V}_{Y}^{T}$. Again,
this decomposition can not be computed directly, however we can go on the
following line of computation.

To express the $\mathbf{U}_{Y}$ we can apply a similar approach to what is
exploited in the computation of the kernel principal component analysis
[Schölkopf et al., 1998]. Recall that the kernel matrix on columns of
$\phi(\mathbf{Y})$ is $\mathbf{K}_{Y}=\phi(\mathbf{Y})^{T}\phi(\mathbf{Y})$.
From the singular value decomposition we can derive that
$\mathbf{K}_{Y}=\mathbf{V}_{Y}\mathbf{S}_{Y}^{2}\mathbf{V}^{T}$. This kernel
has a reasonably small size, $n_{y}\times n_{y}$, thus its eigenvalue
decomposition can be computed, which yields $\mathbf{V}_{Y}$ and the squares
of the singular values of the diagonal elements of $\mathbf{S}_{Y}$. By
combining these expressions we have

$\phi(\mathbf{Y})\mathbf{V}_{Y}\mathbf{S}_{Y}=\mathbf{U}_{Y}\mathbf{S}_{Y}^{2}\
\Rightarrow\
\framebox{$\mathbf{U}_{Y}=\phi(\mathbf{Y})\mathbf{V}_{Y}\mathbf{S}_{Y}^{+}$}$
(17)

with the help of the Moore-Penrose generalized inverse. Our algorithm hinges
on evaluating products between projectors and the variable vectors. We can now
write the products of the $\mathbf{U}_{Y}^{T}$ with an arbitrary vector
represented in $\mathcal{H}$ as

$\begin{array}[]{ll}\framebox{$\mathbf{U}_{Y}^{T}\phi(\mathbf{x})$}=\mathbf{S}_{Y}^{-1}\mathbf{V}_{Y}^{T}\phi(\mathbf{Y})^{T}\phi(\mathbf{x})=\framebox{
$\mathbf{S}_{Y}^{-1}\mathbf{V}_{Y}^{T}\kappa(\mathbf{Y},\mathbf{x})$}.\end{array}$
(18)

Thus the product can be expressed with the help of the kernel on the variables
with complexity $O(n_{y}^{2})$ if the $\mathbf{K}_{Y}$, $\mathbf{V}_{Y}$ and
$\mathbf{S}_{Y}^{-1}$ are precomputed.

###  4.2 The recursive selection procedure

To calculate the projection operator efficiently in each iteration we can
exploit the structure of
$\mathbf{P}_{\mathcal{L}_{Y}\cap\mathcal{L}_{\tilde{X}_{t}^{\perp}}}\phi(\mathbf{x})$
introduced in Proposition 13. To this end, we define an intermediate operator,
projection into the complement subspace of vector
$\mathbf{q}\in\mathbb{R}^{n}$ as:

$\mathbf{Q}(\mathbf{q})=\mathbf{I}-\dfrac{\mathbf{q}\mathbf{q}^{T}}{||\mathbf{q}||^{2}}.$
(19)

Since $\mathbf{Q}(\mathbf{q})$ is a projection, we have
$\mathbf{Q}(\mathbf{q})=\mathbf{Q}(\mathbf{q})\mathbf{Q}(\mathbf{q})$ and
$\mathbf{Q}(\mathbf{q})=\mathbf{Q}(\mathbf{q})^{T}$. It can also be seen that
multiplying $\mathbf{Q}$ with a matrix $\mathbf{A}\in\mathbb{R}^{n\times n}$,

$\mathbf{Q}(\mathbf{q})\mathbf{A}=\left(\mathbf{I}-\dfrac{\mathbf{q}\mathbf{q}^{T}}{||\mathbf{q}||^{2}}\right)\mathbf{A}=\mathbf{A}-\dfrac{\mathbf{q}(\mathbf{q}^{T}\mathbf{A})}{||\mathbf{q}||^{2}},$
(20)

has the complexity of only $O(n^{2})$ since only matrix-vector and outer
product are needed. We are also going to use the following recurrent matrix
products for a fixed $t$

$\tilde{\mathbf{U}}_{t}=\mathbf{U}_{Y}\prod_{s=1}^{t-1}\mathbf{Q}(\mathbf{q}_{s})=\tilde{\mathbf{U}}_{t-1}\mathbf{Q}(\mathbf{q}_{t-1}).$
(21)

Now we can write up the sequence of projections corresponding to the Algorithm
(2):

$\begin{array}[]{@{}l@{\:}l@{}}\text{Let}\hfil\>&\mathbf{U}_{0}=\mathbf{U}_{Y},\
\mathcal{I}_{0}=\emptyset,\\\
\framebox{$\mathbf{P}_{0}$}\hfil\>&=\mathbf{P}_{\phi(Y)}=\framebox{$\mathbf{U}_{0}\mathbf{U}_{0}^{T}$},\\\
\mathbf{q}_{1}\hfil\>&=\mathbf{U}_{0}^{T}\phi(\mathbf{x}_{k_{1*}}),\
k_{1*}=\arg\max_{k\in[n_{x}]\setminus\mathcal{I}_{0}}||\mathbf{P}_{0}\phi(\mathbf{x}_{k})||^{2},\
\mathcal{I}_{1}=\mathcal{I}_{0}\cup{k_{1^{*}}},\\\
\framebox{$\mathbf{P}_{1}$}\hfil\>&=\mathbf{U}_{0}\mathbf{U}_{0}^{T}-\dfrac{\mathbf{U}_{0}\mathbf{q}_{1}\mathbf{q}_{1}^{T}\mathbf{U}_{0}^{T}}{||\mathbf{q}_{1}||^{2}}=\mathbf{U}_{0}\mathbf{Q}(\mathbf{q}_{1})\mathbf{U}_{0}^{T}=\mathbf{U}_{0}\mathbf{Q}(\mathbf{q}_{1})\mathbf{Q}(\mathbf{q}_{1})\mathbf{U}_{0}^{T}=\framebox{$\mathbf{U}_{1}\mathbf{U}_{1}^{T}$},\\\
\mathbf{q}_{2}\hfil\>&=\mathbf{U}_{1}^{T}\phi(\mathbf{x}_{k_{2*}}),\
k_{2*}=\arg\max_{k\in[n_{x}]\setminus\mathcal{I}_{1}}||\mathbf{P}_{1}\phi(\mathbf{x}_{k})||^{2},\
\mathcal{I}_{2}=\mathcal{I}_{1}\cup{k_{2^{*}}},\\\ \vdots\hfil\>\\\
\framebox{$\mathbf{P}_{t}$}\hfil\>&=\mathbf{U}_{t-1}\mathbf{U}_{t-1}^{T}-\dfrac{\mathbf{U}_{t-1}\mathbf{q}_{t}\mathbf{q}_{t}^{T}\mathbf{U}_{t-1}^{T}}{||\mathbf{q}_{t}||^{2}}=\mathbf{U}_{t-1}\mathbf{Q}(\mathbf{q}_{t})\mathbf{U}_{t-1}^{T}\\\
\hfil\>&=\mathbf{U}_{t-1}\mathbf{Q}(\mathbf{q}_{t})\mathbf{Q}(\mathbf{q}_{t})\mathbf{U}_{t-1}^{T}=\framebox{$\mathbf{U}_{t}\mathbf{U}_{t}^{T}$},\\\
\mathbf{q}_{t+1}\hfil\>&=\mathbf{U}_{t}^{T}\phi(\mathbf{x}_{k_{(t+1)*}}),\
k_{(t+1)*}=\arg\max_{k\in[n_{x}]\setminus\mathcal{I}_{t}}||\mathbf{P}_{t}\phi(\mathbf{x}_{k})||^{2},\
\mathcal{I}_{t+1}=\mathcal{I}_{t}\cup{k_{(t+1)^{*}}},\\\
\vdots\hfil\>\end{array}$

######  Proposition 17.

The sequence of projections above correctly computes the projection operators
of Algorithm in Figure 2.

###### Proof.

We apply induction on $t$ to prove the statement. In case of $t=1$ we have by
Proposition 13, that

$\begin{array}[]{@{}l@{\;}l@{}}\mathbf{P}_{1}&=\mathbf{P}_{0}-\dfrac{\mathbf{P}_{0}\phi(\mathbf{x}_{k_{1*}})\phi(\mathbf{x}_{k_{1*}})^{T}\mathbf{P}_{0}}{||\mathbf{P}_{0}\phi(\mathbf{x}_{k_{1*}})||}=\mathbf{U}_{0}\mathbf{U}_{0}^{T}-\dfrac{\mathbf{U}_{0}\mathbf{U}_{0}^{T}\phi(\mathbf{x}_{k_{1*}})\phi(\mathbf{x}_{k_{1*}})^{T}\mathbf{U}_{0}\mathbf{U}_{0}^{T}}{||\mathbf{U}_{0}\mathbf{U}_{0}^{T}\phi(\mathbf{x}_{k_{1*})})||^{2}}\\\
&=\mathbf{U}_{0}\left(\mathbf{I}-\dfrac{\mathbf{q}_{1}\mathbf{q}_{1}^{T}}{||\mathbf{q}_{1}||^{2}}\right)\mathbf{U}_{0}^{T}=\mathbf{U}_{0}\mathbf{Q}(\mathbf{q}_{1})\mathbf{U}_{0}^{T}=\mathbf{U}_{0}\mathbf{Q}(\mathbf{q}_{1})\mathbf{Q}(\mathbf{q}_{1})\mathbf{U}_{0}^{T}=\framebox{$\mathbf{U}_{1}\mathbf{U}_{1}^{T}$}.\\\
\end{array}$ (22)

In transforming
$||\mathbf{U}_{0}\mathbf{U}_{0}^{T}\phi(\mathbf{x})_{t_{1*}}||^{2}$ into
$||\mathbf{q}_{1}||^{2}$ we exploited that $\mathbf{U}_{0}\mathbf{U}_{0}^{T}$
is a projection, hence it is idempotent. Let $t>1$ be arbitrary. Suppose that

$\begin{array}[]{l}\mathbf{P}_{t}=\mathbf{U}_{t-1}\mathbf{U}_{t-1}^{T}-\dfrac{\mathbf{U}_{t-1}\mathbf{q}_{t}\mathbf{q}_{t}^{T}\mathbf{U}_{t-1}^{T}}{||\mathbf{q}_{t}||^{2}}=\framebox{$\mathbf{U}_{t}\mathbf{U}_{t}^{T}$},\\\
\end{array}$ (23)

holds true. Now, computing the projector $t+1$ we obtain

$\begin{array}[]{@{}l@{\;}l@{}}\mathbf{P}_{t+1}&=\mathbf{P}_{t}-\dfrac{\mathbf{P}_{t}\phi(\mathbf{x}_{k_{(t+1)*}})\phi(\mathbf{x}_{k_{(t+1)*}})^{T}\mathbf{P}_{t}}{||\mathbf{P}_{t}\phi(\mathbf{x}_{k_{(t+1)*}})||^{2}}\\\
&=\mathbf{U}_{t}\mathbf{U}_{t}^{T}-\dfrac{\mathbf{U}_{t}\mathbf{U}_{t}^{T}\phi(\mathbf{x}_{k_{(t+1)*}})\phi(\mathbf{x}_{k_{(t+1)*}})^{T}\mathbf{U}_{t}\mathbf{U}_{t}^{T}}{||\mathbf{U}_{t}\mathbf{U}_{t}^{T}\phi(\mathbf{x}_{k_{(t+1)*}})||^{2}}=\mathbf{U}_{t}\left(\mathbf{I}-\dfrac{\mathbf{q}_{t+1}\mathbf{q}_{t+1}}{||\mathbf{q}_{t+1}||^{2}}\right)\mathbf{U}_{t}^{T}\\\
&=\mathbf{U}_{t}\mathbf{Q}(\mathbf{q}_{t+1})\mathbf{U}_{t}^{T}=\mathbf{U}_{t}\mathbf{Q}(\mathbf{q}_{t+1})\mathbf{Q}(\mathbf{q}_{t+1})\mathbf{U}_{t}^{T}=\framebox{$\mathbf{U}_{t+1}\mathbf{U}_{t+1}^{T}$}.\end{array}$

In the norm we again applied that $\mathbf{U}_{t}\mathbf{U}_{t}^{T}$ is
idempotent. ∎

1. Input: Output data matrix $\mathbf{Y}\in\mathbb{R}^{m\times n_{y}}$, input data matrix $\mathbf{X}\in\mathbb{R}^{m\times n_{x}}$, $D\leq n_{y}$ the number of variables to select, $\kappa$ kernel on variables 2. Output: Set $\mathcal{I}_{D}$ of indices of selected variables from $\mathbf{X}$ in the selection order. The algorithm: 1. Compute the cross-kernel matrix $\mathbf{K}_{YX}=[\kappa(y_{i},x_{j})],\ i\in[n_{y}],j\in[n_{x}]$. 2. Compute $\mathbf{K}_{YY}=[\kappa(y_{i},y_{j})],\ i,j\in[n_{x}]$; compute $\mathbf{K}_{YY}=\mathbf{V}_{Y}\mathbf{S}_{Y}^{2}\mathbf{V}_{Y}^{T}$ 3. Express $\mathbf{U}_{Y}^{T}\phi(X)^{T}$ as $\mathbf{R}^{(0)}=\mathbf{S}_{Y}^{+}\mathbf{V}_{Y}^{T}\mathbf{K}_{YX}$. 4. Let $t=0$, $\mathcal{I}_{t}=\emptyset$. 5. For $t\in\\{0,\dots,D-1\\}$ Do (a) If $t>0$ Then $\mathbf{R}^{(t)}=\mathbf{R}^{(t-1)}-\dfrac{\mathbf{q}_{t-1}\mathbf{q}_{t-1}^{T}\mathbf{R}^{(t-1)}}{||\mathbf{q}_{t-1}||^{2}}$ projection update (b) $k_{t}^{*}=\operatorname*{arg\,max}\limits_{j\in[n_{x}]\setminus\mathcal{I}_{t}}\left(\sum_{i=1}^{n_{y}}\left(\mathbf{R}^{(t)}\right)_{i,j}^{\circ 2}\right)$ (c) $\mathbf{q}_{t}=\mathbf{R}^{(t)}[:,k^{*}_{t}]$; $\mathcal{I}_{t}=\mathcal{I}_{t-1}\cup\\{k_{t}^{*}\\}$

Figure 3: Efficient implementation of the kernelized realization of supervised
variable selection by projection algorithm, ProjSe. Note the notation, e.g.
$\mathbf{R}^{(t)}=\mathbf{U}_{t}\phi(X)$.

We can express the main computation step, Step 2.b in Algorithm 2, by
exploiting the kernelized recursive iteration. From the sequential procedure
we can see that a key step of the computation is the calculation of the
vectors $\mathbf{q}_{i}$ via Equation (18),
$\mathbf{U}_{Y}^{T}\phi(\mathbf{x})=\mathbf{S}_{Y}^{-1}\mathbf{V}_{Y}^{T}\kappa(\mathbf{Y},\mathbf{x})$
for an arbitrary $\phi(\mathbf{x})\in\mathcal{H}$. In iteration $t$, we have

$\begin{array}[]{@{}l@{}l@{}}\mathbf{q}_{t+1}=\mathbf{U}_{t}^{T}\phi(\mathbf{x})=\left(\mathbf{U}_{Y}\prod_{s=1}^{t-1}\mathbf{Q}(\mathbf{q}_{s})\right)^{T}\phi(\mathbf{x})=\mathbf{Q}(\mathbf{q}_{t})\cdot\dots\cdot\mathbf{Q}(\mathbf{q}_{1})\underbrace{\mathbf{U}_{Y}^{T}\phi(\mathbf{x})}_{\mathbf{S}_{Y}^{-1}\mathbf{V}_{Y}^{T}\kappa(\mathbf{Y},\mathbf{x})}.\end{array}$
(24)

Taking advantage of the recursive definition of
$\mathbf{U}_{t}^{T}\phi(\mathbf{x})$ we also have that

$\begin{array}[]{ll}\mathbf{U}_{t+1}^{T}\phi(\mathbf{x})&=\mathbf{Q}(\mathbf{q}_{t+1})\mathbf{U}_{t}^{T}\phi(\mathbf{x})=\left(\mathbf{I}-\dfrac{\mathbf{q}_{t+1}\mathbf{q}_{t+1}^{T}}{||\mathbf{q}_{t+1}||^{2}}\right)\mathbf{U}_{t}^{T}\phi(\mathbf{x}),\end{array}$
(25)

where $\mathbf{q}_{t+1}=\mathbf{U}_{t}^{T}\phi(\mathbf{x}_{k_{(t+1)*}})$, thus
all terms relate to those computed in the previous iteration. The computation
of the norm $||\mathbf{q}_{t+1}||^{2}$ can also exploit the recursive nature
of the algorithm. Finally, all the feature representations $\phi(\mathbf{x})$
and $\phi(\mathbf{Y})$ are implicit, and are only expressed via kernel
evaluations since they only appear in inner products.

Based on these statements and Proposition 17 we can present a concrete
implementation of our algorithm in Figure 3. In the first step the kernels are
computed, where $\mathbf{K}_{YX}$ requires $O(mn_{y}n_{x})$, and
$\mathbf{K}_{Y}$ $O(mn_{y}^{2})$ operations in case of for example linear and
Gaussian kernels. For the eigenvalue decomposition of $\mathbf{K}_{Y}$ we need
$O(n_{y}^{3})$ operations, where $D\leq\min(n_{y},n_{x})$. In the algorithm,
the critical step is Step 4.a. Its complexity in step $t$ is
$O(n_{y}(n_{x}-t))$, thus, in general for selecting $D$ variables we need
$O(n_{y}n_{x}D)$ operations. Assuming that $m\gg n_{y},n_{x}$ , the dominating
part is the computation of the kernels, thus the entire complexity is equal to
$O(mn_{y}\max(n_{x},n_{y}))$.

##  5 Experiments

In this section we experimentally validate our approach.333The code for the
algorithm is available https://github.com/aalto-ics-kepaco/ProjSe. 444The
experiments are run on a machine with this parameters: 12th Gen Intel CoreTM
i5-12600K * 10. We first show demonstrations on our algorithm’s scalability on
synthetic data, before moving on to experimenting with real data and analysing
the stability of the feature selection.

###  5.1 Scalability demonstration with synthetic data

This test is implemented via a scheme presented by (27) in Figure 4 and by
(26). The components of the input matrix, $\mathbf{X}$ and the components of a
transformation matrix $\mathbf{W}$ are independently sampled from normal
distribution. Then output matrix is constructed, and finally random noise is
added to the output.

$\begin{array}[]{l@{\ }c@{\ }c@{\ }c@{\ }c}\text{Input}&&\text{Linear
transformation}&&\text{Noise}\\\
\hline\cr&&\mathbf{W}\sim[\mathcal{N}(0,\sigma)]^{n_{x}\times
n_{y}}&&\mathbf{E}\sim[\mathcal{N}(0,\sigma)]^{m\times n_{y}}\\\
&&\Downarrow&&\Downarrow\\\ \mathbf{X}\sim[\mathcal{N}(0,\sigma)]^{m\times
n_{x}}&\Longrightarrow&\mathbf{Y}=\mathbf{X}\mathbf{W}&\Longrightarrow&\tilde{\mathbf{Y}}=\mathbf{Y}+\mathbf{E}.\end{array}$
(26)

We apply ProjSe to this data with various sample sizes. Figure 4 presents the
dependence of the selection time on the sample size, where the maximum sample
size is $10$ million and the number of variables is $10$ – the variable
selection is performed in less than four seconds.

|  $\begin{array}[]{ll}\mathbf{Y}=\mathbf{X}\mathbf{W}+\mathbf{E}\\\
\hline\cr\mathbf{X}\sim[\mathcal{N}(0,\sigma)]^{m\times n_{x}},&m=1,\dots,10\
\text{million},\\\ \mathbf{W}\sim[\mathcal{N}(0,\sigma)]^{n_{x}\times
n_{y}},&n_{x}=100,\\\ \mathbf{E}\sim[\mathcal{N}(0,\sigma)]^{m\times
n_{y}},&n_{y}=100.\\\ \end{array}$ (27)
---|---
Figure 4: The dependence of the variable selection time on the sample size is
shown in seconds on the left, and the random data generation scheme applied is
on right

###  5.2 Feature selection from biological data

In this set of experiments, we compare our approach to [Brouard et al., 2022]
– a kernel-based feature selection method, where kernels are considered
traditionally on data samples instead of features. We experiment with the two
gene expression datasets, ”Carcinoma”, ”Glioma”, considered there for
unsupervised feature selection. While this kind of setting with one-view fat
data is not the one our method was developed for, as the scalability we
propose is first and foremost for large sample sizes, these experiments still
serve for illustrative comparison of feature selection performance.

Figure 5: Clustering results (NMI) on Carcinoma (right) and Glioma (left)
datasets with ProjSe as functions of number of variables chosen, averaged over
20 runs of k-means. ”Full” refers to results when full set of features is
used.  Table 2: ProjSe clustering results (NMI and time) with selected
features (10 or 300) compared to results reported in [Brouard et al., 2022].
Running time of ProjSe for choosing 300 features; running time and variation
of k-means is negligible.

| Carcinoma ($m$=174, $n_{x}$=9182, $C$=11) | Glioma ($m$=50, $n_{x}$=4434, $C$=4)
---|---|---
| NMI(10) | NMI(300) | t (s) | NMI(10) | NMI(300) | t (s)
lapl | 0.36 (0.02) | 0.64 (0.04) |  0.25 (0.04) |  0.50 (0.03) | 0.47 (0.06) |  0.02 (0.00)
NDFS | 0.22 (0.28) | 0.78 (0.03) | 6,162 (305) | 0.20 (0.04) | 0.36 (0.07) | 368 (21)
UKFS |  0.57 (0.03) | 0.75 (0.05) | 326 (52) | 0.26 (0.05) | 0.42 (0.05) | 23.74 (4.03)
ProjSe lin | 0.31 (0.01) | 0.58 (0.03) | 1,146111Running time of ProjSe for choosing 300 features; running time and variation of k-means is negligible. | 0.37 (0.01) |  0.52 (0.03) | 210111Running time of ProjSe for choosing 300 features; running time and variation of k-means is negligible.
ProjSe poly | 0.34 (0.01) | 0.79 (0.03) | 1,239111Running time of ProjSe for choosing 300 features; running time and variation of k-means is negligible. | 0.13 (0.02) | 0.34 (0.07) | 299111Running time of ProjSe for choosing 300 features; running time and variation of k-means is negligible.
ProjSe RBF | 0.33 (0.01) |  0.82 (0.03) | 1,263111Running time of ProjSe for choosing 300 features; running time and variation of k-means is negligible. | 0.38 (0.05) | 0.23 (0.06) | 284111Running time of ProjSe for choosing 300 features; running time and variation of k-means is negligible.

As the data is only available in one view in unsupervised setting, we apply
our algorithm by using the data in both views: as the reference/label view and
as the view the feature selection is performed on. Intuitively, this would
filter out the noise and redundant features in the view. In our method we
consider linear kernel, $k(\mathbf{x},\mathbf{z})=\mathbf{x}^{T}\mathbf{z}$,
polynomial kernel of degree 3,
$k(\mathbf{x},\mathbf{z})=(\mathbf{x}^{T}\mathbf{z})^{3}$, and RBF kernel,
$k(\mathbf{x},\mathbf{z})=\exp(\|\mathbf{x}-\mathbf{z}\|^{2}/(2\sigma^{2}))$
with the kernel parameter $\sigma$ set as mean of pairwise distances. We
assess the performance of the feature selection by measuring the normalised
mutual information (NMI) of k-means clustering results. Here the clusterer has
been given the amount of classes in the data as the number of clusters.

The results are displayed in Table 2, with comparison to selected methods from
[Brouard et al., 2022]: UKFS proposed there, as well as a scoring-based method
”lapl” [He et al., 2005] that performed well on Glioma dataset, and NDFS [Li
et al., 2012], a clustering-based approach that performed well with Carcinoma
dataset. Our method is very competitive with these, sometimes achieving better
performance. As in our method the kernel is calculated on features, the
running time is slower than for UKFS where the kernel on samples is
considered. However notably we are still competitive when compared to the
NDFS.

Additionally, Figure 2 displays more detailed clustering results with respect
to the number of variables chosen by ProjSe. These results also highlight the
differences that can be obtained by applying different kernels on the
features: with Carcinoma dataset the non-linear kernels, RBF and polynomial
kernel of degree 3, are clearly superior, while with Glioma linearity works
the best.

###  5.3 Feature selection from vector-valued output setting

Table 3: The time series classification datasets. Dtaset name | # tr. samples | # test samples | # features | # classes
---|---|---|---|---
Crop | 7200 | 16800 | 46 | 24
NonInvasiveFetalECGThorax1 | 1800 | 1965 | 750 | 42
ShapesAll | 600 | 600 | 512 | 60

Figure 6: Results with time series datasets. The top row reports the kernel
alignment of the input kernel with chosen variables (either RBF or linear) to
the ideal output kernel. The bottom row reports accuracy on test set, again
with both linear and RBF kernels for SVM; comparison is shown to randomly
selected features and to full feature set. The colors differentiate which
kernel is used on features in ProjSe, while the line style indicates if
traditional linear or RBF kernel is used on samples.

We next consider a setting more aligned with our method: supervised feature
selection with vector-valued output as the reference view. Here we consider
three datasets from the UEA & UCR Time Series Classification
Repository555http://www.timeseriesclassification.com/dataset.php, Crop,
NonInvasiveFetalECGThorax1 (”Thorax”), and ShapesAll, as detailed in Table 3.
These datasets are associated with multi-class classification tasks, and we
use the one-hot encoding of the labels as the vector-valued output to perform
the feature selection with ProjSe. As before, we consider linear, polynomial
and RBF kernels. We assess the success of the feature selection task by
performing classification with SVM with the selected features - here we
consider both linear and RBF kernels on the data samples.

The results are displayed in Figure 6, where both kernel alignment
($KA(\mathbf{K},\mathbf{K}^{\prime})=\langle\mathbf{K}_{c},\mathbf{K}^{\prime}_{c}\rangle_{F}/(\|\mathbf{K}_{c}\|_{F}\|\mathbf{K}^{\prime}_{c}\|_{F})$
where $c$ denotes centering) to the linear kernel on the one-hot-encoded
outputs, and accuracy of SVM classification are shown. The different kernels
used in feature selection give slightly different features; however the
performance on the subsequent classification task is mostly dependent on which
kernel is used on the samples. Especially for Thorax and ShapesAll datasets
with higher dimensionality, it can be seen that all the ProjSe results with
linear SVM outperform using the full set of features.

###  5.4 Experiments with two-view data

Table 4: The detailed computation times in seconds for the variable selection method, where 10, 20, 50 and 100 variables are extracted. |  | MediaMill | Cifar 100
---|---|---
| Linear | RBF | Linear | RBF
Computing $K_{yx}$ | 0.015 | 0.019 | 0.221 | 0.470
Computing $K_{yy}$ | 0.007 | 0.015 | 0.006 | 0.025
Eigen decomp. of $K_{yy}$ | 0.003 | 0.001 | 0.001 | 0.001
Centralization of $\mathbf{K}_{xx}$ 111Optional | 0.005 | 0.014 | 1.462 | 2.065
10 variables | 0.035 | 0.045 | 1.848 | 2.640
20 variables | 0.024 | 0.048 | 1.793 | 2.657
50 variables | 0.037 | 0.050 | 1.804 | 2.671
100 variables | 0.049 | 0.061 | 1.833 | 2.706

In our last set of experiments, we consider the following datasets:

  * •

MNIST handwritten digits [LeCun et al., 1998, LeCun, 1998]: This dataset
contains 60000 training and 10000 test examples of handwritten digits in
greyscale. The number of pixels in each image is $28\times 28=784$, resulting
in total to 784 variables. To construct the two sets of variables, the image
columns are split into half similarly as in [Andrew et al., 2013]. Thus both
views comprise of 392 variables.

  * •

MediaMill dataset [Snoek et al., 2006]: This dataset contains 43907 examples
which are extracted from keyframes of video shots. There are two views in this
data: text annotations (101 variables) and visual features (120 variables).

  * •

Cifar100 dataset [Krizhevsky, 2009]: This dataset, chosen to demonstrate the
scalability of ProjSe, contains 50000 training and 10000 test examples of
color images. The number of pixels in each image is $32\times 32=1024$, where
to each pixel 3 colors are assigned. The examples belong to 100 classes, where
each class contains 500 training and 100 test examples. The classes are
represted by indicator vectors.

We perform variable selection independently on both views in the datasets.
After the variable selection is performed on both sides, we compute canonical
correlations between all subset pairs of the extracted variables, starting
from the first ones and incrementally growing them to the entire sets. To
demonstrate the performance of the proposed variable selection algorithm
ProjSe, it is compared to the following methods: large-scale sparse kernel
canonical correlation analysis (GradKCCA), [Uurtio et al., 2019], deep
canonical correlation analysis (DCCA), [Andrew et al., 2013], randomized non-
linear CCA (RCCA), [Lopez-Paz et al., 2014], kernel non-linear orthogonal
iterations (KNOI), [Wang and Livescu, 2015], and CCA through Hilbert-Schmidt
independence criterion (SCCA-HSIC), [Uurtio et al., 2018].

These CCA variants explicitly or implicitly rely on singular value
decomposition, and their performance highly depends on the distribution of the
singular values of the data matrices. Since the data matrices have small
number of dominating singular values, we expect from a variable selection
method that it can capture a relatively small set of variables to reproduce
similar accuracy, measured in canonical correlation. We need to bear in mind
that the CCA-based methods add up all information represented by the full
collection of variables, however the projector based selection only relies a
relatively small subset of those variables.

(a)  MNIST dataset

(b)  MediaMill dataset

Figure 7: Variable selection results w.r.t the number of selected variables.  Table 5: CCA comparison on the MNIST and MediaMill datasets.  |  | MNIST | MediaMill
---|---|---
| $\rho_{\text{TEST}}$ | TIME (s) | $\rho_{\text{TEST}}$ | TIME (s)
Generic CCA | 0.923 | 2.40 | 0.675 | 0.429
GradKCCA | 0.952 $\pm$ 0.001 | 56 $\pm$6 | 0.657 $\pm$ 0.007 | 8 $\pm$4
DCCA | 0.943 $\pm$ 0.003 | 4578 $\pm$203 | 0.633 $\pm$ 0.003 | 1280 $\pm$112
RCCA | 0.949 $\pm$ 0.010 | 78 $\pm$13 | 0.626 $\pm$ 0.005 | 23 $\pm$9
KNOI | 0.950 $\pm$ 0.005 | 878 $\pm$62 | 0.645 $\pm$ 0.003 | 218 $\pm$73
SCCA-HSIC | 0.934 $\pm$ 0.006 | 5611 $\pm$193 | 0.625 $\pm$ 0.002 | 1804$\pm$143
ProjSe 10 var. | 0.847 |  | 0.542 |
ProjSe 20 var. | 0.890 |  | 0.586 |
ProjSe 50 var. | 0.918 |  | 0.631 |
ProjSe 100 var. | 0.935 | 0.45111The algorithm was run once with 100 variables, the smaller amounts were subsampled  | 0.672 | 0.041111The algorithm was run once with 100 variables, the smaller amounts were subsampled

Figure 7 shows the performance of ProjSe on MNIST and MediaMill datasets with
both linear and Gaussian kernels, as functions of the number of selected
variables. The results are measured by canonical correlation between the
subsets of variables selected from the two views. Comparing the results to the
other CCA methods in Table 5 (taken from [Uurtio et al., 2019]), we observe
ProjSe obtaining comparable performance after 20 or 50 selected variables,
while being orders of magnitude faster than the other methods. Since ProjSe is
fully deterministic, there is no variance is reported for it. The heaviest
computation for ProjSe is in the beginning when eigenvalue decomposition of
$\mathbf{K}_{YY}$ is calculated (see Table 4). Thus the running time varies
only minimally when different number of variables is selected. This is also
demonstrated in 4 where the running times for MNIST and Cifar100 datasets are
detailed.

|
---|---
Figure 8: Means of the two stability scores on random sub-samples computed on
the data sets the MNIST (left) and the MediaMill (right), with linear kernel

As a variable selection method, we are interested in evaluating the stability
of the selection process. In order to measure this, we here consider the
stability index [Nogueira et al., 2018] and the average Pearson correlation of
relevance [Hamer and Dupont, 2021]. In both measures, higher values indicate
higher stability; the maximum value in both is 1.

First the number of extracted variables is chosen from $(1,2,\dots,5)$ in
MNIST, and from $(10,20,\dots,50)$ in MediaMill. For each number of selected
variables the subsamples are taken with the following percentage of the entire
training sets: $(10\%,20\%,\dots,50\%)$. Then random subsets are extracted
$10$ times of the size given above. The scores are computed for each number of
selected variables, and for each subsample size and finally averaged on all
random samples. They are shown in Figure 8, where the averages for all pairs
of subsets of variables and for all subsample sizes are presented.

##  6 Conclusion

In this paper we introduced a novel variable selection method for two-view
settings. Our method is deterministic, selecting variables based on
correlation defined with projection operators. The kernelised formulation of
our approach paves way for efficient and highly scalable implementation,
allowing the application of our method to datasets with millions of data
samples. We empirically demonstrated this efficiency and the suitability of
our approach for feature selection task, with both synthetic and real data.

## Declarations

The authors wish to acknowledge the financial support by Academy of Finland
through the grants 334790 (MAGITICS), 339421 (MASF) and 345802 (AIB), as well
as the Global Programme by Finnish Ministry of Education and Culture

## References

  * [Aghazadeh et al., 2018] Aghazadeh, A., Spring, R., LeJeune, D., Dasarathy, G., Shrivastava, A., et al. (2018).  Mission: Ultra large-scale feature selection using count-sketches.  In ICML, pages 80–88. PMLR.
  * [Anderson. and Duffin, 1969] Anderson., J. W. N. and Duffin, R. J. (1969).  Series and parallel addition of matrices.  SIAM J. Appl. Math., 26:576–594.
  * [Andrew et al., 2013] Andrew, G., Arora, R., Bilmes, J., and Livescu, K. (2013).  Deep canonical correlation analysis.  In Dasgupta, S. and McAllester, D., editors, Proceedings of the 30th ICML, volume 28(3) of Proceedings of Machine Learning Research, pages 1247–1255, Atlanta, Georgia, USA. PMLR.
  * [Anette and Nokto, 2018] Anette, K. and Nokto, D. (2018).  A benchmark of prevalent feature selection algorithms on a diverse set of classification problems.  Master’s thesis.
  * [Ben-Israel, 2015] Ben-Israel, A. (2015).  Projectors on intersections of subspaces.  Contemporary Mathematics, 636.
  * [Ben-Israel and Greville, 2003] Ben-Israel, A. and Greville, T. N. (2003).  Generalized inverses: Theory and applications.  Springer New York, NY, second edition.
  * [Bie et al., 2005] Bie, T. D., Cristianini, N., and Rosipal, R. (2005).  Eigenproblems in pattern recognition.  In Handbook of Geometric Computing : Applications in Pattern Recognition, Computer Vision, Neuralcomputing, and Robotics, pages 129–170. Springer-Verlag.
  * [Bommert et al., 2020] Bommert, A., Sun, X., Bischl, B., Rahnenführer, J., and Lang, M. (2020).  Benchmark for filter methods for feature selection in high-dimensional classification data.  Computational Statistics & Data Analysis, 143(C).
  * [Boyd et al., 2018] Boyd, N., Hastie, T., Boyd, S., Recht, B., and Jordan, M. I. (2018).  Saturating splines and feature selection.  JMLR, 18(197):1–32.
  * [Brogat-Motte et al., 2022] Brogat-Motte, L., Rudi, A., Brouard, C., Rousu, J., and d’Alché-Buc, F. (2022).  Vector-valued least-squares regression under output regularity assumptions.  JMLR, 23(344):1–50.
  * [Brouard et al., 2022] Brouard, C., Mariette, J., Flamary, R., and Vialaneix, N. (2022).  Feature selection for kernel methods in systems biology.  NAR genomics and bioinformatics, 4(1):lqac014.
  * [Conway, 1997] Conway, J. B. (1997).  Course in Functional Analysis.  Springer, 2nd edition.
  * [Deshwal et al., 2019] Deshwal, A., Doppa, J. R., and Roth, D. (2019).  Learning and inference for structured prediction: A unifying perspective.  In IJCAI-19.
  * [Fan et al., 2009] Fan, J., Samworth, R., and Wu, Y. (2009).  Ultrahigh dimensional feature selection: Beyond the linear model.  JMLR, 10:2013–2038.
  * [Golub and Loan, 2013] Golub, G. H. and Loan, C. F. V. (2013).  Matrix Computations.  The Johns Hopkins University Press, Baltimore, MD, 4th edition.
  * [Halperin, 1962] Halperin, I. (1962).  The product of projection operators.  Acta Sei. Math. (Szeged), 23:96–99.
  * [Hamer and Dupont, 2021] Hamer, V. and Dupont, P. (2021).  An importance weighted feature selection stability measure.  JMLR, 22(116):1–57.
  * [Hayashi and Nagaoka, 2002] Hayashi, M. and Nagaoka, H. (2002).  A general formula for the classical capacity of a general quantum channel.  In Proceedings IEEE International Symposium on Information Theory,. IEEE.
  * [He et al., 2005] He, X., Cai, D., and Niyogi, P. (2005).  Laplacian score for feature selection.  Advances in neural information processing systems, 18.
  * [Hofmann et al., 2008] Hofmann, T., Schölkopf, B., and Smola, A. J. (2008).  Kernel methods in machine learning.  The annals of statistics, pages 1171–1220.
  * [Hotelling, 1936] Hotelling, H. (1936).  Relations between two sets of variates.  Biometrika, 28(3/4):321.
  * [Jordan et al., 2021] Jordan, M. I., Liu, K., and Ruan, F. (2021).  On the self-penalization phenomenon in feature selection.  arXiv preprint arXiv:2110.05852.
  * [Kreyszig, 1989] Kreyszig, E. (1989).  Introductory Functional Analysis with Applications.  Wiley.
  * [Krizhevsky, 2009] Krizhevsky, A. (2009).  Learning multiple layers of features from tiny images.  Technical report, University of Totonto.
  * [LeCun, 1998] LeCun, Y. (1998).  The mnist database of handwritten digits.  exdb.
  * [LeCun et al., 1998] LeCun, Y., Bottou, L., Bengio, Y., and Haffner, P. (1998).  Gradient-based learning applied to document recognition.  In Proceedings of the IEEE, volume 86(11), pages 2278–2324.
  * [Li et al., 2017] Li, J., Cheng, K., Wang, S., Morstatter, F., Trevino, R. P., Tang, J., and Liu, H. (2017).  Feature selection.  ACM Computing Surveys, 50(6):1–45.
  * [Li et al., 2012] Li, Z., Yang, Y., Liu, J., Zhou, X., and Lu, H. (2012).  Unsupervised feature selection using nonnegative spectral analysis.  In Proceedings of the AAAI, volume 26, pages 1026–1032.
  * [Lopez-Paz et al., 2014] Lopez-Paz, D., Sra, S., Smola, A., Ghahramani, Z., and Schoelkopf, B. (2014).  Randomized nonlinear component analysis.  In Xing, E. P. and Jebara, T., editors, Proceedings of the 31st ICML, volume 32(2) of Proceedings of Machine Learning Research, pages 1359–1367, Bejing, China. PMLR.
  * [Micchelli and Pontil, 2005] Micchelli, C. A. and Pontil, M. (2005).  On learning vector-valued functions.  Neural computation, 17(1):177–204.
  * [Minh et al., 2016a] Minh, H. Q., Bazzani, L., and Murino, V. (2016a).  A unifying framework in vector-valued reproducing kernel hilbert spaces for manifold regularization and co-regularized multi-view learning.  volume 17, pages 769–840. JMLR. org.
  * [Minh et al., 2016b] Minh, H. Q., Biagio, M. S., Bazzani, L., and Murino, V. (2016b).  Approximate log-hilbert-schmidt distances between covariance operators for image classification.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 5195–5203.
  * [Muandet et al., 2017] Muandet, K., Fukumizu, K., Sriperumbudur, B., Schölkopf, B., et al. (2017).  Kernel mean embedding of distributions: A review and beyond.  Foundations and Trends® in Machine Learning, 10(1-2):1–141.
  * [Nielsen and Chuang, 2000] Nielsen, M. and Chuang, I. (2000).  Quantum Computation and Quantum Information.  Cambridge University Press.
  * [Nogueira et al., 2018] Nogueira, S., Sechidis, K., and Brown, G. (2018).  On the stability of feature selection algorithms.  JMLR, 18(174):1–54.
  * [Rasmussen and Williams, 2005] Rasmussen, C. E. and Williams, C. K. I. (2005).  Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning).  The MIT Press.
  * [Schölkopf et al., 1998] Schölkopf, B., Smola, A., and Müller, K.-R. (1998).  Nonlinear Component Analysis as a Kernel Eigenvalue Problem.  Neural Computation, 10(5):1299–1319.
  * [Snoek et al., 2006] Snoek, C. G. M., Worring, M., van Gemert, J. C., Geusebroek, J.-M., and Smeulders, A. W. M. (2006).  The challenge problem for automated detection of 101 semantic concepts in multimedia.  In Proceedings of the 14th ACM International Conference on Multimedia, MM ’06, page 421–430.
  * [Song et al., 2012] Song, L., Smola, A., Gretton, A., Bedo, J., and Borgwardt, K. (2012).  Feature selection via dependence maximization.  JMLR, 13(1):1393–1434.
  * [Tjur, 1984] Tjur, T. (1984).  Analysis of variance models in orthogonal designs.  International Statistical Review, Revue Internationale de Statistique, 52:33–65.
  * [Uurtio et al., 2018] Uurtio, V., Bhadra, S., and Rousu, J. (2018).  Sparse non-linear cca through hilbert-schmidt independence criterion.  In 2018 IEEE International Conference on Data Mining (ICDM), pages 1278–1283, Los Alamitos, CA, USA. IEEE Computer Society.
  * [Uurtio et al., 2019] Uurtio, V., Bhadra, S., and Rousu, J. (2019).  Large-scale sparse kernel canonical correlation analysis.  In Chaudhuri, K. and Salakhutdinov, R., editors, Proceedings of the 36th ICML, volume 97 of Proceedings of Machine Learning Research, pages 6383–6391. PMLR.
  * [Valcárcel et al., 2022] Valcárcel, L. V., San José-Enériz, E., Cendoya, X., Rubio, A., Agirre, X., Prósper, F., and Planes, F. J. (2022).  Boso: A novel feature selection algorithm for linear regression with high-dimensional data.  PLOS Computational Biology, 18(5):1–29.
  * [von Neumann, 1950] von Neumann, J. (1950).  Functional Operators. Vol II: The Geometry of Orthogonal Subspaces, volume 39.  Annals of Math. Studies, Princeton University Press, Princeton.
  * [Wang and Livescu, 2015] Wang, W. and Livescu, K. (2015).  Large-scale approximate kernel canonical correlation analysis.  CoRR, abs/1511.04773.
  * [Xu et al., 2013] Xu, C., Tao, D., and Xu, C. (2013).  A survey on multi-view learning.  arXiv preprint arXiv:1304.5634.
  * [Zebari et al., 2020] Zebari, R., Abdulazeez, A., Zeebaree, D., Zebari, D., and Saeed, J. (2020).  A comprehensive review of dimensionality reduction techniques for feature selection and feature extraction.  J. Appl. Sci. Technol. Trends, 1(2):56–70.


# Accurate and precise measurement of the wavelength of
4d${}^{4}D_{7/2}\rightarrow$5p${}^{4}P^{\circ}_{5/2}$ transition of Kr II

Y. Dancheva  The author to whom correspondence may be addressed:
<EMAIL_ADDRESS>P. Coniglio  D. Pagano  A. Garde  F.
Scortecci  Aerospazio Tecnologie Srl, via dei Tessili 99, 53040 Rapolano Terme
(SI), Italy

###### Abstract

Electric propulsion requires exhaustive ground test campaigns to obtain an
accurate characterization of the propulsion devices, or thrusters, used by the
spacecraft. Among the many plasma parameters, accurately measured during the
tests, that of the ion velocity is key, and can be measured using non-
intrusive tools such as Laser-Induced Fluorescence (LIF) diagnostics. The ion
velocity is inferred by Doppler shift measurements that presupposes a precise
and accurate knowledge of the wavelength of excitation of the ions at rest.
Today electric propulsion is moving towards the use of Krypton as a
propellant, due to the dramatic increase in cost of the more advantageous
Xenon gas propellant, commonly used until now. The utilization of Krypton
implies that LIF diagnostic tool used be adapted accordingly. In this work an
accurate measurement of the wavelength of
4d${}^{4}D_{7/2}\rightarrow$5p${}^{4}P^{\circ}_{5/2}$ transition of Kr II (one
of the most advantageous transitions for LIF diagnostics) is performed in the
plasma of a hollow cathode acting as a neutralizer. These measurements will
significantly improve accuracy in the determination of the ion velocity.

##  I Introduction

Traditionally, Xenon has been the favorite propellant gas for electric
propulsion applications as it provides an optimal compromise between
performance and ease of handling. Although Xenon has numerous technical
advantages, its relative scarcity and resulting increase in price are now
posing considerable budgetary constrains for its use, especially for deep
space exploration Tirila_2023 ; Nakles_2011 . Therefore, a number of
alternative propellants are being examined intensively for use in
electrostatic spacecraft propulsion thrusters Duchemin_2009 . A more economic
alternative to Xe needs to be identified, while maintaining thruster
performance levels, and without overlooking performance optimization aspects.
Among various types of electric thrusters, the Hall Effect Thruster (HET) is
one of the most widely used owing to its simple structure, high reliability,
and long life Mazouffre_2016 ; Levchenko_2018 .

Like Xenon, Krypton is a noble gas and can be integrated easily into existing
spacecraft propellant management systems with minimal modifications. As this
gas has similar ionization potential it should not dramatically affect
thruster efficiency, and the lower atomic mass could possibly produce a 25%
increase in specific impulse due to the increased propellant exit velocity of
lighter ions Linnell_2007 ; Hargus_2012 . The higher specific impulse also
provides advantages in space missions such as maintaining spacecraft orbit
Kurzyna_2018 ; Lim2019 . A lower anode efficiency of approximately 5-15$\%$ is
expected with respect to Xenon propellant at the same operating conditions
Hargus_2012 ; Kurzyna_2018 ; Jorns_2022 . Krypton is also about ten times more
common in the atmosphere than Xenon, and hence less expensive. Currently,
Krypton is the most widely-used alternative propellant in HET both for space
applications (for example the Starlink constellation), as well as ground
testing Tirila_2023 . Table 1 summarizes the properties of Xenon and Krypton
specifically relevant to electrostatic spacecraft propulsion.

Table 1: Kr and Xe properties comparison Property | Xe | Kr
---|---|---
Atomic mass | 131.293 amu | 83.798 amu
First ionization energy | 12.1 eV | 14.0 eV
Atmospheric |  |
concentration | 87 ppb | 1000 ppb
Stable isotopes | 9 | 6
Odd isotopes | 2 | 1
Critical pressure | 57.65 atm | 54.3 atm
Critical temperature | 290 K | 209 K
Boiling point (1 atm) | 161 K | 120 K

The characterization and qualification of ion thrusters requires extensive,
long-duration test campaigns in space simulators. A valuable tool for
characterizing thruster performance and comparing this with already-developed
models is Light-Induced Fluorescence (LIF) spectroscopy, which measures the
velocity of ions employing the Doppler effect. Many excellent works on LIF
characterization of ion thrusters operating with Xe propellant are available
in the literature, but only few of them consider operation with Krypton.
Advantageous transitions for LIF diagnostics are characterized by large
intensity and narrow intrinsic width (transitions involving long-living
levels). The Xe transition at 834 nm
(5d${}^{2}[4]_{7/2}\rightarrow$6p${}^{2}[3]^{\circ}_{5/2}$) is used for LIF
diagnostics by the majority of the electric propulsion research community. As
to Kr, such favoured transitions are at 820 nm (the
4d${}^{4}F_{7/2}\rightarrow$5p${}^{2}D^{\circ}_{5/2}$ transition) Lejeune_2012
and at 729 nm (the 4d${}^{4}D_{7/2}\rightarrow$5p${}^{4}P^{\circ}_{5/2}$
transition) Hargus_2011 ; Buorgeois_2011 . An accurate knowledge of the
features of the transition used for LIF spectroscopy is required for better
determination of the ion velocity.

The air wavelength of the Kr II
4d${}^{4}D_{7/2}\rightarrow$5p${}^{4}P^{\circ}_{5/2}$ transition has been
measured to be 728.982 nm with an accuracy of 0.7 pm Dzierżęga_2001 ;
Saloman_2007 . Such accuracy would introduce an error of about 300 m/s in the
determination of the Kr ion velocity. In this work, the vacuum/air wavelength
of the Kr II transition 5d${}^{4}D_{7/2}\rightarrow$5p${}^{4}P^{\circ}_{5/2}$
is measured with an improved accuracy (34 fm). For this purpose, the LIF
spectroscopy is applied in the plasma of a hollow cathode discharge, normally
used as a neutralizer in electric propulsion.

##  II Set-up

The LIF is a particularly powerful optical-diagnostic tool providing valuable,
species selective information with excellent spatial resolution. A highly
coherent, single-frequency, laser source that can be tuned over a broad
spectral interval permits to establish a versatile laser-plasma interaction
tool. Given its remote detection nature, LIF diagnostic is minimally invasive.
The plasma is interrogated locally with an electromagnetic tool: all the
optical devices used to illuminate the plasma and to detect its fluorescence
are placed far from the interaction point, can have dimensions in the mm
scale, and the interaction volume can have sizes as small as few millimeters.
The ion velocity distribution function (IVDF) can be measured along different
directions with excellent spatial and time resolution. The measurement
procedure is based on detecting the intensity of the fluorescence emitted by
the plasma in the interaction point in response to the illumination by means
of narrow-band radiation characterized by an accurately and precisely known
and scanned frequency.

###  II.1 The test facility

The measurements shown in this work are conducted in a non-magnetic stainless
steel vacuum chamber. Vacuum is obtained by a single stage cryogenic panel.
The base pressure of the vacuum chamber is as low as $10^{-7}$mbar and
increases up to about $10^{-5}$mbar during the cathode operation. The hollow
cathode is operated in diode regime (see Fig.1) at currents of about 3 A with
an anode positioned 20 mm away from it. The measurements are performed
positioning the LIF equipment outside the vacuum chamber illuminating the
detection region and detecting the emitted fluorescence through an optical
access (view port).

Figure 1: Schematics of the laser excitation and fluorescence collection
arrangement. FAS - fluorescence analyzing system, FCO - fluorescence
collection objective, VP - view port of the vacuum chamber, SMF, MMF - single-
mode and multi-mode fibers and AOM - acousto-optical modulator.

###  II.2 The LIF set-up

The light source used to excite the
4d${}^{4}D_{7/2}\rightarrow$5p${}^{4}P^{\circ}_{5/2}$ Kr II transition is a
tunable diode laser in a primary-and-secondary configuration (the primary
laser is an extended-cavity diode laser and the secondary laser is a tapered
amplifier). The laser linewidth is better than 1 MHz and the power available
is in excess.

Figure 2: Schematic of the LIF set-up: Proportional-Integral-Derivative (PID)
controller; acousto-optical modulator (AOM); photo-multiplier tube (PMT); and
data acquisition card (DAQ).

A schematic of the LIF set-up is given in Fig.2 (a more detailed description
can be found in Refs.Dancheva_2022 ; Dancheva_2023 ). The laser wavelength is
locked and scanned with the help of a high accuracy ($\pm~10~$MHz) wavelength
meter and a proportional-integral-derivative (PID) controller. The wavelength
meter is calibrated periodically using a diode laser that is frequency
stabilized to the absorption profile of the Caesium D2 line and typically
drifts below 2 MHz/day.

The detection region is near to the cathode exit plane at the cathode tip
center, where the laser beam and the fluorescence collection objective (FCO)
view regions are crossing (see Fig.1). The laser beam is aligned at a small
angle with respect to the normal to cathode exit plane ($\alpha=6^{\circ}$ in
Fig.1) to excite the plasma at the tip center. The spatial resolution is
provided by the laser intensity amplitude modulation at 3 kHz by a fibre-
coupled, acousto-optical modulator (AOM). The AOM output is coupled to a
single-mode optical fiber (see Fig.2) that provides a laser beam with waist
diameter of 0.3 mm using a fiber collimator positioned at 0.6 m away from the
cathode.

The view spot of the FCO is about 2.5 mm in diameter at the region of
detection. The fluorescence signal at air wavelength of 473.9 nm
(corresponding to the transition
5p${}^{4}P^{\circ}_{5/2}\rightarrow$5s${}^{4}P_{5/2}$) is selected using a
grating monochromator. Subsequently, the signal at 473.9 nm is detected using
a photo-multiplier tube (PMT). The LIF signal demodulation is performed by a
lock-in amplifier, whose input is the PMT current, converted to voltage by a
trans-impedance amplifier. The poor signal-to-noise ratio of the LIF spectra
is counteracted by using a long lock-in integration time, implying a slow scan
of the laser wavelength and averaging over many measurements. Indeed, the
lock-in amplifier integration time is 300 ms, which with a third order low-
pass filter at the output determines a bandwidth of about 0.3 Hz .

The laser vacuum wavelength is locked and scanned by the wavelength meter
using a PID control (the schematics of the set-up are shown in Fig.2). A scan
of 20 pm (about 11 GHz) is performed backwards and forwards using the
wavelength meter in an overall time of 100 s. The bidirectional scanning
permits to account for possible shifts of the spectral line due to integration
time.

The wavelength meter used in this work is temperature and pressure stabilized.
A variation less than 50 MHz is expected for temperature variations of about
10∘C, that are well compensated for rates lower than
$\nicefrac{{2$^{\circ}$C}}{{hour}}$, which is higher than the typical
temperature excursion in the laboratory.

##  III The Kr II 4d4D${}_{7/2}\rightarrow$5p4P${}^{\circ}_{5/2}$ line.

In this analysis the Kr isotopes with concentration higher than 0.3$\%$ are
considered and listed in Table 2. A theoretical line-shape is calculated and
used to fit the experimental data.

Table 2: Isotope shift and hyperfine coefficients | Relative shift |
---|---|---
Isotope | with respect to | Abundance $\alpha$
|  84Kr (MHz) | ($\%$)
78 | 1185.7 | 0.35
80 | 768.0 | 2.27
82 | 372.6 | 11.56
83 | 175.0 | 11.55
84 | 0 | 56.9
86 | -365.2 | 17.37
Electronic | A coefficient | B coefficient
state | (MHz) | (MHz)
4d4D7/2 | -43.513 | -294.921
5p4P5/2 | -167.2 | +91

Only one isotope (83Kr) has a non zero nuclear spin (I=9/2) and thus its
hyperfine structure should be included in the line-shape analysis. The
hyperfine constants of both the states of the 83Kr II transition have been
measured previously Scholl_1986 ; Schuessler_1992  and are given in Table 2.

Table 3: Kr hyperfine structure components |  | Relative | Relative |
---|---|---|---|---
|  | shift | shift $\delta_{i}$ | Intensity
F | F′ | (with respect | (with respect |  $\beta$($\%$)
|  | to J$\rightarrow$J′) (MHz) | to 84Kr) (MHz) |
8 | 7 | -1099.2 | -924.2 | 0.21
7 | 7 | -1559.6 | -1384.6 | 0.034
7 | 6 | -431.7 | -256.7 | 0.15
6 | 7 | -1901.1 | -1726.1 | 0.0031
6 | 6 | -773.2 | -598.2 | 0.055
6 | 5 | 223.2 | 398.2 | 0.10
5 | 6 | -1020.2 | -845.2 | 0.0083
5 | 5 | -23.8 | 151.2 | 0.064
5 | 4 | 827.3 | 1002.3 | 0.065
4 | 5 | -197.5 | -22.5 | 0.015
4 | 4 | 653.7 | 828.7 | 0.063
3 | 4 | 535.8 | 710.8 | 0.022
3 | 3 | 1230.4 | 1405.4 | 0.052
3 | 2 | 1759.3 | 1934.3 | 0.013
2 | 3 | 1154.3 | 1329.3 | 0.03
2 | 2 | 1683.2 | 1858.2 | 0.033
1 | 2 | 1638.3 | 1813.3 | 0.038

As can be seen from Fig.3 the lower 83Kr II state splits in 8 hyperfine
sublevels and the upper one in 6 with total angular momenta as denoted on the
energy schematic. The electric-dipole allowed 18 transitions are shown in the
plot. A detailed description about the hyperfine structure and the relevant
selection rules can be found elsewhere Arimondo_1977 ; Hargus_2010 ;
Hargus_2001 . The spectrum of the
4d${}^{4}D_{7/2}\rightarrow$5p${}^{4}P^{\circ}_{5/2}$ transition at Doppler
linewidth of 35 MHz (to highlight the line structure) and 1900 MHz (to show
the line-shape at ion temperatures near to the measured one) are reported in
Fig.3.

Figure 3: The transition line-shape at two different temperature/Doppler
widths. The hyperfine structure of 83Kr is shown in the insert.

##  IV Results

Fig.4 shows the spectrum of the
4d${}^{4}D_{7/2}\rightarrow$5p${}^{4}P^{\circ}_{5/2}$ recorded when scanning
the laser frequency in both directions. Each trace is obtained as an average
over 36 laser wavelength scans and the cathode is kept at constant operational
parameters during all the measurements.

Figure 4: Kr II line at about 729 nm when scanning the laser wavelength in
both directions and measuring the vacuum a) and the air wavelength c). For the
sake of convenience, the measurements are also presented as a function of the
frequency - vacuum case, and as a function of $\nu n(p,T,\lambda)$ \- air
case, where n is the index of refraction of the dry air at the wavelength of
interest. Each measurement is an average over 36 laser scans. The residual is
shown for each case in plots b) and d).

The spectrum is registered measuring both the vacuum and the air wavelength
(at standard dry air: T=15∘C and p=760 mmHg). The experimental data are fitted
using an orthogonal distance regression algorithm Pallavi_2022  considering
the wavelength meter accuracy and the LIF signal noise outside the spectral
line as errors.

The fitting function is given by:

$\begin{split}G(\nu)=A\sum_{i}B_{i}\exp^{-4\ln{2}\big{(}\frac{\nu-(\nu_{0}\pm\delta_{i})}{\Delta\nu_{D}}\big{)}^{2}}+D,\\\
~\Delta\nu_{D}=\frac{\nu_{0}}{c}\sqrt{8\ln{2}\frac{k_{B}T}{m}},\end{split}$
(1)

where the $\delta_{i}$ is the relative shift for each spectral component
(given in Table 3), $B_{i}$ are the relative intensities
($B_{i}$=$\alpha_{i}\times\beta_{i}$), $k_{B}$ is the Boltzmann constant, $T$
is the ion temperature, $m$ is the ion mass, $A$ is a scaling factor, $c$ is
the speed of light, and $D$ accounts for possible data offset is used to fit
the experimental data (see Fig.4). Among the four fitting parameters ($A$,
$\nu_{0}$, $T$, and $D$) the line central frequency $\nu_{0}$ and the ion
temperature $T$ are of interest. The residuals, given in Fig.4, are indicative
for the fit goodness.

Figure 4 shows the measured $\nu_{0}$ when scanning the laser wavelength in
both directions. A shift, occurring due to the finite lock-in settling time of
about 100 MHz is estimated from the fitted $\nu_{0}$ values. As this
instrumental shift occurs symmetrically in both cases, the exact value of the
$\nu_{0}$ has to be estimated as the average of the two results.

The $\nu_{0}$ value is temperature dependent through the hyperfine and isotope
structure contribution. A negative systematic error in the determination of
$\nu_{0}$ is introduced at different Doppler widths (ion temperatures). The
calculated shift is given in Fig.5.

Figure 5: Kr line center shift as a function of the ion temperature (Doppler
width).

Considering a coverage factor of 2, for a confidence level of approximately
95$\%$, the measured central wavelength $\lambda_{0}$ and frequency $\nu_{0}$
are given in Table 4. Taking into account a Doppler width of 1900 MHz the
corrected $\lambda_{0}$ values are also provided.

Table 4: Measured central wavelength and frequency. HFS - hyperfine structure. | $\lambda_{0}$ | Uncertainty
---|---|---
Condition: | (nm) | (fm)
vacuum | 729.1794(87) | 34
air | 728.9786(02) | 34
| $\nu_{0}$ | Uncertainty
Condition: | (THz) | (MHz)
vacuum/air | 411.1367(14) | 20
| $\lambda_{0}$ | Uncertainty
Condition: | (nm) | (fm)
vacuum (corrected for |  |
HFS and isotope shift) | 729.1795(04) | 34
air (corrected for |  |
HFS and isotope shifts) | 728.9786(19) | 34

The fitted ion temperature results in the range from 3300 K to 3500 K and can
be related to the cathode wall temperature which is about 1.5 times lower.

##  V Conclusions

The plasma of a hollow cathode has been used for accurate measurement of the
wavelength of the Kr II 4d${}^{4}D_{7/2}\rightarrow$5p${}^{4}P^{\circ}_{5/2}$
transition. By performing long-lasting measurements an improved signal to
noise is obtained as to make the residual noise contribute negligibly to the
error of the measurements. The theoretical spectrum of the transition is
calculated and used to fit the experimental data. Both the frequency and the
vacuum/air wavelength centers are given. Such measurements can be used, for
example, to perform velocity distribution measurements that would find
application in the field of the electric propulsion.

The systematic error due to the presence of hyperfine and isotope structure in
determination of the line centre is calculated and taken into account as well
as its variation with the ion temperature. In the field of electric
propulsion, the achieved accuracy in the characterization of the Kr II line at
729 nm will determine an uncertainty of the Kr ion velocity measurements of
only 14 m/s.

##  VI Acknowledgements

This work has received funding from the European Union’s Horizon 2020 research
and innovation programme under grant agreement No 101004140. The authors would
like to thank Patricia Nugent for revising the English of the manuscript.

##  VII Data availability statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * (1) V.-G. Tirila, A. Demairé, and C. N. Ryan, “Review of alternative propellants in Hall thrusters,” Acta Astronautica, vol. 212, pp. 284–306, 2023.
  * (2) M. R. Nakles, W. A. Hargus Jr., J. J. Delgadoz, and R. L. Coreyx, “A performance comparison of xenon and krypton propellant on an SPT-100 Hall thruster,” International Electric Propulsion Conference, Weisbaden, Germany, pp. IEPC–2011–003, 2011.
  * (3) O. Duchemin, D. Valentian, and N. Cornu, “Cryostorage of propellants for electric propulsion,” 45th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, 2009.
  * (4) S.Mazouffre, “Electric propulsion for satellites and spacecraft: established technologies and novel approaches,” Plasma Sources Sci. Technol., vol. 25, p. 033002, 2016.
  * (5) I. Levchenko, K. Bazaka, Y. Ding, Y. Raitses, S. Mazouffre, T. Henning, P. J. Klar, S. Shinohara, J. Schein, L. Garrigues, M. Kim, D. Lev, F. Taccogna, R. W. Boswell, C. Charles, H. Koizumi, Y. Shen, C. Scharlemann, M. Keidar, and S. Xu, “Space micropropulsion systems for Cubesats and small satellites: From proximate targets to furthermost frontiers,” Applied Physics Reviews, vol. 5, p. 011104, 02 2018.
  * (6) J. A. Linnell, “An evaluation of krypton propellant in hall thrusters,” Ph. D. Thesis, University of Michigan, 2007.
  * (7) W. A. Hargus Jr., G. M. Azarnia, and M. R. Nakles, “Kr II laser-induced fluorescence for measuring plasma acceleration,” Review of Scientific Instruments, vol. 83, p. 103111, 10 2012.
  * (8) J. Kurzyna, M. Jakubczak, A. Szelecka, and K. Dannenmayer, “Performance tests of IPPLM’s krypton Hall thruster,” Laser and Particle Beams, vol. 36, pp. 105–114, 2018.
  * (9) J. Lim, I. Levchenko, S. Huang, L. Xu, R. Sim, J. Yee, G. Potrivitu, Y. Sun, K. Bazaka, X. Wen, et al., “Plasma parameters and discharge characteristics of lab-based krypton-propelled miniaturized hall thruster,” Plasma Sources Science and Technology, vol. 28, no. 6, p. 064003, 2019.
  * (10) L.Su, R.Marks, and B.Jorns, “Investigation into the efficiency gap between Krypton and Xenon operation on a magnetically shielded Hall thruster,” International Electric Propulsion Conference 2022, June 19-23, 2022, Boston, MA, 2022.
  * (11) A. Lejeune, G. Bourgeois, and S. Mazouffre, “Kr II and Xe II axial velocity distribution functions in a cross-field ion source,” Physics of Plasmas, vol. 19, p. 073501, 07 2012.
  * (12) W. Hargus JR., G. Azarnia, and M. Nakles, “Demonstration of laser-induced fluorescence on a krypton Hall effect thruster,” IEPC - Electric Rocket Society, pp. Paper No. 2011–018, 2011.
  * (13) G. Bourgeois, A. Lejune, and S. Mazouffre, “Ion velocity evolution with channel width, magnetic topology and propellant in a 200 W Hall thruster,” IEPC - Electric Rocket Society, pp. Paper No. 2011–123, 2010\.
  * (14) K. Dzierżęga, U. Griesmann, G. Nave, and Łukasz Bratasz, “Absolute transition rates for transitions from 5p levels in kr ii,” Physica Scripta, vol. 63, p. 209, mar 2001.
  * (15) E. B. Saloman, “Energy Levels and Observed Spectral Lines of Krypton, Kr I through Kr XXXVI,” Journal of Physical and Chemical Reference Data, vol. 36, pp. 215–386, 02 2007.
  * (16) Y. Dancheva, D. Pagano, S. Scaranzin, R. Mercatelli, M. Presi, F. Scortecci, and G. Castellini, “Non-intrusive tools for electric propulsion diagnostics,” CEAS Space Journal, 2022.
  * (17) Y. Dancheva, P. Coniglio, M. Da Valle, and F. Scortecci, “Ion dynamic characterization using phase-resolved laser-induced fluorescence spectroscopy in a Hall effect thruster,” Review of Scientific Instruments, vol. 94, p. 073503, 07 2023.
  * (18) T. J. Scholl, T. D. Gaily, R. A. Holt, and S. D. Rosner, “Fast-ion-beam laser and laser-rf double-resonance measurements of hyperfine structure in 83KrII,” Phys. Rev. A, vol. 33, pp. 2396–2400, Apr 1986.
  * (19) H. A. Schuessler, A. Alousi, M. Idrees, Y. F. Li, F. Buchinger, R. M. Evans, and C. F. Fischer, “Isotope shifts and hyperfine-structure-splitting constants of the 4d-5p transition of Kr II at $\lambda$=729 nm,” Phys. Rev. A, vol. 45, pp. 6459–6467, May 1992.
  * (20) E. Arimondo, M. Inguscio, and P. Violino, “Experimental determinations of the hyperfine structure in the alkali atoms,” Rev. Mod. Phys., vol. 49, pp. 31–75, Jan 1977.
  * (21) W. Hargus, “A preliminary study of krypton laser-induced fluorescence,” 46th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, 25-28 July 2010, Nashvill TN.
  * (22) W.A.Hargus and M.A.Cappelli, “Laser-induced fluorescence measurements of velocity within a Hall discharge,” Appl. Phys. B, vol. 72, pp. 961–9, 2001.
  * (23) P. Pallavi, S. Joshi, D. Singh, M. Kaur, and H. Lee, “Comprehensive review of orthogonal regression and its applications in different domains,” Archives of Computational Methods in Engineering, vol. 29, p. 4027–4047, 2022\.