diff --git a/src/interface/orthnull.jl b/src/interface/orthnull.jl index b20a841d..ad60648a 100644 --- a/src/interface/orthnull.jl +++ b/src/interface/orthnull.jl @@ -443,7 +443,9 @@ left_orth_alg(alg::LeftOrthAlgorithm) = alg left_orth_alg(alg::QRAlgorithms) = LeftOrthViaQR(alg) left_orth_alg(alg::PolarAlgorithms) = LeftOrthViaPolar(alg) left_orth_alg(alg::SVDAlgorithms) = LeftOrthViaSVD(alg) +left_orth_alg(alg::DiagonalAlgorithm) = LeftOrthViaQR(alg) left_orth_alg(alg::TruncatedAlgorithm{<:SVDAlgorithms}) = LeftOrthViaSVD(alg) +left_orth_alg(alg::TruncatedAlgorithm{DiagonalAlgorithm}) = LeftOrthViaSVD(alg) """ right_orth_alg(alg::AbstractAlgorithm) -> RightOrthAlgorithm @@ -478,7 +480,9 @@ right_orth_alg(alg::RightOrthAlgorithm) = alg right_orth_alg(alg::LQAlgorithms) = RightOrthViaLQ(alg) right_orth_alg(alg::PolarAlgorithms) = RightOrthViaPolar(alg) right_orth_alg(alg::SVDAlgorithms) = RightOrthViaSVD(alg) +right_orth_alg(alg::DiagonalAlgorithm) = RightOrthViaLQ(alg) right_orth_alg(alg::TruncatedAlgorithm{<:SVDAlgorithms}) = RightOrthViaSVD(alg) +right_orth_alg(alg::TruncatedAlgorithm{DiagonalAlgorithm}) = RightOrthViaSVD(alg) """ left_null_alg(alg::AbstractAlgorithm) -> LeftNullAlgorithm diff --git a/test/orthnull.jl b/test/orthnull.jl index dec12946..26c81747 100644 --- a/test/orthnull.jl +++ b/test/orthnull.jl @@ -19,16 +19,18 @@ for T in (BLASFloats..., GenericFloats...), n in (37, m, 63) if T ∈ BLASFloats if CUDA.functional() TestSuite.test_orthnull(CuMatrix{T}, (m, n); test_nullity = false) - n == m && TestSuite.test_orthnull(Diagonal{T, CuVector{T}}, m; test_orthnull = false) + n == m && TestSuite.test_orthnull(Diagonal{T, CuVector{T}}, m) end if AMDGPU.functional() TestSuite.test_orthnull(ROCMatrix{T}, (m, n); test_nullity = false) - n == m && TestSuite.test_orthnull(Diagonal{T, ROCVector{T}}, m; test_orthnull = false) + n == m && TestSuite.test_orthnull(Diagonal{T, ROCVector{T}}, m) + end + if AMDGPU.functional() end end if !is_buildkite TestSuite.test_orthnull(T, (m, n)) AT = Diagonal{T, Vector{T}} - TestSuite.test_orthnull(AT, m; test_orthnull = false) + TestSuite.test_orthnull(AT, m) end end diff --git a/test/testsuite/TestSuite.jl b/test/testsuite/TestSuite.jl index 361f125b..95fec775 100644 --- a/test/testsuite/TestSuite.jl +++ b/test/testsuite/TestSuite.jl @@ -73,9 +73,11 @@ is_pivoted(alg::MatrixAlgebraKit.LQViaTransposedQR) = is_pivoted(alg.qr_alg) isleftcomplete(V, N) = V * V' + N * N' ≈ I isleftcomplete(V::AnyCuMatrix, N::AnyCuMatrix) = isleftcomplete(collect(V), collect(N)) isleftcomplete(V::AnyROCMatrix, N::AnyROCMatrix) = isleftcomplete(collect(V), collect(N)) +isleftcomplete(V::Diagonal{TV, <:AnyROCVector}, N::Diagonal{TN, <:AnyROCVector}) where {TV, TN} = isleftcomplete(Diagonal(collect(V.diag)), Diagonal(collect(N.diag))) isrightcomplete(Vᴴ, Nᴴ) = Vᴴ' * Vᴴ + Nᴴ' * Nᴴ ≈ I isrightcomplete(V::AnyCuMatrix, N::AnyCuMatrix) = isrightcomplete(collect(V), collect(N)) isrightcomplete(V::AnyROCMatrix, N::AnyROCMatrix) = isrightcomplete(collect(V), collect(N)) +isrightcomplete(V::Diagonal{TV, <:AnyROCVector}, N::Diagonal{TN, <:AnyROCVector}) where {TV, TN} = isrightcomplete(Diagonal(collect(V.diag)), Diagonal(collect(N.diag))) instantiate_unitary(T, A, sz) = qr_compact(randn!(similar(A, eltype(T), sz, sz)))[1] # AMDGPU can't generate ComplexF32 random numbers diff --git a/test/testsuite/orthnull.jl b/test/testsuite/orthnull.jl new file mode 100644 index 00000000..16b86c82 --- /dev/null +++ b/test/testsuite/orthnull.jl @@ -0,0 +1,324 @@ +using TestExtras +using LinearAlgebra + +include("../linearmap.jl") + +_left_orth_svd(x; kwargs...) = left_orth(x; alg = :svd, kwargs...) +_left_orth_svd!(x, VC; kwargs...) = left_orth!(x, VC; alg = :svd, kwargs...) +_left_orth_qr(x; kwargs...) = left_orth(x; alg = :qr, kwargs...) +_left_orth_qr!(x, VC; kwargs...) = left_orth!(x, VC; alg = :qr, kwargs...) +_left_orth_polar(x; kwargs...) = left_orth(x; alg = :polar, kwargs...) +_left_orth_polar!(x, VC; kwargs...) = left_orth!(x, VC; alg = :polar, kwargs...) + +_right_orth_svd(x; kwargs...) = right_orth(x; alg = :svd, kwargs...) +_right_orth_svd!(x, CVᴴ; kwargs...) = right_orth!(x, CVᴴ; alg = :svd, kwargs...) +_right_orth_lq(x; kwargs...) = right_orth(x; alg = :lq, kwargs...) +_right_orth_lq!(x, CVᴴ; kwargs...) = right_orth!(x, CVᴴ; alg = :lq, kwargs...) +_right_orth_polar(x; kwargs...) = right_orth(x; alg = :polar, kwargs...) +_right_orth_polar!(x, CVᴴ; kwargs...) = right_orth!(x, CVᴴ; alg = :polar, kwargs...) + +function test_orthnull(T::Type, sz; test_nullity = true, kwargs...) + summary_str = testargs_summary(T, sz) + return @testset "orthnull $summary_str" begin + test_left_orthnull(T, sz; kwargs...) + test_nullity && test_left_nullity(T, sz; kwargs...) + test_right_orthnull(T, sz; kwargs...) + test_nullity && test_right_nullity(T, sz; kwargs...) + end +end + +function test_left_nullity( + T::Type, sz; + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "left_nullity $summary_str" begin + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + m, n = size(A) + minmn = min(m, n) + if m > n + nullity = 5 + V, C = @testinferred left_orth(A) + N = @testinferred left_null(A; trunc = (; maxnullity = nullity)) + @test V isa typeof(A) && size(V) == (m, minmn) + @test C isa typeof(A) && size(C) == (minmn, n) + @test eltype(N) == eltype(A) && size(N) == (m, nullity) + @test V * C ≈ A + @test isisometric(V) + @test LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T)) + @test isisometric(N) + end + + rtol = eps(real(eltype(T))) + for (trunc_orth, trunc_null) in ( + ((; rtol = rtol), (; rtol = rtol)), + (trunctol(; rtol), trunctol(; rtol, keep_below = true)), + ) + V, C = left_orth(A) + N = left_null(A) + V2, C2 = @testinferred left_orth!(copy!(Ac, A), (V, C); trunc = trunc_orth) + N2 = @testinferred left_null!(copy!(Ac, A), N; trunc = trunc_null) + @test V2 * C2 ≈ A + @test isisometric(V2) + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T)) + @test isisometric(N2) + @test isleftcomplete(V2, N2) + end + + alg = :svd + V2, C2 = @testinferred _left_orth_svd!(copy!(Ac, A), (V, C); trunc = (; atol)) + N2 = @testinferred left_null!(copy!(Ac, A), N; alg, trunc = (; atol)) + @test V2 * C2 ≈ A + @test isisometric(V2) + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T)) + @test isisometric(N2) + @test isleftcomplete(V2, N2) + + V2, C2 = @testinferred _left_orth_svd!(copy!(Ac, A), (V, C); trunc = (; rtol)) + N2 = @testinferred left_null!(copy!(Ac, A), N; alg, trunc = (; rtol)) + @test V2 * C2 ≈ A + @test isisometric(V2) + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T)) + @test isisometric(N2) + @test isleftcomplete(V2, N2) + + # doesn't work on AMD... + atol = eps(real(eltype(T))) + V2, C2 = @testinferred left_orth!(copy!(Ac, A), (V, C); trunc = (; atol = atol)) + N2 = @testinferred left_null!(copy!(Ac, A), N; trunc = (; atol = atol)) + @test V2 * C2 ≈ A + @test isisometric(V2) + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N2) + @test isleftcomplete(V2, N2) + + end +end + +function test_left_orthnull( + T::Type, sz; + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "left_orth! and left_null! $summary_str" begin + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + V, C = @testinferred left_orth(A) + N = @testinferred left_null(A) + m, n = size(A) + minmn = min(m, n) + @test V isa typeof(A) && size(V) == (m, minmn) + @test C isa typeof(A) && size(C) == (minmn, n) + @test eltype(N) == eltype(A) && size(N) == (m, m - minmn) + @test V * C ≈ A + @test isisometric(V) + @test LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N) + @test isleftcomplete(V, N) + + M = LinearMap(A) + VM, CM = @testinferred _left_orth_svd(M) + @test parent(VM) * parent(CM) ≈ A + + # passing a kind and some kwargs + V, C = @testinferred _left_orth_qr(A; positive = true) + N = @testinferred left_null(A; alg = :qr, positive = true) + @test V isa typeof(A) && size(V) == (m, minmn) + @test C isa typeof(A) && size(C) == (minmn, n) + @test eltype(N) == eltype(A) && size(N) == (m, m - minmn) + @test V * C ≈ A + @test isisometric(V) + @test LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N) + @test isleftcomplete(V, N) + + # passing an algorithm + V, C = @testinferred left_orth(A; alg = MatrixAlgebraKit.default_qr_algorithm(A)) + N = @testinferred left_null(A; alg = :qr, positive = true) + @test V isa typeof(A) && size(V) == (m, minmn) + @test C isa typeof(A) && size(C) == (minmn, n) + @test eltype(N) == eltype(A) && size(N) == (m, m - minmn) + @test V * C ≈ A + @test isisometric(V) + @test LinearAlgebra.norm(A' * N) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N) + @test isleftcomplete(V, N) + + V2, C2 = @testinferred left_orth!(copy!(Ac, A), (V, C)) + N2 = @testinferred left_null!(copy!(Ac, A), N) + @test V2 * C2 ≈ A + @test isisometric(V2) + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N2) + @test isleftcomplete(V2, N2) + + for alg in (:qr, :polar, :svd) # explicit kind kwarg + m < n && alg === :polar && continue + if alg == :svd + V2, C2 = @testinferred _left_orth_svd!(copy!(Ac, A), (V, C)) + elseif alg == :qr + V2, C2 = @testinferred _left_orth_qr!(copy!(Ac, A), (V, C)) + elseif alg == :polar + V2, C2 = @testinferred _left_orth_polar!(copy!(Ac, A), (V, C)) + end + @test V2 * C2 ≈ A + @test isisometric(V2) + if alg != :polar + N2 = @testinferred left_null!(copy!(Ac, A), N; alg) + @test LinearAlgebra.norm(A' * N2) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(N2) + @test isleftcomplete(V2, N2) + end + + # with kind and tol kwargs + if alg != :svd + @test_throws ArgumentError left_orth!(copy!(Ac, A), (V, C); alg, trunc = (; atol)) + @test_throws ArgumentError left_orth!(copy!(Ac, A), (V, C); alg, trunc = (; rtol)) + alg == :polar && continue + @test_throws ArgumentError left_null!(copy!(Ac, A), N; alg, trunc = (; atol)) + @test_throws ArgumentError left_null!(copy!(Ac, A), N; alg, trunc = (; rtol)) + end + end + end +end + +function test_right_nullity( + T::Type, sz; + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "right_nullity $summary_str" begin + A = instantiate_matrix(T, sz) + Ac = deepcopy(A) + m, n = size(A) + minmn = min(m, n) + + C, Vᴴ = @testinferred right_orth(A) + Nᴴ = @testinferred right_null(A) + atol = eps(real(eltype(T))) + C2, Vᴴ2 = @testinferred right_orth!(copy!(Ac, A), (C, Vᴴ); trunc = (; atol)) + Nᴴ2 = @testinferred right_null!(copy!(Ac, A), Nᴴ; trunc = (; atol)) + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T)) + @test isisometric(Nᴴ; side = :right) + @test isrightcomplete(Vᴴ2, Nᴴ2) + + rtol = eps(real(eltype(T))) + C2, Vᴴ2 = @testinferred right_orth!(copy!(Ac, A), (C, Vᴴ); trunc = (; rtol)) + Nᴴ2 = @testinferred right_null!(copy!(Ac, A), Nᴴ; trunc = (; rtol)) + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T)) + @test isisometric(Nᴴ2; side = :right) + @test isrightcomplete(Vᴴ2, Nᴴ2) + + alg = :svd + C2, Vᴴ2 = @testinferred _right_orth_svd!(copy!(Ac, A), (C, Vᴴ); trunc = (; atol)) + Nᴴ2 = @testinferred right_null!(copy!(Ac, A), Nᴴ; alg = alg, trunc = (; atol)) + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(Nᴴ2; side = :right) + @test isrightcomplete(Vᴴ2, Nᴴ2) + + C2, Vᴴ2 = @testinferred _right_orth_svd!(copy!(Ac, A), (C, Vᴴ); trunc = (; rtol)) + Nᴴ2 = @testinferred right_null!(copy!(Ac, A), Nᴴ; alg = alg, trunc = (; rtol)) + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(T) + @test isisometric(Nᴴ2; side = :right) + @test isrightcomplete(Vᴴ2, Nᴴ2) + end +end + +function test_right_orthnull( + T::Type, sz; + atol::Real = 0, rtol::Real = precision(T), + kwargs... + ) + summary_str = testargs_summary(T, sz) + return @testset "right_orth! and right_null! $summary_str" begin + A = instantiate_matrix(T, sz) + m, n = size(A) + minmn = min(m, n) + Ac = deepcopy(A) + C, Vᴴ = @testinferred right_orth(A) + Nᴴ = @testinferred right_null(A) + @test C isa typeof(A) && size(C) == (m, minmn) + @test Vᴴ isa typeof(A) && size(Vᴴ) == (minmn, n) + @test eltype(Nᴴ) == eltype(A) && size(Nᴴ) == (n - minmn, n) + @test C * Vᴴ ≈ A + @test isisometric(Vᴴ; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T)) + @test isisometric(Nᴴ; side = :right) + @test isrightcomplete(Vᴴ, Nᴴ) + + M = LinearMap(A) + CM, VMᴴ = @testinferred _right_orth_svd(M) + @test parent(CM) * parent(VMᴴ) ≈ A + + # passing a kind and some kwargs + C, Vᴴ = @testinferred _right_orth_lq(A; positive = true) + Nᴴ = @testinferred right_null(A; alg = :lq, positive = true) + @test C isa typeof(A) && size(C) == (m, minmn) + @test Vᴴ isa typeof(A) && size(Vᴴ) == (minmn, n) + @test eltype(Nᴴ) == eltype(A) && size(Nᴴ) == (n - minmn, n) + @test C * Vᴴ ≈ A + @test isisometric(Vᴴ; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T)) + @test isisometric(Nᴴ; side = :right) + @test isrightcomplete(Vᴴ, Nᴴ) + + # passing an algorithm + C, Vᴴ = @testinferred right_orth(A; alg = MatrixAlgebraKit.default_lq_algorithm(A)) + @test C isa typeof(A) && size(C) == (m, minmn) + @test Vᴴ isa typeof(A) && size(Vᴴ) == (minmn, n) + Nᴴ = @testinferred right_null(A; alg = :lq, positive = true) + @test eltype(Nᴴ) == eltype(A) && size(Nᴴ) == (n - minmn, n) + @test C * Vᴴ ≈ A + @test isisometric(Vᴴ; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T)) + @test isisometric(Nᴴ; side = :right) + @test isrightcomplete(Vᴴ, Nᴴ) + + C2, Vᴴ2 = @testinferred right_orth!(copy!(Ac, A), (C, Vᴴ)) + Nᴴ2 = @testinferred right_null!(copy!(Ac, A), Nᴴ) + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T)) + @test isisometric(Nᴴ; side = :right) + @test isrightcomplete(Vᴴ2, Nᴴ2) + + for alg in (:lq, :polar, :svd) + n < m && alg == :polar && continue + if alg == :lq + C2, Vᴴ2 = @testinferred _right_orth_lq!(copy!(Ac, A), (C, Vᴴ)) + elseif alg == :polar + C2, Vᴴ2 = @testinferred _right_orth_polar!(copy!(Ac, A), (C, Vᴴ)) + elseif alg == :svd + C2, Vᴴ2 = @testinferred _right_orth_svd!(copy!(Ac, A), (C, Vᴴ)) + end + @test C2 * Vᴴ2 ≈ A + @test isisometric(Vᴴ2; side = :right) + if alg != :polar + Nᴴ2 = @testinferred right_null!(copy!(Ac, A), Nᴴ; alg = alg) + @test LinearAlgebra.norm(A * adjoint(Nᴴ2)) ≈ 0 atol = MatrixAlgebraKit.defaulttol(eltype(T)) + @test isisometric(Nᴴ2; side = :right) + @test isrightcomplete(Vᴴ2, Nᴴ2) + end + + if alg != :svd + @test_throws ArgumentError right_orth!(copy!(Ac, A), (C, Vᴴ); alg, trunc = (; atol)) + @test_throws ArgumentError right_orth!(copy!(Ac, A), (C, Vᴴ); alg, trunc = (; rtol)) + alg == :polar && continue + @test_throws ArgumentError right_null!(copy!(Ac, A), Nᴴ; alg, trunc = (; atol)) + @test_throws ArgumentError right_null!(copy!(Ac, A), Nᴴ; alg, trunc = (; rtol)) + end + end + end +end