Two-by-Two Matrix Decomposition via Jacobi Rotation

Hard

The Singular Value Decomposition (SVD) is one of the most powerful and versatile tools in numerical linear algebra, enabling us to decompose any matrix into a product of three special matrices that reveal its fundamental geometric properties. For a real 2×2 matrix A, the SVD expresses it as: $$A = U \cdot \Sigma \cdot V^T$$ Where: U is a 2×2 orthogonal matrix whose columns are the left singular vectors Σ (Sigma) is a diagonal matrix containing the singular values (σ₁ and σ₂) V^T is the transpose of the 2×2 orthogonal matrix V, whose columns are the right singular vectors The Jacobi Rotation Approach: One elegant method to compute the SVD of a 2×2 matrix involves a Jacobi rotation—a planar rotation designed to diagonalize symmetric matrices. The algorithm works as follows: Form the symmetric product: Compute B = AᵀA (this is always symmetric and positive semi-definite) Compute the Jacobi rotation angle: For the symmetric matrix B with elements b₁₁, b₁₂, b₂₂, calculate the rotation angle θ that diagonalizes B: If b₁₂ = 0, the matrix is already diagonal (θ = 0) Otherwise: θ = 0.5 × arctan(2 × b₁₂ / (b₁₁ - b₂₂)) Construct the rotation matrix V: This is an orthogonal matrix built from cos(θ) and sin(θ) Compute singular values: The eigenvalues of B are the squared singular values of A Derive U: Calculate U = A × V × Σ⁻¹ (where Σ⁻¹ contains the reciprocals of singular values) This single-rotation approach provides an approximate SVD that is exact for symmetric matrices and offers a good approximation for general 2×2 matrices. Your Task: Implement a Python function that computes the approximate SVD of a real 2×2 matrix using a single Jacobi rotation step. You must not use any high-level SVD library functions (like numpy.linalg.svd). Instead, build the decomposition from fundamental matrix operations.

Example 1

Input

A = np.array([[2.0, 1.0], [1.0, 2.0]])

Output

U ≈ [[0.7071, -0.7071], [0.7071, 0.7071]]
S = [3.0, 1.0]
Vt ≈ [[0.7071, 0.7071], [-0.7071, 0.7071]]

Explanation

This is a symmetric matrix with eigenvalues 3 and 1. For symmetric matrices, the SVD coincides with the eigendecomposition—the singular values equal the absolute values of the eigenvalues, and the singular vectors are the eigenvectors. The decomposition satisfies A = U × diag(S) × Vt: [[0.7071, -0.7071], [0.7071, 0.7071]] × [[3, 0], [0, 1]] × [[0.7071, 0.7071], [-0.7071, 0.7071]] = [[2, 1], [1, 2]]

Example 2

Input

A = np.array([[3.0, 0.0], [0.0, 2.0]])

Output

U = [[1.0, 0.0], [0.0, 1.0]]
S = [3.0, 2.0]
Vt = [[1.0, 0.0], [-0.0, 1.0]]

Explanation

For a diagonal matrix, the SVD is trivial—the matrix is already in 'decomposed' form. The singular values are simply the diagonal entries (3 and 2), and both U and V are identity matrices (or have sign adjustments). No rotation is needed since the matrix is already aligned with the coordinate axes.

Example 3

Input

A = np.array([[1.0, 2.0], [3.0, 4.0]])

Output

U ≈ [[0.3826, 0.9239], [0.9239, -0.3826]]
S ≈ [5.465, 0.366]
Vt ≈ [[0.8174, 0.576], [-0.576, 0.8174]]

Explanation

For this non-symmetric matrix, the Jacobi rotation method computes the SVD by first analyzing AᵀA. The larger singular value (≈5.465) represents the maximum 'stretching' the matrix applies to any unit vector, while the smaller singular value (≈0.366) represents the minimum stretching. The ratio σ₁/σ₂ ≈ 14.9 indicates this matrix is moderately ill-conditioned.

Constraints

Input A is always a 2×2 NumPy array with real (float) values
-10⁶ ≤ A[i][j] ≤ 10⁶
You may NOT use numpy.linalg.svd or any other high-level SVD routine
You may use basic NumPy operations: matrix multiplication (@), transpose (.T), element-wise math, np.sqrt, np.arctan, np.cos, np.sin, etc.
Use only a single Jacobi rotation step (no iterative refinement)
The returned singular values should be non-negative (use absolute values if needed)
Round output values to 4 decimal places for consistency