A Consistently Oriented Basis for Eigenanalysis: Improved Directional Statistics

Extending eigenvector orientation from π to 2π angular intervals via a modified arctan2 algorithm, with applications to random matrix theory and eigenvector stabilization.

Jay Damask · arXiv:2402.08139 · February 2024 · thucyd Python package

§1  Overview

A model-free algorithm to consistently orient the eigenvectors returned from an svd or eig software call was detailed in an earlier article, and the Python package thucyd is now used in industry. To consistently orient a matrix V of eigenvectors is to replace V with \mathcal{V} such that the latter is only a pure rotation R away from the identity matrix.

Problem Statement

Given an eigenvector matrix V from SVD/eigendecomposition, produce \mathcal{V} = VS that is only a pure rotation away from I. The original algorithm restricted angles to [-\pi/2, \pi/2]; this paper extends the primary rotation at each subspace to a full 2\pi interval.

The original algorithm sorts eigenvectors by their associated eigenvalues in descending order, then decides—one by one, top to bottom—whether the sign of an eigenvector needs to be flipped. In the context of running linear regressions on the eigenbasis \mathcal{V} over an evolving data stream, the signs of regression weights \boldsymbol{\beta} are stabilized by regressing on \mathcal{V} rather than V.

The key insight of this article is that only the first rotation upon entry into a new subspace of V requires a 2\pi rotation range, while the remaining angles in the subspace are inherently bound to a \pi interval. This leads to a modified arctan2 algorithm that avoids the edge cases that break a naive arctan2 approach.

Tradeoff: Sign flips in the original algorithm can be replaced by rotations on a 2\pi range (although not in a one-to-one manner). The revised algorithm is better suited to problems that need full angular disambiguation, while the original is better suited to regression where a sign flip in V has knock-on effects to \boldsymbol{\beta}.

§2  Review of the Original Algorithm

The central equation for consistent orientation of an eigenvector matrix V is:

(1) $$R^T V S = I$$

where R, S, and I are rotation, reflection, and identity matrices respectively. The reflection matrix is defined as S = \prod S_k where:

(2) $$S_k = \operatorname{diag}(1, 1, \ldots, s_k, \ldots, 1), \quad s_k = \pm 1$$

The VS product produces the orientated eigenvector matrix:

(3) $$\mathcal{V} = VS = R$$

The algorithm visits N - 1 subspaces, and the pattern that develops is:

(4) $$V \to VS_1 \to R_1^T V S_1 \to R_1^T V S_1 S_2 \to \cdots$$

leading to the expansion:

(5) $$R_n^T \cdots R_2^T R_1^T V S_1 S_2 \cdots S_n = I$$

The partially oriented eigenvector matrix \mathcal{V}_k is defined as:

(6) $$\mathcal{V}_k := R_k^T \cdots R_1^T V S_1 \cdots S_k$$
Key structural property: Reflections lie solely to the right of V whereas rotations lie solely to the left. Section 7.1 of the original paper showed that substituting Householder reflections H for rotations R leads to scattering of the carefully constructed reflections in S, nullifying interpretability.

Subspace Reduction Pattern

Writing V in terms of a working matrix W:

(7–8) $$V = \begin{pmatrix} W_{(n)} \end{pmatrix}, \quad \mathcal{V}_1 = \begin{pmatrix} 1 \\ & W_{(n-1)} \end{pmatrix}, \quad \mathcal{V}_2 = \begin{pmatrix} 1 \\ & 1 \\ && W_{(n-2)} \end{pmatrix}$$

Each rotation R_k^T leads to partial diagonalization of subspace W_{(n-k-1)} such that a 1 appears on the kth diagonal entry.

System of Transcendental Equations

Rotations R_k^T are composed from a cascade of n - k - 1 Givens rotations, each operating within a two-dimensional plane. For R_1 \in \mathbb{R}^4:

$$R_1(\theta_{1,2}, \theta_{1,3}, \theta_{1,4}) = R_{1,2}(\theta_{1,2})\; R_{1,3}(\theta_{1,3})\; R_{1,4}(\theta_{1,4})$$

Multiplying through the Givens cascade and equating to the target yields the governing transcendental system:

(10) $$\begin{pmatrix} c_2\, c_3\, c_4 \\ s_2\, c_3\, c_4 \\ s_3\, c_4 \\ s_4 \end{pmatrix} = \begin{pmatrix} a_1 \\ a_2 \\ a_3 \\ a_4 \end{pmatrix}$$

where \{s_k | c_k\} = \{\sin|\cos\}(\theta_{1,k}) and a_i are symbolic substitutes for w_{4,\cdot,1}.

The arcsin Solution

The original solution takes entries pairwise from the bottom up:

(11) $$\begin{aligned} \theta_{1,4} &= \arcsin(a_4) \\ \theta_{1,3} &= \arcsin(a_3 / \cos\theta_{1,4}) \\ \theta_{1,2} &= \arcsin\!\bigl(a_2 / (\cos\theta_{1,4}\cos\theta_{1,3})\bigr) \end{aligned}$$

In total, there are N(N-1)/2 angles embedded in V \in \mathbb{R}^N. The issue is that all angles are minor, restricted to \theta \in [-\pi/2, \pi/2], meaning apparent pointing directions are limited to a hyper-hemisphere rather than a full hypersphere.

The Wrap-Around Problem

In the presence of noise in an evolving system, vector wobble can lead to angular wrap-around such that \mathbf{v}(t_2) and \mathbf{v}(t_1) appear on opposite sides of the hyper-hemisphere even when these vectors are in proximity on a hypersphere. The original article proposed an arctan2 method to span 2\pi, but edge cases were shown to break the algorithm.

§3  Handedness in High Dimensions

This section revisits equation (1) by considering the handedness of an orthogonal basis holistically—in contrast with the "opportunistic" approach of the original article.

Group-Theoretic Foundation

Eigenvector matrix V and rotation matrix R are both unitary:

(12) $$V^T V = I \quad \text{and} \quad R^T R = I$$

But they differ in determinant:

(13) $$\det(V) = \pm 1 \quad \text{and} \quad \det(R) = 1$$
Key Observation

V belongs to the orthogonal group O(N), while R belongs to the special orthogonal group SO(N) \subset O(N). The problem is one of squeezing V into SO(N).

Taking determinants of (1):

(14) $$\det(V)\det(S) = \det(R) = 1$$

In light of this, no reflections are inherently required when \det(V) = 1, and only a single reflection is necessary otherwise. The original algorithm's S = \operatorname{diag}(\pm 1, \pm 1, \ldots, \pm 1) often introduces more reflections than necessary.

Replacing Reflection with Major-Angle Rotation

The key forward step is recognizing that a reflection in a reducible subspace can be replaced by a rotation through a major angle (an angle exceeding \pi/2).

(a) π₁ π₂ v₁ v₂ (b) π₁ π₂ v₁∥π₁ v₂ reflect v₂ (c) v₂ (d) v₂∥π₂

Original method: v₂ is reflected in step (b)→(c) to point toward π₂, then rotated to align. The reflection creates s_2 = -1. Resulting \mathcal{S} = (1,-1,1).

(a') π₁ π₂ v₁ v₂ (b') v₂ major-angle rotation (c') v₃↓ irreducible: must reflect (d') v₂∥π₂

New method: v₂ is rotated through a major angle in step (b') rather than reflected. The reflection is deferred to the irreducible subspace (c'→d'). Resulting \mathcal{S}_{\tan} = (1,1,-1).

Figure 1. Two methods to orient a left-handed basis in ℝ³ onto the right-handed identity matrix I. Top: arcsin method with opportunistic reflection. Bottom: modified arctan2 method deferring reflection to the irreducible subspace.

The Modified arctan2 Solution

The modified arctan2 solution to (10), for four dimensions, is:

(17) $$\begin{aligned} \theta_{1,2} &= \operatorname{arctan2}(a_2, a_1) \\ \theta_{1,3} &= \operatorname{arctan2}(a_3, |a_2 \csc\theta_2|) \\ \theta_{1,4} &= \operatorname{arctan2}(a_4, |a_3 \csc\theta_3|) \end{aligned}$$

To avoid overflow, this is rewritten as:

(18) $$\theta_{\cdot,k} = \operatorname{arctan2}\bigl(a_k\,|\sin\theta_{k-1}|,\; |a_{k-1}|\bigr)$$
Critical Difference from arcsin

Only the first angle is calculated from all four quadrants; subsequent angles enforce the non-negative nature of the vector as projected onto the \pi_1 axis. As a side effect, the catastrophic edge case from the original article is avoided because a signed zero cannot appear in the x argument.

Handling Sparse Vectors

When the column vector \mathbf{a} is sparse, with a_{k} = 0, the indexing requires adjustment to skip the zero entry:

(19) $$\theta_{\cdot,k} = \operatorname{arctan2}\bigl(a_k\,|\sin\theta_{k-j}|,\; |a_{k-j}|\bigr), \quad j \geq 1$$

where a_{k-j} points to the first nonzero entry behind a_k. When a_2 = 0:

(20) $$\theta_{\cdot,k} = \operatorname{arctan2}(a_k, |a_1|)$$

The Revised Reflection Matrix

The central equation remains R_{\tan}^T V S_{\tan} = I in structure, but now:

(21) $$S_{\tan} = \operatorname{diag}(1, 1, \ldots, \pm 1)$$

Only the last entry can be negative (i.e., only the irreducible subspace may need reflection). The angles \boldsymbol{\theta}_{\tan} can now contain major angles for the first rotation of each new subspace.

Relationship Between \mathcal{S}_{\sin} and \mathcal{S}_{\tan}

Untwisting the reflection vectors

A rotation through a major angle in subspace k flips the sign of both s_k and s_{k+1}. Starting from \mathcal{S}_{\sin}, pairs of adjacent negative entries can be "untwisted" via major-angle rotations.

Example 1: \mathcal{S}_{\sin} = (1,-1,1) \to \mathcal{S}_{\tan} = (1,1,-1) — major-angle rotation in subspace 2 flips entries 2 and 3.

Example 2: \mathcal{S}_{\sin} = (-1,-1,-1,1) \to (1,1,-1,1) \to \mathcal{S}_{\tan} = (1,1,1,-1) — two successive untwistings.

When \det(V) = +1 and \mathcal{S}_{\sin} has an even number of negative signs, no reflection is needed for the arctan2 method.

§4  Reference Implementation

The Python package thucyd is freely available from PyPI and Conda-Forge (version 0.2.5+). The interface exposes two functions:

orient_eigenvectors
Implements both arcsin and arctan2 methods. The method keyword defaults to arcsin for backward compatibility but arctan2 is now preferred.
generate_oriented_eigenvectors
Same as in the original article—reconstructs eigenvectors from angle matrices.
OrientToFirstOrthant Flag

When the first eigenvector is expected to point into the first orthant (all elements positive), but the sign from eig/svd is negative, setting this flag records a reflection at the top-level dimension (s_1 = -1). The general case becomes:

(23) $$S_{\tan} = \operatorname{diag}(\pm 1, 1, \ldots, 1, \pm 1)$$

With s_1 = -1, rotation R_{1,2}^T is always minor.

Pseudocode: Algorithm 1

def OrientEigenvectors(V, E, OrientToFirstOrthant): Vsort, Esort, SortIdx SortEigenvectors(V, E) Vwork copy(Vsort), N dim(V) if OrientToFirstOrthant then SignFlipVec[0] (Vwork[0,0] >= 0) ? 1 : −1 Vwork[:,0] Vwork[:,0] * SignFlipVec[0] for i = 0 to N−2 do Vwork, AnglesMtx[i,:] ReduceDimensionByOne(N, i, Vwork) SignFlipVec[−1] (Vwork[−1,−1] >= 0) ? 1 : −1 Vor Vsort · diag(SignFlipVec), Eor Esort return Vor, Eor, AnglesMtx, SignFlipVec, SortIdx ────────────────────────────────────────── def SolveRotationAnglesInSubDim(N, i, Vcol): j i + 1 AnglesWork zeros(0:N) AnglesWork[j] arctan2(Vcol[j], Vcol[j−1]) # full 2π for j = i+2 to N do AnglesWork[j] arctan2( Vcol[j] |sin(AnglesWork[j−1])|, # y |Vcol[j−1]| # x (always ≥ 0) ) return AnglesWork
Functions that differ from the original algorithm are: OrientEigenvectors (top-level logic for arctan2), SolveRotationAnglesInSubDim (uses arctan2 instead of arcsin), and the final sign flip logic. Consult the thucyd source for sparsity handling.

§5  Empirical Example & Random Matrix Theory

Data Description

Chicago Mercantile Exchange (CME) futures quote and trade data was obtained for seven FX pairs for May 2023 (23 business days). The FX pairs were:

PairDescription
EURUSDEuro / US Dollar
USDJPYUS Dollar / Japanese Yen
GBPUSDBritish Pound / US Dollar
USDCHFUS Dollar / Swiss Franc
AUDUSDAustralian Dollar / US Dollar
NZDUSDNew Zealand Dollar / US Dollar
CADUSDCanadian Dollar / US Dollar
Data Processing Pipeline
  • June 2023 contract (most liquid) was the focus
  • Daily datasets sliced to 00:00–16:00 Chicago time
  • Mid-price generated from best bid/offer quotes
  • Common panel created via time arbitration (14 columns: 7 quotes + 7 trades)
  • Causal linear filters applied to a common timescale
  • Log quotes → price change; signed trades → directionality
  • Downsampled to remove autocorrelation (500–1000 records/day)
  • Columns mapped through Student-t copula to zero-mean, unit-variance Gaussian

SVD and Eigendecomposition

Panels analyzed via SVD:

(24) $$P = U \Lambda^{1/2} V^T$$

The eigendecomposition of the correlation matrix is:

(25) $$\operatorname{Corr}(P) = V \Lambda V^T$$

The Angle Matrix Structure

For N = 7, the embedded angles are organized into an upper triangular matrix:

(27) $$\boldsymbol{\theta} = \begin{pmatrix} 0 & \theta_{1,2} & \theta_{1,3} & \theta_{1,4} & \theta_{1,5} & \theta_{1,6} & \theta_{1,7} \\ & 0 & \theta_{2,3} & \theta_{2,4} & \theta_{2,5} & \theta_{2,6} & \theta_{2,7} \\ && 0 & \theta_{3,4} & \theta_{3,5} & \theta_{3,6} & \theta_{3,7} \\ &&& 0 & \theta_{4,5} & \theta_{4,6} & \theta_{4,7} \\ &&&& 0 & \theta_{5,6} & \theta_{5,7} \\ &&&&& 0 & \theta_{6,7} \\ &&&&&& 0 \end{pmatrix}$$

The diagonal color stripes highlight modes 1 (top row) through 6 (bottom). In each mode row, only the first angle (green) spans the full 2\pi interval; remaining angles are minor.

Participation Score

The inverse participation ratio (IPR) quantifies how localized an eigenvector is:

(28) $$\text{IPR} = \sum_{i=1}^{N} v_i^4, \quad 1/N \leq \text{IPR} \leq 1$$

To score higher participation above lower, the participation score is:

(29) $$\text{PS} = \frac{1}{N \times \text{IPR}}, \quad 1/N \leq \text{PS} \leq 1$$

When all elements participate equally, \text{PS} = 1.

Eigenvector Angular Distributions

Figure 4. Hemispherical angles (arcsin method) for quote eigenvectors over 23 days. All angles are minor (\pi interval). Radii indicate participation scores. Wrap-around ambiguity is visible in modes 3–6.
Figure 5. Spherical angles (modified arctan2) for quote eigenvectors. First rotation angle spans the full circle (2\pi). Mode 1 is highly directed; modes 4–6 scatter randomly. The red circle highlights a CPI-release outlier.

Key Empirical Findings

arcsin → arctan2 Transition
  • Quote Mode 1: Highly directed into the first orthant with both methods.
  • Quote Modes 2–3: Somewhat directed, but the arctan2 method reveals a CPI-release outlier (May 10, 2023) that is invisible in the arcsin plots.
  • Quote Modes 4–6: The arctan2 method makes it evident that these modes scatter randomly across the full 2\pi, confirming they are noise-dominated.
  • Trade Mode 1: Clearly directed; trade modes 2–6 scatter uniformly.
CPI Outlier: The US Consumer Price Index update on May 10, 2023 caused an anomalous EURUSD response. Removal of this pair from the panel restored the outlying point to the north-facing cluster. The arctan2 method was required to reveal this outlier.

Random Matrix Theory Connection

The Marčenko–Pastur (MP) distribution for eigenvalues of a Gaussian random matrix is:

(26) $$\rho(\lambda) = \frac{\sqrt{(\lambda_+ - \lambda)(\lambda - \lambda_-)}}{2\pi q \lambda}, \quad \lambda_\pm = (1 \pm \sqrt{q})^2$$

where q = N/T > 0. Eigenvalues within the MP distribution are indistinguishable from sample noise. The rescaled distribution accounts for informative ("spike") eigenvalues:

(30–31) $$\bar{\lambda} = \text{Mean}(\lambda_2, \lambda_3, \ldots, \lambda_N), \quad \rho_{\bar\lambda}(\lambda) = \frac{1}{\bar\lambda}\,\rho\!\left(\lambda/\bar\lambda\right)$$

Eigenvalue Spectra vs. Marčenko–Pastur

Figure 7. Empirical eigenvalues and participation scores overlaid with averaged rescaled MP distributions. Quotes: modes 1–3 fall outside MP → informative. Trades: only mode 1 falls outside → the rest is noise. This matches the directional analysis from the angular plots.
Cross-Validation of Eigenvalue and Eigenvector Analysis

The eigenvalue analysis (which modes fall outside the MP distribution) provides a very good match to the conclusions from eigenvector pointing-direction analysis. However, the eigenvalue analysis involves a rescaling that can be somewhat self-fulfilling—adjusting q' = (N-1)/T and rescaling by \bar{\lambda}. The directional analysis introduces genuinely new information to cross-validate the traditional RMT approach.

§6  Eigenvector Stabilization

Informative modes are at least somewhat directed, so their average direction is meaningful. Noninformative modes scatter randomly, so their average can be replaced by a fixed direction—exchanging variance for bias. These are called dynamic and static stabilization, respectively.

Dynamic Stabilization

The angular inclination of a collection of wobbling vectors is correctly computed by stacking the vectors end-to-end and measuring the resultant vector (not by averaging angles directly). In general:

(32) $$M_h[n] = (h * \mathcal{V})[n]$$

followed by renormalization:

(33) $$M[n] = M_h[n] \, \|M_h[n]\|^{-1}$$

Then re-orthogonalize via the orientation algorithm:

(34–35) $$M[n] \longrightarrow \boldsymbol{\theta}_M[n] \longrightarrow \tilde{\mathcal{V}}[n]$$

Eigenvalues are also filtered: \bar{\Lambda}[n] = (h * \Lambda)[n]. A five-point causal filter h[n] = [1,2,3,2,1]/9 was applied (average delay: 2 samples/days).

Filtering Results

Filtered pointing directions show reduced scatter. Participation scores tighten and slightly increase. Quote mode 2 retains a curiously low participation score but vectors remain directed. Increasing filter length adds stability at the expense of delay.

Static Stabilization

For noninformative modes, the angles are simply set to zero. For quotes (3 informative modes):

(37) $$\boldsymbol{\theta}_{\text{modal}} = \begin{pmatrix} 0 & \theta_{1,2} & \theta_{1,3} & \theta_{1,4} & \theta_{1,5} & \theta_{1,6} & \theta_{1,7} \\ & 0 & \theta_{2,3} & \theta_{2,4} & \theta_{2,5} & \theta_{2,6} & \theta_{2,7} \\ && 0 & \theta_{3,4} & \theta_{3,5} & \theta_{3,6} & \theta_{3,7} \\ &&& 0 & 0 & 0 & 0 \\ &&&& 0 & 0 & 0 \\ &&&&& 0 & 0 \\ &&&&&& 0 \end{pmatrix}$$

Geometrically, \mathcal{V}_{\text{modal}} = R_1 R_2 R_3 I since R_{4,5,6} = I. This is not the same as PCA (which discards modes); all modes remain present in \mathcal{V}.

(38) $$\mathcal{V}_{\text{modal}} = R_1 R_2 R_3 R_4 R_5 R_6\, I = R_1 R_2 R_3\, I$$

Correlation Matrix Reconstruction

A correlation matrix is calculated as:

(39) $$\overline{\operatorname{Corr}(P)} = \operatorname{diag}(\mathbf{s}^{-1})\left(\tilde{\mathcal{V}}\bar{\Lambda}\tilde{\mathcal{V}}^T\right) \operatorname{diag}(\mathbf{s}^{-1})$$

where \mathbf{s}^2 = \operatorname{diag}\!\left(\tilde{\mathcal{V}}\bar{\Lambda}\tilde{\mathcal{V}}^T\right). The renormalization corrects for small deviations from unity along the diagonal introduced by the filtering.

Correlation Stabilization Progression

High structure
Original
V Λ Vᵀ
Less structure
Averaged directions
V̄ Λ̄ V̄ᵀ
Least structure
Informative modes only
V̄_modal Λ̄ V̄_modalᵀ
Figure 10. Quote and trade correlation matrices progress from original (high structure) through dynamically stabilized (less structure) to dynamically + statically stabilized (least structure, most persistent over time).
Figure 11. Pairwise correlation scatter. Blue: original. Gray: after directional averaging. Orange: after retaining only informative modes. The dispersion progressively decreases.

Connection to Ledoit–Wolf Shrinkage

Ledoit–Wolf shrinkage blends an empirical matrix with an identity matrix:

(40) $$\Sigma_{\text{shr}} = \alpha\hat{\Sigma} + (1 - \alpha)I, \quad \alpha \in [0,1]$$

However, this estimator applies to matrices with no informative modes. Informative modes cannot simply be dropped; but they can be rotated away:

(41) $$\hat{\Sigma}_R = R_3^T R_2^T R_1^T \mathcal{V} = \begin{pmatrix} 1 \\ & 1 \\ && 1 \\ &&& \cdots \\ &&& \cdots \\ &&& \cdots \\ &&& \cdots \end{pmatrix}$$

which is admissible into (40), after which the shrunk matrix is rotated back. Only the noise-related eigenvalues should be included when calculating \alpha. Tying Ledoit–Wolf shrinkage to static stabilization: the effective shrinkage weight for the latter is \alpha = 0.

§7  Conclusions

Summary of Contributions
  1. Modified arctan2 method: Can defer the rectification of reflections embedded in V until the last (irreducible) subspace. N-1 angles within V \in \mathbb{R}^{N \times N} span a 2\pi interval rather than \pi.
  2. Angular disambiguation: Avoids wrap-around on the \pi interval, as demonstrated empirically with FX data where an outlier was only detectable with the new method.
  3. RMT cross-validation: Directional analysis of eigenvectors introduces new information to cross-validate the traditional eigenvalue-based RMT approach to distinguishing informative from noisy modes.
  4. Eigenvector stabilization: Dynamic (temporal filtering) and static (zeroing noisy-mode angles) stabilization methods reduce eigenvector wobble and correlation matrix dispersion.
  5. Shrinkage connection: Static stabilization relates to Ledoit–Wolf shrinkage; informative modes can be rotated away before shrinkage is applied.
Method Selection Guide
Use CasePreferred MethodReason
Regression on eigenbasisarcsinSign flips in V propagate to β; stability of signs matters
Directional statisticsModified arctan2Full 2π disambiguation of pointing directions
Informative vs. noise mode detectionModified arctan2Scatter patterns distinguish directed from random modes
Correlation matrix cleaningModified arctan2 + stabilizationCombines directional analysis with filtering & zeroing
Handedness insight: What is learned is that handedness, which is either left- or right-handed in 3D, does not generalize well in higher dimensions. What can be said is whether \det(V) = +1 (pure rotation, even number of reflections) or \det(V) = -1 (odd number of embedded reflections plus rotation).