On the Adversarial Robustness of Subspace Learning

TL;DR

This paper investigates the vulnerability of subspace learning algorithms to powerful adversaries capable of fully observing and modifying data matrices, revealing optimal attack strategies based on data singular values.

Contribution

It introduces the first characterization of optimal adversarial attacks on subspace learning, considering both rank-one and general scenarios, dependent on data singular values and adversary's energy.

Findings

Optimal attack strategies depend on data singular values.

Adversaries can significantly increase subspace distance.

Numerical experiments validate attack effectiveness.

Abstract

In this paper, we study the adversarial robustness of subspace learning problems. Different from the assumptions made in existing work on robust subspace learning where data samples are contaminated by gross sparse outliers or small dense noises, we consider a more powerful adversary who can first observe the data matrix and then intentionally modify the whole data matrix. We first characterize the optimal rank-one attack strategy that maximizes the subspace distance between the subspace learned from the original data matrix and that learned from the modified data matrix. We then generalize the study to the scenario without the rank constraint and characterize the corresponding optimal attack strategy. Our analysis shows that the optimal strategies depend on the singular values of the original data matrix and the adversary's energy budget. Finally, we provide numerical experiments and practical applications to demonstrate the efficiency of the attack strategies.

Equations

$$\cos ({\theta }_i)=$$

$$\Vert {\mathbf{u}}_i\Vert =\Vert {\mathbf{v}}_i\Vert =1,$$

$${\mathbf{u}}_j^{\top }{\mathbf{u}}_i={\mathbf{v}}_j^{\top }{\mathbf{v}}_i=0,\forall j=1,2,\cdots ,i-1.$$

$$\underset{\mathbf{\Delta X}\in {\mathbb{R}}^{d\times n}}{\max }:$$

$$\hat{\mathbf{X}}=\mathbf{X}+\mathbf{\Delta X},$$

$$\Vert \mathbf{\Delta X}{\Vert }_{\text{F}}\le \eta .$$

$$\underset{\mathbf{a}\in {\mathbb{R}}^d,\mathbf{b}\in {\mathbb{R}}^n}{\max }:$$

$$\hat{\mathbf{X}}=\mathbf{X}+\mathbf{\Delta X},$$

$$\mathbf{\Delta X}=\mathbf{a}{\mathbf{b}}^{\top },$$

$$\tilde{\mathbf{U}}=\mathbf{U}+(\alpha \mathbf{Uw}+\beta \mathbf{s}){\mathbf{w}}^{\top },$$

$${\mathbf{a}}_{u^{\perp }}=(\mathbf{I}-\mathbf{U}{\mathbf{U}}^{\top })\mathbf{a},$$

$$\tilde{\mathbf{w}}=-{\mathbf{W}}^{-\top }\mathbf{b},$$

$$\omega =(1-{\mathbf{a}}^{\top }\mathbf{U}\tilde{\mathbf{w}})/\Vert {\mathbf{a}}_{u^{\perp }}\Vert ,$$

$$\alpha =|\omega |/\Vert \mathbf{g}\Vert -1,$$

$$\underset{\mathbf{a},\mathbf{b}}{\min }:$$

$$\Vert \mathbf{a}{\mathbf{b}}^{\top }{\Vert }_{\text{F}}=\Vert \mathbf{a}\Vert \Vert \mathbf{b}\Vert \le \eta ,$$

$$\frac{|\omega |}{\Vert \mathbf{g}\Vert }=\frac{|1+{\mathbf{a}}_u^{\top }{\mathbf{W}}^{-\top }\mathbf{b}|}{\Vert [\Vert {\mathbf{a}}_{u^{\perp }}\Vert {\mathbf{W}}^{-\top }\mathbf{b},1+{\mathbf{a}}_u^{\top }{\mathbf{W}}^{-\top }\mathbf{b}]\Vert },$$

$$\frac{|\omega |}{\Vert \mathbf{g}\Vert }=\frac{|1+{\mathbf{a}}_u^{\top }{\mathbf{\Sigma }}^{-1}\mathbf{b}|}{\Vert [\Vert {\mathbf{a}}_{u^{\perp }}\Vert {\mathbf{\Sigma }}^{-1}\mathbf{b},1+{\mathbf{a}}_u^{\top }{\mathbf{\Sigma }}^{-1}\mathbf{b}]\Vert }.$$

$$|{\mathbf{a}}_u^{\top }{\mathbf{\Sigma }}^{-1}\mathbf{b}|$$

$$\stackrel{(a)}{\le }\frac{1}{{\sigma }_n}\Vert \mathbf{a}\Vert \Vert \mathbf{b}\Vert =\frac{1}{{\sigma }_n}\Vert \mathbf{a}{\mathbf{b}}^{\top }{\Vert }_{\text{F}}\stackrel{(b)}{\le }\frac{\eta }{{\sigma }_n},$$

$${\mathbf{a}}_u=[0,0,\cdots ,-{\sigma }_n{]}^{\top },\mathbf{b}=[0,0,\cdots ,1{]}^{\top },$$

$$\mathbf{a}=-{\sigma }_n{\mathbf{u}}_n+\hat{a}{\mathbf{u}}_q,\mathbf{b}={\mathbf{v}}_n,$$

$$\frac{\Vert {\mathbf{a}}_{u^{\perp }}{\Vert }^2\Vert {\mathbf{\Sigma }}^{-1}\mathbf{b}{\Vert }^2}{(1+{\mathbf{a}}_u^{\top }{\mathbf{\Sigma }}^{-1}\mathbf{b}{)}^2}.$$

$${\mathbf{a}}_u=[0,0,\cdots ,-\Vert {\mathbf{a}}_u\Vert {]}^{\top },\mathbf{b}=[0,0,\cdots ,1{]}^{\top }.$$

$$\underset{c,h}{\max }:$$

$$(c^2+h^2)\le {\eta }^2,$$

$$\mathbf{a}=-{\eta }^2/{\sigma }_n{\mathbf{u}}_n\pm \eta \sqrt{1-{\eta }^2/{\sigma }_n^2}{\mathbf{u}}_q,\mathbf{b}={\mathbf{v}}_n.$$

$${\theta }^{\ast }=\{\begin{array}{ll}\pi /2,& \text{mbox}if\eta >{\sigma }_n\\ \arcsin (\eta /{\sigma }_n),& \text{mbox}if\eta \le {\sigma }_n\end{array}.$$

$$\underset{\mathbf{a}\in {\mathbb{R}}^d,\mathbf{b}\in {\mathbb{R}}^n}{\max }:$$

$$\hat{\mathbf{X}}=\mathbf{X}+\mathbf{a}{\mathbf{b}}^{\top },$$

$$\Vert \mathbf{a}\Vert \Vert \mathbf{b}\Vert \le \eta .$$

$$\underset{\mathbf{a}\in {\mathbb{R}}^{k+1},\mathbf{b}\in {\mathbb{R}}^{k+1}}{\max }:$$

$$\mathbf{Y}=\tilde{\mathbf{\Sigma }}+\mathbf{a}{\mathbf{b}}^{\top },$$

$$\Vert \mathbf{a}\Vert \Vert \mathbf{b}\Vert \le \eta ,$$

$$\mathbf{a}=[0,\cdots ,0,a_k,a_{k+1}{]}^{\top },\mathbf{b}=[0,0,\cdots ,0,1,0{]}^{\top },$$

$$\mathbf{Y}{\mathbf{Y}}^{\top }=\left[\begin{array}{cc}{\mathbf{\Lambda }}_{k-1}^2& 0\\ 0& \mathbf{c}{\mathbf{c}}^{\top }\end{array}\right],$$

$$\hat{\mathbf{U}}=\left[\begin{array}{cc}{\mathbf{I}}_{k-1}& 0\\ 0& \mathbf{z}\end{array}\right],$$

$${\left[\begin{array}{c}{\mathbf{I}}_k\\ 0\end{array}\right]}^{\top }\cdot \left[\begin{array}{cc}{\mathbf{I}}_{k-1}& 0\\ 0& \mathbf{z}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{I}}_{k-1}& 0\\ 0& z_1\end{array}\right].$$

Full text

000F. Li, L. Lai and S. Cui are with the Department of Electrical and Computer Engineering, University of California, Davis, CA. Email:{fli,lflai,sgcui}@ucdavis.edu. The work of F. Li and S. Cui was supported by the National Science Foundation with Grants CNS-1824553, DMS-1622433, AST-1547436, and ECCS-1659025. The work of L. Lai was supported by the National Science Foundation under Grants CCF-1717943, CNS-1824553 and CCF-1908258. This paper was presented in part at IEEE International Conference on Acoustics, Speech, and Signal Processing, Brighton, UK, May. 2019[1].

I Introduction

Subspace learning has a wide range of applications, such as surveillance video analysis, recommendation system, anomaly detection, etc[2, 3, 4, 5, 6, 7, 8, 9]. Among a large variety of subspace learning algorithms, principal component analysis (PCA) is one of the most widely used algorithms. In this paper, we will use PCA as the subspace learning algorithm. PCA computes a small number of principal components, which are orthogonal to each other and represent the majority of the variability of the data samples, and treats the span of these principal components as the desired low-dimensional subspace. Furthermore, many works have proposed robust PCA that can mitigate the impact of certain percentages of outliers and small dense random noise[10, 11, 12, 13].

In this paper, we investigate the adversarial robustness of subspace learning algorithms. Particularly, we examine the robustness of subspace learning algorithms against not only random noise or unintentional corrupted data as considered in existing works, but also malicious data produced by powerful adversaries who can modify the whole data set. Our study is motivated by the fact that subspace learning and many other machine learning algorithms are increasingly being used in safety critical and security related applications, such as autonomous vehicle system[14], voice recognition[15], medical image processing[16], etc. In these applications, there might exist powerful adversaries who can modify the data with the goal of maneuvering the machine learning algorithms to make the wrong decision or leave backdoor in the system[17]. To ensure the security and safety of these systems, it is important to understand the impact of these adversarial attacks on the performance of machine learning algorithms.

In our problem, given the original data matrix, we learn a low-dimensional subspace via PCA. However, there is an adversary who can observe the whole data matrix and then carefully design a modification matrix to change the original data. The goal of the adversary is to modify the original data so as to maximize the subspace distance between the subspace learned from the original data and that learned from the modified data. In this paper, we use Asimov distance[18], defined as the largest principal angle between two subspaces, to measure the subspace distance. Asimov distance has a close relationship with the chordal 2-norm distance and the Finsler distance, which are used in the analysis of optimization on Manifolds [19, 20]. Additionally, Asimov distance is closely related to the projection 2-norm distance and the gap distance, which are used in the control theory to describe the stability and robustness of a system [21, 22, 23]. As the Asimov distance depends on the modification matrix in a complex manner, to characterize the optimal attack strategy that maximizes the Asimov distance, we need to solve a complicated non-convex optimization problem.

Towards this goal, we first solve the optimization problem with an additional rank-one constraint on the modification matrix. We note that a rank-one modification is already powerful enough to capture many common modifications such as changing one data sample, inserting one adversarial data, deleting one feature, etc. Furthermore, the techniques and insights obtained from this special case are useful for the general case without the rank-one constraint. For the rank-one attack case, we study two different scenarios depending on whether the dimension of the selected subspace is equal to the rank of the data matrix or not. Our study reveals that the optimal attack strategy depends on the energy budget and the singular values of the data matrix. Specifically, in the scenario where the dimension of the selected subspace is the same as the rank of the data matrix, we show that the optimal rank-one strategy depends solely on the energy budget and the smallest singular value of the data matrix. In the scenario where the dimension of the selected subspace is less than the rank of the original data matrix, the optimal strategy depends not only on the energy budget but also on the $k$th and ($k+1$)th singular values, where $k$ is the dimension of the selected subspace.

Relying on the insights gained from the rank-one case, we then extend our study to the more general case where no rank constraint is imposed. Compared with the case with the rank-one constraint, the attacker now has more degrees of freedom to modify the data, which makes the characterization of the optimal attack strategy significantly more challenging. To solve this optimization problem, we first prove that, under the basis of the principal components of the original data matrix, the optimal attack matrix has only few non-zero entries at particular locations. This result greatly reduces the complexity of our problem. With the help of this result, we then simplify our problem to an optimization problem with the objective function being ratio of two quadratic functions. To solve this non-convex problem, we further convert our optimization problem to a feasibility problem and find the close-form solution to this problem. Our result shows that the optimal strategy depends on the energy budget and the $k$th and ($k+1$)th singular values of the data matrix. Our analysis shows that, compared with the optimal rank-one strategy, this strategy leads to a larger subspace distance.

Our study is related to the recent works on adversary machine learning. For example, [24] studies how to change the data to manipulate the result of the regression learning system. [25] investigates the optimal modification strategy to maximize the inference errors in a multivariate estimation system. In an interesting related work [26], the authors study how to design an adversarial data sample and add it to the data matrix in order to maximize the Asimov distance between the subspace estimated by PCA from the contaminated data matrix and that from the original data matrix. [26] focuses on the case where the original data matrix is low-rank and the dimension of the selected subspace is equal to the rank of the data matrix. By contrast, we consider a more powerful adversarial setting, where the data matrix is not constrained to being low-rank, the dimension of the selected subspace does not necessarily equal the rank of the data matrix, and the adversary can modify the whole data matrix instead of only adding one data sample.

The remainder of the paper is organized as follows. In Section II, we describe the precise problem formulation. In Section III, we investigate the optimal rank-one attack strategy. We generalize our results to the case without the rank constraint in Section IV. In Section V, we provide numerical experiments with both synthesized data and real data to illustrate results obtained in this paper. Finally, we offer concluding remarks in Section VI.

II Problem formulation

In this section, we introduce the problem formulation. Given a data matrix $\mathbf{X}=[\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{n}]$ with each $\mathbf{x}_{i}\in\mathbb{R}^{d}$, our goal is to learn a low-dimension subspace via PCA. In the data matrix $\mathbf{X}$, we assume that all the preprocessing steps (such as data centering and standardization) have been done. In this paper, we consider an adversarial setup in which an adversary will first observe $\mathbf{X}$ and then carefully design a modification (attack) matrix $\mathbf{\Delta X}$ to change $\mathbf{X}$ to $\hat{\mathbf{X}}=\mathbf{X}+\mathbf{\Delta X}$. We denote function $g_{k}(\cdot)$ as the PCA operation that computes the $k$ leading principal components. Furthermore, let $\mathbb{X}=\textrm{span}(g_{k}(\mathbf{X}))$ be a $k$-dimensional subspace learned from $\mathbf{X}$ and $\hat{\mathbb{X}}=\textrm{span}(g_{k}(\hat{\mathbf{X}}))$ be a $k$-dimensional subspace learned from the modified data $\hat{\mathbf{X}}$. The goal of the adversary is to design the modification matrix $\mathbf{\Delta X}$ so as to make the distance between $\mathbb{X}$ and $\hat{\mathbb{X}}$ as large as possible. To measure such a distance, we use the largest principal angle between $\mathbb{X}$ and $\hat{\mathbb{X}}$ as defined below [18]. The largest principal angle plays an important role in the subspace classification problem[27]. It is closely related to the projection 2-norm which is widely used in engineering applications[28, 18, 29]. The projection 2-norm also provides a way to measure the discrepancy of the projections of a vector on two distinct subspaces. It is useful in the robustness analysis of the principal component regression (PCR), as one is actually projecting the response value vector onto the selected feature subspace in PCR. We will provide an example to illustrate it in Section V using real data.

Definition 1**.**

Let $\mathbb{X}$ and $\hat{\mathbb{X}}$ be two $k$-dimensional subspaces in $\mathbb{R}^{d}$. The principal angles $\{\theta_{i}\}_{i=1}^{k}$ are defined recursively:

[TABLE]

In this paper, we will use $\|\cdot\|$ to denote the $\ell_{2}$ norm and \theta\big{(}g_{k}(\mathbf{X}),g_{k}(\hat{\mathbf{X}})\big{)} or simply $\theta$ to denote the Asimov distance between the subspace $\mathbb{X}$ estimated from $\mathbf{X}$ and the subspace $\hat{\mathbb{X}}$ estimated from $\hat{\mathbf{X}}$. Given an orthonormal basis $\mathbf{U}_{\mathbb{X}}$ of $\mathbb{X}$ and an orthonormal basis $\mathbf{U}_{\hat{\mathbb{X}}}$ of $\hat{\mathbb{X}}$, $\{\cos(\theta_{1}),\cdots,\cos(\theta_{k})\}$ are the singular values of $\mathbf{U}_{\mathbb{X}}^{\top}\mathbf{U}_{\hat{\mathbb{X}}}$ [18]. Hence, the Asimov distance is determined by the smallest singular value of $\mathbf{U}_{\mathbb{X}}^{\top}\mathbf{U}_{\hat{\mathbb{X}}}$. It is easy to see that, if no constraint is imposed on $\mathbf{\Delta X}$, $\hat{\mathbf{X}}$ can be arbitrary and $\theta$ can be easily made to be $\pi/2$. Therefore, we impose an energy constraint on $\mathbf{\Delta X}$. In particular, we assume that the energy of $\mathbf{\Delta X}$ is less than or equal to $\eta$. In this paper, we use the Frobenius norm $\|\mathbf{\Delta X}\|_{\text{F}}$ to measure the energy. Hence, the goal of this attacker is to solve the following optimization problem:

[TABLE]

Even though (1) is a complicated non-convex optimization problem, we will fully characterize the optimal solution to (1) for any given $\eta$. This characterization will enable us to investigate the impact of this optimal attack with respect to the energy budget $\eta$.

III Optimal rank-one adversarial strategy

In this section, we will solve (1) for the special case where the modification matrix $\mathbf{\Delta X}$ is limited to being rank-one. The techniques and insights obtained from this special case will be useful for the general case considered in Section IV.

With this additional rank-one constraint, $\mathbf{\Delta X}$ can be written as $\mathbf{a}\mathbf{b}^{\top}$ for some $\mathbf{a}\in\mathbb{R}^{d}$ and $\mathbf{b}\in\mathbb{R}^{n}$, and the optimization problem (1) becomes

[TABLE]

It is easy to see that, for any feasible solution $(\tilde{\mathbf{a}},\tilde{\mathbf{b}})$ with $||\tilde{\mathbf{b}}||\neq 1$, we can construct another feasible solution $(||\tilde{\mathbf{b}}||\tilde{\mathbf{a}},\tilde{\mathbf{b}}/||\tilde{\mathbf{b}}||)$ that gives the same objective function value. Hence, without loss of optimality, we will fix the norm of $\mathbf{b}$ to be $1$ throughout this section.

Based on the value of $k$, i.e., the dimension of the subspace we select, we will first present the solution to the case when $k=\text{rank}(\mathbf{X})$, and then generalize the result to the case when $k<\text{rank}(\mathbf{X})$.

III-A Case with $k=\text{rank}(\mathbf{X})$

In this subsection, we consider the case when the dimension of the subspace selected is equal to the rank of the data matrix. In this case, the span of $\mathbf{X}$ equals the span of $g_{k}(\mathbf{X})$. Furthermore, we divide this case into two scenarios where the data matrix is full-rank and the data matrix is low-rank.

III-A1 Full-Rank Case

In the full column rank case, $\text{rank}(\mathbf{X})=n$, where $n\leq d$. This case arises when the number of samples is limited, for example, at the beginning of online PCA. In this case, the span of $\hat{\mathbf{X}}$ is equal to the span of $g_{k}(\hat{\mathbf{X}})$, and hence we can write \theta\big{(}g_{k}(\mathbf{X}),g_{k}(\hat{\mathbf{X}})\big{)} as $\theta(\mathbf{X},\hat{\mathbf{X}})$. In the following, we first find the expression of $\theta(\mathbf{X},\hat{\mathbf{X}})$ for any given $\hat{\mathbf{X}}=\mathbf{X}+\mathbf{a}\mathbf{b}^{T}$. Using this expression, we then characterize the optimal attack matrix $\mathbf{\Delta}\mathbf{X}$.

Suppose the compact SVD of $\mathbf{X}$ is $\mathbf{X}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}=\mathbf{U}\mathbf{W},$ where $\mathbf{\Sigma}=\text{diag}(\sigma_{1},\sigma_{2},\cdots,\sigma_{n})$. One set of orthonormal bases for the column space of $\mathbf{X}$ is $\mathbf{U}$. We can also use SVD to find a set of orthonormal bases $\tilde{\mathbf{U}}$ of $\text{span}(\hat{\mathbf{X}})$.

Since $\hat{\mathbf{X}}=\mathbf{X}+\mathbf{a}\mathbf{b}^{\top}$, $\tilde{\mathbf{U}}$ can be directly expressed as a function of $\mathbf{U}$ [30]:

[TABLE]

where

[TABLE]

and $\mathbf{W}^{-\top}=\mathbf{(}\mathbf{W}^{-1})^{\top}$. Hence, we have $\mathbf{U}^{\top}\tilde{\mathbf{U}}=\mathbf{U}^{\top}\left(\mathbf{U}+(\alpha\mathbf{U}\mathbf{w}+\beta\mathbf{s})\mathbf{w}^{\top}\right)=\mathbf{I}+\alpha\mathbf{w}\mathbf{w}^{\top}.$ The singular values of $\mathbf{I}+\alpha\mathbf{w}\mathbf{w}^{\top}$ are $\{1,1,\cdots,1+\alpha\mathbf{w}^{\top}\mathbf{w}\}$. Since $\mathbf{w}^{\top}\mathbf{w}=1,1+\alpha=|\omega|/\|\mathbf{g}\|$, the smallest singular value of $\mathbf{U}^{\top}\tilde{\mathbf{U}}$ is $\cos(\theta)=|\omega|/\|\mathbf{g}\|$. Our objective is to maximize $\theta$, which is equivalent to minimizing the smallest singular value of $\mathbf{U}^{\top}\tilde{\mathbf{U}}$. Hence, the optimization problem (2) is simplified as

[TABLE]

where we use the identity $\|\mathbf{a}\|\|\mathbf{b}\|=\|\mathbf{a}\cdot\mathbf{b}^{\top}\|_{\text{F}}$. Expanding the objective function, we have

[TABLE]

where $\mathbf{a}_{u}=\mathbf{U}^{\top}\mathbf{a}$.

Since $\mathbf{W}=\mathbf{\Sigma}\mathbf{V}^{\top}$, we have $\mathbf{W}^{-\top}\mathbf{b}=\mathbf{\Sigma}^{-1}\mathbf{V}^{\top}\mathbf{b}$. As $\mathbf{V}$ is a unitary matrix, changing the coordinate $\mathbf{b}\Leftarrow\mathbf{V}^{\top}\mathbf{b}$ does not change the constraint. The value $\mathbf{a}_{u}^{\top}\mathbf{W}^{-\top}\mathbf{b}$ in the original coordinate is the same as $\mathbf{a}_{u}^{\top}\mathbf{\Sigma}^{-1}\mathbf{b}$ in the new coordinate. In the following, we will use this new coordinate system and the cost function in (3) can be written as

[TABLE]

The objective function (4) is zero if and only if the numerator is zero. Using the matrix norm inequality[31], we have

[TABLE]

where $\|\mathbf{\Sigma}^{-1}\|_{2}$ is the induced 2-norm of matrix $\mathbf{\Sigma}^{-1}$, in (a) we use $\|\mathbf{a}_{u}\|\leq\|\mathbf{a}\|$, and (b) is due to the energy constraint. From the inequalities, we conclude that when $\eta<\sigma_{n}$, we can not make the numerator to be zero. We now consider two different cases depending on whether we can make the numerator to be zero or not.

Case 1: When $\eta>\sigma_{n}$, if we set

[TABLE]

and any $\|\mathbf{a}_{u^{\perp}}\|^{2}=\hat{a}^{2}$ with $0<\hat{a}^{2}<\eta^{2}-\sigma_{n}^{2}$, the numerator will be zero. Since $\mathbf{a}=\mathbf{U}\mathbf{a}_{u}+(\mathbf{I}-\mathbf{U}\mathbf{U}^{\top})\mathbf{a}_{u^{\perp}}$, the attacker can make the Asimov distance to be $\pi/2$ by setting:

[TABLE]

where $\mathbf{u}_{q}$ is any vector orthogonal to the column space of $\mathbf{X}$ and $0<\hat{a}^{2}<\eta^{2}-\sigma_{n}^{2}$.

Case 2: When $\eta\leq\sigma_{n}$, the value of $1+\mathbf{a}_{u}^{\top}\mathbf{\Sigma}^{-1}\mathbf{b}$ can not reach zero. In this case, it is easy to check that minimizing (4) is equivalent to maximizing

[TABLE]

As $\|\mathbf{b}\|=1$, $\|\mathbf{\Sigma}^{-1}\mathbf{b}\|^{2}$ is maximized when $\mathbf{b}=[0,0,\cdots,1]^{\top}$. Furthermore, for any fixed norm of $\mathbf{a}_{u}$, $(1+\mathbf{a}_{u}^{\top}\mathbf{\Sigma}^{-1}\mathbf{b})^{2}$ is minimized when $\mathbf{a}_{u}=[0,0,\cdots,-\|\mathbf{a}_{u}\|]^{\top}$, $\mathbf{b}=[0,0,\cdots,1]^{\top}$. Hence, for fixed norms of $\mathbf{a}_{u}$, $\mathbf{a}_{u^{\perp}}$, the objective function (6) is maximized when

[TABLE]

Let $c=\|\mathbf{a}_{u^{\perp}}\|,h=\|\mathbf{a}_{u}\|$. Using the optimal form of $\mathbf{a}_{u}$ and $\mathbf{b}$ in (7), the objective function (6) can be simplified to

[TABLE]

It is easy to check that the objective function is maximized when $c^{2}+h^{2}=\eta^{2}$. Hence, we have $c^{2}=\eta^{2}-h^{2}$. Inserting this value of $c$ into the objective function and setting the derivative with respect to $h$ to be [math], we get a unique solution $h=\eta^{2}/\sigma_{n}$. At this value of $h$, the second derivative is $\frac{-2\sigma_{n}^{2}}{(\sigma_{n}^{2}-\eta^{2})^{3}}$, which is negative. It indicates that $h=\eta^{2}/\sigma_{n}$ is indeed the maximum point. Hence $c=\pm\eta\sqrt{1-\eta^{2}/\sigma_{n}^{2}}$. This implies that the optimal solution to problem (2) for Case 2 is

[TABLE]

Summarizing the discussion above, we have the following proposition regarding the optimal value of problem (2) in the full-rank case.

Proposition 1**.**

In the full rank case, the optimal value of (2) is

[TABLE]

III-A2 Low-Rank Case

We now consider the case where $\mathbf{X}$ is not full rank. Let $k<\text{min}(d,n)$ be the rank of $\mathbf{X}$. In this subsection, with a slight abuse of notation, we write the full SVD of $\mathbf{X}$ as $\mathbf{X}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$. The optimal attack matrix could be found by solving

[TABLE]

We can further simplify this optimization problem as

[TABLE]

where $\tilde{\mathbf{\Sigma}}=\text{diag}(\sigma_{1},\sigma_{2},\cdots,\sigma_{k},0)$ and $\{\sigma_{1},\sigma_{2},\cdots,\sigma_{k}\}$ are singular values of $\mathbf{X}$. Detailed proof of the equivalence between (9) and (10) can be found in Appendix A. Here, we describe the main idea of the proof. The main step of the simplification is to left multiply the unitary matrix $\mathbf{U}^{\top}$ and right multiply the unitary matrix $\mathbf{V}$ on both $\mathbf{X}$ and $\hat{\mathbf{X}}$. Note that multiplying a unitary matrix does not change the column space and its singular values. In addition, a rank-one modification can only add at most one principal component orthogonal to its original column subspace. Hence, by changing the coordinates, $\mathbf{a}$ and $\mathbf{b}$ are $k+1$ dimensional vectors.

To solve problem (10), we divide it into two cases based on the value of the energy budget.

Case 1: When $\eta>\sigma_{k}$, it is simple to verify that the solution

$\mathbf{a}=[0,0,\cdots,\eta]^{\top}$, $\mathbf{b}=[0,0,\cdots,1]^{\top}$ leads to the maximal Asimov distance, which is $\pi/2$.

Case 2: When $\eta\leq\sigma_{k}$, the following theorem characterizes the form of optimal $\mathbf{a}$ and $\mathbf{b}$.

Theorem 1**.**

There exists an optimal solution to problem (10) in the following form

[TABLE]

with $a_{k}^{2}+a_{k+1}^{2}=\eta^{2}$.

Proof.

Please see Appendix B. ∎

In the following, we will find the optimal values of $a_{k}$ and $a_{k+1}$. Since $\|\mathbf{a}\|^{2}=\eta^{2}$ and $\mathbf{a}$ is in the form of (11), we can write $\mathbf{a}=\eta[0,0,\cdots,\cos(\alpha),\sin(\alpha)]^{\top}$, where $\alpha\in[0,2\pi)$. To compute the $k$ leading principal components of $\mathbf{Y}$, we can perform the eigenvalue decomposition of $\mathbf{Y}\mathbf{Y}^{\top}$,

[TABLE]

where $\mathbf{c}=[\sigma_{k}+\eta\cos{\alpha},\eta\sin(\alpha)]^{\top}$, $\mathbf{\Lambda}_{k-1}=\text{diag}(\sigma_{1},\sigma_{2},\cdots,\sigma_{k-1})$. Suppose the compact SVD of $\mathbf{Y}\mathbf{Y}^{\top}$ is $\mathbf{Y}\mathbf{Y}^{\top}=\hat{\mathbf{U}}\hat{\mathbf{\Sigma}}\hat{\mathbf{V}}^{\top},$ where

[TABLE]

and $\mathbf{z}\in\mathbb{R}^{2}$ is the eigenvector of $\mathbf{c}\mathbf{c}^{\top}$ corresponding to its nonzero eigenvalue. Since one orthonormal basis of $\text{span}(\tilde{\mathbf{\Sigma}})$ is $[\mathbf{I}_{k},\mathbf{0}]^{\top}$, the Asimov distance is determined by the singular values of

[TABLE]

Hence, the Asimov distance is $\arccos(|z_{1}|)$. Since $\mathbf{c}$ is the eigenvector of $\mathbf{c}\mathbf{c}^{\top}$ corresponding to its nonzero eigenvalue, we have $|z_{1}|=\frac{|c_{1}|}{\|\mathbf{c}\|}$. Our objective function is reduced to

[TABLE]

It is simple to show that the optimal solution to (12) is

[TABLE]

or

[TABLE]

Substitute the optimal solution of $\alpha^{*}$ in (13) or (14) into the objective of problem (12), we have $\sin(\theta^{*})=\eta/\sigma_{k}.$ Hence, the optimal solution to problem (10) is

[TABLE]

which indicates that the optimal solution to problem (9) is

[TABLE]

where $\mathbf{u}_{q}$ is any vector orthogonal to the column space of $\mathbf{X}$. The corresponding optimal subspace distance is $\theta^{*}=\arcsin(\eta/\sigma_{k}).$ In summary, we have

Proposition 2**.**

The optimal Asimov distance in the low-rank case is

[TABLE]

The result is similar to the full column rank case characterized in Proposition 1.

III-B Case with $k<\text{rank}(\mathbf{X})$

In this section, we consider the more practical but much more challenging case with $k<\text{rank}(\mathbf{X})$.

Given the data matrix $\mathbf{X}\in\mathbb{R}^{d\times n}$, without loss of generality, we assume $d\leq n$ and $\text{rank}(\mathbf{X})=d$. Assume the full SVD of $\mathbf{X}$ is $\mathbf{X}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}$, where $\mathbf{U}\in\mathbb{R}^{d\times d}$, $\mathbf{\Sigma}\in\mathbb{R}^{d\times n}$, $\mathbf{V}\in\mathbb{R}^{n\times n}$, and the singular values of $\mathbf{X}$ are $\{\sigma_{1},\sigma_{2},\cdots,\sigma_{k},\cdots,\sigma_{d}\}$. Recall that we denote $g_{k}(\cdot)$ as the PCA operation that computes the $k$ leading principal components. In this scenario, as the original data matrix is not low-rank, we will perform PCA both on the original data matrix and on the modified data matrix. Hence, the optimal rank-one modification matrix can be found by solving the following optimization problem

[TABLE]

By diagonalizing the data matrix and using similar arguments in Appendix A, (16) can be further simplified as

[TABLE]

where $g_{k}(\mathbf{\Sigma})=[\mathbf{I}_{k},\mathbf{0}]^{\top}\in\mathbb{R}^{d\times k}$. Here we also perform variable change $\mathbf{a}\Leftarrow\mathbf{U}^{\top}\mathbf{a}$ and $\mathbf{b}\Leftarrow\mathbf{V}^{\top}\mathbf{b}$. To solve this optimization problem, we divide it into two cases depending on the energy budget and the difference between $\sigma_{k}$ and $\sigma_{k+1}$.

Case 1: When $\eta\geq\sigma_{k}-\sigma_{k+1}$, we have one simple solution $\mathbf{a}=[0,0,\cdots,0,\eta,0,\cdots,0]^{\top}$, where $\eta$ is in the ($k+1$)th coordinate, and $\mathbf{b}=[0,0,\cdots,0,1,0,\cdots,0]^{\top}$, where element $1$ is in the ($k+1$)th coordinate. Clearly, this setting of $\mathbf{a}$ and $\mathbf{b}$ leads to the maximal subspace distance, which is $\pi/2$.

Case 2: When $\eta<\sigma_{k}-\sigma_{k+1}$, the following theorem gives the form of the optimal solution.

Theorem 2**.**

The optimal solution to problem (17) should be in the form of

[TABLE]

where $a_{k}^{2}+a_{k+1}^{2}=\eta^{2}$ and $b_{k}^{2}+b_{k+1}^{2}=1$.

Proof.

Please see Appendix C for details. ∎

As the optimal solution of $\mathbf{a}$ and $\mathbf{b}$ are in the form of (18) and (19), we can parametrize $\mathbf{a}$ and $\mathbf{b}$ with parameters $\alpha$ and $\beta$ using $\mathbf{a}=\eta[0,0,\cdots,\cos(\alpha),\sin(\alpha),0,\cdots,0]^{\top}$ and $\mathbf{b}=[0,0,\cdots,\cos(\beta),\sin(\beta),0,\cdots,0]^{\top}$ respectively.

As a result, the modified data matrix $\mathbf{Y}$ can be written as

[TABLE]

where $\mathbf{\Sigma}_{1}=\text{diag}(\sigma_{1},\sigma_{2},\cdots,\sigma_{k-1})$, $\mathbf{\Sigma}_{3}=\text{diag}(\sigma_{k+2},\cdots,\sigma_{d})$, and

[TABLE]

Since $\mathbf{Y}$ has the pseudo block diagonal form, the singular values and principal components of $\mathbf{Y}$ are determined by the SVD of $\mathbf{\Sigma}_{1}$, $\mathbf{\Sigma_{2}}$, and $\mathbf{\Sigma}_{3}$. For notation convenience, we denote $\mathbf{\Sigma_{2}}=\mathbf{D}+\eta\bar{\mathbf{a}}\bar{\mathbf{b}}^{\top}$, where $\mathbf{D}=\text{diag}(\sigma_{k},\sigma_{k+1})$, $\bar{\mathbf{a}}=[\cos{\alpha},\sin{\alpha}]^{\top}$, and $\bar{\mathbf{b}}=[\cos{\beta},\sin{\beta}]^{\top}$. Let $\xi_{1}$ and $\xi_{2}$ be the two singular values of $\mathbf{\Sigma}_{2}$ and denote their corresponding left singular vectors as

[TABLE]

The following lemma characterizes the form of $k$-dimensional subspace learned by PCA from $\mathbf{Y}$.

Lemma 1**.**

[TABLE]

Proof.

According to the perturbation theory [32], the singular values of $\mathbf{\Sigma}_{2}$ must satisfy

[TABLE]

It indicates that $\xi_{1}>\sigma_{k}$, $\xi_{2}>\sigma_{k}$ and $\xi_{1}<\sigma_{k+1}$, $\xi_{2}<\sigma_{k+1}$ will not happen. Hence, we will select the eigenvector corresponding to singular value $\xi_{1}$ as one of the leading $k$ principal components, which completes the proof. ∎

Since one set of orthonormal bases for $g_{k}(\mathbf{\Sigma})$ is $[\mathbf{I}_{k},\mathbf{0}]^{\top}$, then the subspace distance \theta\big{(}g_{k}(\mathbf{\Sigma}),g_{k}(\mathbf{Y})\big{)} is determined by the singular values of

[TABLE]

Hence, the subspace distance is $\arccos(|\cos\varphi|)$ and our optimization problem can be equivalently formulated as

[TABLE]

Let $\mathbf{Z}=\mathbf{\Sigma}_{2}\mathbf{\Sigma}_{2}^{\top}$, we can compute $\mathbf{W}$ through eigenvalue decomposition of $\mathbf{Z}$. According to the equality $\mathbf{\Sigma}_{2}\mathbf{\Sigma}_{2}^{\top}=\mathbf{W}\cdot\text{diag}(\xi_{1}^{2},\xi_{2}^{2})\cdot\mathbf{W}^{\top}$, we have

[TABLE]

From this equation, we obtain

[TABLE]

Then we can compute $\varphi$ through

[TABLE]

where $\text{atan2}(\cdot,\cdot)$ is the four-quadrant inverse tangent function, $a_{x}=Z_{1,1}-Z_{2,2}$, and $a_{y}=Z_{1,2}+Z_{2,1}$. In our case, the specific expressions of $a_{x}$ and $a_{y}$ are

[TABLE]

Let us write $a_{x}$ and $a_{y}$ as a function of $\alpha$ and $\beta$: $a_{x}=a_{x}(\alpha,\beta)$ and $a_{y}=a_{y}(\alpha,\beta)$. To further restrict the domains of $\alpha$ and $\beta$, we analyze the properties of the angle $\varphi$ in (23) as a function of $\alpha$ and $\beta$. First, we have $a_{x}(\alpha,\beta)=a_{x}(\pi+\alpha,\pi+\beta)$ and $a_{y}(\alpha,\beta)=a_{y}(\pi+\alpha,\pi+\beta)$. So $\varphi(\alpha,\beta)=\varphi(\pi+\alpha,\pi+\beta)$. This property indicates that we only need to consider the function value in the domain $\alpha\in[0,\pi],\beta\in[-\pi,\pi]$. Second, $a_{x}(\alpha,\beta)=a_{x}(\pi-\alpha,\pi-\beta)$ and $a_{y}(\alpha,\beta)=-a_{y}(\pi-\alpha,\pi-\beta)$, and then we have $\varphi(\alpha,\beta)=-\varphi(\pi-\alpha,\pi-\beta)$. Since $\cos(\varphi)$ is an even function, we only need to consider the function with domain $\alpha\in[0,\pi/2],\beta\in[-\pi,\pi]$. Note that $\mathbf{\Sigma}_{2}$ is in the form of (20), the variance in the direction of $\mathbf{e}_{k}$ is $v_{k}=\cos(\alpha)^{2}+\sigma_{k}^{2}+2\cos(\alpha)\cos(\beta)$, and the variance in the direction of $\mathbf{e}_{k+1}$ is $v_{k+1}=\sin(\alpha)^{2}+\sigma_{k+1}^{2}+2\sin(\alpha)\sin(\beta)$. To maximize the subspace distance, we should make $v_{k}$ small and make $v_{k+1}$ large. Apparently, the sign of $\cos(\alpha)\cos(\beta)$ should be negative and the sign of $\sin(\alpha)\sin(\beta)$ should be positive. Hence the optimal $\alpha$ and $\beta$ should satisfy $\alpha\in[0,\pi/2]$ and $\beta\in[\pi/2,\pi]$. As a result, the optimization problem (22) can be written as

[TABLE]

The following theorem characterizes the optimal solution to problem (25).

Theorem 3**.**

The optimal solution to problem (25) is

[TABLE]

where $H=\sigma_{k}^{4}+\sigma_{k+1}^{4}+\eta^{4}-2\sigma_{k}^{2}\sigma_{k+1}^{2}-2\sigma_{k}^{2}\eta^{2}-2\sigma_{k+1}^{2}\eta^{2}$.

Proof.

Please see Appendix D. ∎

Accordingly, the optimal solution to problem (16) is

[TABLE]

Furthermore, the optimal subspace distance $\theta^{*}$ can be computed according to (24) and (23). Moreover, according to the properties of the function $\varphi(\alpha,\beta)$ we have discussed before, there are other three optimal solutions

[TABLE]

which lead to the same optimal objective value.

IV Optimal adversarial strategy without the rank constraint

Using the insights gained from Section III, we now characterize the optimal attack strategy in the general case without the rank-one constraint by solving (1). We will directly consider the general case with $k\leq\text{rank}(\mathbf{X})$.

Following the similar transformation from (9) to (10), we can simplify the optimization problem (1) as

[TABLE]

where without loss of generality we assume $d\leq n$, the full SVD of the data matrix is $\mathbf{X}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^{\top}$, the singular values of the data matrix are $\{\sigma_{1},\sigma_{2},\cdots,\sigma_{d}\}$, and $\mathbf{B}=\mathbf{U}^{\top}\mathbf{\Delta X}\mathbf{V}$. To identify the optimal modification matrix $\mathbf{B}$ in problem (29), we divide it into two cases.

Case 1: When $\eta\geq\frac{\sigma_{k}-\sigma_{k+1}}{\sqrt{2}}$, by setting $b_{k,k}=-\eta/\sqrt{2}$, $b_{k+1,k+1}=\eta/\sqrt{2}$, and all other entries of $\mathbf{B}$ to zero, where $b_{i,j}$ is the element in the $i$th row and $j$th column of $\mathbf{B}$, this will lead to the maximal subspace distance, $\pi/2$.

Case 2: When $\eta<\frac{\sigma_{k}-\sigma_{k+1}}{\sqrt{2}}$, the following theorem states the form of the optimal $\mathbf{B}$.

Theorem 4**.**

The optimal $\mathbf{B}$ to problem (29) has only four possible non-zero entries: $b_{k,k},b_{k,k+1},b_{k+1,k}$ and $b_{k+1,k+1}$.

Proof.

Please see Appendix E. ∎

This characterization reduces the complexity of problem (29). Using this optimal form of $\mathbf{B}$ and following similar steps leading to (23), we can write the subspace distance as

[TABLE]

where

[TABLE]

It is easy to see that we can change the sign of $b_{y}$ by changing the signs of $b_{k,k+1}$ and $b_{k+1,k}$. We also have $b_{x}>0$, as

[TABLE]

Using these two facts and the fact that $\text{atan2}(b_{y},b_{x})$ is an odd function of $b_{y}$ when $b_{x}>0$, we know that maximizing $\theta$ in (30) is equivalent to maximizing $b_{y}/b_{x}$. Hence, our optimization problem can be written as

[TABLE]

where $\mathbf{u}\triangleq\bar{\mathbf{b}}+\boldsymbol{\sigma}$ with $\bar{\mathbf{b}}=[b_{k,k},b_{k+1,k},b_{k,k+1},b_{k+1,k+1}]^{\top}$ and $\boldsymbol{\sigma}=[\sigma_{k},0,0,\sigma_{k+1}]^{\top}$,

[TABLE]

The objective function is the ratio of two quadratic functions. It is a non-convex problem in general. In the following, we transform this problem into a feasibility problem and obtain the closed-form solution analytically.

Let $\lambda$ denote the value of the objective function in (31). We can rewrite the optimization problem (31) as

[TABLE]

The first constraint can be written as $\mathbf{u}^{\top}(\mathbf{A}_{1}-\lambda\mathbf{A}_{2})\mathbf{u}=0$, where

[TABLE]

To further simplify the constraint, we perform eigenvalue decomposition on $\mathbf{Q}=\mathbf{P}\mathbf{\Lambda}\mathbf{P}^{\top},$ where $\mathbf{\Lambda}=\text{diag}(\sqrt{\lambda^{2}+1},-\sqrt{\lambda^{2}+1})$ and

[TABLE]

with $t=1/\sqrt{(\sqrt{\lambda^{2}+1}+\lambda)^{2}+1}$.

We further perform variable change $\mathbf{v}\triangleq\text{diag}(\mathbf{P}^{\top},\mathbf{P}^{\top})\mathbf{u}$. Thus, the constraint (32) is equivalent to $\mathbf{v}^{\top}\mathbf{\Lambda}\mathbf{v}=0$, which indicates $v_{1}^{2}+v_{3}^{2}=v_{2}^{2}+v_{4}^{2}$. With this, the optimization problem is simplified as

[TABLE]

where $\bar{\boldsymbol{\sigma}}=\text{diag}(\mathbf{P}^{\top},\mathbf{P}^{\top})\boldsymbol{\sigma}=[p_{1,1}\sigma_{k},\,p_{1,2}\sigma_{k},\,p_{2,1}\sigma_{k+1},\,p_{2,2}\sigma_{k+1}]^{\top}.$ Note that $p_{1,2}=-p_{2,1}$ and $p_{2,2}=p_{1,1}$, we have $\bar{\boldsymbol{\sigma}}=[p_{1,1}\sigma_{k},\,-p_{2,1}\sigma_{k},\,p_{2,1}\sigma_{k+1},\,p_{1,1}\sigma_{k+1}]^{\top}.$

Now, problem (34) can be solved by checking the feasibility of (35) and (36) given a particular $\lambda$. Given $\lambda$, the feasibility of problem(34) is equivalent to the feasibility of

[TABLE]

Note that $\bar{\boldsymbol{\sigma}}$ depends on $\lambda$, we denote the left hand side of inequality (37) as $f(\mathbf{v},\lambda)=\|\mathbf{v}-\bar{\boldsymbol{\sigma}}\|^{2}$ and parametrize $\mathbf{v}$ as

[TABLE]

It is easy to verify that the minimum point of $f(\mathbf{v},\lambda)$ in terms of $\mathbf{v}$ is obtained at the following stationary point

[TABLE]

Plug the optimal $r$, $\alpha$, $\beta$ of (39) into $f(\mathbf{v},\lambda)$, and we have

[TABLE]

According to inequality (37), inequality $f(\lambda)\leq\eta^{2}$ now is equivalent to

[TABLE]

Denote the right hand of the above inequality as $c\triangleq(\sigma_{k}^{2}+\sigma_{k+1}^{2})/2-\eta^{2}$. Since $\eta<(\sigma_{k}-\sigma_{k+1})/\sqrt{2}$, we have $c>\sigma_{k}\sigma_{k+1}$. Furthermore, we notice that $p_{1,1}^{2}=1-p_{2,1}^{2}$. Plug it into inequality (40), and we have

[TABLE]

Let

[TABLE]

and since $\sigma_{k}\sigma_{k+1}<c\leq(\sigma_{k}^{2}+\sigma_{k+1}^{2})/2$, we have $0<w\leq\frac{(\sigma_{k}^{2}+\sigma_{k+1}^{2})^{2}/4-\sigma_{k}^{2}\sigma_{k+1}^{2}}{(\sigma_{k}^{2}-\sigma_{k+1}^{2})^{2}}=1/4$. Denote the left hand of inequality (41) as $h(p_{2,1})$, and we have

[TABLE]

Moreover, since $1/\sqrt{2}<p_{2,1}<1$, we must have

[TABLE]

where $p_{2,1}^{H}=\sqrt{(1+\sqrt{1-4w})/2}$ is the largest root of $h(p_{2,1})=0$. Pluging the expressions of $p_{2,1}$ and $p_{2,1}^{H}$ into (43), we can get

[TABLE]

Simplifying this inequality leads to $\lambda\leq\frac{e^{2}-1}{2e}$, where

[TABLE]

Thus we can conclude that

[TABLE]

Accordingly, the optimal subspace distance in (1) is

[TABLE]

In summary, given energy budget $\eta$, we first compute $w$ according to (42) and compute $e$ according to (44), from which we can get $\lambda_{\max}$ and $\theta^{*}$ using (45) and (46). Having obtained the optimal $\lambda_{\max}$, we can compute $\mathbf{P}$ in (33) and compute $\mathbf{v}$ using (39) and (38), and sequentially compute $\mathbf{u}$ and $\bar{\mathbf{b}}$. Finally, if the optimal solution of problem (29) is $\mathbf{B}^{*}$ with non-zero entries $\bar{\mathbf{b}}^{*}=[b_{k,k}^{*},b_{k+1,k}^{*},b_{k,k+1}^{*},b_{k+1,k+1}^{*}]^{\top}$, we also have another paired feasible optimal solution with non-zero entries being $[b_{k,k}^{*},-b_{k+1,k}^{*},-b_{k,k+1}^{*},b_{k+1,k+1}^{*}]^{\top}$, which leads to the same optimal value. Accordingly, the optimal solution to problem (1) is $\mathbf{\Delta X}^{*}=\mathbf{U}\mathbf{B}^{*}\mathbf{V}^{\top}$.

V Numerical experiments and applications

In this section, we provide numerical examples to illustrate the results obtained in this paper. We will also apply the results to principal component regression[33] to illustrate potential applications in practice.

V-A Numerical experiments

In this subsection, we illustrate the results with synthesized data.

In the first experiment, we employ different attack strategies in a low-rank data matrix. In this simulation, we set $d=5$, $n=5$, and $k=3$. We generate the original data matrix as $\mathbf{X}=\mathbf{A}\mathbf{B}^{\top}$, where $\mathbf{A}\in\mathbb{R}^{d\times k}$, $\mathbf{B}\in\mathbb{R}^{n\times k}$, and each entry of $\mathbf{A}$ and $\mathbf{B}$ is i.i.d. generated according to a standard normal distribution. First, we conduct our optimal rank-one attack strategy. In this strategy, we use the result from the analysis of the optimal rank-one modification matrix to design $\mathbf{a},\mathbf{b}$ and add the attack matrix $\mathbf{\Delta X}=\mathbf{a}\mathbf{b}^{\top}$ to the original data matrix $\mathbf{X}$. We then perform SVD on $\hat{\mathbf{X}}$ and select the $k$ leading principal components. Finally, we compute the distance between the selected subspace and the original subspace. We also conduct a test using a random rank-one attack strategy, in which we randomly generate $\mathbf{a},\mathbf{b}$ with each entry of $\mathbf{a},\mathbf{b}$ being i.i.d. generated according to the standard normal distribution. Then we normalize the energy of $\mathbf{a}\mathbf{b}^{\top}$to be $\eta^{2}$. For each $\eta$, we repeatedly generate 100000 pairs of $\mathbf{a}$ and $\mathbf{b}$ and compute their corresponding subspace distances. In addition, we compare it with the strategy where the modification matrix is free of rank constraint. Although our analysis is deliberately designed for general data matrices, we set the ($k+1$)th singular value to be zero so that it can be applied to the low-rank data matrix. We design the modification matrix $\mathbf{\Delta X}$ according to our analysis in this paper and calculate the subspace distance between the original subspace and that after modification. Moreover, we conduct another random attack strategy in which we randomly generate the modification matrix without any rank constraint. Each entry of the modification matrix is i.i.d. generated according to a standard normal distribution. After that, we normalize its Frobenious norm equal to $\eta$. We repeat this attack 100000 times for each $\eta$ and record its corresponding subspace distance. Furthermore, we also compare it with the strategy described in [26], which adds one adversarial data sample into the data set.

Fig. 1 demonstrates the subspace distances obtained by the five strategies. In this figure, r1-opt represents the rank-one optimal attack obtained in this paper, r1-rnd represents the maximal subspace distance obtained among the 100000 times random rank-one attacks, wr-opt stands for our optimal attack without the rank constraint, wr-rnd is the maximal subspace distance among the 100000 random attacks without the rank constraint, and ad-pca is the algorithm described in [26]. The $x$ axis is the ratio between $\eta$ and the smallest singular value of the original data matrix. From the figure, we can see our optimal strategies are much better than the ad-pca strategy. It is because our strategies can modify the data matrix, and thus have higher degrees of freedom to manipulate the data. The optimal strategies designed in this paper also have a larger subspace distance compared with their corresponding random attack strategies. In the region where $\eta/\sigma_{k}\in[0,1/\sqrt{2}]$, both of our two optimal strategies provide the same subspace distances, which can be verified by setting $\sigma_{k+1}=0$, computing $\theta^{*}$ in equation (46) and comparing it with the value in equation (15). When $\eta/\sigma_{k}>1/\sqrt{2}$, the optimal attack without the rank constraint leads to the largest subspace distance, $\pi/2$, which is much larger than the distance obtained by the optimal rank-one attack strategy. That means, without the rank constraint, it indeed provides a larger subspace distance.

In the second numerical experiment, we test these strategies except the ad-pca in the general data matrix in which the data matrix is not low-rank. In this experiment, we set $d=5$, $n=5$, and $k=3$. We randomly generate the data matrix $\mathbf{X}\in\mathbb{R}^{d\times n}$ with each entries i.i.d generated according to a standard normal distribution. We also design the optimal rank-one attack matrix and the optimal modification matrix without the rank constraint according to the analysis provided in this paper. In addition, we do random attacks 100000 times using the randomly generated modification matrix with the rank-one constraint and without the rank constraint respectively.

Fig. 2 shows the subspace distances obtained through different strategies over different energy budgets. In this figure, the $x$ axis is the ratio between $\eta$ and $\sigma_{k}-\sigma_{k+1}$. For the two random attack strategies, we demonstrate the maximal subspace distances achieved by the 100000 times random attacks. As the figure shows, both of the two random strategies have smaller subspace distances compared with their perspective optimal strategies. Different from the low-rank case, the strategy without the rank constraint provides larger subspace distances consistently over all the energy budgets.

V-B Applications

In this subsection, we use real data to illustrate the results obtained in this paper.

In particular, we illustrate the impact of adversarial attack on PCR, which is widely used in statistical learning especially when collinearity exists in the data. Ordinary regression will increase the standard error of the coefficients when there are high correlations or even collinearities between features. This happens particularly when the number of features is much larger than the number of data samples. PCR deals with this issue by performing PCA on the feature matrix and only selecting the leading $k$ principal components as the predictors, and thus dramatically decreases the number of predictors. The regression process of PCR can be seen as projecting the response values onto the subspace spanned by the leading $k$ principal components. So, the accuracy of the subspace will significantly influence the regression results. More details of PCR can be found in [33].

In this experiment, our task is to use the gasoline spectral intensity to predict its octane rating. We use the gasoline spectral data set [34], which comprises spectral intensities of 60 samples of gasoline at 401 wavelengths, and their octane ratings. Fig. 3 shows the spectral intensities of the data set. This figure indicates that the correlation of intensity among different wavelengths is very high. To complete the regression task, we can use PCR.

In this experiment, we randomly select $80$ percent of the data as the training set and the remaining $20$ percent as the test set. We choose $4$ principal components as our predictors and perform regression based on these principal components. We also record the r-squared values both in the training phase and the test phase. The r-squared value is defined as $r^{2}=1-\frac{\|\mathbf{y}-\hat{\mathbf{y}}\|^{2}}{\|\mathbf{y}-\bar{\mathbf{y}}\|^{2}}$, where $r^{2}$ is the r-squared value, $\mathbf{y}$ is the response values, $\hat{\mathbf{y}}$ is the predicted values, $\|\mathbf{y}-\bar{\mathbf{y}}\|^{2}$ represents the total variance of the response values, and $\bar{\mathbf{y}}=\textrm{mean}(\mathbf{y})\cdot\mathbf{1}$ stands for the mean vector of the response values. R-squared value measures how well the model fits the data and larger r-squared value indicates better regression. Firstly, we perform regular PCR without attack and let na-train and na-test denote the r-squared values of the training and test respectively. We then attack the feature matrix using the optimal rank-one strategy proposed in this paper with different energies and denote r1-train and r1-test as its r-squared values in the training and test processes. Finally, we also carry out the optimal attack without the rank constraint and denote wr-train, wr-test as the r-squared values in the training and test procedures.

Fig. 4 illustrates the r-squared values with different attack strategies under different energy budgets. As shown in this figure, with the increasing of the energy budget, r-squared values of training and test decrease for both attack strategies. This figure also indicates that the strategy with no rank constraint is more efficient than the rank-one strategy considering its smaller r-squared values. Furthermore, the r-squared value of the strategy without the rank constraint has a tremendous drop at the point $\eta/(\sigma_{4}-\sigma_{5})=1/\sqrt{2}$, which is consistent with our analysis that beyond this particular point the maximal subspace distance is $\pi/2$.

VI Conclusion

In this paper, we have investigated the adversarial robustness of PCA problem. We have characterized the optimal rank-one adversarial modification strategy and the optimal strategy without the rank constraint to modify the data. Our analysis has showed that both of the two strategies depend on the singular values of the data matrix and the adversary’s energy budget. We have also performed numerical simulations and investigated the impact of this attack on PCR. Both the numerical experiments and the PCR application illustrate that adversarial attacks degrade the performance of subspace learning significantly. In the future, it is of interest to investigate the defense strategy to mitigate the effects of this attack.

Appendix A Poof of the equivalence between problem (9) and problem (10)

Before giving the proof, we first examine the unitary invariant property of the Asimov distance, which is helpful in our subsequent proof.

Proposition 3**.**

Let $\mathbf{P}$ and $\mathbf{T}$ be unitary matrices, and then for the Asimov distance function $\theta(\cdot,\cdot)$, we have

[TABLE]

Proof.

First, we show $\theta(\mathbf{X}_{1},\mathbf{X}_{2})=\theta(\mathbf{P}\mathbf{X}_{1}\mathbf{T}^{\top},\mathbf{P}\mathbf{X}_{2}\mathbf{T}^{\top})$. Suppose the thin QR decompositions of $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$ are $\mathbf{X}_{1}=\mathbf{Q}_{1}\mathbf{R}_{1},\quad\mathbf{X}_{2}=\mathbf{Q}_{2}\mathbf{R}_{2}$, and then the subspace distance between the two subspaces spanned by the columns of $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$ is determined by the singular values of $\mathbf{Q}_{1}^{\top}\mathbf{Q}_{2}$. Since $(\mathbf{P}\mathbf{Q}_{1})^{\top}(\mathbf{P}\mathbf{Q}_{2})=\mathbf{Q}_{1}^{\top}\mathbf{Q}_{2}$ and right multiplying an unitary matrix does not change the singular values and the column subspace of a matrix, we have $\theta(\mathbf{X}_{1},\mathbf{X}_{2})=\theta(\mathbf{P}\mathbf{X}_{1}\mathbf{T}^{\top},\mathbf{P}\mathbf{X}_{2}\mathbf{T}^{\top})$.

Second, suppose the full SVD of $\mathbf{X}_{2}$ is $\mathbf{X}_{2}=\mathbf{U}_{2}\mathbf{\Sigma}_{2}\mathbf{V}_{2}^{\top}$, where $\mathbf{U}_{2}=[\mathbf{u}_{21},\mathbf{u}_{22},\cdots,\mathbf{u}_{2d}]$. Then

[TABLE]

which can be verified by checking that $\mathbf{P}\mathbf{U}_{2}\mathbf{\Sigma}_{2}\mathbf{V}_{2}^{\top}$ is a valid SVD of $\mathbf{P}\mathbf{X}_{2}$. It completes the proof. ∎

With the help of this proposition, let $\mathbf{P}=\mathbf{U}^{\top}$, $\mathbf{T}=\mathbf{V}^{\top}$, right multiply $\mathbf{P}$ and left multiply $\mathbf{T}^{\top}$ on both $\mathbf{X}$ and $\mathbf{\hat{X}}$, and we can simplify problem (9) as the following

[TABLE]

where we assume $n>d$, $\mathbf{\Sigma}=[\textrm{diag}(\sigma_{1},\sigma_{2},\cdots,\sigma_{k},\mathbf{0}),\mathbf{0}]\in\mathbb{R}^{d\times n}$. Also, from problem (9) to problem (47), we do variable change $\mathbf{a}\Leftarrow\mathbf{U}^{\top}\mathbf{a},\mathbf{b}\Leftarrow\mathbf{V}^{\top}\mathbf{b}$.

To further simplify this optimization problem, we split $\mathbf{a}$ and $\mathbf{b}$ into $\mathbf{a}=[\mathbf{a}_{1}^{\top},\mathbf{a}_{2}^{\top}]^{\top},\mathbf{b}=[\mathbf{b}_{1}^{\top},\mathbf{b}_{2}^{\top}]^{\top}$, where $\mathbf{a}_{1}\in\mathbb{R}^{k}$, $\mathbf{a}_{2}\in\mathbb{R}^{d-k}$, $\mathbf{b}_{1}\in\mathbb{R}^{k}$, and $\mathbf{b}_{2}\in\mathbb{R}^{n-k}$. In addition, utilizing the Householder transformation[31], we construct an orthogonal matrix

[TABLE]

where

[TABLE]

Similarly, we can construct another Householder transformation matrix $\mathbf{H}_{2}$ for $\mathbf{b}_{2}$ and the corresponding orthogonal matrix $\mathbf{M}_{2}=\text{diag}(\mathbf{I}_{k},\mathbf{H}_{2})$. Left multiplying $\mathbf{M}_{1}^{\top}$ and right multiplying $\mathbf{M}_{2}$ on $\tilde{\mathbf{Y}}$, we have

[TABLE]

where $s_{2}=\pm 1$.

Let $\mathbf{a}\Leftarrow[\mathbf{a}_{1}^{\top},\,s_{1}\|\mathbf{a}_{2}\|]^{\top}$ and $\mathbf{b}\Leftarrow[\mathbf{b}_{1}^{\top},\,s_{2}\|\mathbf{b}_{2}\|]^{\top}$. Utilizing Proposition 3, it is clear that problem (10) and problem (47) are equivalent.

Appendix B Proof of Theorem 1

The proof follows similar steps to those in [26]. In problem (10), $\tilde{\mathbf{\Sigma}}$ is a diagonal matrix with diagonal elements $\{\sigma_{1},\sigma_{2},\cdots,\sigma_{k},0\}$. The subspace spanned by $g_{k}(\mathbf{Y})$ is a $k$ dimensional subspace in $\mathbb{R}^{k+1}$. We denote this subspace as $\mathbb{Q}$, denote $\mathbb{P}$ as the subspace spanned by $\tilde{\mathbf{\Sigma}}$ and further denote their intersection as $\mathbb{T}=\mathbb{P}\cap\mathbb{Q}$. Note that $\mathbb{P}$ is not equal to $\mathbb{Q}$ (otherwise the Asimov distance will be zero), so we have $\textrm{dim}(\mathbb{P}\cup\mathbb{Q})=k+1$. Since $\textrm{dim}(\mathbb{P})+\textrm{dim}(\mathbb{Q})-\textrm{dim}(\mathbb{T})=\textrm{dim}(\mathbb{P}\cup\mathbb{Q}),$ we have $\textrm{dim}(\mathbb{T})=k-1$. Let $\mathbf{T}$ be an orthonormal basis of $\mathbb{T}$. Let $[\mathbf{T},\mathbf{p}]$ be an orthonormal basis of $\mathbb{P}$ and let $[\mathbf{T},\mathbf{q}]$ be an orthonormal basis of $\mathbb{Q}$. By the definition of Asimov distance, the subspace distance between $\mathbb{P}$ and $\mathbb{Q}$ is the angle between $\mathbf{p}$ and $\mathbf{q}$.

Firstly, it is easy to see that $a_{k+1}\neq 0$. Otherwise, $\mathbb{Q}$ will be equal to $\mathbb{P}$, which means that their Asimov distance is zero.

Secondly, it is easy to see $\mathbf{q}\in\textrm{span}[\mathbf{T},\mathbf{p},\mathbf{e}_{k+1}]$, where $\mathbf{e}_{k+1}$ is an ordinary basis vector that only has element $1$ in the ($k+1$)th coordinate. Since $\mathbf{T}$ is orthogonal to $\mathbf{q}$, we have $\mathbf{q}\in\textrm{span}[\mathbf{p},\mathbf{e}_{k+1}]$. It is easy to see that the larger variance in the direction of $\mathbf{p}$ is, the closer $\mathbf{p}$ and $\mathbf{q}$ will be. Then we should select $\mathbf{p}$ as the direction with the smallest variance in $\mathbf{X}$. Since we are assuming that $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{k}$, $\mathbf{p}$ should be $\mathbf{e}_{k}$.

Thirdly, for a fixed direction of $\mathbf{a}$, let $\hat{\mathbf{a}}$ be the projection of $\mathbf{a}$ onto $\text{span}[\mathbf{e}_{k},\mathbf{e}_{k+1}]$. Clearly, $\mathbf{q}$ will be closer to $\mathbf{a}$ as $\hat{\mathbf{a}}$ grows. As a result, the angle between $\mathbf{q}$ and $\mathbf{p}$ will be larger. This also implies that the length of $\mathbf{a}$ should be maximized: $\|\mathbf{a}\|=\eta$. Hence, the Asimov distance is maximized when $\mathbf{a}=\hat{\mathbf{a}}$ and $\|\mathbf{a}\|=\eta$, implying that $\mathbf{a}$ only has nonzero elements in its $k$th and $k+1$th coordinates.

Finally, for a fixed $\mathbf{a}$ in the form of (11), the projected variance of $\mathbf{Y}$ on the direction of $\mathbf{e}_{k}$ is $v_{1}=\sum_{i\neq k}(a_{k}b_{i})^{2}+(a_{k}b_{k}+\sigma_{k})^{2}=a_{k}^{2}+\sigma_{k}^{2}+2a_{k}b_{k}\sigma_{k}$ and the projected variance of $\mathbf{Y}$ on the direction of $\mathbf{e}_{k+1}$ is $v_{2}=\sum_{i}(a_{k+1}b_{i})^{2}=a_{k+1}^{2}$. To maximize the Asimov distance, we need to make $v_{1}$ small and $v_{2}$ large. Apparently, for fixed $\mathbf{a}$, $v_{1}$ is minimized when $b_{k}=-\text{sign}(a_{k})$, which implies $b_{i}=0,\forall i\neq k$. To avoid the sign ambiguity, we set $b_{k}=1$.

Appendix C Proof of Theorem 2

This proof follows the similar steps in the proof of the low-rank case. Denote $\mathbb{P}$ as the subspace spanned by $g_{k}(\mathbf{\Sigma})$ and $\mathbb{Q}$ as the subspace spanned by $g_{k}(\mathbf{Y})$, and denote their intersection as $\mathbb{T}=\mathbb{P}\cap\mathbb{Q}$. We further denote $\mathbf{T}$ as an orthonormal basis of $\mathbb{T}$, $[\mathbf{T},\mathbf{p}]$ as an orthonormal basis of $\mathbb{P}$, and $[\mathbf{T},\mathbf{q}]$ as an orthonormal basis of $\mathbb{Q}$. From the definition of Asimov distance, the subspace distance between $\mathbb{P}$ and $\mathbb{Q}$ is the subspace distance between the span of $\mathbf{p}$ and the span of $\mathbf{q}$.

First, it is apparent that $\mathbf{q}\in\text{span}[\mathbf{T},\mathbf{p},\mathbf{e}_{k+1},\mathbf{e}_{k+2},\cdots,\mathbf{e}_{d}]$. Since $\mathbf{q}\perp\mathbf{T}$, we have $\mathbf{q}\in\text{span}[\mathbf{p},\mathbf{e}]$, where $\mathbf{e}\in\text{span}[\mathbf{e}_{k+1},\cdots,\mathbf{e}_{d}]$. It is easy to see that the subspace distance between the span of $\mathbf{q}$ and the span of $\mathbf{p}$ will be large if the variance of $\mathbf{\Sigma}$ in the span of $\mathbf{p}$ is large and the variance of $\mathbf{\Sigma}$ in the span of $\mathbf{q}$ is small. So we should select $\mathbf{p}$ as the direction in $\text{span}[\mathbf{e}_{1},\cdots,\mathbf{e}_{k}]$ that has the smallest variance of $\mathbf{\Sigma}$ and select $\mathbf{e}$ as the direction among $\text{span}[\mathbf{e}_{k+1},\cdots,\mathbf{e}_{n}]$ that has the largest variance of $\mathbf{\Sigma}$. Since $\mathbf{e}\in\text{span}[\mathbf{e}_{1},\mathbf{e}_{2},\cdots,\mathbf{e}_{d}]$ and $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{k}\geq\sigma_{k+1}\geq\cdots\geq\sigma_{d}$, $\mathbf{p}$ should be $\mathbf{e}_{k}$ and $\mathbf{e}$ should be $\mathbf{e}_{k+1}$. So, we have $\mathbf{q}\in\text{span}[\mathbf{e}_{k},\mathbf{e}_{k+1}]$.

Second, for a fixed direction of $\mathbf{a}$, let $\hat{\mathbf{a}}$ be the projection of $\mathbf{a}$ onto $\text{span}[\mathbf{e}_{k},\mathbf{e}_{k+1}]$. It is easy to see that $\mathbf{q}$ will be closer to $\mathbf{a}$ as $\hat{\mathbf{a}}$ grows, and as a result, the angle between $\mathbf{q}$ and $\mathbf{p}$ will be larger. This implies the length of $\mathbf{a}$ should be maximized, which indicates $\|\mathbf{a}\|=\eta$ and the distance is maximized when $\mathbf{a}=\hat{\mathbf{a}}$. It also indicates $a_{i}=0$ if $i\neq k,\,k+1$.

Finally, for a fixed $\mathbf{a}$ in the form of (18), the projected variance of $\mathbf{Y}$ in the direction of $\mathbf{e}_{k}$ is $v_{k}=\sum_{i\neq k}(a_{k}b_{i})^{2}+(a_{k}b_{k}+\sigma_{k})^{2}=a_{k}^{2}+\sigma_{k}^{2}+2a_{k}b_{k}\sigma_{k}$ and the projected variance of $\mathbf{Y}$ in the direction of $\mathbf{e}_{k+1}$ is $v_{k+1}=\sum_{i\neq k+1}(a_{k+1}b_{i})^{2}+(\sigma_{k+1}+a_{k+1}b_{k+1})^{2}=a_{k+1}^{2}+\sigma_{k+1}^{2}+2a_{k+1}b_{k+1}$. To maximize the Asimov distance, we should make $v_{k}$ small and make $v_{k+1}$ large. With the constraint that $\|\mathbf{b}\|=1$, we should have $b_{k}^{2}+b_{k+1}^{2}=1$, which implies $b_{i}=0$ for all $i\neq k$ and $i\neq(k+1)$.

As shown above, the optimal $\mathbf{a}$ and $\mathbf{b}$ should be in the form of (18) and (19), which completes our proof.

Appendix D Proof of Theorem 3

The optimal solution to problem (25) either locates at the boundary, or at the stationary points.

We first characterize the stationary points. At the stationary points, the value $(\alpha^{*},\beta^{*})$ satisfies the necessary conditions

[TABLE]

Since $\sin\varphi^{*}\neq 0$, we have

[TABLE]

in which

[TABLE]

Eliminating $\sin(\alpha)\sin(\beta)$ from (50), we have

[TABLE]

where $C=2\eta\sigma_{k}(\sigma_{k}^{2}-\sigma_{k+1}^{2})(\sigma_{k+1}^{2}+\eta^{2}-\sigma_{k}^{2})$, $D=(\sigma_{k}^{2}-\sigma_{k+1}^{2})\left(-(\sigma_{k}^{2}-\sigma_{k+1}^{2})^{2}-2\eta^{2}(\sigma_{k}^{2}+\sigma_{k+1}^{2})+3\eta^{4}\right)$ and $F=\eta\sigma_{k}\left(\sigma_{k}^{4}+\sigma_{k+1}^{4}+\eta^{4}-2\sigma_{k}^{2}\sigma_{k+1}^{2}-2\sigma_{k}^{2}\eta^{2}-2\sigma_{k+1}^{2}\eta^{2}\right).$

Further, we rewrite the first equation of (50) as

[TABLE]

where $c=\sigma_{k}(\sigma_{k+1}^{2}+\eta^{2}-\sigma_{k}^{2})$, $d=\sigma_{k+1}(\sigma_{k}^{2}+\eta^{2}-\sigma_{k+1}^{2}),$ and $e=2\eta\sigma_{k}\sigma_{k+1}.$

Combining (51) and (52) and eliminating $\cos^{2}(\alpha)$ and $\cos^{2}(\beta)$, we have

[TABLE]

The left side of the equation is a quadratic function with respect to $r=\cos(\alpha)\cos(\beta)$. The two roots are:

[TABLE]

Note that $\eta\in[0,\sigma_{k}-\sigma_{k+1})$, so we have $r_{1}\in(-\frac{\sigma_{k}}{\sigma_{k}+\sigma_{k+1}},0]$, $r_{2}\in(-\infty,-\frac{\sigma_{k}}{\sigma_{k}-\sigma_{k+1}})$. Since $|\cos(\alpha)\cos(\beta)|\leq 1$, $\frac{\sigma_{k}}{\sigma_{k}+\sigma_{k+1}}<1$, and $\frac{\sigma_{k}}{\sigma_{k}-\sigma_{k+1}}>1$, we should only retain the first root $r_{1}$. Substitute $\cos(\alpha)\cos(\beta)=r_{1}=-\frac{\eta\sigma_{k}}{\sigma_{k}^{2}-\sigma_{k+1}^{2}}$ into (51), and we have $C\cos^{4}(\alpha)+(Dr_{1}+F)\cos^{2}(\alpha)+Cr_{1}^{2}=0.$ The left side of the equation is a quadratic function with respect to $s=\cos^{2}(\alpha)$, so we can easily find its roots. Let us denote $s_{1}$ and $s_{2}$ as the two roots:

[TABLE]

where $H=\sigma_{k}^{4}+\sigma_{k+1}^{4}+\eta^{4}-2\sigma_{k}^{2}\sigma_{k+1}^{2}-2\sigma_{k}^{2}\eta^{2}-2\sigma_{k+1}^{2}\eta^{2}$. We need to check that $H$ is positive. Viewing $H$ as a function of $\eta$ and taking derivative, we have $H^{\prime}(\eta)=2\eta(2\eta^{2}-2(\sigma_{k}^{2}+\sigma_{k+1}^{2}))<0$. Since $\eta^{2}\in[0,(\sigma_{k}-\sigma_{k+1})^{2})$, we have $H(\eta)\in(0,(\sigma_{k}^{2}-\sigma_{k+1}^{2})^{2}]$.

As $\cos(\alpha)^{2}\leq 1$, we need to check whether $s_{1},s_{2}\in[0,1]$.

Firstly, as $H$ is a decreasing function of $\eta$ in the considered range, $s_{1}$ is a increasing function of $\eta$. Therefore, we have $\min(s_{1})=s_{1}(\eta)|_{\eta=0}=0$ and $\max(s_{1})=s_{1}(\eta)|_{\eta=\sigma_{k}-\sigma_{k+1}}=\frac{\sigma_{k}}{\sigma_{k}+\sigma_{k+1}}<1$. Hence, $s_{1}$ is a valid solution.

Secondly, it is easy to check that $s_{2}$ is a decreasing function of $\eta$. So, we have $\max(s_{2})=s_{2}(\eta)|_{\eta=0}=1$ and $\min(s_{2})=s_{2}(\eta)|_{\eta=\sigma_{k}-\sigma_{k+1}}=\frac{\sigma_{k}}{\sigma_{k}+\sigma_{k+1}}<1$, which means $s_{2}$ is also a valid solution. Hence, we have two stationary points

[TABLE]

Since there are two sets of solutions in (53), we should determine which one is better. The variance of $\mathbf{Y}$ in the direction of $\mathbf{e}_{k}$ is $v_{k}=\cos^{2}(\alpha)+\sigma_{k}^{2}+2\cos(\alpha)\cos(\beta)$ and the variance of $\mathbf{Y}$ in the direction of $\mathbf{e}_{k+1}$ is $v_{k+1}=\sin^{2}(\alpha)+\sigma_{k+1}^{2}+2\sin(\alpha)\sin(\beta)$. Both of the two sets of solutions in (53) lead to $\cos(\alpha)\cos(\beta)=-\frac{\eta\sigma_{k}}{\sigma_{k}^{2}-\sigma_{k+1}^{2}}$ and $\sin(\alpha)\sin(\beta)=\frac{\eta\sigma_{k+1}}{\sigma_{k}^{2}-\sigma_{k+1}^{2}}$. For fixed $\cos(\alpha)\cos(\beta)$ and $\sin(\alpha)\sin(\beta)$, the smaller $\cos^{2}(\alpha)$ is, the smaller $v_{k}$ will be, and the larger the subspace distance will be. Hence, we conclude the stationary point that satisfies

[TABLE]

leads to a larger subspace distance.

Finally, it is easy to compute the objective values of problem (25) at the boundary points. Comparing these values with the objective values induced by the point in equation (54), we can readily conclude the point in equation (54) gives a larger objective value. In summary, given that $\alpha\in[0,\pi/2]$ and $\beta\in[\pi/2,\pi]$, the optimal $\alpha$ and $\beta$ are shown in (26).

Appendix E Proof of theorem 4

The proof has two main steps. In the first step, we show that non-zero entries of $\mathbf{B}$ are in the $k$th and ($k+1$)th rows. In the second step, we will further prove the entries except in the $k$th and ($k+1$)th columns should be zero.

In the first step, we follow similar proof procedures in Theorem 2. We use $\mathbb{P}$ to denote the subspace spanned by $g_{k}(\mathbf{\Sigma})$ and $\mathbb{Q}$ to denote the subspace spanned by $g_{k}(\mathbf{Y})$. We also use $\mathbb{T}$ to represent the intersection of the two subspaces and further denote $\mathbf{T}$ as one set of orthonormal bases of $\mathbb{T}$, $[\mathbf{T},\mathbf{p}]$ as one set of orthonormal bases of $\mathbb{P}$ and $[\mathbf{T},\mathbf{q}]$ as one set of orthonormal bases of $\mathbb{Q}$. So, the subspace distance between $\mathbb{P}$ and $\mathbb{Q}$ is the subspace distance between the subspace spanned by $\mathbf{p}$ and that spanned by $\mathbf{q}$. Following the same arguments in Theorem 2, by setting all the entries of $\mathbf{B}$ to be zero except the $k$th and ($k+1$)th rows, we can guarantee achieving the maximal subspace distance and further we have $\mathbf{q}\in\text{span}[\mathbf{e}_{k},\mathbf{e}_{k+1}]$ and $\mathbf{p}=\mathbf{e}_{k}$.

In the second step, since the non-zero elements of $\mathbf{B}$ only locate in the $k$th and ($k+1$)th rows and $\mathbf{q}\in\text{span}[\mathbf{e}_{k},\mathbf{e}_{k+1}]$, it indicates $\mathbf{q}$ is the direction with the maximal variance on the span of $\mathbf{e}_{k}$ and $\mathbf{e}_{k+1}$. Assuming $\mathbf{q}=[0,\cdots,\cos(\gamma),\sin(\gamma),\cdots,0]^{\top}$ with $\cos(\gamma)$ and $\sin(\gamma)$ being in the $k$th and ($k+1$)th coordinates respectively and according to the definition of principal components, we can find $\gamma$ by solving the optimization problem

[TABLE]

Plug $\mathbf{q}=[0,\cdots,\cos(\gamma),\sin(\gamma),\cdots,0]^{\top}$ into the objective function, and we have

[TABLE]

where $b_{x1}=\|\mathbf{b}_{k}+\mathbf{e}_{k}\sigma_{k}\|^{2}$, $b_{x2}=\|\mathbf{b}_{k+1}+\mathbf{e}_{k+1}\sigma_{k+1}\|^{2}$, $b_{y}=2(\mathbf{b}_{k}+\mathbf{e}_{k}\sigma_{k})^{\top}(\mathbf{b}_{k+1}+\mathbf{e}_{k+1}\sigma_{k+1}),$ with $\mathbf{b}_{k}$ and $\mathbf{b}_{k+1}$ being the transpose of the $k$th and ($k+1$)th rows of $\mathbf{B}$ respectively and $\mathbf{e}_{k}\in\mathbb{R}^{n}$, $\mathbf{e}_{k+1}\in\mathbb{R}^{n}$ being the standard bases.

We can solve (56) by computing the first principal component of the middle matrix of the right hand of (56). Using the result from equation (23), we have $\gamma=0.5\text{atan2}(b_{y},b_{x}),$ where $b_{x}=b_{x1}-b_{x2}.$ Since the subspace distance is the distance between $\mathbf{q}$ and $\mathbf{e}_{k}$, it is apparent that the subspace distance is $|\gamma|$. To maximize $|\gamma|$, we first determine the sign of $b_{y}$ or $b_{x}$. We have

[TABLE]

where inequality (57) is the result of the energy constraint that $\eta\geq\|\mathbf{B}\|_{\textrm{F}}=\sqrt{\|\mathbf{b}_{k}\|^{2}+\|\mathbf{b}_{k+1}\|^{2}}\geq\frac{1}{\sqrt{2}}(\|\mathbf{b}_{k}\|+\|\mathbf{b}_{k+1}\|)$, and inequality (58) is due to the assumption that $\eta<\frac{\sigma_{k}-\sigma_{k+1}}{\sqrt{2}}$. In summary, $b_{x}$ is positive. Using the property of atan2 function, when $b_{x}>0$, maximizing $|\gamma|$ is equivalent to maximizing $|b_{y}/b_{x}|$. Thus, we can formulate our problem as

[TABLE]

In the objective function,

[TABLE]

where $\mathbf{b}_{1}=[b_{k,1},\,b_{k,2},\,\cdots,\,b_{k,k-1},\,b_{k,k+2},\,\cdots,\,b_{k,n}]^{\top}$ and

$\mathbf{b}_{2}=[b_{k+1,1},\,b_{k+1,2},\,\cdots,\,b_{k+1,k-1},\,b_{k+1,k+2},\,\cdots,\,b_{k+1,n}]^{\top}$ which are the vectors obtained by deleting the $k$th and ($k+1$)th elements of $\mathbf{b}_{k}$ and $\mathbf{b}_{k+1}$ respectively. We can change the sign of $b_{y}/b_{x}$ by changing the signs of $\mathbf{b}_{1},b_{k+1,k}$, and $b_{k,k+1}$. Since both of the values $b_{y}/b_{x}$ and $-b_{y}/b_{x}$ are obtainable, we can remove the absolute value operation. Thus, our objective can be further simplified to maximize $b_{y}/b_{x}$. To complete the proof of Theorem 4, we should further demonstrate that when the optimality of our objective function is obtained, $\mathbf{b}_{1}$ and $\mathbf{b}_{2}$ should be vectors with all their entries being zero. To prove that, we examine the objective function further

[TABLE]

Inequality (60) implies that the optimal value is determined by the norms of $\mathbf{b}_{1}$ and $\mathbf{b}_{2}$ instead of their specific values. Inequality (61) is true as

[TABLE]

Inequality (62) is due to $-(\sqrt{b_{k+1,k+1}^{2}+\|\mathbf{b}_{2}\|^{2}}+\sigma_{k+1})^{2}\leq-\|\mathbf{b_{2}}\|^{2}-(b_{k+1,k+1}+\sigma_{k+1})^{2}$. The equalities in (61) and (62) hold when $\|\mathbf{b}_{1}\|=0$ and $\|\mathbf{b}_{2}\|=0$. This means that, for any feasible solution $(\mathbf{b}_{1},\mathbf{b}_{2},b_{k,k},b_{k,k+1},b_{k+1,k},b_{k+1,k+1})$ in (59), there is another corresponding feasible solution $(\mathbf{0},\mathbf{0},b_{k,k},\sqrt{b_{k,k+1}^{2}+\|\mathbf{b}_{1}\|^{2}},b_{k+1,k},\sqrt{b_{k+1,k+1}^{2}+\|\mathbf{b}_{2}\|^{2}})$, which has a larger objective value. In conclusion, $\mathbf{b}_{1}$ and $\mathbf{b}_{2}$ should be zero vectors when the optimality of (59) is obtained. This completes our proof.