From 3b07931c4f2e0f46b80039258ac1233e11b4ca9f Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 23 Jan 2026 12:21:43 +0100 Subject: [PATCH 01/10] Add some tests for randomized SVD --- src/implementations/svd.jl | 12 +----------- test/svd.jl | 3 +++ test/testsuite/svd.jl | 14 ++++++++++++++ 3 files changed, 18 insertions(+), 11 deletions(-) diff --git a/src/implementations/svd.jl b/src/implementations/svd.jl index 9e7e6600..89fc6642 100644 --- a/src/implementations/svd.jl +++ b/src/implementations/svd.jl @@ -404,23 +404,13 @@ function svd_trunc_no_error!(A::AbstractMatrix, USVᴴ, alg::TruncatedAlgorithm{ end function svd_trunc!(A::AbstractMatrix, USVᴴ, alg::TruncatedAlgorithm{<:GPU_Randomized}) - U, S, Vᴴ = USVᴴ - check_input(svd_trunc!, A, (U, S, Vᴴ), alg.alg) - _gpu_Xgesvdr!(A, diagview(S), U, Vᴴ; alg.alg.kwargs...) - - # TODO: make sure that truncation is based on maxrank, otherwise this might be wrong - (Utr, Str, Vᴴtr), _ = truncate(svd_trunc!, (U, S, Vᴴ), alg.trunc) - + Utr, Str, Vᴴtr = svd_trunc_no_error!(A, USVᴴ, alg) # normal `truncation_error!` does not work here since `S` is not the full singular value spectrum normS = norm(diagview(Str)) normA = norm(A) # equivalent to sqrt(normA^2 - normS^2) # but may be more accurate ϵ = sqrt((normA + normS) * (normA - normS)) - - do_gauge_fix = get(alg.alg.kwargs, :fixgauge, default_fixgauge())::Bool - do_gauge_fix && gaugefix!(svd_trunc!, Utr, Vᴴtr) - return Utr, Str, Vᴴtr, ϵ end diff --git a/test/svd.jl b/test/svd.jl index affe2942..4c1714ad 100644 --- a/test/svd.jl +++ b/test/svd.jl @@ -25,6 +25,9 @@ for T in (BLASFloats..., GenericFloats...), m in (0, 54), n in (0, 37, m, 63) CUSOLVER_Jacobi(), ) TestSuite.test_svd_algs(CuMatrix{T}, (m, n), CUDA_SVD_ALGS) + k = 5 + p = min(m, n) - 2 + TestSuite.test_randomized_svd(CuMatrix{T}, (m, n), (CUSOLVER_Randomized(; k, p, niters = 20),)) if n == m TestSuite.test_svd(Diagonal{T, CuVector{T}}, m) TestSuite.test_svd_algs(Diagonal{T, CuVector{T}}, m, (DiagonalAlgorithm(),)) diff --git a/test/testsuite/svd.jl b/test/testsuite/svd.jl index 800018b2..23f05a41 100644 --- a/test/testsuite/svd.jl +++ b/test/testsuite/svd.jl @@ -312,3 +312,17 @@ function test_svd_trunc_algs( end end end + +function test_randomized_svd(T::Type, sz, algs; kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "randomized svd_trunc! algorithm $alg $summary_str" for alg in algs + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + m, n = size(A) + minmn = min(m, n) + S₀ = collect(svd_vals(A)) + U1, S1, V1ᴴ, ϵ1 = @testinferred svd_trunc(A; alg) + @test length(diagview(S1)) == alg.alg.k + @test collect(diagview(S1)) ≈ S₀[1:alg.alg.k] + end +end From 7401ab37a14677a1d2dbea4c66709152a9e79ab6 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 23 Jan 2026 14:50:26 +0100 Subject: [PATCH 02/10] Don't run randomized SVD on empty arrays --- test/svd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/svd.jl b/test/svd.jl index 4c1714ad..e8781f31 100644 --- a/test/svd.jl +++ b/test/svd.jl @@ -27,7 +27,7 @@ for T in (BLASFloats..., GenericFloats...), m in (0, 54), n in (0, 37, m, 63) TestSuite.test_svd_algs(CuMatrix{T}, (m, n), CUDA_SVD_ALGS) k = 5 p = min(m, n) - 2 - TestSuite.test_randomized_svd(CuMatrix{T}, (m, n), (CUSOLVER_Randomized(; k, p, niters = 20),)) + min(m, n) > k && TestSuite.test_randomized_svd(CuMatrix{T}, (m, n), (CUSOLVER_Randomized(; k, p, niters = 20),)) if n == m TestSuite.test_svd(Diagonal{T, CuVector{T}}, m) TestSuite.test_svd_algs(Diagonal{T, CuVector{T}}, m, (DiagonalAlgorithm(),)) From 42c448cd1497ef61c399ed5bae1153a6186c20f6 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 23 Jan 2026 15:32:22 +0100 Subject: [PATCH 03/10] Fix again --- test/svd.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/svd.jl b/test/svd.jl index e8781f31..5de9cfab 100644 --- a/test/svd.jl +++ b/test/svd.jl @@ -26,8 +26,8 @@ for T in (BLASFloats..., GenericFloats...), m in (0, 54), n in (0, 37, m, 63) ) TestSuite.test_svd_algs(CuMatrix{T}, (m, n), CUDA_SVD_ALGS) k = 5 - p = min(m, n) - 2 - min(m, n) > k && TestSuite.test_randomized_svd(CuMatrix{T}, (m, n), (CUSOLVER_Randomized(; k, p, niters = 20),)) + p = min(m, n) - k - 2 + min(m, n) > k + 2 && TestSuite.test_randomized_svd(CuMatrix{T}, (m, n), (CUSOLVER_Randomized(; k, p, niters = 20),)) if n == m TestSuite.test_svd(Diagonal{T, CuVector{T}}, m) TestSuite.test_svd_algs(Diagonal{T, CuVector{T}}, m, (DiagonalAlgorithm(),)) From 5e1f7821419963dc4591ef6e568e9444ea78b4f8 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 23 Jan 2026 15:41:42 +0100 Subject: [PATCH 04/10] One more fix --- test/testsuite/svd.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/testsuite/svd.jl b/test/testsuite/svd.jl index 23f05a41..06683689 100644 --- a/test/testsuite/svd.jl +++ b/test/testsuite/svd.jl @@ -322,7 +322,7 @@ function test_randomized_svd(T::Type, sz, algs; kwargs...) minmn = min(m, n) S₀ = collect(svd_vals(A)) U1, S1, V1ᴴ, ϵ1 = @testinferred svd_trunc(A; alg) - @test length(diagview(S1)) == alg.alg.k - @test collect(diagview(S1)) ≈ S₀[1:alg.alg.k] + @test length(diagview(S1)) == alg.k + @test collect(diagview(S1)) ≈ S₀[1:alg.k] end end From eb86e19aba91daf7dc8f4218c780eb51c6f83796 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 23 Jan 2026 15:50:55 +0100 Subject: [PATCH 05/10] Another fix --- test/testsuite/svd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/testsuite/svd.jl b/test/testsuite/svd.jl index 06683689..1ed5a1e4 100644 --- a/test/testsuite/svd.jl +++ b/test/testsuite/svd.jl @@ -323,6 +323,6 @@ function test_randomized_svd(T::Type, sz, algs; kwargs...) S₀ = collect(svd_vals(A)) U1, S1, V1ᴴ, ϵ1 = @testinferred svd_trunc(A; alg) @test length(diagview(S1)) == alg.k - @test collect(diagview(S1)) ≈ S₀[1:alg.k] + @test collect(diagview(S1))[1:alg.k] ≈ S₀[1:alg.k] end end From 86eeb484d6b040a20dc9c4cd3c922a26e6b7cab9 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 23 Jan 2026 16:01:24 +0100 Subject: [PATCH 06/10] And one more --- src/implementations/svd.jl | 3 ++- test/testsuite/svd.jl | 1 - 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/implementations/svd.jl b/src/implementations/svd.jl index 89fc6642..01a151d6 100644 --- a/src/implementations/svd.jl +++ b/src/implementations/svd.jl @@ -398,7 +398,8 @@ function svd_trunc_no_error!(A::AbstractMatrix, USVᴴ, alg::TruncatedAlgorithm{ (Utr, Str, Vᴴtr), _ = truncate(svd_trunc!, (U, S, Vᴴ), alg.trunc) do_gauge_fix = get(alg.alg.kwargs, :fixgauge, default_fixgauge())::Bool - do_gauge_fix && gaugefix!(svd_trunc!, Utr, Vᴴtr) + # the output matrices here are the same size as for svd_full! + do_gauge_fix && gaugefix!(svd_full!, Utr, Vᴴtr) return Utr, Str, Vᴴtr end diff --git a/test/testsuite/svd.jl b/test/testsuite/svd.jl index 1ed5a1e4..b65bd49e 100644 --- a/test/testsuite/svd.jl +++ b/test/testsuite/svd.jl @@ -322,7 +322,6 @@ function test_randomized_svd(T::Type, sz, algs; kwargs...) minmn = min(m, n) S₀ = collect(svd_vals(A)) U1, S1, V1ᴴ, ϵ1 = @testinferred svd_trunc(A; alg) - @test length(diagview(S1)) == alg.k @test collect(diagview(S1))[1:alg.k] ≈ S₀[1:alg.k] end end From f9cd6d5125ccdc165befe6d5d52e906d8242fe7b Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 23 Jan 2026 16:17:12 +0100 Subject: [PATCH 07/10] Not sure about this --- src/implementations/svd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/implementations/svd.jl b/src/implementations/svd.jl index 01a151d6..f2b54997 100644 --- a/src/implementations/svd.jl +++ b/src/implementations/svd.jl @@ -411,7 +411,7 @@ function svd_trunc!(A::AbstractMatrix, USVᴴ, alg::TruncatedAlgorithm{<:GPU_Ran normA = norm(A) # equivalent to sqrt(normA^2 - normS^2) # but may be more accurate - ϵ = sqrt((normA + normS) * (normA - normS)) + ϵ = sqrt((normA + normS) * abs(normA - normS)) return Utr, Str, Vᴴtr, ϵ end From 0da38370edaf5fbcce5657b93f7ec088457e2bfb Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 23 Jan 2026 16:21:54 +0100 Subject: [PATCH 08/10] Fix trunc arg whoops --- src/implementations/svd.jl | 2 +- test/svd.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/implementations/svd.jl b/src/implementations/svd.jl index f2b54997..7dc19c8d 100644 --- a/src/implementations/svd.jl +++ b/src/implementations/svd.jl @@ -399,7 +399,7 @@ function svd_trunc_no_error!(A::AbstractMatrix, USVᴴ, alg::TruncatedAlgorithm{ do_gauge_fix = get(alg.alg.kwargs, :fixgauge, default_fixgauge())::Bool # the output matrices here are the same size as for svd_full! - do_gauge_fix && gaugefix!(svd_full!, Utr, Vᴴtr) + do_gauge_fix && gaugefix!(svd_trunc!, Utr, Vᴴtr) return Utr, Str, Vᴴtr end diff --git a/test/svd.jl b/test/svd.jl index 5de9cfab..295ba5ba 100644 --- a/test/svd.jl +++ b/test/svd.jl @@ -27,7 +27,7 @@ for T in (BLASFloats..., GenericFloats...), m in (0, 54), n in (0, 37, m, 63) TestSuite.test_svd_algs(CuMatrix{T}, (m, n), CUDA_SVD_ALGS) k = 5 p = min(m, n) - k - 2 - min(m, n) > k + 2 && TestSuite.test_randomized_svd(CuMatrix{T}, (m, n), (CUSOLVER_Randomized(; k, p, niters = 20),)) + min(m, n) > k + 2 && TestSuite.test_randomized_svd(CuMatrix{T}, (m, n), (TruncatedAlgorithm(CUSOLVER_Randomized(; k, p, niters = 20), truncrank(k)),)) if n == m TestSuite.test_svd(Diagonal{T, CuVector{T}}, m) TestSuite.test_svd_algs(Diagonal{T, CuVector{T}}, m, (DiagonalAlgorithm(),)) From bc586e589f0e65a83ff87e31a0fdbf5155a5aacb Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 23 Jan 2026 16:45:19 +0100 Subject: [PATCH 09/10] Another typo --- test/svd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/svd.jl b/test/svd.jl index 295ba5ba..bc9e8c17 100644 --- a/test/svd.jl +++ b/test/svd.jl @@ -27,7 +27,7 @@ for T in (BLASFloats..., GenericFloats...), m in (0, 54), n in (0, 37, m, 63) TestSuite.test_svd_algs(CuMatrix{T}, (m, n), CUDA_SVD_ALGS) k = 5 p = min(m, n) - k - 2 - min(m, n) > k + 2 && TestSuite.test_randomized_svd(CuMatrix{T}, (m, n), (TruncatedAlgorithm(CUSOLVER_Randomized(; k, p, niters = 20), truncrank(k)),)) + min(m, n) > k + 2 && TestSuite.test_randomized_svd(CuMatrix{T}, (m, n), (MatrixAlgebraKit.TruncatedAlgorithm(CUSOLVER_Randomized(; k, p, niters = 20), truncrank(k)),)) if n == m TestSuite.test_svd(Diagonal{T, CuVector{T}}, m) TestSuite.test_svd_algs(Diagonal{T, CuVector{T}}, m, (DiagonalAlgorithm(),)) From 3afd7bccc9a642c448f29684941f9b535dde3164 Mon Sep 17 00:00:00 2001 From: Katharine Hyatt Date: Fri, 23 Jan 2026 17:42:52 +0100 Subject: [PATCH 10/10] Keyword args --- test/testsuite/svd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/testsuite/svd.jl b/test/testsuite/svd.jl index b65bd49e..6d5799ba 100644 --- a/test/testsuite/svd.jl +++ b/test/testsuite/svd.jl @@ -322,6 +322,6 @@ function test_randomized_svd(T::Type, sz, algs; kwargs...) minmn = min(m, n) S₀ = collect(svd_vals(A)) U1, S1, V1ᴴ, ϵ1 = @testinferred svd_trunc(A; alg) - @test collect(diagview(S1))[1:alg.k] ≈ S₀[1:alg.k] + @test collect(diagview(S1))[1:alg.alg.k] ≈ S₀[1:alg.alg.k] end end