Extending eigenvector orientation from π to 2π angular intervals via a modified arctan2 algorithm, with applications to random matrix theory and eigenvector stabilization.
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.
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.
The central equation for consistent orientation of an eigenvector matrix V is:
where R, S, and I are rotation, reflection, and identity matrices respectively. The reflection matrix is defined as S = \prod S_k where:
The VS product produces the orientated eigenvector matrix:
The algorithm visits N - 1 subspaces, and the pattern that develops is:
leading to the expansion:
The partially oriented eigenvector matrix \mathcal{V}_k is defined as:
Writing V in terms of a working matrix W:
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.
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:
Multiplying through the Givens cascade and equating to the target yields the governing transcendental system:
where \{s_k | c_k\} = \{\sin|\cos\}(\theta_{1,k}) and a_i are symbolic substitutes for w_{4,\cdot,1}.
The original solution takes entries pairwise from the bottom up:
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.
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.
This section revisits equation (1) by considering the handedness of an orthogonal basis holistically—in contrast with the "opportunistic" approach of the original article.
Eigenvector matrix V and rotation matrix R are both unitary:
But they differ in determinant:
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):
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.
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).
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).
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).
The modified arctan2 solution to (10), for four dimensions, is:
To avoid overflow, this is rewritten as:
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.
When the column vector \mathbf{a} is sparse, with a_{k} = 0, the indexing requires adjustment to skip the zero entry:
where a_{k-j} points to the first nonzero entry behind a_k. When a_2 = 0:
The central equation remains R_{\tan}^T V S_{\tan} = I in structure, but now:
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.
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.
The Python package thucyd is freely available from PyPI and Conda-Forge (version 0.2.5+). The interface exposes two functions:
orient_eigenvectorsarcsin and arctan2 methods. The method keyword defaults to arcsin for backward compatibility but arctan2 is now preferred.generate_oriented_eigenvectorsWhen 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:
With s_1 = -1, rotation R_{1,2}^T is always minor.
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.
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:
| Pair | Description |
|---|---|
| EURUSD | Euro / US Dollar |
| USDJPY | US Dollar / Japanese Yen |
| GBPUSD | British Pound / US Dollar |
| USDCHF | US Dollar / Swiss Franc |
| AUDUSD | Australian Dollar / US Dollar |
| NZDUSD | New Zealand Dollar / US Dollar |
| CADUSD | Canadian Dollar / US Dollar |
Panels analyzed via SVD:
The eigendecomposition of the correlation matrix is:
For N = 7, the embedded angles are organized into an upper triangular matrix:
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.
The inverse participation ratio (IPR) quantifies how localized an eigenvector is:
To score higher participation above lower, the participation score is:
When all elements participate equally, \text{PS} = 1.
The Marčenko–Pastur (MP) distribution for eigenvalues of a Gaussian random matrix is:
where q = N/T > 0. Eigenvalues within the MP distribution are indistinguishable from sample noise. The rescaled distribution accounts for informative ("spike") eigenvalues:
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.
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.
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:
followed by renormalization:
Then re-orthogonalize via the orientation algorithm:
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).
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.
For noninformative modes, the angles are simply set to zero. For quotes (3 informative modes):
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}.
A correlation matrix is calculated as:
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.
Ledoit–Wolf shrinkage blends an empirical matrix with an identity matrix:
However, this estimator applies to matrices with no informative modes. Informative modes cannot simply be dropped; but they can be rotated away:
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.
| Use Case | Preferred Method | Reason |
|---|---|---|
| Regression on eigenbasis | arcsin | Sign flips in V propagate to β; stability of signs matters |
| Directional statistics | Modified arctan2 | Full 2π disambiguation of pointing directions |
| Informative vs. noise mode detection | Modified arctan2 | Scatter patterns distinguish directed from random modes |
| Correlation matrix cleaning | Modified arctan2 + stabilization | Combines directional analysis with filtering & zeroing |