Z-Eigenvector Residual Recovery for Dominant Fourth-Order Structure Missed by PCA

1. Introduction

Principal component analysis (PCA) is one of the standard methods for dimension reduction, approximation, and denoising. Its origins go back to the least-squares geometric construction of Pearson and the variance-based formulation of Hotelling [1,2]; see also [3]. For centered data, PCA orders directions according to the quadratic form determined by the covariance matrix. It is therefore intrinsically a second-order method. A direction with a large second moment is favored, regardless of whether that direction represents the structural component one wishes to preserve.

Higher-order criteria provide a different way to rank directions. Projection pursuit searches for low-dimensional projections that maximize a chosen measure of non-Gaussian or otherwise interesting structure [4,5]. Fourth-order cumulants and related contrast functions also play a central role in independent component analysis and blind source separation [6,7]. Those theories have objectives and normalizations that differ from the present one, but they illustrate the general principle that second-order information need not reveal all distinguished directions in a data set.

In this paper, we introduce a Z-eigenvector residual recovery mechanism based on successive orthogonal maximization. Symmetric tensor eigenvalues and their variational interpretation were developed in [9,10,11]. The maximization problem is also closely related to symmetric best rank-one approximation [12], and shifted higher-order power methods [13]. General background on tensor decompositions may be found in [8].

The method presented here is more specific than a general tensor decomposition. We first compute a PCA subspace B, remove its contribution from the data, and then maximize the residual quartic form over $B^{\perp }$. The resulting direction is used to augment the original PCA subspace. This residual construction is motivated by a conceivable scenario: a structural component may occur less frequently, or may contribute less covariance energy, than a nuisance component while still producing significantly large fourth moment. In that case, second- and fourth-order methods rank the two directions in opposite order.

The paper is organized as follows. Section 2 introduces fourth-order Z-eigenvectors and the residual recovery procedure. Section 3 proves the exact second- versus fourth-order separation result. Section 4 establishes the perturbation bound and the controlled reconstruction-error comparisons. Section 5 presents a three-dimensional example, and the final section summarizes the conclusions and limitations.

2. Fourth-Order Z-Eigenvectors and Residual Recovery

We first recall the definition of a Z-eigenpair and then describe the residual recovery procedure used throughout the paper. Let $\mathcal{C}=\left(c_{i_1\dots i_m}\right)$ be a real symmetric tensor of order $m\ge 2$ and dimension $n\ge 2$. This means,

$$c_{i_1\dots i_m}=c_{i_{\pi \left(1\right)}\dots i_{\pi \left(m\right)}},1\le i_1,\dots ,i_m\le n.$$

for every permutation $\pi$ of $\{1,\dots ,m\}$. The homogeneous polynomial associated with $\mathcal{C}$ is

$$f_{\mathcal{C}}\left(x\right)=\mathcal{C}{x}^m:={\displaystyle \sum _{i_1,\dots ,i_m=1}^n}{c}_{i_1\dots i_m}{x}_{i_1}\dots x_{i_m}.$$

For $x\in {\mathbb{R}}^n$, the vector contraction $\mathcal{C}{x}^{m-1}\in {\mathbb{R}}^n$ is defined componentwise via

$${\left(\mathcal{C}{x}^{m-1}\right)}_i:={\displaystyle \sum _{i_2,\dots ,i_m=1}^n}{c}_{i{i}_2\dots i_m}{x}_{i_2}\dots x_{i_m}.$$

We use the standard definition of a Z-eigenpair from [9,10,11].

Definition 1

(Z-eigenpair). A pair $(\lambda ,v)\in \mathbb{R}\times ({\mathbb{R}}^n\setminus \left\{\mathbf{0}\right\})$ is called a Z-eigenpair of $\mathcal{C}$ if

$$\mathcal{C}{v}^{m-1}=\lambda v,{\parallel v\parallel }_2=1.$$

The scalar λ is a Z-eigenvalue of $\mathcal{C}$ and v is its associated Z-eigenvector.

We shall henceforth focus on real symmetric tensors of order four. By a Lagrange multiplier argument, every unit vector attaining ${\displaystyle \underset{{\parallel u\parallel }_2=1}{max}\mathcal{C}{u}^4}$ is a Z-eigenvector of $\mathcal{C}$. We shall refer to a vector in ${\displaystyle \underset{{\parallel u\parallel }_2=1}{argmax}\mathcal{C}{u}^4}$ as a dominant fourth-order direction of $\mathcal{C}$.

For the remainder of this section, let $z_1,\dots ,z_K\in {\mathbb{R}}^n$ be centered data vectors so that

$${\displaystyle \frac{1}{K}}{\displaystyle \sum _{q=1}^K}{z}_q=0.$$

We define the associated covariance matrix and fourth-order symmetric tensor by

$$\widehat{\Sigma }:={\displaystyle \frac{1}{K}}{\displaystyle \sum _{q=1}^K}{z}_q{z}_q^{\top }\mathrm{and}\mathcal{C}:={\displaystyle \frac{1}{K}}{\displaystyle \sum _{q=1}^K}{z}_q^{\otimes 4}.$$

For every $u\in {\mathbb{R}}^n$, we then have

$$\mathcal{C}{u}^4={\displaystyle \frac{1}{K}}{\displaystyle \sum _{q=1}^K}{\langle z_q,u\rangle }^4\mathrm{and}\mathcal{C}{u}^3={\displaystyle \frac{1}{K}}{\displaystyle \sum _{q=1}^K}{\langle z_q,u\rangle }^3{z}_q.$$

Thus, the tensor-vector contraction can be evaluated directly from the data without explicitly forming or storing all $n^4$ entries of $\mathcal{C}$.

Our objective is not to compute a complete low-rank decomposition of the fourth-order tensor. Instead, we use a dominant fourth-order direction to augment a PCA subspace and recover residual structure that is weak in covariance. The construction is as follows.

1.

Fix an integer $1\le s<n$. Let the matrix

$$U_s=[u_1,\dots ,u_s]\in {\mathbb{R}}^{n\times s}$$

have orthonormal columns consisting of eigenvectors associated with the s largest eigenvalues of $\widehat{\Sigma }$. We define

$$B:=range\left(U_s\right)\mathrm{and}{P}_B:=U_s{U}_s^{\top }.$$

2.

We also define the residual data vectors by

$$r_q:=(I-P_B)z_q,q=1,\dots ,K,$$

and form the residual fourth-order tensor

$${\mathcal{C}}_{\mathrm{res}}:={\displaystyle \frac{1}{K}}{\displaystyle \sum _{q=1}^K}{r}_q^{\otimes 4}.$$

Since each $r_q$ belongs to $B^{\perp }$, we have that

$${\mathcal{C}}_{\mathrm{res}}{u}^3\in B^{\perp }\mathrm{for}\mathrm{every}u\in {\mathbb{R}}^n.$$

3.

If ${\mathcal{C}}_{\mathrm{res}}=0$, then no residual direction is added.

If ${\mathcal{C}}_{\mathrm{res}}\ne 0$, then we select a unit vector

$$\widehat{\mathbf{v}}\in \underset{\begin{array}{c}u\in B^{\perp }\\ {\parallel u\parallel }_2=1\end{array}}{argmax}{\mathcal{C}}_{\mathrm{res}}{u}^4.$$

It follows immediately $\widehat{\mathbf{v}}$ corresponds to a residual direction attaining the largest fourth-order response. Since ${\mathcal{C}}_{\mathrm{res}}{u}^3\in B^{\perp }$ for every $u\in {\mathbb{R}}^n$, applying the Lagrange multiplier argument for the maximization problem on the unit sphere in $B^{\perp }$ gives

$${\mathcal{C}}_{\mathrm{res}}{\widehat{\mathbf{v}}}^3=\widehat{\lambda }\widehat{\mathbf{v}}\mathrm{and}\widehat{\lambda }={\mathcal{C}}_{\mathrm{res}}{\widehat{\mathbf{v}}}^4.$$

Thus, $\widehat{\mathbf{v}}$ is a Z-eigenvector of ${\mathcal{C}}_{\mathrm{res}}$.

4.

In the case ${\mathcal{C}}_{\mathrm{res}}\ne 0$, let $y\in {\mathbb{R}}^n$ be a target vector to which the constructed projection is applied. The PCA plus residual-recovery estimator is then defined by

$${\widehat{y}}_{\mathrm{rec}}:=P_{B}y+{\widehat{\mathbf{v}}\widehat{\mathbf{v}}}^{\top }(I-P_B)y=\left(P_B+{\widehat{\mathbf{v}}\widehat{\mathbf{v}}}^{\top }\right)y.$$

The second equality follows from $\widehat{\mathbf{v}}\in B^{\perp }$. Thus, $P_B+{\widehat{\mathbf{v}}\widehat{\mathbf{v}}}^{\top }$ is an orthogonal projector of rank $s+1$. The first term retains the component of y in the leading covariance subspace, while the second term restores its component along the dominant fourth-order residual direction.

In numerical implementations, the contraction ${\mathcal{C}}_{\mathrm{res}}{u}^3$ can be evaluated from

$${\mathcal{C}}_{\mathrm{res}}{u}^3={\displaystyle \frac{1}{K}}{\displaystyle \sum _{q=1}^K}{\langle r_q,u\rangle }^3{r}_q.$$

A stationary Z-eigenvector may be sought using the shifted symmetric higher-order power method (SS-HOPM) of [13], preferably with multiple initializations. In general, a numerical iteration may converge to a local stationary point; our results below nevertheless concern global maximizers of the stated quartic forms.

In the next two sections, we first analyze an exact model in which covariance and fourth-order maximization select different subspaces. We then establish stability under controlled perturbations of the residual fourth-order tensor.

3. A Fourth-Order-Dominant Regime in Which Z-Eigenvector Recovery Outperforms PCA

We now present a model in which the distinction between second- and fourth-order is transparent. Let $\{v_1,\dots ,v_k;w_1,\dots ,w_r\}$ be an orthonormal set in ${\mathbb{R}}^n$. Let $S=span\{v_1,\dots ,v_k\}$ denote the structural subspace and $N=span\{w_1,\dots ,w_r\}$ represent the nuisance subspace.

Let $Z\in {\mathbb{R}}^n$ be a centered random vector supported on these coordinate directions as follows. For each structural direction $v_j$, with probability $p_j/2$, we assign $Z=a_j{v}_j$, and with probability $p_j/2$, we assign $Z=-a_j{v}_j$. For each nuisance direction $w_{\ell }$, with probability $q_{\ell }/2$, we assign $Z=b_{\ell }{w}_{\ell }$, and with probability $q_{\ell }/2$, we assign $Z=-b_{\ell }{w}_{\ell }$. Here $p_j,q_{\ell }>0$ and $a_j,b_{\ell }>0$ are chosen such that

$${\displaystyle \sum _{j=1}^k}{p}_j+{\displaystyle \sum _{\ell =1}^r}{q}_{\ell }=1.$$

Consequently, the covariance matrix and fourth-order tensor are

$${\Sigma }_2={\displaystyle \sum _{j=1}^k}{p}_j{a}_j^2{v}_j{v}_j^{\top }+{\displaystyle \sum _{\ell =1}^r}{q}_{\ell }{b}_{\ell }^2{w}_{\ell }{w}_{\ell }^{\top }$$

(1)

and

$${\mathcal{C}}_4={\displaystyle \sum _{j=1}^k}{p}_j{a}_j^4{v}_j^{\otimes 4}+{\displaystyle \sum _{\ell =1}^r}{q}_{\ell }{b}_{\ell }^4{w}_{\ell }^{\otimes 4}.$$

(2)

Let ${\nu }_1,\dots ,{\nu }_k$ be chosen recursively so that

$${\nu }_t\in \underset{\begin{array}{c}{\parallel u\parallel }_2=1\\ u\perp {\nu }_1,\dots ,{\nu }_{t-1}\end{array}}{argmax}{\mathcal{C}}_4{u}^4,t=1,\dots ,k,$$

(3)

where the orthogonality condition is vacuous when $t=1$.

It is worth noting for the special fourth-order tensor ${\mathcal{C}}_4$ constructed in (2), the successive orthogonal maximizers are indeed coordinate directions and therefore are automatically Z-eigenvectors of ${\mathcal{C}}_4$. For a general real symmetric tensor, however, an orthogonally constrained maximizer need not be a Z-eigenvector of the original tensor.

Theorem 1

(Fourth-order advantage). Assume $r\ge k$ and

$$\begin{array}{cc}\hfill & \underset{1\le \ell \le r}{min}{q}_{\ell }{b}_{\ell }^2>\underset{1\le j\le k}{max}{p}_j{a}_j^2,\mathrm{while}\hfill \end{array}$$

(b)

$$\begin{array}{cc}\hfill & \underset{1\le j\le k}{min}{p}_j{a}_j^4>\underset{1\le \ell \le r}{max}{q}_{\ell }{b}_{\ell }^4.\hfill \end{array}$$

(b)

We then have:

1.

The leading k covariance eigenvectors of ${\Sigma }_2$ are nuisance directions in N, whereas the successive orthogonal maximizers of ${\mathcal{C}}_4$, $\{{\nu }_1,\dots ,{\nu }_k\}$ in (3), coincide with the structural directions $\{v_1,\dots ,v_k\}$, up to signs and ordering.

2.

Let $P_S$ denote the orthogonal projector onto S and let $P_{\mathrm{PCA}}$ denote an orthogonal projector onto a leading k-dimensional eigenspace of ${\Sigma }_2$. Let $y=x+\epsilon ,$ where $x\in S$ and $\epsilon \in {\mathbb{R}}^n$ is an arbitrary error vector. Then

$$\parallel P_S{y-x\parallel }_2={\parallel P_S\epsilon \parallel }_2,$$

(5)

whereas

$$\parallel P_{\mathrm{PCA}}{y-x\parallel }_2^2={\parallel x\parallel }_2^2+{\parallel P_{\mathrm{PCA}}\epsilon \parallel }_2^2.$$

(6)

Consequently,

$$\parallel P_{\mathrm{PCA}}{y-x\parallel }_2^2-\parallel P_S{y-x\parallel }_2^2={\parallel x\parallel }_2^2+\parallel P_{\mathrm{PCA}}{\epsilon \parallel }_2^2-{\parallel P_S\epsilon \parallel }_2^2.$$

(7)

Thus, the fourth-order structural projection has smaller squared reconstruction error whenever

$$\parallel P_S{\epsilon \parallel }_2^2-\parallel P_{\mathrm{PCA}}{\epsilon \parallel }_2^2<{\parallel x\parallel }_2^2.$$

(8)

In particular, the sufficient condition

$${\parallel \epsilon \parallel }_2<{\parallel x\parallel }_2$$

(9)

guarantees strict improvement. In the noiseless case $\epsilon =0$, the improvement is exactly ${\parallel x\parallel }_2^2$.

Proof. 

Equation (1) shows that the covariance eigenvalues in the structural directions are $p_j{a}_j^2$, while the covariance eigenvalues in the nuisance directions are $q_{\ell }{b}_{\ell }^2$. Assumption (a) implies that the leading k covariance directions are nuisance directions.

On the other hand, for the fourth-order tensor ${\mathcal{C}}_4$, we set

$${\gamma }_j=p_j{a}_j^4\mathrm{and}{\eta }_{\ell }=q_{\ell }{b}_{\ell }^4.$$

If $u={\sum }_{j=1}^k{\alpha }_j{v}_j+{\sum }_{\ell =1}^r{\beta }_{\ell }{w}_{\ell }+h$, where h is orthogonal to the span of the displayed directions, then

$${\mathcal{C}}_4{u}^4={\displaystyle \sum _{j=1}^k}{\gamma }_j{\alpha }_j^4+{\displaystyle \sum _{\ell =1}^r}{\eta }_{\ell }{\beta }_{\ell }^4\le M\left({\displaystyle \sum _{j=1}^k}{\alpha }_j^4+{\displaystyle \sum _{\ell =1}^r}{\beta }_{\ell }^4\right)\le M{\left({\displaystyle \sum _{j=1}^k}{\alpha }_j^2+{\displaystyle \sum _{\ell =1}^r}{\beta }_{\ell }^2\right)}^2\le M,$$

(10)

where

$$M=max\{\underset{j}{max}{\gamma }_j,\underset{\ell }{max}{\eta }_{\ell }\}.$$

Equalities hold throughout (10) only if $h=0$ and exactly one of the coefficients ${\alpha }_j,{\beta }_{\ell }$ has absolute value one. By assumption (b), since every ${\gamma }_j$ is larger than every ${\eta }_{\ell }$, a maximizing coordinate must be structural. After one such direction is selected, its coefficient is zero on the orthogonal complement. Repeating the same argument selects each of the remaining structural directions before any nuisance direction is chosen. Thus, the recursively selected fourth-order directions are exactly $\{v_1,\dots ,v_k\}$, up to signs and ordering. Assertion (1) follows.

Since $x\in S$, we have $P_{S}x=x$. It follows that

$$P_{S}y-x=P_S(x+\epsilon )-x=P_S\epsilon ,$$

which proves (5).

By assertion (1), the range of $P_{\mathrm{PCA}}$ is contained in the nuisance subspace N, which is orthogonal to S. Hence,

$$P_{\mathrm{PCA}}x=0\mathrm{and}{P}_{\mathrm{PCA}}y-x=P_{\mathrm{PCA}}\epsilon -x.$$

Since $P_{\mathrm{PCA}}\epsilon \in N$ and $x\in S$, $P_{\mathrm{PCA}}\epsilon \perp x$, whence

$$\parallel P_{\mathrm{PCA}}{y-x\parallel }_2^2=\parallel P_{\mathrm{PCA}}{\epsilon \parallel }_2^2+{\parallel x\parallel }_2^2,$$

which proves (6).

Combining (5) and (6) gives (7). Condition (8) is therefore equivalent to strict improvement by the structural projection. Furthermore, since

$$\parallel P_S{\epsilon \parallel }_2^2\le {\parallel \epsilon \parallel }_2^2\mathrm{and}{\parallel P_{\mathrm{PCA}}\epsilon \parallel }_2^2\ge 0,$$

if ${\parallel \epsilon \parallel }_2<{\parallel x\parallel }_2$, then

$$\parallel P_S{y-x\parallel }_2^2\le {\parallel \epsilon \parallel }_2^2{<\parallel x\parallel }_2^2\le {\parallel P_{\mathrm{PCA}}y-x\parallel }_2^2.$$

When $\epsilon =0$, the two squared errors are respectively 0 and ${\parallel x\parallel }_2^2$, which completes the proof. □

Remark 1.

The assumptions (a) and (b) have a straightforward interpretation. The nuisance directions dominate second-order terms, so PCA selects them. The structural directions dominate the fourth-order terms, so successive fourth-order maximization chooses them. For one structural direction with probability p and amplitude a, and one nuisance direction with probability q and amplitude b, it is possible to have $p{a}^2<q{b}^2$ while $p{a}^4>q{b}^4.$ Equivalently,

$${\left({\displaystyle \frac{q}{p}}\right)}^{1/4}<{\displaystyle \frac{a}{b}}<{\left({\displaystyle \frac{q}{p}}\right)}^{1/2},$$

which is possible whenever $p<q$. This is the rare, high-amplitude regime in which fourth-order information can reveal structure that covariance missed.

In practice, after a prescribed baseline PCA subspace B has been removed, one can apply Theorem 1 to $B^{\perp }$. In that situation, both structural and nuisance directions in the theorem are residual directions, so Theorem 1 explains why fourth-order maximization may recover a residual structural component missed by an additional PCA direction.

4. Stable Residual Recovery Under Tensor Perturbations

The previous theorem treats a highly structured fourth-order tensor ${\mathcal{C}}_4$. However, in real-world applications, the residual fourth-order tensor will generally contain sampling error, cross terms, or small variations from the idealized model. In this section, we endeavor to show that residual Z-eigenvector recovery is stable under a bounded tensor perturbation.

Let $B\subseteq {\mathbb{R}}^n$ be a fixed s-dimensional PCA subspace, and let $P_B$ denote the orthogonal projector onto B. Let the perturbation tensor $\mathcal{E}$ be a fourth-order n-dimensional real symmetric tensor. We then define its semi-norm on $B^{\perp }$ by

$${\parallel \mathcal{E}\parallel }_{B^{\perp }}:=max\left\{\right|\mathcal{E}{u}^4|:u\in B^{\perp }{,\parallel u\parallel }_2=1\}.$$

Suppose $\mathbf{v},w_1,\dots ,w_r\in B^{\perp }$ are orthonormal. The vector $\mathbf{v}$ is regarded as the residual structural direction, while $w_1,\dots ,w_r$ are residual nuisance directions. Consider the unperturbed residual tensor

$${\mathcal{C}}_0:=\gamma {\mathbf{v}}^{\otimes 4}+{\displaystyle \sum _{\ell =1}^r}{\eta }_{\ell }{w}_{\ell }^{\otimes 4},\gamma >0,{\eta }_{\ell }\ge 0.$$

We set

$${\eta }_{max}:=\underset{1\le \ell \le r}{max}{\eta }_{\ell }\mathrm{and}\Delta :=\gamma -{\eta }_{max}.$$

The quantity $\Delta$ measures the fourth-order spectral gap between the structural direction and the strongest residual nuisance direction.

Theorem 2

(Stability of residual Z-eigenvector recovery). Assume $\Delta >0$ and the observed residual tensor has the form ${\mathcal{C}}_{\mathrm{res}}={\mathcal{C}}_0+\mathcal{E},$ where the perturbation tensor $\mathcal{E}$ satisfies ${\parallel \mathcal{E}\parallel }_{B^{\perp }}=\rho \ge 0$. Let $\widehat{\mathbf{v}}\in argmax\{{\mathcal{C}}_{\mathrm{res}}{u}^4:u\in B^{\perp }{,\parallel u\parallel }_2=1\}.$ Then

$$1-|\langle \widehat{\mathbf{v}},\mathbf{v}\rangle {|}^2\le min\left\{1,{\displaystyle \frac{2\rho }{\Delta }}\right\}.$$

Equivalently, if θ is the acute angle between $\widehat{\mathbf{v}}$ and $\mathbf{v}$, then

$${sin}^2\theta \le min\left\{1,{\displaystyle \frac{2\rho }{\Delta }}\right\}.$$

In particular, if $\rho <\Delta /2$, then the recovered direction has nonzero alignment with $\mathbf{v}$.

Proof. 

We write

$$\widehat{\mathbf{v}}=\alpha \mathbf{v}+{\displaystyle \sum _{\ell =1}^r}{\beta }_{\ell }{w}_{\ell }+h,$$

where $h\in B^{\perp }\cap span{\{\mathbf{v},w_1,\dots ,w_r\}}^{\perp }.$ It is straightforward to calculate

$$\delta :=1-|\langle \widehat{\mathbf{v}},\mathbf{v}\rangle {|}^2=1-{cos}^2\theta ={sin}^2\theta .$$

Since $\parallel \widehat{\mathbf{v}}{\parallel }_2=1$, we have ${\alpha }^2=1-\delta$ and ${\sum }_{\ell =1}^r{\beta }_{\ell }^2\le \delta .$ Consequently,

$${\displaystyle \sum _{\ell =1}^r}{\beta }_{\ell }^4\le {\left({\displaystyle \sum _{\ell =1}^r}{\beta }_{\ell }^2\right)}^2\le {\delta }^2.$$

Hence,

$${\mathcal{C}}_0{\widehat{\mathbf{v}}}^4=\gamma {\alpha }^4+{\displaystyle \sum _{\ell =1}^r}{\eta }_{\ell }{\beta }_{\ell }^4\le \gamma {(1-\delta )}^2+{\eta }_{max}{\delta }^2.$$

(11)

Since $\widehat{\mathbf{v}}$ is a global maximizer of the perturbed quartic form on the unit sphere in $B^{\perp }$, it is clear ${\mathcal{C}}_{\mathrm{res}}{\widehat{\mathbf{v}}}^4\ge {\mathcal{C}}_{\mathrm{res}}{\mathbf{v}}^4.$

On one hand, the perturbation bound gives

$${\mathcal{C}}_{\mathrm{res}}{\mathbf{v}}^4=\gamma +\mathcal{E}{\mathbf{v}}^4\ge \gamma -\rho .$$

On the other hand, (11) gives ${\mathcal{C}}_{\mathrm{res}}{\widehat{\mathbf{v}}}^4\le \gamma {(1-\delta )}^2+{\eta }_{max}{\delta }^2+\rho .$ Thus,

$$\gamma -\rho \le \gamma {(1-\delta )}^2+{\eta }_{max}{\delta }^2+\rho .$$

As a consequence,

$$\delta \left(2\gamma -(\gamma +{\eta }_{max})\delta \right)\le 2\rho .$$

(12)

Since $0\le \delta \le 1$,

$$2\gamma -(\gamma +{\eta }_{max})\delta \ge 2\gamma -(\gamma +{\eta }_{max})=\gamma -{\eta }_{max}=\Delta .$$

It follows from (12) that $\delta \Delta \le 2\rho ,$ and hence, $\delta \le {\displaystyle \frac{2\rho }{\Delta }}.$ By combining this with the trivial bound $\delta \le 1$, the assertion follows. □

Remark 2.

In the special case $\rho =0$, it implies $|\langle \widehat{\mathbf{v}},\mathbf{v}\rangle |=1$, that is, $\widehat{\mathbf{v}}=\pm \mathbf{v}$, hence, the above theorem provides exact residual recovery. More generally, the direction error is controlled by the ratio of the perturbation tensor’s semi-norm to the fourth-order spectral gap. Thus, bounded perturbations do not bear major impact on the recovery when the fourth-order spectral gap remains sufficiently large.

We next establish some useful reconstruction-error bounds.

Corollary 1

(Controlled reconstruction-error estimate).

Assume the hypotheses of Theorem 2. Let $x=x_B+\tau \mathbf{v},$ $x_B\in B,$ and let $y=x+\xi ,$ where $\xi \in {\mathbb{R}}^n$ is an arbitrary error vector. Let

$$Q_{\mathrm{rec}}:=P_B+{\widehat{\mathbf{v}}\widehat{\mathbf{v}}}^{\top }\mathrm{and}{\widehat{y}}_{\mathrm{rec}}:=Q_{\mathrm{rec}}y.$$

Then the following hold:

1.

$$\parallel {\widehat{y}}_{\mathrm{rec}}{-x\parallel }_2^2={\tau }^2\left(1-|\langle \widehat{v},v\rangle {|}^2\right)+{\parallel Q_{\mathrm{rec}}\xi \parallel }_2^2.$$

(13)

2.

$$\parallel {\widehat{y}}_{\mathrm{rec}}{-x\parallel }_2^2\le {\tau }^{2}min\left\{1,{\displaystyle \frac{2\rho }{\Delta }}\right\}+{\parallel Q_{\mathrm{rec}}\xi \parallel }_2^2.$$

(14)

Suppose additionally that an equal-rank PCA projector has the form

$$P_{\mathrm{PCA}}=P_B+w{w}^{\top },$$

where

$$w\in B^{\perp },{\parallel w\parallel }_2=1,w\perp \mathbf{v}.$$

Then

$$\parallel P_{\mathrm{PCA}}{y-x\parallel }_2^2={\tau }^2+{\parallel P_{\mathrm{PCA}}\xi \parallel }_2^2.$$

(15)

Moreover,

$$\parallel P_{\mathrm{PCA}}{y-x\parallel }_2^2-\parallel {\widehat{y}}_{\mathrm{rec}}{-x\parallel }_2^2={\tau }^2|\langle \widehat{\mathbf{v}},\mathbf{v}\rangle {|}^2+|\langle w,\xi \rangle {|}^2-{\left|\langle \widehat{\mathbf{v}},\xi \rangle \right|}^2.$$

(16)

In particular, if $\rho <\Delta /2$, then

$$\parallel P_{\mathrm{PCA}}{y-x\parallel }_2^2-\parallel {\widehat{y}}_{\mathrm{rec}}{-x\parallel }_2^2\ge \left(1-{\displaystyle \frac{2\rho }{\Delta }}\right){\tau }^2+{\left|\langle w,\xi \rangle \right|}^2-{\left|\langle \widehat{\mathbf{v}},\xi \rangle \right|}^2.$$

(17)

Thus, a sufficient condition for residual recovery to have smaller squared error than equal-rank PCA is given by

$$|\langle \widehat{\mathbf{v}},\xi \rangle {|}^2-{\left|\langle w,\xi \rangle \right|}^2<\left(1-{\displaystyle \frac{2\rho }{\Delta }}\right){\tau }^2.$$

(18)

The residual-recovery estimator can also be compared with the unaugmented rank-s PCA estimate

$${\widehat{y}}_B:=P_{B}y.$$

In this case,

$$\parallel {\widehat{y}}_B{-x\parallel }_2^2={\tau }^2+{\parallel P_B\xi \parallel }_2^2,$$

(19)

and

$$\parallel {\widehat{y}}_B{-x\parallel }_2^2-\parallel {\widehat{y}}_{\mathrm{rec}}{-x\parallel }_2^2={\tau }^2|\langle \widehat{\mathbf{v}},\mathbf{v}\rangle {|}^2-{\left|\langle \widehat{\mathbf{v}},\xi \rangle \right|}^2.$$

(20)

Therefore, whenever $\rho <\Delta /2$, a sufficient condition for the augmented estimator to improve upon rank-s PCA is

$$|\langle \widehat{\mathbf{v}},\xi \rangle {|}^2<\left(1-{\displaystyle \frac{2\rho }{\Delta }}\right){\tau }^2.$$

(21)

Proof. 

Since $\widehat{\mathbf{v}}\in B^{\perp }$, $Q_{\mathrm{rec}}=P_B+{\widehat{\mathbf{v}}\widehat{\mathbf{v}}}^{\top }$ defines the orthogonal projector onto $B\oplus span\left\{\widehat{\mathbf{v}}\right\}.$ Since $x_B\in B$, we have:

$$Q_{\mathrm{rec}}x-x=Q_{\mathrm{rec}}(x_B+\tau \mathbf{v})-(x_B+\tau \mathbf{v})=\tau \left({\widehat{\mathbf{v}}\widehat{\mathbf{v}}}^{\top }\mathbf{v}-\mathbf{v}\right).$$

(22)

The vector in (22) is clearly orthogonal to B. It is also orthogonal to $\widehat{\mathbf{v}}$, since

$$\langle {\widehat{\mathbf{v}}\widehat{\mathbf{v}}}^{\top }\mathbf{v}-\mathbf{v},\widehat{\mathbf{v}}\rangle =\langle \widehat{\mathbf{v}},\mathbf{v}\rangle -\langle \mathbf{v},\widehat{\mathbf{v}}\rangle =0.$$

Hence,

$$Q_{\mathrm{rec}}x-x\in range{\left(Q_{\mathrm{rec}}\right)}^{\perp }.$$

On the other hand,

$$Q_{\mathrm{rec}}\xi \in range\left(Q_{\mathrm{rec}}\right).$$

The two vectors are therefore orthogonal. Since

$${\widehat{y}}_{\mathrm{rec}}-x=Q_{\mathrm{rec}}x-x+Q_{\mathrm{rec}}\xi ,$$

it is immediate that

$$\parallel {\widehat{y}}_{\mathrm{rec}}{-x\parallel }_2^2=\parallel Q_{\mathrm{rec}}{x-x\parallel }_2^2+{\parallel Q_{\mathrm{rec}}\xi \parallel }_2^2.$$

Furthermore,

$$\parallel Q_{\mathrm{rec}}{x-x\parallel }_2^2={\tau }^2{\parallel {\widehat{\mathbf{v}}\widehat{\mathbf{v}}}^{\top }\mathbf{v}-\mathbf{v}\parallel }_2^2={\tau }^2\left(1-|\langle \widehat{\mathbf{v}},\mathbf{v}\rangle {|}^2\right).$$

This proves (13). The estimate (14) follows directly from Theorem 2.

Now consider the equal-rank PCA projector $P_{\mathrm{PCA}}=P_B+w{w}^{\top }.$ Since $w\perp \mathbf{v}$ and $\mathbf{v}\in B^{\perp }$, we have $P_{\mathrm{PCA}}x=x_B.$ Thus, $P_{\mathrm{PCA}}x-x=-\tau \mathbf{v}.$ We note the vector $\mathbf{v}$ is orthogonal to the range of $P_{\mathrm{PCA}}$, while $P_{\mathrm{PCA}}\xi$ belongs to that range. Therefore, the Pythagorean Theorem yields

$$\parallel P_{\mathrm{PCA}}{y-x\parallel }_2^2={\tau }^2+{\parallel P_{\mathrm{PCA}}\xi \parallel }_2^2,$$

which proves (15).

Because $w,\widehat{\mathbf{v}}\in B^{\perp }$, we also have

$$\parallel P_{\mathrm{PCA}}{\xi \parallel }_2^2=\parallel P_B{\xi \parallel }_2^2+{\left|\langle w,\xi \rangle \right|}^2$$

and

$$\parallel Q_{\mathrm{rec}}{\xi \parallel }_2^2=\parallel P_B{\xi \parallel }_2^2+{\left|\langle \widehat{\mathbf{v}},\xi \rangle \right|}^2.$$

Subtracting (13) from (15) yields (16). If $\rho <\Delta /2$, Theorem 2 shows

$$|\langle \widehat{\mathbf{v}},\mathbf{v}\rangle {|}^2\ge 1-{\displaystyle \frac{2\rho }{\Delta }},$$

which proves (17) and the sufficient condition (18).

Lastly, since $\mathbf{v}\in B^{\perp }$, $P_{B}x-x=-\tau \mathbf{v},$ which is orthogonal to $P_B\xi$. Hence,

$$\parallel {\widehat{y}}_B{-x\parallel }_2^2={\tau }^2+{\parallel P_B\xi \parallel }_2^2.$$

Subtracting (13) gives (20). The sufficient condition (21) again follows from Theorem 2. □

5. A Three-Dimensional Example

We illustrate the residual recovery procedure with a finite data set in ${\mathbb{R}}^3$. Let $e_1,e_2,e_3$ denote the standard basis vectors in ${\mathbb{R}}^3$, and consider the eight centered data vectors

$$\begin{array}{cccc}\hfill z_1& =6e_1,\hfill & \hfill z_2& =-6e_1,\hfill \\ \hfill z_3& =4e_2,\hfill & \hfill z_4& =-4e_2,\hfill \\ \hfill z_5& =3e_3,\hfill & \hfill z_6& =-3e_3,\hfill \\ \hfill z_7& =3e_3,\hfill & \hfill z_8& =-3e_3.\hfill \end{array}$$

Their mean is zero. The relative frequencies and the corresponding second- and fourth-order coefficients are carefully chosen to make the calculations more friendly. We summarize the parameters below:

The covariance matrix is

$$\widehat{\Sigma }={\displaystyle \frac{1}{8}}{\displaystyle \sum _{q=1}^8}{z}_q{z}_q^{\top }=\left(\begin{array}{ccc}9& 0& 0\\ 0& 4& 0\\ 0& 0& {\displaystyle \frac{9}{2}}\end{array}\right).$$

Thus, rank-one PCA selects the subspace

$$B=span\left\{e_1\right\}\mathrm{and}{P}_B=e_1{e}_1^{\top }.$$

Since $9>{\displaystyle \frac{9}{2}}>4,$ rank-two PCA selects $span\{e_1,e_3\}.$ In particular, after the first PCA direction is removed, covariance prefers the nuisance direction $e_3$ to the structural direction $e_2$.

The residual vectors are

$$\begin{array}{cc}& r_1=0,r_2=0,r_3=4e_2,r_4=-4e_2,\hfill \\ & r_5=3e_3,r_6=-3e_3,r_7=3e_3,\mathrm{and}{r}_8=-3e_3.\hfill \end{array}$$

Consequently, the residual fourth-order tensor is

$${\mathcal{C}}_{\mathrm{res}}=64e_2^{\otimes 4}+{\displaystyle \frac{81}{2}}{e}_3^{\otimes 4}.$$

For a unit vector

$$u=u_2{e}_2+u_3{e}_3\in B^{\perp },u_2^2+u_3^2=1,$$

the residual quartic form is

$${\mathcal{C}}_{\mathrm{res}}{u}^4=64u_2^4+{\displaystyle \frac{81}{2}}{u}_3^4.$$

Since $81/2<64$, we have

$${\mathcal{C}}_{\mathrm{res}}{u}^4\le 64\left(u_2^4+u_3^4\right)\le 64{\left(u_2^2+u_3^2\right)}^2=64.$$

Equality holds precisely when $u=\pm e_2$. Hence, the residual fourth-order method selects $\widehat{v}=\pm e_2.$ We now compare the resulting projections. Let $y_0={(6,4,1)}^{\top }$ be the observed vector, and let $x_0={(6,4,0)}^{\top }$ be the corresponding projected vector. The third coordinate of $y_0$ is therefore regarded as a nuisance component.

The rank-one PCA estimate is

$${\widehat{y}}_{\mathrm{PCA},1}=e_1{e}_1^{\top }{y}_0={(6,0,0)}^{\top },$$

and its squared error is

$${∥{\widehat{y}}_{\mathrm{PCA},1}-x_0∥}_2^2=16.$$

Rank-two PCA selects $span\{e_1,e_3\}$ and therefore gives

$${\widehat{y}}_{\mathrm{PCA},2}=\left(e_1{e}_1^{\top }+e_3{e}_3^{\top }\right)y_0={(6,0,1)}^{\top }.$$

Its squared error is

$${∥{\widehat{y}}_{\mathrm{PCA},2}-x_0∥}_2^2=4^2+1^2=17.$$

Thus, adding the second PCA direction admits the nuisance coordinate without recovering the missing structural component.

By contrast, the PCA plus residual-recovery estimator uses the projector

$$P_{\mathrm{rec}}=e_1{e}_1^{\top }+e_2{e}_2^{\top }.$$

It follows that

$${\widehat{y}}_{\mathrm{rec}}=P_{\mathrm{rec}}{y}_0={(6,4,0)}^{\top }=x_0,$$

and hence

$${∥{\widehat{y}}_{\mathrm{rec}}-x_0∥}_2^2=0.$$

This example exhibits the separation analyzed in the preceding sections. The nuisance direction $e_3$ has a larger second-order coefficient than the structural direction $e_2$: ${\displaystyle \frac{9}{2}}>4,$ so residual PCA prefers $e_3$. On the other hand, $e_2$ has the larger fourth-order coefficient: $64>{\displaystyle \frac{81}{2}},$ so residual fourth-order maximization selects $e_2$.

In the notation of Theorem 2, we have

$$\mathbf{v}=e_2,w_1=e_3,\gamma =64,\mathrm{and}{\eta }_1={\displaystyle \frac{81}{2}}.$$

The fourth-order spectral gap is therefore

$$\Delta =\gamma -{\eta }_1=64-{\displaystyle \frac{81}{2}}={\displaystyle \frac{47}{2}}.$$

Because the example is unperturbed, $\rho =0$, and the perturbation Theorem 2 gives exact recovery:

$$\widehat{\mathbf{v}}=\pm \mathbf{v}.$$

Thus, even at the same output rank, the fourth-order residual method recovers the structural direction missed by PCA while excluding the variance-dominant nuisance direction.

6. Conclusion

We have introduced a residual dimension-reduction procedure that combines a prescribed PCA subspace with a direction obtained by maximizing a quartic form on its orthogonal complement. The analysis identifies an orthogonal setting in which the quadratic and quartic criteria order the relevant directions differently. In this setting, PCA selects directions with the largest quadratic coefficients, whereas successive quartic maximization recovers structural directions with smaller quadratic coefficients but larger fourth-order tensor coefficients.

Under this model, the proposed procedure recovers a structural direction that may be removed by an equal-rank PCA projection. The improvement is therefore not merely a consequence of enlarging the projection subspace. For an arbitrary perturbation vector, we derived exact bounds for the squared reconstruction errors and explicit conditions under which the residual fourth-order projection improves upon both equal-rank PCA and the original unaugmented PCA projection.

Our residual recovery approach may be useful to supplement PCA- or SVD-based filtering in imaging and spatiotemporal signal processing. In particular, a dominant low-rank component may represent smooth background or clutter, while a weak feature of interest may hide in the residual.