math.NASep 26, 2026

Fast Differentiable SVD on GPU via Polar Decomposition

Authors: Uliana Parkina, Askar Tsyganov, Sergei Kudriashov, Sergey Samsonov, Maxim Rakhuba

Organizations: HSE University

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

Abstract

We present a fully GPU-oriented SVD pipeline based on polar decomposition, motivated by iterative methods that rely solely on matrix multiplications, such as the Newton-Schulz iteration. We show that this approach enables up to a 2×2\times speedup compared to standard implementations. Furthermore, we derive a numerically stable backward pass for the polar decomposition and leverage it to obtain a fully differentiable SVD. Our methods are released as open-source implementations in both PyTorch and JAX: https://github.com/fallnlove/cans_svd.

Figures & tables

Figure 1: Comparison of different parts of the SVD computation by their share of the total runtime for a random matrix M∈Rn×nM\in\mathbb{R}^{n\times n} with κ2(M)=10\kappa_{2}(M)=10 on an NVIDIA B200 GPU. We report the shares of polar-factor computation (in blue), symmetric EVD (EIGH) (in pink), and other routines (in green) for different polar-decomposition-based SVD algorithms.
Figure 2: Comparison of SVD algorithms on random matrices. The total wall-clock time is averaged over 100100 trials and reported with the standard deviation. Note: The algorithms have different accuracy levels, see Tables 1 and 2 for a more comprehensive analysis. The experiments are conducted on an NVIDIA B200 GPU.
Method101010210^{2}10310^{3}10410^{4}10610^{6}
CANS SVD (fp32)4.9_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}4.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}4.3_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-6}4.4_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}5.3_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-6}
CANS SVD (tf32)9.7_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-4}8.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-4}8.5_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-4}8.7_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-4}1.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-3}
QDWH SVD3.4_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}3.3_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}3.7_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}5.7_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}6.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}
CUDA POLAR4.3_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}5.5_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}7.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{1}7.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{1}7.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{1}
CUDA QR1.8_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-5}2.7_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-5}1.5_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-5}1.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-5}1.1_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-5}
CUDA JACOBI9.7_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.4}\cdot 10^{-4}8.4_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.4}\cdot 10^{-4}6.8_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.8}\cdot 10^{-4}6.1_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.6}\cdot 10^{-4}5.9_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.7}\cdot 10^{-4}
Table 1: Comparison of the accuracy of SVD algorithms on a random square matrix M∈R4096×4096M\in\mathbb{R}^{4096\times 4096} for different condition numbers. The reported metric is the relative reconstruction error ∥M−UΣV⊤∥F/∥M∥F\|M-U\Sigma V^{\top}\|_{F}/\|M\|_{F} , with standard deviation. Results are averaged over 100100 sampled matrices. Red indicates cases with significantly large errors.
Method101010210^{2}10310^{3}10410^{4}10610^{6}
CANS SVD (fp32)3.7_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}3.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}3.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.9}\cdot 10^{-6}3.3_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}2.4_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}
CANS SVD (tf32)6.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-4}6.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-4}2.8_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}6.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-4}6.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-4}
QDWH SVD3.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}3.5_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}3.4_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}3.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}3.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}
CUDA POLAR4.1_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}4.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}8.9_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}8.1_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}7.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}
CUDA QR9.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.2}\cdot 10^{-6}1.4_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-5}1.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-5}1.4_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-5}9.7_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.4}\cdot 10^{-6}
CUDA JACOBI2.1_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-3}2.1_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-3}2.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-3}2.3_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-3}2.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-3}
Table 2: Comparison of the accuracy of SVD algorithms on a random square matrix M∈R4096×4096M\in\mathbb{R}^{4096\times 4096} for different condition numbers. The reported metric is the relative orthogonality error max⁡{∥U⊤U−I∥F/∥I∥F,∥V⊤V−I∥F/∥I∥F}\max\{\|U^{\top}U-I\|_{F}/\|I\|_{F},\|V^{\top}V-I\|_{F}/\|I\|_{F}\} , with standard deviation. Results are averaged over 100100 sampled matrices.
Method \ Spec.gap1110−110^{-1}10−210^{-2}10−310^{-3}10−410^{-4}10−510^{-5}10−610^{-6}
Analytical Townsend (2016)66.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.6}67.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm 1.0}100.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm 41.0}690.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm 580.0}5500.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm 4600.0}NaNNaN
Taylor Wang et al. (2021)49.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.3}49.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.3}49.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.3}49.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.3}49.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.3}49.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.3}49.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.3}
Tikhonov (ours)66.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.6}67.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm 1.0}100.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm 41.0}690.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm 580.0}5500.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm 4600.0}66.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.6}66.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.6}
Table 3: Comparison of the gradient norm of SVD backward pass algorithms on a random square matrix M∈R4096×4096M\in\mathbb{R}^{4096\times 4096} for different spectral gaps. The reported metric is the Frobenius norm of the gradient ∥M‾∥F\|\overline{M}\|_{F} . Results are averaged over 100100 sampled matrices.
RatioMethodPPLVal Losssec. / iter.
100%LLaMA-2-7B{\color[rgb]{0.5,0.5,0.5}{5.53}}{\color[rgb]{0.5,0.5,0.5}{6.25}}{\color[rgb]{0.5,0.5,0.5}{-}}
 Dobi-SVD6.957.1775.09
80%Our ( ε=0.1\varepsilon=0.1 )6.907.1955.62
 Our ( ε=1.0\varepsilon=1.0 )6.877.1155.62
 Our ( ε=10\varepsilon=10 )6.937.1755.62
 Dobi-SVD7.477.8375.09
Table 4: Metric values of various compression methods. Experiments were conducted using the LLaMA-2-7B model on the WikiText-2 dataset.
Appendix figures & tables13 assets

Supplementary material from the paper’s appendix.

Appendix

Figure 3: Comparison of computation time for matrix multiplications in tf32 precision (left plot) and QR decompositions (right plot) across different GPUs on random square matrices of various sizes. The total wall-clock time is averaged over 100100 trials.
Figure 4: Comparison illustrating how many matrix multiplications can be executed within the runtime of a single QR decomposition for square matrices of different sizes across different GPUs. In the left plot, matrix multiplications are performed in tf32 precision. In the right plot, matrix multiplications are performed in fp32 precision. The total wall-clock time is averaged over 100100 trials.
Figure 5: Comparison of SVD algorithms on a random low rank matrix M∈Rn×nM\in\mathbb{R}^{n\times n} with rank⁡(M)=⌊n/10⌋\operatorname{rank}(M)=\lfloor n/10\rfloor . Total wall-clock time was averaged over 100100 trials and reported with the standard deviation. Note: we do not report CUDA POLAR for this setting, since it has substantially large reconstruction error. And algorithms have different accuracy, see Tables 5 and 6 for a more comprehensive analysis. The experiments are conducted on an NVIDIA B200 GPU.
Method116464102410244095409540964096
CANS SVD (fp32)1.9_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}2.7_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-6}4.3_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}8.4_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.4}\cdot 10^{-6}4.9_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}
CANS SVD (tf32)5.3_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.2}\cdot 10^{-4}7.5_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.2}\cdot 10^{-4}1.1_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-3}9.8_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-4}9.7_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-4}
QDWH SVD1.7_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-6}1.9_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}2.7_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}3.5_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}3.4_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}
CUDA POLAR8.3_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm 22.7}\cdot 10^{0}7.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{1}7.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{1}7.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{1}4.3_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}
CUDA QR7.4_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.5}\cdot 10^{-7}1.3_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.2}\cdot 10^{-6}6.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.4}\cdot 10^{-6}1.8_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-5}1.8_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-5}
CUDA JACOBI3.5_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.6}\cdot 10^{-5}9.5_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.2}\cdot 10^{-5}5.1_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.4}\cdot 10^{-4}1.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.2}\cdot 10^{-4}9.7_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.4}\cdot 10^{-4}
Appendix
Table 5: Comparison of SVD algorithms accuracy on a random square matrix M∈R4096×4096M\in\mathbb{R}^{4096\times 4096} for different ranks. Reported relative reconstruction error ∥M−UΣV⊤∥F/∥M∥F\|M-U\Sigma V^{\top}\|_{F}/\|M\|_{F} , with standard deviation. Metrics were averaged over 100100 sampled matrices. Red indicates cases with significantly large error.
Method116464102410244095409540964096
CANS SVD (fp32)1.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}1.5_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}1.8_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}3.1_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm 1.2}\cdot 10^{-6}3.7_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}
CANS SVD (tf32)1.9_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.2}\cdot 10^{-6}1.9_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}2.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}3.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}6.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-4}
QDWH SVD1.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}1.5_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}1.8_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}3.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}3.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}
CUDA POLAR7.5_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm 1.4}\cdot 10^{-6}3.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}6.9_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}1.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-5}4.1_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}
CUDA QR7.9_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.2}\cdot 10^{-6}7.7_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.2}\cdot 10^{-6}6.4_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.2}\cdot 10^{-6}9.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.3}\cdot 10^{-6}9.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.2}\cdot 10^{-6}
CUDA JACOBI2.5_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.2}\cdot 10^{-6}3.8_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-4}1.7_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-3}2.1_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-3}2.1_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-3}
Appendix
Table 6: Comparison of SVD algorithms accuracy on a random square matrix M∈R4096×4096M\in\mathbb{R}^{4096\times 4096} for different ranks. Reported relative orthogonality max⁡{∥U⊤U−I∥F/∥I∥F,∥V⊤V−I∥F/∥I∥F}\max\{\|U^{\top}U-I\|_{F}/\|I\|_{F},\|V^{\top}V-I\|_{F}/\|I\|_{F}\} , with standard deviation. Metrics were averaged over 100100 sampled matrices.
Method101010210^{2}10310^{3}10410^{4}10610^{6}
CANS SVD (fp32)4.9_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}4.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}4.3_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-6}4.4_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}5.3_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-6}
Extended5.8_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}5.6_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-6}5.4_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.1}\cdot 10^{-6}1.1_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-5}9.0_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}
SciPy5.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}5.1_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}5.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}5.2_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}5.1_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}
Gram4.4_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-6}2.3_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-5}1.1_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-4}3.9_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-4}6.9_{\color[rgb]{0.5,0.5,0.5}\scriptscriptstyle\pm.0}\cdot 10^{-4}
Appendix
Table 7: Comparison of the accuracy of SVD algorithms on a random square matrix M∈R4096×4096M\in\mathbb{R}^{4096\times 4096} for different condition numbers. The reported metric is the relative reconstruction error ∥M−UΣV⊤∥F/∥M∥F\|M-U\Sigma V^{\top}\|_{F}/\|M\|_{F} , with standard deviation. Results are averaged over 100100 sampled matrices.
Figure 6: Comparison of SVD algorithms on a square random matrices M∈Rn×nM\in\mathbb{R}^{n\times n} with κ2(M)=10\kappa_{2}(M)=10 . Total wall-clock time was averaged over 100100 trials and reported with the standard deviation. The experiments are conducted on an NVIDIA B200 GPU.
Figure 7: Comparison of different methods for computing polar decomposition on a square random matrices M∈Rn×nM\in\mathbb{R}^{n\times n} with κ2(M)=106\kappa_{2}(M)=10^{6} . The reported metric is relative orthogonality error ∥W⊤W−I∥F/∥I∥F\|W^{\top}W-I\|_{F}/\|I\|_{F} .
Figure 8: Pareto front for different CUDA JACOBI tolerance settings on square random matrices M∈R4096×4096M\in\mathbb{R}^{4096\times 4096} with κ2(M)=102\kappa_{2}(M)=10^{2} .
κ2(M)↓ δ→\kappa_{2}(M)\downarrow\ \delta\rightarrow0.90.90.990.990.9990.9990.99990.9999
10210^{2}3.30⋅10−63.30\cdot 10^{-6}4.26⋅10−64.26\cdot 10^{-6}7.25⋅10−67.25\cdot 10^{-6}1.84⋅10−51.84\cdot 10^{-5}
10410^{4}3.77⋅10−63.77\cdot 10^{-6}4.41⋅10−64.41\cdot 10^{-6}1.05⋅10−51.05\cdot 10^{-5}2.51⋅10−52.51\cdot 10^{-5}
10610^{6}4.00⋅10−64.00\cdot 10^{-6}5.38⋅10−65.38\cdot 10^{-6}1.69⋅10−51.69\cdot 10^{-5}4.88⋅10−54.88\cdot 10^{-5}
Appendix
Table 8: Ablation on different δ\delta for preprocessing on a random matrix M∈R4096×4096M\in\mathbb{R}^{4096\times 4096} for different condition numbers. The reported metric is the relative reconstruction error of CANS SVD (fp32).
κ2(M)↓ δ→\kappa_{2}(M)\downarrow\ \delta\rightarrow0.90.90.990.990.9990.9990.99990.9999
10107.20⋅10−47.20\cdot 10^{-4}9.73⋅10−49.73\cdot 10^{-4}2.98⋅10−32.98\cdot 10^{-3}9.15⋅10−39.15\cdot 10^{-3}
10210^{2}6.57⋅10−46.57\cdot 10^{-4}8.24⋅10−48.24\cdot 10^{-4}1.84⋅10−31.84\cdot 10^{-3}4.33⋅10−34.33\cdot 10^{-3}
10310^{3}7.04⋅10−47.04\cdot 10^{-4}8.52⋅10−48.52\cdot 10^{-4}1.91⋅10−31.91\cdot 10^{-3}3.63⋅10−33.63\cdot 10^{-3}
Appendix
Table 9: Ablation on different δ\delta for preprocessing on a random matrix M∈R4096×4096M\in\mathbb{R}^{4096\times 4096} for different condition numbers. The reported metric is the relative reconstruction error of CANS SVD (tf32).
Figure 9: Performance comparison of different precisions (tf32 and bf16) for matrix multiplication. The left plot shows the runtime of the matrix multiplications themselves. The right plot shows how the SVD runtime changes when the multiplication precision is varied. The total wall-clock time is averaged over 100100 trials. Experiments are conducted on an NVIDIA B200 GPU.
Figure 10: Comparison of different methods for solving the symmetric Lyapunov equation under different scenarios. The total wall-clock time was averaged over 100100 trials and reported with the standard deviation. The experiments were conducted on an NVIDIA B200 GPU.

Explore similar work

CardsList
Jul 24, 2026

Singular value soft-thresholding via the polar decomposition

Singular value soft-thresholding can be computed via a reduction to the matrix polar decomposition, which allows one to exploit GPU-friendly algorithms for computing the polar decomposition. Empirically, there is a significant speed-up on GPUs compared to the standard approach using the SVD. We leave the investigation of robustness to future work, but note that due to the discontinuous nature of the sign function, the reduction to the polar decomposition is likely only suitable for low-accuracy applications.

Stephen Becker

University of Colorado Boulder, USA

Topic pending
AskPDF
Jul 28, 2026

Quasi-SVD: Learning a Lie-constrained matrix factorisation for real-time imaging

Singular Value Decomposition (SVD) underlies matrix factorisation tasks across many fields, with imaging applications demanding real-time processing. Yet SVD algorithms are inherently sequential, constraining real-time GPU throughput and limit online deployment in imaging pipelines. This study introduces a fully parallelized matrix factorization framework for GPUs by enforcing matrix orthogonality on left singular vectors via Lie-parametrised algebra and recovering the remaining components through soft constraints. This asymmetric constraint design enables an efficient parallel and provably valid decomposition, achieves high reconstruction fidelity and substantially accelerates computation relative to the exact SVD, with real-time throughput exceeding standard video frame rates. Performance is evaluated on multiple imaging tasks spanning complementary computational regimes: (1) spatio-temporal background subtraction for ultrasound localisation microscopy, requiring high-dimensional matrix separation, (2) Mueller matrix polarimetry for neurosurgical tissue characterisation, requiring massive batch processing of small matrices, and (3) an MNIST denoising benchmark at an intermediate scale with known ground truth. Across regimes and instruments, the proposed framework demonstrates robust domain transfer at various matrix scales, sufficient for live image-guided workflows that classical solvers cannot currently support in these settings. By prioritising downstream reconstruction fidelity over exact spectral recovery, the proposed SVD framework makes structured matrix factorisation practical for real-time processing.

Christopher Hahne

University of Bern 3008 Bern, Switzerland

Lie Group MethodsComputational ImagingSingular Value DecompositionMatrix Factorization
AskPDF
Apr 7, 2026

Training Without Orthogonalization, Inference With SVD: A Gradient Analysis of Rotation Representations

Recent work has shown that removing orthogonalization during training and applying it only at inference improves rotation estimation in deep learning, with empirical evidence favoring 9D representations with SVD projection. However, the theoretical understanding of why SVD orthogonalization specifically harms training, and why it should be preferred over Gram-Schmidt at inference, remains incomplete. We provide a detailed gradient analysis of SVD orthogonalization specialized to 3×33 \times 3 matrices and SO(3)SO(3) projection. Our central result derives the exact spectrum of the SVD backward pass Jacobian: it has rank 33 (matching the dimension of SO(3)SO(3)) with nonzero singular values 2/(si+sj)2/(s_i + s_j) and condition number κ=(s1+s2)/(s2+s3)κ= (s_1 + s_2)/(s_2 + s_3), creating quantifiable gradient distortion that is most severe when the predicted matrix is far from SO(3)SO(3) (e.g., early in training when s3≈0s_3 \approx 0). We further show that even stabilized SVD gradients introduce gradient direction error, whereas removing SVD from the training loop avoids this tradeoff entirely. We also prove that the 6D Gram-Schmidt Jacobian has an asymmetric spectrum: its parameters receive unequal gradient signal, explaining why 9D parameterization is preferable. Together, these results provide the theoretical foundation for training with direct 9D regression and applying SVD projection only at inference.

Chris Choy

NVIDIA

3D Representation LearningRepresentation LearningSingular Value Decomposition
AskPDF