Update null and docs

Author: mtfishman

src/implementations/orthnull.jl

@@ -156,59 +156,70 @@ end
 
 # Implementation of null functions
 # --------------------------------
-function left_null!(A::AbstractMatrix, N; kwargs...)
+function null_truncation_strategy(; atol=nothing, rtol=nothing, maxrank=nothing)
+    if isnothing(maxrank) && isnothing(atol) && isnothing(rtol)
+        return NoTruncation()
+    end
+    atol = @something atol 0
+    rtol = @something rtol 0
+    trunc = TruncationKeepBelow(atol, rtol)
+    return !isnothing(maxrank) ? trunc & truncrank(maxrank; rev=false) : trunc
+end
+
+function left_null!(A::AbstractMatrix, N; trunc=nothing,
+                    kind=isnothing(trunc) ? :qr : :svd, alg_qr=(; positive=true),
+                    alg_svd=(;))
     check_input(left_null!, A, N)
-    atol = get(kwargs, :atol, 0)
-    rtol = get(kwargs, :rtol, 0)
-    kind = get(kwargs, :kind, iszero(atol) && iszero(rtol) ? :qrpos : :svd)
-    if !(iszero(atol) && iszero(rtol)) && kind != :svd
-        throw(ArgumentError("nonzero tolerance not supported for left_orth with kind=$kind"))
+    if !isnothing(trunc) && kind != :svd
+        throw(ArgumentError("truncation not supported for left_null with kind=$kind"))
     end
     if kind == :qr
-        alg = get(kwargs, :alg, select_algorithm(qr_null!, A))
-        return qr_null!(A, N, alg)
-    elseif kind == :qrpos
-        alg = get(kwargs, :alg, select_algorithm(qr_null!, A; positive=true))
-        return qr_null!(A, N, alg)
-    elseif kind == :svd && iszero(atol) && iszero(rtol)
-        alg = get(kwargs, :alg, select_algorithm(svd_full!, A))
-        U, _, _ = svd_full!(A, alg)
+        alg_qr′ = algorithm_or_select_algorithm(qr_null!, A, alg_qr)
+        return qr_null!(A, N, alg_qr′)
+    elseif kind == :svd && isnothing(trunc)
+        alg_svd′ = algorithm_or_select_algorithm(svd_full!, A, alg_svd)
+        U, _, _ = svd_full!(A, alg_svd′)
+        (m, n) = size(A)
+        return copy!(N, view(U, 1:m, (n + 1):m))
+    elseif kind == :svd && isnothing(trunc)
+        alg_svd′ = algorithm_or_select_algorithm(svd_full!, A, alg_svd)
+        U, _, _ = svd_full!(A, alg_svd′)
         (m, n) = size(A)
         return copy!(N, view(U, 1:m, (n + 1):m))
     elseif kind == :svd
-        alg = get(kwargs, :alg, select_algorithm(svd_full!, A))
-        U, S, _ = svd_full!(A, alg)
-        trunc = TruncationKeepBelow(atol, rtol)
-        return truncate!(left_null!, (U, S), trunc)
+        alg_svd′ = algorithm_or_select_algorithm(svd_full!, A, alg_svd)
+        U, S, _ = svd_full!(A, alg_svd′)
+        trunc′ = trunc isa TruncationStrategy ? trunc :
+                 trunc isa NamedTuple ? null_truncation_strategy(; trunc...) :
+                 throw(ArgumentError("Unknown truncation strategy: $trunc"))
+        return truncate!(left_null!, (U, S), trunc′)
     else
         throw(ArgumentError("`left_null!` received unknown value `kind = $kind`"))
     end
 end
 
-function right_null!(A::AbstractMatrix, Nᴴ; kwargs...)
+function right_null!(A::AbstractMatrix, Nᴴ; trunc=nothing,
+                     kind=isnothing(trunc) ? :lq : :svd, alg_lq=(; positive=true),
+                     alg_svd=(;))
     check_input(right_null!, A, Nᴴ)
-    atol = get(kwargs, :atol, 0)
-    rtol = get(kwargs, :rtol, 0)
-    kind = get(kwargs, :kind, iszero(atol) && iszero(rtol) ? :lqpos : :svd)
-    if !(iszero(atol) && iszero(rtol)) && kind != :svd
-        throw(ArgumentError("nonzero tolerance not supported for left_orth with kind=$kind"))
+    if !isnothing(trunc) && kind != :svd
+        throw(ArgumentError("truncation not supported for right_null with kind=$kind"))
     end
     if kind == :lq
-        alg = get(kwargs, :alg, select_algorithm(lq_null!, A))
-        return lq_null!(A, Nᴴ, alg)
-    elseif kind == :lqpos
-        alg = get(kwargs, :alg, select_algorithm(lq_null!, A; positive=true))
-        return lq_null!(A, Nᴴ, alg)
-    elseif kind == :svd && iszero(atol) && iszero(rtol)
-        alg = get(kwargs, :alg, select_algorithm(svd_full!, A))
-        _, _, Vᴴ = svd_full!(A, alg)
+        alg_lq′ = algorithm_or_select_algorithm(lq_null!, A, alg_lq)
+        return lq_null!(A, Nᴴ, alg_lq′)
+    elseif kind == :svd && isnothing(trunc)
+        alg_svd′ = algorithm_or_select_algorithm(svd_full!, A, alg_svd)
+        _, _, Vᴴ = svd_full!(A, alg_svd′)
         (m, n) = size(A)
         return copy!(Nᴴ, view(Vᴴ, (m + 1):n, 1:n))
     elseif kind == :svd
-        alg = get(kwargs, :alg, select_algorithm(svd_full!, A))
-        _, S, Vᴴ = svd_full!(A, alg)
-        trunc = TruncationKeepBelow(atol, rtol)
-        return truncate!(right_null!, (S, Vᴴ), trunc)
+        alg_svd′ = algorithm_or_select_algorithm(svd_full!, A, alg_svd)
+        _, S, Vᴴ = svd_full!(A, alg_svd′)
+        trunc′ = trunc isa TruncationStrategy ? trunc :
+                 trunc isa NamedTuple ? null_truncation_strategy(; trunc...) :
+                 throw(ArgumentError("Unknown truncation strategy: $trunc"))
+        return truncate!(right_null!, (S, Vᴴ), trunc′)
     else
         throw(ArgumentError("`right_null!` received unknown value `kind = $kind`"))
     end

src/implementations/truncation.jl

@@ -69,11 +69,11 @@ TruncationKeepBelow(atol::Real, rtol::Real) = TruncationKeepBelow(promote(atol,
 
 # TODO: better names for these functions of the above types
 """
-    truncrank(howmany::Int, by=abs, rev=true)
+    truncrank(howmany::Int; by=abs, rev=true)
 
 Truncation strategy to keep the first `howmany` values when sorted according to `by` or the last `howmany` if `rev` is true.
 """
-truncrank(howmany::Int, by=abs, rev=true) = TruncationKeepSorted(howmany, by, rev)
+truncrank(howmany::Int; by=abs, rev=true) = TruncationKeepSorted(howmany, by, rev)
 
 """
     trunctol(atol::Real)

src/interface/orthnull.jl

@@ -28,6 +28,9 @@ The keyword argument `kind` can be used to specify the specific orthogonal decom
 that should be used to factor `A`, whereas `trunc` can be used to control the
 precision in determining the rank of `A` via its singular values.
 
+`trunc` can either be a truncation strategy object or a NamedTuple with fields
+`atol`, `rtol`, and `maxrank`.
+
 This is a high-level wrapper and will use one of the decompositions
 `qr_compact!`, `svd_compact!`/`svd_trunc!`, and `left_polar!` to compute the orthogonal basis `V`, as controlled
 by the keyword arguments.
@@ -75,15 +78,18 @@ function left_orth(A::AbstractMatrix; kwargs...)
 end
 
 """
-    right_orth(A; [kind::Symbol, atol::Real=0, rtol::Real=0, alg]) -> C, Vᴴ
-    right_orth!(A, [CVᴴ]; [kind::Symbol, atol::Real=0, rtol::Real=0, alg]) -> C, Vᴴ
+    right_orth(A; [kind::Symbol, trunc, alg_lq, alg_polar, alg_svd]) -> C, Vᴴ
+    right_orth!(A, [CVᴴ]; [kind::Symbol, trunc, alg_lq, alg_polar, alg_svd]) -> C, Vᴴ
 
 Compute an orthonormal basis `V = adjoint(Vᴴ)` for the coimage of the matrix `A`, i.e.
 for the image of `adjoint(A)`, as well as a matrix `C` such that `A = C * Vᴴ`.
 The keyword argument `kind` can be used to specify the specific orthogonal decomposition
 that should be used to factor `A`, whereas `trunc` can be used to control the
 precision in determining the rank of `A` via its singular values.
 
+`trunc` can either be a truncation strategy object or a NamedTuple with fields
+`atol`, `rtol`, and `maxrank`.
+
 This is a high-level wrapper and will use call one of the decompositions
 `lq_compact!`, `svd_compact!`/`svd_trunc!`, and `right_polar!` to compute the
 orthogonal basis `V`, as controlled by the keyword arguments.
@@ -109,7 +115,7 @@ When `kind` is not provided, the default value is `:lq` when `isnothing(trunc)`
 and `:svd` otherwise. Finally, finer control is obtained by providing an explicit algorithm
 for backend factorizations through the `alg_lq`, `alg_polar`, and `alg_svd` keyword arguments,
 which will only be used if the corresponding factorization is called based on the other inputs.
-If NamedTuples are passed as `alg_lq`, `alg_polar`, or `alg_svd`, a default algorithm is chosen
+If `alg_lq`, `alg_polar`, or `alg_svd` are NamedTuples, a default algorithm is chosen
 with `select_algorithm` and the NamedTuple is passed as keyword arguments to that algorithm.
 `alg_lq` defaults to `(; positive=true)` so that by default a positive QR decomposition will
 be used.
@@ -133,36 +139,38 @@ end
 # Null functions
 # --------------
 """
-    left_null(A; [kind::Symbol, atol::Real=0, rtol::Real=0, alg]) -> N
-    left_null!(A, [N]; [kind::Symbol, atol::Real=0, rtol::Real=0, alg]) -> N
+    left_null(A; [kind::Symbol, trunc, alg_qr, alg_svd]) -> N
+    left_null!(A, [N]; [kind::Symbol, alg_qr, alg_svd]) -> N
 
 Compute an orthonormal basis `N` for the cokernel of the matrix `A` of size `(m, n)`, i.e.
 the nullspace of `adjoint(A)`, such that `adjoint(A)*N ≈ 0` and `N'*N ≈ I`.
 The keyword argument `kind` can be used to specify the specific orthogonal decomposition
-that should be used to factor `A`, whereas `atol` and `rtol` can be used to control the
-precision in determining the rank of `A` via its singular values.
+that should be used to factor `A`, whereas `trunc` can be used to control the
+the rank of `A` via its singular values.
+
+`trunc` can either be a truncation strategy object or a NamedTuple with fields
+`atol`, `rtol`, and `maxrank`.
 
 This is a high-level wrapper and will use one of the decompositions `qr!` or `svd!`
 to compute the orthogonal basis `N`, as controlled by the keyword arguments.
 
 When `kind` is provided, its possible values are
 
-*   `kind == :qrpos`: `N` is computed using the positive QR decomposition.
-    This requires `iszero(atol) && iszero(rtol)` and `left_null!(A, [N], kind=:qrpos)` is equivalent to
+*   `kind == :qr`: `N` is computed using the QR decomposition.
+    This requires `isnothing(trunc)` and `left_null!(A, [N], kind=:qr)` is equivalent to
     `qr_null!(A, [N], alg)` with a default value `alg = select_algorithm(qr_compact!, A; positive=true)`
 
-*   `kind == :qr`: `N` is computed using the (nonpositive) QR decomposition.
-    This requires `iszero(atol) && iszero(rtol)` and `left_null!(A, [N], kind=:qr)` is equivalent to
-    `qr_null!(A, [N], alg)` with a default value `alg = select_algorithm(qr_compact!, A)`
-
 *   `kind == :svd`: `N` is computed using the singular value decomposition and will contain 
-    the left singular vectors corresponding to the singular values that
-    are smaller than `max(atol, rtol * σ₁)`, where `σ₁` is the largest singular value of `A`.
+    the left singular vectors corresponding to either the zero singular values if `trunc`
+    isn't specified or the singular values specified by `trunc`.
 
-When `kind` is not provided, the default value is `:qrpos` when `iszero(atol) && iszero(rtol)`
+When `kind` is not provided, the default value is `:qr` when `isnothing(trunc)`
 and `:svd` otherwise. Finally, finer control is obtained by providing an explicit algorithm
-using the `alg` keyword argument, which should be compatible with the chosen or default value
-of `kind`.
+using the `alg_qr` and `alg_svd` keyword arguments, which will only be used by the corresponding
+factorization backend. If `alg_qr` or `alg_svd` are NamedTuples, a default algorithm is chosen
+with `select_algorithm` and the NamedTuple is passed as keyword arguments to that algorithm.
+`alg_qr` defaults to `(; positive=true)` so that by default a positive QR decomposition will
+be used.
 
 !!! note
     The bang method `left_null!` optionally accepts the output structure and possibly destroys
@@ -181,33 +189,32 @@ function left_null(A::AbstractMatrix; kwargs...)
 end
 
 """
-    right_null(A; [kind::Symbol, atol::Real=0, rtol::Real=0, alg]) -> Nᴴ
-    right_null!(A, [Nᴴ]; [kind::Symbol, atol::Real=0, rtol::Real=0, alg]) -> Nᴴ
+    right_null(A; [kind::Symbol, alg_lq, alg_svd]) -> Nᴴ
+    right_null!(A, [Nᴴ]; [kind::Symbol, alg_lq, alg_svd]) -> Nᴴ
 
 Compute an orthonormal basis `N = adjoint(Nᴴ)` for the kernel or nullspace of the matrix `A`
 of size `(m, n)`, such that `A*adjoint(Nᴴ) ≈ 0` and `Nᴴ*adjoint(Nᴴ) ≈ I`.
 The keyword argument `kind` can be used to specify the specific orthogonal decomposition
-that should be used to factor `A`, whereas `atol` and `rtol` can be used to control the
-precision in determining the rank of `A` via its singular values.
+that should be used to factor `A`, whereas `trunc` can be used to control the
+the rank of `A` via its singular values.
+
+`trunc` can either be a truncation strategy object or a NamedTuple with fields
+`atol`, `rtol`, and `maxrank`.
 
 This is a high-level wrapper and will use one of the decompositions `lq!` or `svd!`
 to compute the orthogonal basis `Nᴴ`, as controlled by the keyword arguments.
 
 When `kind` is provided, its possible values are
 
-*   `kind == :lqpos`: `Nᴴ` is computed using the positive LQ decomposition.
-    This requires `iszero(atol) && iszero(rtol)` and `right_null!(A, [Nᴴ], kind=:lqpos)` is equivalent to
-    `lq_null!(A, [Nᴴ], alg)` with a default value `alg = select_algorithm(lq_compact!, A; positive=true)`
-
 *   `kind == :lq`: `Nᴴ` is computed using the (nonpositive) LQ decomposition.
-    This requires `iszero(atol) && iszero(rtol)` and `right_null!(A, [Nᴴ], kind=:lq)` is equivalent to
-    `lq_null!(A, [Nᴴ], alg)` with a default value `alg = select_algorithm(lq_compact!, A)`
+    This requires `isnothing(trunc)` and `right_null!(A, [Nᴴ], kind=:lq)` is equivalent to
+    `lq_null!(A, [Nᴴ], alg)` with a default value `alg = select_algorithm(lq_compact!, A; positive=true)`
 
 *   `kind == :svd`: `N` is computed using the singular value decomposition and will contain 
     the left singular vectors corresponding to the singular values that
     are smaller than `max(atol, rtol * σ₁)`, where `σ₁` is the largest singular value of `A`.
 
-When `kind` is not provided, the default value is `:lqpos` when `iszero(atol) && iszero(rtol)`
+When `kind` is not provided, the default value is `:lq` when `isnothing(trunc)`
 and `:svd` otherwise. Finally, finer control is obtained by providing an explicit algorithm
 using the `alg` keyword argument, which should be compatible with the chosen or default value
 of `kind`.

test/orthnull.jl

@@ -36,7 +36,7 @@ using LinearAlgebra: LinearAlgebra, I
 
         atol = eps(real(T))
         V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); trunc=(; atol=atol))
-        N2 = @constinferred left_null!(copy!(Ac, A), N; atol=atol)
+        N2 = @constinferred left_null!(copy!(Ac, A), N; trunc=(; atol=atol))
         @test V2 !== V
         @test C2 !== C
         @test N2 !== C
@@ -48,7 +48,7 @@ using LinearAlgebra: LinearAlgebra, I
 
         rtol = eps(real(T))
         V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); trunc=(; rtol=rtol))
-        N2 = @constinferred left_null!(copy!(Ac, A), N; rtol=rtol)
+        N2 = @constinferred left_null!(copy!(Ac, A), N; trunc=(; rtol=rtol))
         @test V2 !== V
         @test C2 !== C
         @test N2 !== C
@@ -77,7 +77,8 @@ using LinearAlgebra: LinearAlgebra, I
             if kind == :svd
                 V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); kind=kind,
                                                    trunc=(; atol=atol))
-                N2 = @constinferred left_null!(copy!(Ac, A), N; kind=kind, atol=atol)
+                N2 = @constinferred left_null!(copy!(Ac, A), N; kind=kind,
+                                               trunc=(; atol=atol))
                 @test V2 !== V
                 @test C2 !== C
                 @test N2 !== C
@@ -89,7 +90,8 @@ using LinearAlgebra: LinearAlgebra, I
 
                 V2, C2 = @constinferred left_orth!(copy!(Ac, A), (V, C); kind=kind,
                                                    trunc=(; rtol=rtol))
-                N2 = @constinferred left_null!(copy!(Ac, A), N; kind=kind, rtol=rtol)
+                N2 = @constinferred left_null!(copy!(Ac, A), N; kind=kind,
+                                               trunc=(; rtol=rtol))
                 @test V2 !== V
                 @test C2 !== C
                 @test N2 !== C
@@ -103,8 +105,10 @@ using LinearAlgebra: LinearAlgebra, I
                                                       trunc=(; atol=atol))
                 @test_throws ArgumentError left_orth!(copy!(Ac, A), (V, C); kind=kind,
                                                       trunc=(; rtol=rtol))
-                @test_throws ArgumentError left_null!(copy!(Ac, A), N; kind=kind, atol=atol)
-                @test_throws ArgumentError left_null!(copy!(Ac, A), N; kind=kind, rtol=rtol)
+                @test_throws ArgumentError left_null!(copy!(Ac, A), N; kind=kind,
+                                                      trunc=(; atol=atol))
+                @test_throws ArgumentError left_null!(copy!(Ac, A), N; kind=kind,
+                                                      trunc=(; rtol=rtol))
             end
         end
     end
@@ -142,7 +146,7 @@ end
 
         atol = eps(real(T))
         C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); trunc=(; atol=atol))
-        Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; atol=atol)
+        Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; trunc=(; atol=atol))
         @test C2 !== C
         @test Vᴴ2 !== Vᴴ
         @test Nᴴ2 !== Nᴴ
@@ -154,7 +158,7 @@ end
 
         rtol = eps(real(T))
         C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); trunc=(; rtol=rtol))
-        Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; rtol=rtol)
+        Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; trunc=(; rtol=rtol))
         @test C2 !== C
         @test Vᴴ2 !== Vᴴ
         @test Nᴴ2 !== Nᴴ
@@ -182,7 +186,8 @@ end
             if kind == :svd
                 C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); kind=kind,
                                                      trunc=(; atol=atol))
-                Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; kind=kind, atol=atol)
+                Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; kind=kind,
+                                                 trunc=(; atol=atol))
                 @test C2 !== C
                 @test Vᴴ2 !== Vᴴ
                 @test Nᴴ2 !== Nᴴ
@@ -194,7 +199,8 @@ end
 
                 C2, Vᴴ2 = @constinferred right_orth!(copy!(Ac, A), (C, Vᴴ); kind=kind,
                                                      trunc=(; rtol=rtol))
-                Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; kind=kind, rtol=rtol)
+                Nᴴ2 = @constinferred right_null!(copy!(Ac, A), Nᴴ; kind=kind,
+                                                 trunc=(; rtol=rtol))
                 @test C2 !== C
                 @test Vᴴ2 !== Vᴴ
                 @test Nᴴ2 !== Nᴴ
@@ -209,9 +215,9 @@ end
                 @test_throws ArgumentError right_orth!(copy!(Ac, A), (C, Vᴴ); kind=kind,
                                                        trunc=(; rtol=rtol))
                 @test_throws ArgumentError right_null!(copy!(Ac, A), Nᴴ; kind=kind,
-                                                       atol=atol)
+                                                       trunc=(; atol=atol))
                 @test_throws ArgumentError right_null!(copy!(Ac, A), Nᴴ; kind=kind,
-                                                       rtol=rtol)
+                                                       trunc=(; rtol=rtol))
             end
         end
     end