Add SVD algorithm kwarg to nullspace

src/dense.jl

@@ -1838,7 +1838,7 @@ end
 ## Basis for null space
 
 """
-    nullspace(M; atol::Real=0, rtol::Real=atol>0 ? 0 : n*ϵ)
+    nullspace(M; atol::Real=0, rtol::Real=atol>0 ? 0 : n*ϵ, alg::Algorithm=default_svd_algorithm(A))
     nullspace(M, rtol::Real) = nullspace(M; rtol=rtol) # to be deprecated in Julia 2.0
 
 Computes a basis for the nullspace of `M` by including the singular
@@ -1849,6 +1849,13 @@ By default, the relative tolerance `rtol` is `n*ϵ`, where `n`
 is the size of the smallest dimension of `M`, and `ϵ` is the [`eps`](@ref) of
 the element type of `M`.
 
+The desired algorithm, `alg` is passed through to `svd`. The available algorithms will
+be the same as that of ['svd'](@ref).
+
+!!! compat "Julia 1.3"
+    The `alg` keyword argument requires Julia 1.3 or later.
+
+
 # Examples
 ```jldoctest
 julia> M = [1 0 0; 0 1 0; 0 0 0]
@@ -1876,10 +1883,10 @@ julia> nullspace(M, atol=0.95)
  1.0
 ```
 """
-function nullspace(A::AbstractVecOrMat; atol::Real=0, rtol::Real = (min(size(A, 1), size(A, 2))*eps(real(float(oneunit(eltype(A))))))*iszero(atol))
+function nullspace(A::AbstractVecOrMat; atol::Real=0, rtol::Real = (min(size(A, 1), size(A, 2))*eps(real(float(oneunit(eltype(A))))))*iszero(atol), alg::Algorithm = default_svd_alg(A))
     m, n = size(A, 1), size(A, 2)
     (m == 0 || n == 0) && return Matrix{eigtype(eltype(A))}(I, n, n)
-    SVD = svd(A; full=true)
+    SVD = svd(A; full=true, alg)
     tol = max(atol, SVD.S[1]*rtol)
     indstart = sum(s -> s .> tol, SVD.S) + 1
     return copy((@view SVD.Vt[indstart:end,:])')

src/svd.jl

@@ -639,3 +639,53 @@ AbstractMatrix(F::SVD) = (F.U * Diagonal(F.S)) * F.Vt
 AbstractArray(F::SVD) = AbstractMatrix(F)
 Matrix(F::SVD) = Array(AbstractArray(F))
 Array(F::SVD) = Matrix(F)
+
+"""
+    nullspace(S::SVD{<:Any, T}; atol::Real=0, rtol::Real=min(n,m)*ϵ) where {T}
+
+Returns a basis for the nullspace of `S` where `S` is the SVD of an m \\times n matrix. The same functionality as ['nullspace'](@ref) but using an SVD directly rather than creating one from a matrix. Numerical rank is computed using ['rank'](@ref). 
+
+By default, the relative tolerance `rtol` is `n*ϵ`, where `n`
+is the size of the smallest dimension of `M`, and `ϵ` is the [`eps`](@ref) of
+the element type of `M`.
+
+# Examples
+```jldoctest
+julia> A = [1 1; 2 0; -1 1]*[1 2 3; 3 2 1]
+3×3 Matrix{Int64}:
+ 4  4   4
+ 2  4   6
+ 2  0  -2
+
+julia> F = svd(A)
+SVD{Float64, Float64, Matrix{Float64}, Vector{Float64}}
+U factor:
+3×3 Matrix{Float64}:
+ -0.673105    0.462165  -0.57735
+ -0.736799   -0.351843   0.57735
+  0.0636942   0.814008   0.57735
+singular values:
+3-element Vector{Float64}:
+ 10.027069108300799
+  3.384949792442987
+  5.768888059150692e-16
+Vt factor:
+3×3 Matrix{Float64}:
+ -0.402773  -0.562439  -0.722106
+  0.819212   0.130367  -0.558477
+ -0.408248   0.816497  -0.408248
+
+julia> nullspace(F)
+3×1 Matrix{Float64}:
+ -0.4082482904638629
+  0.816496580927726
+ -0.4082482904638631
+```
+"""
+
+function nullspace(F::SVD; atol::Real=0, rtol::Real = (min(size(F.U, 1), size(F.V, 1))*eps(real(float(oneunit(eltype(F))))))*iszero(atol))
+    r = rank(F; atol, rtol)
+    indstart = r + 1 # nullspace starts after the numerical rank + 1
+    return copy((@view F.Vt[indstart:end,:])')
+end
+

test/dense.jl

@@ -136,6 +136,13 @@ bimg  = randn(n,2)/2
                 @test size(@inferred nullspace(transpose(a[:, 1]))) == (n, n - 1)
                 @test size(@inferred nullspace(transpose(b[1, :]))) == (2, 1)
             end
+
+            @testset "Test SVD algorithm kwarg, issue #1571" begin
+                a_null_QR = nullspace(a, alg=LinearAlgebra.QRIteration())
+                @test norm(a * a_null_QR, Inf) ≈ zero(eltya) atol=n*ε
+                a_null_div = nullspace(a, alg=LinearAlgebra.DivideAndConquer())
+                @test norm(a * a_null_div, Inf) ≈ zero(eltya) atol=n*ε
+            end
         end
     end # for eltyb
 

test/svd.jl

@@ -327,4 +327,15 @@ end
     @test rank(svd([1.0 2.0 3.0; 4.0 5.0 6.0 ; 7.0 8.0 9.0])) == 2
 end
 
+@testset "nullspcae svd" begin
+    # Test that we can get the nullspace
+    A1 = [1.0 1.0; 2.0 0.0; -1.0 1.0] * [1.0 1.0 1.0; 1.0 -1.0 1.0]
+    N1 = nullspace(svd(A1))
+    @test norm(A1*N1) < 3*eps(eltype(A1))
+
+    A2 = [1.0 1.0; -1.0 2.0; 3.0 3.0; 2.0 4.0] * [1.0 1.0 1.0 1.0; 1.0 -1.0 1.0 -1.0]
+    null2 = nullspace(svd(A2))
+    @test norm(A2 * null2) < 8 * maximum(A2) * eps(eltype(A2))
+end
+
 end # module TestSVD