Drop RandomizedLinAlg (REQUIRE-only release breaks Aqua persistent-tasks) and implement RSVD locally

Project.toml

@@ -8,7 +8,6 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-RandomizedLinAlg = "0448d7d9-159c-5637-8537-fd72090fea46"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
@@ -19,7 +18,6 @@ DocStringExtensions = "0.9.5"
 LinearAlgebra = "1"
 ModelingToolkit = "11.42"
 PrecompileTools = "1.2.1"
-RandomizedLinAlg = "0.1"
 SciMLBase = "3.48"
 SparseArrays = "1"
 SymbolicIndexingInterface = "0.3.49"

src/DataReduction/POD.jl

@@ -1,5 +1,4 @@
 import TSVD as TruncatedSVD
-using RandomizedLinAlg: rsvd
 
 function matricize(VoV::Vector{Vector{T}})::Matrix{T} where {T}
     return reduce(hcat, VoV)
@@ -24,7 +23,14 @@ function _rsvd(data::Vector{Vector{T}}, n::Int, p::Int) where {T}
     return _rsvd(mat_data, n, p)
 end
 
-_rsvd(data, n::Int, p::Int) = rsvd(data, n, p)
+# Halko, Martinsson & Tropp (2011), Algorithms 4.1 and 5.1.
+function _rsvd(A, n::Int, p::Int)
+    l = n + p
+    l > size(A, 1) && throw(ArgumentError("Cannot find $l linearly independent vectors of $(size(A, 1)) x $(size(A, 2)) matrix"))
+    Q = Matrix(qr!(A * randn(size(A, 2), l)).Q)
+    F = svd!(Q' * A)
+    return (Q * F.U)[:, 1:n], F.S[1:n], F.V[:, 1:n]
+end
 
 """
     POD(snapshots; min_renergy = 1.0, min_nmodes = 1,
@@ -215,8 +221,8 @@ randomized SVD.
 - `nothing`: `pod` is updated in place.
 
 # Throws
-- An exception propagated by `RandomizedLinAlg.rsvd` if the snapshot matrix is not
-  accepted by the randomized SVD backend.
+- `ArgumentError`: if `pod.nmodes + alg.p` exceeds the number of rows of the snapshot
+  matrix.
 
 # Examples
 ```jldoctest

src/ModelOrderReduction.jl

@@ -8,7 +8,7 @@ using ModelingToolkit: ModelingToolkit, @parameters, @variables, Differential, E
     System, SymbolicUtils, Symbolics, complete, mtkcompile, substitute
 using SciMLBase: SciMLBase
 using SymbolicIndexingInterface: SymbolicIndexingInterface
-using LinearAlgebra: LinearAlgebra, mul!, qr, svd, lyap, Symmetric, Diagonal,
+using LinearAlgebra: LinearAlgebra, mul!, qr, qr!, svd, svd!, lyap, Symmetric, Diagonal,
     eigen, eigvals, I, factorize, norm, ColumnNorm
 using SparseArrays: SparseArrays, SparseMatrixCSC, findnz, sparse
 

src/Types.jl

@@ -90,7 +90,7 @@ end
 Randomized singular value decomposition backend for [`reduce!`](@ref).
 
 # Arguments
-- `p::Int = 0`: number of oversampling vectors used by `RandomizedLinAlg.rsvd`.
+- `p::Int = 0`: number of oversampling vectors used by the randomized range finder.
 
 # Fields
 - `p::Int`: the oversampling parameter.

test/DataReduction.jl

@@ -1,5 +1,7 @@
 using Test, ModelOrderReduction
+using LinearAlgebra: svdvals
 using OrdinaryDiffEq
+using Random: randn, seed!
 
 function lorenz_prob()
     function lorenz!(du, u, p, t)
@@ -65,3 +67,14 @@ end
     @test size(reducer.rbasis, 1) == size(solution, 1)
     @test reducer.renergy > 0.7
 end
+
+@testset "POD - RSVD" begin
+    seed!(1)
+    A = randn(40, 3) * randn(3, 25)
+    pod = POD(A, 3)
+    reduce!(pod, RSVD(2))
+
+    @test pod.spectrum ≈ svdvals(A)[1:3] rtol = 1.0e-10
+    @test pod.rbasis * (pod.rbasis' * A) ≈ A rtol = 1.0e-10
+    @test_throws ArgumentError reduce!(POD(A[1:3, :], 3), RSVD(1))
+end

test/Project.toml

@@ -4,6 +4,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MethodOfLines = "94925ecb-adb7-4558-8ed8-f975c56a0bf4"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
@@ -15,6 +16,7 @@ LinearAlgebra = "1, < 0.0.1"
 MethodOfLines = "1"
 ModelingToolkit = "11.42"
 OrdinaryDiffEq = "7"
+Random = "1"
 SafeTestsets = "0.1"
 SciMLBase = "3"
 SciMLTesting = "2.4"