diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl index bde22d2c2d58a..6eab8b46d145d 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -895,7 +895,6 @@ function pinv(A::StridedMatrix{T}) where T tol = eps(real(float(one(T))))*maximum(size(A)) return pinv(A, tol) end -pinv(a::StridedVector) = pinv(reshape(a, length(a), 1)) function pinv(x::Number) xi = inv(x) return ifelse(isfinite(xi), xi, zero(xi)) diff --git a/base/linalg/generic.jl b/base/linalg/generic.jl index b08354ca30ecd..4af5f66639381 100644 --- a/base/linalg/generic.jl +++ b/base/linalg/generic.jl @@ -794,6 +794,19 @@ function inv(A::AbstractMatrix{T}) where T A_ldiv_B!(factorize(convert(AbstractMatrix{S}, A)), eye(S0, checksquare(A))) end +function pinv(v::AbstractVector{T}, tol::Real=real(zero(T))) where T + res = similar(v, typeof(zero(T) / (abs2(one(T)) + abs2(one(T)))))' + den = sum(abs2, v) + # as tol is the threshold relative to the maximum singular value, for a vector with + # single singular value σ=√den, σ ≦ tol*σ is equivalent to den=0 ∨ tol≥1 + if iszero(den) || tol >= one(tol) + fill!(res, zero(eltype(res))) + else + res .= v' ./ den + end + return res +end + """ \\(A, B) @@ -841,10 +854,11 @@ function (\)(A::AbstractMatrix, B::AbstractVecOrMat) return qrfact(A,Val(true)) \ B end -(\)(a::AbstractVector, b::AbstractArray) = reshape(a, length(a), 1) \ b +(\)(a::AbstractVector, b::AbstractArray) = pinv(a) * b (/)(A::AbstractVecOrMat, B::AbstractVecOrMat) = (B' \ A')' # \(A::StridedMatrix,x::Number) = inv(A)*x Should be added at some point when the old elementwise version has been deprecated long enough # /(x::Number,A::StridedMatrix) = x*inv(A) +/(x::Number, v::AbstractVector) = x*pinv(v) cond(x::Number) = x == 0 ? Inf : 1.0 cond(x::Number, p) = cond(x) diff --git a/base/linalg/rowvector.jl b/base/linalg/rowvector.jl index b7b3270a507ef..37af261e3ea35 100644 --- a/base/linalg/rowvector.jl +++ b/base/linalg/rowvector.jl @@ -227,8 +227,12 @@ Ac_mul_B(::RowVector, ::AbstractVector) = throw(DimensionMismatch("Cannot multip Ac_mul_B(vec::AbstractVector, rowvec::RowVector) = throw(DimensionMismatch("Cannot multiply two transposed vectors")) @inline Ac_mul_B(vec1::AbstractVector, vec2::AbstractVector) = adjoint(vec1)*vec2 +# Pseudo-inverse +pinv(v::RowVector, tol::Real=0) = pinv(v', tol)' + # Left Division # +\(rowvec1::RowVector, rowvec2::RowVector) = pinv(rowvec1) * rowvec2 \(mat::AbstractMatrix, rowvec::RowVector) = throw(DimensionMismatch("Cannot left-divide transposed vector by matrix")) At_ldiv_B(mat::AbstractMatrix, rowvec::RowVector) = throw(DimensionMismatch("Cannot left-divide transposed vector by matrix")) Ac_ldiv_B(mat::AbstractMatrix, rowvec::RowVector) = throw(DimensionMismatch("Cannot left-divide transposed vector by matrix")) diff --git a/test/linalg/dense.jl b/test/linalg/dense.jl index 7eca457de6f2b..1a15a481d407e 100644 --- a/test/linalg/dense.jl +++ b/test/linalg/dense.jl @@ -76,7 +76,7 @@ bimg = randn(n,2)/2 @test_throws DimensionMismatch b'\b @test_throws DimensionMismatch b\b' @test norm(a*x - b, 1)/norm(b) < ε*κ*n*2 # Ad hoc, revisit! - @test zeros(eltya,n)\ones(eltya,n) ≈ zeros(eltya,n,1)\ones(eltya,n,1) + @test zeros(eltya,n)\ones(eltya,n) ≈ (zeros(eltya,n,1)\ones(eltya,n,1))[1,1] end @testset "Test nullspace" begin @@ -613,6 +613,43 @@ end end end +function test_rdiv_pinv_consistency(a, b) + @test (a*b)/b ≈ a*(b/b) ≈ (a*b)*pinv(b) ≈ a*(b*pinv(b)) + @test typeof((a*b)/b) == typeof(a*(b/b)) == typeof((a*b)*pinv(b)) == typeof(a*(b*pinv(b))) +end +function test_ldiv_pinv_consistency(a, b) + @test a\(a*b) ≈ (a\a)*b ≈ (pinv(a)*a)*b ≈ pinv(a)*(a*b) + @test typeof(a\(a*b)) == typeof((a\a)*b) == typeof((pinv(a)*a)*b) == typeof(pinv(a)*(a*b)) +end +function test_div_pinv_consistency(a, b) + test_rdiv_pinv_consistency(a, b) + test_ldiv_pinv_consistency(a, b) +end + +@testset "/ and \\ consistency with pinv for vectors" begin + @testset "Tests for type $elty" for elty in (Float32, Float64, Complex64, Complex128) + c = rand(elty, 5) + r = rand(elty, 5)' + cm = rand(elty, 5, 1) + rm = rand(elty, 1, 5) + @testset "inner prodcuts" begin + test_div_pinv_consistency(r, c) + test_div_pinv_consistency(rm, c) + test_div_pinv_consistency(r, cm) + test_div_pinv_consistency(rm, cm) + end + @testset "outer prodcuts" begin + test_div_pinv_consistency(c, r) + test_div_pinv_consistency(cm, rm) + end + @testset "matrix/vector" begin + m = rand(5, 5) + test_ldiv_pinv_consistency(m, c) + test_rdiv_pinv_consistency(r, m) + end + end +end + @testset "test ops on Numbers for $elty" for elty in [Float32,Float64,Complex64,Complex128] a = rand(elty) @test expm(a) == exp(a) diff --git a/test/linalg/pinv.jl b/test/linalg/pinv.jl index f95f5472b5f8b..0dd73997c441a 100644 --- a/test/linalg/pinv.jl +++ b/test/linalg/pinv.jl @@ -129,6 +129,18 @@ end a = onediag_sparse(eltya, m) test_pinv(a, m, m, default_tol, default_tol, default_tol) end + @testset "Vector" begin + a = rand(eltya, m) + apinv = @inferred pinv(a) + @test pinv(hcat(a)) ≈ apinv + @test apinv isa RowVector{eltya} + end + @testset "RowVector" begin + a = rand(eltya, m)' + apinv = @inferred pinv(a) + @test pinv(vcat(a)) ≈ apinv + @test apinv isa Vector{eltya} + end end end @@ -141,6 +153,10 @@ end @test a[1] ≈ 0.0 @test a[2] ≈ 0.0 + a = pinv([zero(eltya); zero(eltya)]') + @test a[1] ≈ 0.0 + @test a[2] ≈ 0.0 + a = pinv(Diagonal([zero(eltya); zero(eltya)])) @test a.diag[1] ≈ 0.0 @test a.diag[2] ≈ 0.0