Add some (matrix) functions for UniformScaling (#28872)

Co-authored-by: Daniel Karrasch <Daniel.Karrasch@gmx.de>

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,23 +216,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.λ)
 \(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))
 

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,24 +152,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 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))