From 549a09c2fbe47efa562aedc7f0a12780d8864f2f Mon Sep 17 00:00:00 2001 From: Zepp Date: Wed, 19 Feb 2020 22:24:18 +0800 Subject: [PATCH] Add some (matrix) functions for UniformScaling (#28872) Co-authored-by: Daniel Karrasch --- stdlib/LinearAlgebra/src/uniformscaling.jl | 41 +++++++++- stdlib/LinearAlgebra/test/uniformscaling.jl | 85 ++++++++++++++++++++- 2 files changed, 122 insertions(+), 4 deletions(-) diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl index 862abaed29adf2..7fd8426d175078 100644 --- a/stdlib/LinearAlgebra/src/uniformscaling.jl +++ b/stdlib/LinearAlgebra/src/uniformscaling.jl @@ -88,6 +88,8 @@ end copy(J::UniformScaling) = UniformScaling(J.λ) conj(J::UniformScaling) = UniformScaling(conj(J.λ)) +real(J::UniformScaling) = UniformScaling(real(J.λ)) +imag(J::UniformScaling) = UniformScaling(imag(J.λ)) transpose(J::UniformScaling) = J adjoint(J::UniformScaling) = UniformScaling(conj(J.λ)) @@ -106,11 +108,15 @@ issymmetric(::UniformScaling) = true ishermitian(J::UniformScaling) = isreal(J.λ) isposdef(J::UniformScaling) = isposdef(J.λ) +isposdef(J::UniformScaling) = isposdef(J.λ) + (+)(J::UniformScaling, x::Number) = J.λ + x (+)(x::Number, J::UniformScaling) = x + J.λ (-)(J::UniformScaling, x::Number) = J.λ - x (-)(x::Number, J::UniformScaling) = x - J.λ +(^)(J::UniformScaling, x::Number) = UniformScaling(J.λ ^ x) + (+)(J1::UniformScaling, J2::UniformScaling) = UniformScaling(J1.λ+J2.λ) (+)(B::BitArray{2}, J::UniformScaling) = Array(B) + J (+)(J::UniformScaling, B::BitArray{2}) = J + Array(B) @@ -122,6 +128,21 @@ isposdef(J::UniformScaling) = isposdef(J.λ) (-)(J::UniformScaling, B::BitArray{2}) = J - Array(B) (-)(A::AbstractMatrix, J::UniformScaling) = A + (-J) +# matrix functions +for f in ( :exp, :log, + :expm1, :log1p, + :sqrt, :cbrt, + :sin, :cos, :tan, + :asin, :acos, :atan, + :csc, :sec, :cot, + :acsc, :asec, :acot, + :sinh, :cosh, :tanh, + :asinh, :acosh, :atanh, + :csch, :sech, :coth, + :acsch, :asech, :acoth ) + @eval Base.$f(J::UniformScaling) = UniformScaling($f(J.λ)) +end + # Unit{Lower/Upper}Triangular matrices become {Lower/Upper}Triangular under # addition with a UniformScaling for (t1, t2) in ((:UnitUpperTriangular, :UpperTriangular), @@ -181,6 +202,10 @@ end inv(J::UniformScaling) = UniformScaling(inv(J.λ)) opnorm(J::UniformScaling, p::Real=2) = opnorm(J.λ, p) +pinv(J::UniformScaling) = ifelse(iszero(J.λ), + UniformScaling(zero(inv(J.λ))), # type stability + UniformScaling(inv(J.λ))) + function det(J::UniformScaling{T}) where T if isone(J.λ) one(T) @@ -191,10 +216,19 @@ function det(J::UniformScaling{T}) where T end end +function tr(J::UniformScaling{T}) where T + if iszero(J.λ) + zero(T) + else + throw(ArgumentError("Trace of UniformScaling is only well-defined when λ = 0")) + end +end + *(J1::UniformScaling, J2::UniformScaling) = UniformScaling(J1.λ*J2.λ) *(B::BitArray{2}, J::UniformScaling) = *(Array(B), J::UniformScaling) *(J::UniformScaling, B::BitArray{2}) = *(J::UniformScaling, Array(B)) *(A::AbstractMatrix, J::UniformScaling) = A*J.λ +*(v::AbstractVector, J::UniformScaling) = reshape(v, length(v), 1) * J *(J::UniformScaling, A::AbstractVecOrMat) = J.λ*A *(x::Number, J::UniformScaling) = UniformScaling(x*J.λ) *(J::UniformScaling, x::Number) = UniformScaling(J.λ*x) @@ -202,12 +236,11 @@ end /(J1::UniformScaling, J2::UniformScaling) = J2.λ == 0 ? throw(SingularException(1)) : UniformScaling(J1.λ/J2.λ) /(J::UniformScaling, A::AbstractMatrix) = lmul!(J.λ, inv(A)) /(A::AbstractMatrix, J::UniformScaling) = J.λ == 0 ? throw(SingularException(1)) : A/J.λ +/(v::AbstractVector, J::UniformScaling) = reshape(v, length(v), 1) / J /(J::UniformScaling, x::Number) = UniformScaling(J.λ/x) \(J1::UniformScaling, J2::UniformScaling) = J1.λ == 0 ? throw(SingularException(1)) : UniformScaling(J1.λ\J2.λ) -\(A::Union{Bidiagonal{T},AbstractTriangular{T}}, J::UniformScaling) where {T<:Number} = - rmul!(inv(A), J.λ) \(J::UniformScaling, A::AbstractVecOrMat) = J.λ == 0 ? throw(SingularException(1)) : J.λ\A \(A::AbstractMatrix, J::UniformScaling) = rmul!(inv(A), J.λ) \(F::Factorization, J::UniformScaling) = F \ J(size(F,1)) @@ -229,6 +262,10 @@ Broadcast.broadcasted(::typeof(*), J::UniformScaling,x::Number) = UniformScaling Broadcast.broadcasted(::typeof(/), J::UniformScaling,x::Number) = UniformScaling(J.λ/x) +Broadcast.broadcasted(::typeof(\), x::Number,J::UniformScaling) = UniformScaling(x\J.λ) + +Broadcast.broadcasted(::typeof(^), J::UniformScaling,x::Number) = UniformScaling(J.λ^x) + (^)(J::UniformScaling, x::Number) = UniformScaling((J.λ)^x) Base.literal_pow(::typeof(^), J::UniformScaling, x::Val) = UniformScaling(Base.literal_pow(^, J.λ, x)) diff --git a/stdlib/LinearAlgebra/test/uniformscaling.jl b/stdlib/LinearAlgebra/test/uniformscaling.jl index 419627322c318f..1564198443f247 100644 --- a/stdlib/LinearAlgebra/test/uniformscaling.jl +++ b/stdlib/LinearAlgebra/test/uniformscaling.jl @@ -27,6 +27,53 @@ Random.seed!(123) @test opnorm(UniformScaling(1+im)) ≈ sqrt(2) end +@testset "sqrt, exp, log, and trigonometric functions" begin + # convert to a dense matrix with random size + M(J) = (N = rand(1:10); Matrix(J, N, N)) + + # on complex plane + J = UniformScaling(randn(ComplexF64)) + for f in ( exp, log, + sqrt, + sin, cos, tan, + asin, acos, atan, + csc, sec, cot, + acsc, asec, acot, + sinh, cosh, tanh, + asinh, acosh, atanh, + csch, sech, coth, + acsch, asech, acoth ) + @test f(J) ≈ f(M(J)) + end + + # on real axis + for (λ, fs) in ( + # functions defined for x ∈ ℝ + (()->randn(), (exp, + sin, cos, tan, + csc, sec, cot, + atan, acot, + sinh, cosh, tanh, + csch, sech, coth, + asinh, acsch)), + # functions defined for x ≥ 0 + (()->abs(randn()), (log, sqrt)), + # functions defined for -1 ≤ x ≤ 1 + (()->2rand()-1, (asin, acos, atanh)), + # functions defined for x ≤ -1 or x ≥ 1 + (()->1/(2rand()-1), (acsc, asec, acoth)), + # functions defined for 0 ≤ x ≤ 1 + (()->rand(), (asech,)), + # functions defined for x ≥ 1 + (()->1/rand(), (acosh,)) + ) + for f in fs + J = UniformScaling(λ()) + @test f(J) ≈ f(M(J)) + end + end +end + @testset "conjugation of UniformScaling" begin @test conj(UniformScaling(1))::UniformScaling{Int} == UniformScaling(1) @test conj(UniformScaling(1.0))::UniformScaling{Float64} == UniformScaling(1.0) @@ -42,6 +89,10 @@ end @test issymmetric(UniformScaling(complex(1.0,1.0))) @test ishermitian(I) @test !ishermitian(UniformScaling(complex(1.0,1.0))) + @test isposdef(UniformScaling(rand())) + @test !isposdef(UniformScaling(-rand())) + @test !isposdef(UniformScaling(randn(ComplexF64))) + @test !isposdef(UniformScaling(NaN)) @test isposdef(I) @test !isposdef(-I) @test isposdef(UniformScaling(complex(1.0, 0.0))) @@ -64,8 +115,10 @@ end @test I - α == 1 - α @test α .* UniformScaling(1.0) == UniformScaling(1.0) .* α @test UniformScaling(α)./α == UniformScaling(1.0) + @test α.\UniformScaling(α) == UniformScaling(1.0) @test α * UniformScaling(1.0) == UniformScaling(1.0) * α @test UniformScaling(α)/α == UniformScaling(1.0) + @test (2I)^α == (2I).^α == (2^α)I β = rand() @test (α*I)^2 == UniformScaling(α^2) @@ -77,7 +130,11 @@ end @test (α * I) .^ β == UniformScaling(α^β) end -@testset "det and logdet" begin +@testset "tr, det and logdet" begin + for T in (Int, Float64, ComplexF64, Bool) + @test tr(UniformScaling(zero(T))) === zero(T) + end + @test_throws ArgumentError tr(UniformScaling(1)) @test det(I) === true @test det(1.0I) === 1.0 @test det(0I) === 0 @@ -95,17 +152,28 @@ end let λ = complex(randn(),randn()) J = UniformScaling(λ) - @testset "transpose, conj, inv" begin + @testset "transpose, conj, inv, pinv, cond" begin @test ndims(J) == 2 @test transpose(J) == J @test J * [1 0; 0 1] == conj(*(adjoint(J), [1 0; 0 1])) # ctranpose (and A(c)_mul_B) @test I + I === UniformScaling(2) # + @test inv(I) == I @test inv(J) == UniformScaling(inv(λ)) + @test pinv(J) == UniformScaling(inv(λ)) + @test @inferred(pinv(0.0I)) == 0.0I + @test @inferred(pinv(0I)) == 0.0I + @test @inferred(pinv(false*I)) == 0.0I + @test @inferred(pinv(0im*I)) == 0im*I @test cond(I) == 1 @test cond(J) == (λ ≠ zero(λ) ? one(real(λ)) : oftype(real(λ), Inf)) end + @testset "real, imag, reim" begin + @test real(J) == UniformScaling(real(λ)) + @test imag(J) == UniformScaling(imag(λ)) + @test reim(J) == (UniformScaling(real(λ)), UniformScaling(imag(λ))) + end + @testset "copyto!" begin A = Matrix{Int}(undef, (3,3)) @test copyto!(A, I) == one(A) @@ -113,6 +181,19 @@ let @test copyto!(B, J) == [λ zero(λ)] end + @testset "binary ops with vectors" begin + v = complex.(randn(3), randn(3)) + # As shown in #20423@GitHub, vector acts like x1 matrix when participating in linear algebra + @test v * J ≈ v * λ + @test v' * J ≈ v' * λ + @test J * v ≈ λ * v + @test J * v' ≈ λ * v' + @test v / J ≈ v / λ + @test v' / J ≈ v' / λ + @test J \ v ≈ λ \ v + @test J \ v' ≈ λ \ v' + end + @testset "binary ops with matrices" begin B = bitrand(2, 2) @test B + I == B + Matrix(I, size(B))