Add some (matrix) functions for UniformScaling

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,20 @@ 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,
+           :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 )
+    @eval $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),
@@ -183,6 +203,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)
@@ -193,23 +217,31 @@ 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)
 
 /(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.λ)
 
@@ -227,6 +259,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)
+
 ==(J1::UniformScaling,J2::UniformScaling) = (J1.λ == J2.λ)
 
 ## equality comparison with UniformScaling

stdlib/LinearAlgebra/test/uniformscaling.jl

@@ -23,6 +23,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)
@@ -38,6 +85,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)))
@@ -60,11 +111,17 @@ 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
 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
@@ -82,24 +139,48 @@ 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)
         B = Matrix{ComplexF64}(undef, (1,2))
         @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 all(v  * J .== v  * λ)
+        @test all(v' * J .== v' * λ)
+        @test all(J * v  .== λ * v )
+        @test all(J * v' .== λ * v')
+        @test all(v  / J .== v  / λ)
+        @test all(v' / J .== v' / λ)
+        @test all(J \ v  .== λ \ v )
+        @test all(J \ v' .== λ \ v')
+    end
+
     @testset "binary ops with matrices" begin
         B = bitrand(2, 2)
         @test B + I == B + Matrix(I, size(B))