Deprecate sqrtm in favor of sqrt. (#23504)

NEWS.md

@@ -355,6 +355,8 @@ Deprecated or removed
     full path if you need access to executables or libraries in the `JULIA_HOME` directory, e.g.
     `joinpath(JULIA_HOME, "7z.exe")` for `7z.exe` ([#21540]).
 
+  * `sqrtm` has been deprecated in favor of `sqrt` ([#23504]).
+
   * `expm` has been deprecated in favor of `exp` ([#23233]).
 
   * `logm` has been deprecated in favor of `log` ([#23505]).

base/deprecated.jl

@@ -253,7 +253,7 @@ for f in (
         # base/math.jl
         :cbrt, :sinh, :cosh, :tanh, :atan, :asinh, :exp2,
         :expm1, :exp10, :sin, :cos, :tan, :asin, :acos, :acosh, :atanh,
-        :log2, :log10, :lgamma, #=:log1p,=# :sqrt,
+        :log2, :log10, :lgamma, #=:log1p,=#
         # base/floatfuncs.jl
         :abs, :abs2, :angle, :isnan, :isinf, :isfinite,
         # base/complex.jl
@@ -1620,6 +1620,9 @@ function Tridiagonal(dl::AbstractVector{Tl}, d::AbstractVector{Td}, du::Abstract
     Tridiagonal(map(v->convert(Vector{promote_type(Tl,Td,Tu)}, v), (dl, d, du))...)
 end
 
+# deprecate sqrtm in favor of sqrt
+@deprecate sqrtm sqrt
+
 # deprecate expm in favor of exp
 @deprecate expm! exp!
 @deprecate expm exp

base/exports.jl

@@ -609,7 +609,6 @@ export
     schur,
     schurfact!,
     schurfact,
-    sqrtm,
     svd,
     svdfact!,
     svdfact,

base/linalg/dense.jl

@@ -374,8 +374,8 @@ function schurpow(A::AbstractMatrix, p)
         retmat = A ^ floor(p)
         # Real part
         if p - floor(p) == 0.5
-            # special case: A^0.5 === sqrtm(A)
-            retmat = retmat * sqrtm(A)
+            # special case: A^0.5 === sqrt(A)
+            retmat = retmat * sqrt(A)
         else
             retmat = retmat * powm!(UpperTriangular(float.(A)), real(p - floor(p)))
         end
@@ -385,8 +385,8 @@ function schurpow(A::AbstractMatrix, p)
         R = S ^ floor(p)
         # Real part
         if p - floor(p) == 0.5
-            # special case: A^0.5 === sqrtm(A)
-            R = R * sqrtm(S)
+            # special case: A^0.5 === sqrt(A)
+            R = R * sqrt(S)
         else
             R = R * powm!(UpperTriangular(float.(S)), real(p - floor(p)))
         end
@@ -618,7 +618,7 @@ function log(A::StridedMatrix{T}) where T
 end
 
 """
-    sqrtm(A)
+    sqrt(A::AbstractMatrix)
 
 If `A` has no negative real eigenvalues, compute the principal matrix square root of `A`,
 that is the unique matrix ``X`` with eigenvalues having positive real part such that
@@ -642,40 +642,38 @@ julia> A = [4 0; 0 4]
  4  0
  0  4
 
-julia> sqrtm(A)
+julia> sqrt(A)
 2×2 Array{Float64,2}:
  2.0  0.0
  0.0  2.0
 ```
 """
-function sqrtm(A::StridedMatrix{<:Real})
+function sqrt(A::StridedMatrix{<:Real})
     if issymmetric(A)
-        return full(sqrtm(Symmetric(A)))
+        return full(sqrt(Symmetric(A)))
     end
     n = checksquare(A)
     if istriu(A)
-        return full(sqrtm(UpperTriangular(A)))
+        return full(sqrt(UpperTriangular(A)))
     else
         SchurF = schurfact(complex(A))
-        R = full(sqrtm(UpperTriangular(SchurF[:T])))
+        R = full(sqrt(UpperTriangular(SchurF[:T])))
         return SchurF[:vectors] * R * SchurF[:vectors]'
     end
 end
-function sqrtm(A::StridedMatrix{<:Complex})
+function sqrt(A::StridedMatrix{<:Complex})
     if ishermitian(A)
-        return full(sqrtm(Hermitian(A)))
+        return full(sqrt(Hermitian(A)))
     end
     n = checksquare(A)
     if istriu(A)
-        return full(sqrtm(UpperTriangular(A)))
+        return full(sqrt(UpperTriangular(A)))
     else
         SchurF = schurfact(A)
-        R = full(sqrtm(UpperTriangular(SchurF[:T])))
+        R = full(sqrt(UpperTriangular(SchurF[:T])))
         return SchurF[:vectors] * R * SchurF[:vectors]'
     end
 end
-sqrtm(a::Number) = (b = sqrt(complex(a)); imag(b) == 0 ? real(b) : b)
-sqrtm(a::Complex) = sqrt(a)
 
 function inv(A::StridedMatrix{T}) where T
     checksquare(A)

base/linalg/diagonal.jl

@@ -328,8 +328,7 @@ eye(::Type{Diagonal{T}}, n::Int) where {T} = Diagonal(ones(T,n))
 # Matrix functions
 exp(D::Diagonal) = Diagonal(exp.(D.diag))
 log(D::Diagonal) = Diagonal(log.(D.diag))
-sqrtm(D::Diagonal) = Diagonal(sqrt.(D.diag))
-sqrtm(D::Diagonal{<:AbstractMatrix}) = Diagonal(sqrtm.(D.diag))
+sqrt(D::Diagonal) = Diagonal(sqrt.(D.diag))
 
 #Linear solver
 function A_ldiv_B!(D::Diagonal, B::StridedVecOrMat)

base/linalg/linalg.jl

@@ -9,7 +9,7 @@ import Base: USE_BLAS64, abs, big, broadcast, ceil, conj, convert, copy, copy!,
     adjoint, eltype, exp, eye, findmax, findmin, fill!, floor, full, getindex,
     hcat, imag, indices, inv, isapprox, isone, IndexStyle, kron, length, log, map,
     ndims, oneunit, parent, power_by_squaring, print_matrix, promote_rule, real, round,
-    setindex!, show, similar, size, transpose, trunc, typed_hcat
+    setindex!, show, similar, size, sqrt, transpose, trunc, typed_hcat
 using Base: hvcat_fill, iszero, IndexLinear, _length, promote_op, promote_typeof,
     @propagate_inbounds, @pure, reduce, typed_vcat
 # We use `_length` because of non-1 indices; releases after julia 0.5
@@ -124,7 +124,6 @@ export
     schur,
     schurfact!,
     schurfact,
-    sqrtm,
     svd,
     svdfact!,
     svdfact,

base/linalg/symmetric.jl

@@ -607,40 +607,7 @@ function exp(A::Hermitian{T}) where T
     end
 end
 
-for (funm, func) in ([:sqrtm,:sqrt],)
-    @eval begin
-        function ($funm)(A::Symmetric{T}) where T<:Real
-            F = eigfact(A)
-            if all(λ -> λ ≥ 0, F.values)
-                retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors'
-            else
-                retmat = (F.vectors * Diagonal(($func).(complex.(F.values)))) * F.vectors'
-            end
-            return Symmetric(retmat)
-        end
-
-        function ($funm)(A::Hermitian{T}) where T
-            n = checksquare(A)
-            F = eigfact(A)
-            if all(λ -> λ ≥ 0, F.values)
-                retmat = (F.vectors * Diagonal(($func).(F.values))) * F.vectors'
-                if T <: Real
-                    return Hermitian(retmat)
-                else
-                    for i = 1:n
-                        retmat[i,i] = real(retmat[i,i])
-                    end
-                    return Hermitian(retmat)
-                end
-            else
-                retmat = (F.vectors * Diagonal(($func).(complex(F.values)))) * F.vectors'
-                return retmat
-            end
-        end
-    end
-end
-
-for func in (:log, #=:sqrtm=#)
+for func in (:log, :sqrt)
     @eval begin
         function ($func)(A::Symmetric{T}) where T<:Real
             F = eigfact(A)

base/linalg/triangular.jl

@@ -1815,7 +1815,7 @@ function log(A0::UpperTriangular{T}) where T<:Union{Float64,Complex{Float64}}
     end
     s0 = s
     for k = 1:min(s, maxsqrt)
-        A = sqrtm(A)
+        A = sqrt(A)
     end
 
     AmI = A - I
@@ -1871,7 +1871,7 @@ function log(A0::UpperTriangular{T}) where T<:Union{Float64,Complex{Float64}}
             m = tmax
             break
         end
-        A = sqrtm(A)
+        A = sqrt(A)
         AmI = A - I
         s = s + 1
     end
@@ -2015,7 +2015,7 @@ function invsquaring(A0::UpperTriangular, theta)
     end
     s0 = s
     for k = 1:min(s, maxsqrt)
-        A = sqrtm(A)
+        A = sqrt(A)
     end
 
     AmI = A - I
@@ -2073,7 +2073,7 @@ function invsquaring(A0::UpperTriangular, theta)
                 m = tmax
                 break
             end
-            A = sqrtm(A)
+            A = sqrt(A)
             AmI = A - I
             s = s + 1
         end
@@ -2119,7 +2119,7 @@ unw(x::Number) = ceil((imag(x) - pi) / (2 * pi))
 
 # End of auxiliary functions for matrix logarithm and matrix power
 
-function sqrtm(A::UpperTriangular)
+function sqrt(A::UpperTriangular)
     realmatrix = false
     if isreal(A)
         realmatrix = true
@@ -2131,9 +2131,9 @@ function sqrtm(A::UpperTriangular)
             end
         end
     end
-    sqrtm(A,Val(realmatrix))
+    sqrt(A,Val(realmatrix))
 end
-function sqrtm(A::UpperTriangular{T},::Val{realmatrix}) where {T,realmatrix}
+function sqrt(A::UpperTriangular{T},::Val{realmatrix}) where {T,realmatrix}
     B = A.data
     n = checksquare(B)
     t = realmatrix ? typeof(sqrt(zero(T))) : typeof(sqrt(complex(zero(T))))
@@ -2151,7 +2151,7 @@ function sqrtm(A::UpperTriangular{T},::Val{realmatrix}) where {T,realmatrix}
     end
     return UpperTriangular(R)
 end
-function sqrtm(A::UnitUpperTriangular{T}) where T
+function sqrt(A::UnitUpperTriangular{T}) where T
     B = A.data
     n = checksquare(B)
     t = typeof(sqrt(zero(T)))
@@ -2169,8 +2169,8 @@ function sqrtm(A::UnitUpperTriangular{T}) where T
     end
     return UnitUpperTriangular(R)
 end
-sqrtm(A::LowerTriangular) = sqrtm(A.').'
-sqrtm(A::UnitLowerTriangular) = sqrtm(A.').'
+sqrt(A::LowerTriangular) = sqrt(A.').'
+sqrt(A::UnitLowerTriangular) = sqrt(A.').'
 
 # Generic eigensystems
 eigvals(A::AbstractTriangular) = diag(A)

doc/src/manual/linear-algebra.md

@@ -177,8 +177,8 @@ as well as whether hooks to various optimized methods for them in LAPACK are ava
 
 | Matrix type               | `+` | `-` | `*` | `\` | Other functions with optimized methods                              |
 |:------------------------- |:--- |:--- |:--- |:--- |:------------------------------------------------------------------- |
-| [`Symmetric`](@ref)       |     |     |     | MV  | [`inv()`](@ref), [`sqrtm()`](@ref), [`exp()`](@ref)                |
-| [`Hermitian`](@ref)       |     |     |     | MV  | [`inv()`](@ref), [`sqrtm()`](@ref), [`exp()`](@ref)                |
+| [`Symmetric`](@ref)       |     |     |     | MV  | [`inv()`](@ref), [`sqrt()`](@ref), [`exp()`](@ref)                |
+| [`Hermitian`](@ref)       |     |     |     | MV  | [`inv()`](@ref), [`sqrt()`](@ref), [`exp()`](@ref)                |
 | [`UpperTriangular`](@ref) |     |     | MV  | MV  | [`inv()`](@ref), [`det()`](@ref)                                    |
 | [`LowerTriangular`](@ref) |     |     | MV  | MV  | [`inv()`](@ref), [`det()`](@ref)                                    |
 | [`SymTridiagonal`](@ref)  | M   | M   | MS  | MV  | [`eigmax()`](@ref), [`eigmin()`](@ref)                              |

doc/src/stdlib/linalg.md

@@ -92,7 +92,6 @@ Base.repmat
 Base.kron
 Base.SparseArrays.blkdiag
 Base.LinAlg.linreg
-Base.LinAlg.sqrtm
 Base.LinAlg.lyap
 Base.LinAlg.sylvester
 Base.LinAlg.issuccess

test/linalg/dense.jl

@@ -115,10 +115,10 @@ bimg  = randn(n,2)/2
         end
 
         @testset "Matrix square root" begin
-            asq = sqrtm(a)
+            asq = sqrt(a)
             @test asq*asq ≈ a
             asym = a'+a # symmetric indefinite
-            asymsq = sqrtm(asym)
+            asymsq = sqrt(asym)
             @test asymsq*asymsq ≈ asym
         end
 
@@ -370,10 +370,10 @@ end
 
 @testset "issue #2246" begin
     A = [1 2 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 0]
-    Asq = sqrtm(A)
+    Asq = sqrt(A)
     @test Asq*Asq ≈ A
     A2 = view(A, 1:2, 1:2)
-    A2sq = sqrtm(A2)
+    A2sq = sqrt(A2)
     @test A2sq*A2sq ≈ A2
 
     N = 3
@@ -569,12 +569,12 @@ end
     end
 
     for A in (Aa, Ab, Ac, Ad, Ah, ADi)
-        @test A^(1/2) ≈ sqrtm(A)
-        @test A^(-1/2) ≈ inv(sqrtm(A))
-        @test A^(3/4) ≈ sqrtm(A) * sqrtm(sqrtm(A))
-        @test A^(-3/4) ≈ inv(A) * sqrtm(sqrtm(A))
-        @test A^(17/8) ≈ A^2 * sqrtm(sqrtm(sqrtm(A)))
-        @test A^(-17/8) ≈ inv(A^2 * sqrtm(sqrtm(sqrtm(A))))
+        @test A^(1/2) ≈ sqrt(A)
+        @test A^(-1/2) ≈ inv(sqrt(A))
+        @test A^(3/4) ≈ sqrt(A) * sqrt(sqrt(A))
+        @test A^(-3/4) ≈ inv(A) * sqrt(sqrt(A))
+        @test A^(17/8) ≈ A^2 * sqrt(sqrt(sqrt(A)))
+        @test A^(-17/8) ≈ inv(A^2 * sqrt(sqrt(sqrt(A))))
         @test (A^0.2)^5 ≈ A
         @test (A^(2/3))*(A^(1/3)) ≈ A
         @test (A^im)^(-im) ≈ A
@@ -680,7 +680,6 @@ end
 @testset "test ops on Numbers for $elty" for elty in [Float32,Float64,Complex64,Complex128]
     a = rand(elty)
     @test isposdef(one(elty))
-    @test sqrtm(a) == sqrt(a)
     @test lyap(one(elty),a) == -a/2
 end
 

test/linalg/diagonal.jl

@@ -66,7 +66,7 @@ srand(1)
             @test log(Diagonal(abs.(D.diag))) ≈ log(abs.(DM)) atol=n^3*eps(relty)
         end
         if elty <: BlasComplex
-            for func in (logdet, sqrtm)
+            for func in (logdet, sqrt)
                 @test func(D) ≈ func(DM) atol=n^2*eps(relty)*2
             end
         end
@@ -383,7 +383,7 @@ end
 
     @test exp(D) == Diagonal([exp([1 2; 3 4]), exp([1 2; 3 4])])
     @test log(D) == Diagonal([log([1 2; 3 4]), log([1 2; 3 4])])
-    @test sqrtm(D) == Diagonal([sqrtm([1 2; 3 4]), sqrtm([1 2; 3 4])])
+    @test sqrt(D) == Diagonal([sqrt([1 2; 3 4]), sqrt([1 2; 3 4])])
 end
 
 @testset "multiplication with Symmetric/Hermitian" begin

test/linalg/triangular.jl

@@ -237,7 +237,7 @@ for elty1 in (Float32, Float64, BigFloat, Complex64, Complex128, Complex{BigFloa
         @test ladb ≈ fladb atol=sqrt(eps(real(float(one(elty1)))))*n*n
 
         # Matrix square root
-        @test sqrtm(A1) |> t -> t*t ≈ A1
+        @test sqrt(A1) |> t -> t*t ≈ A1
 
         # naivesub errors
         @test_throws DimensionMismatch naivesub!(A1,ones(elty1,n+1))
@@ -391,10 +391,10 @@ end
 # Matrix square root
 Atn = UpperTriangular([-1 1 2; 0 -2 2; 0 0 -3])
 Atp = UpperTriangular([1 1 2; 0 2 2; 0 0 3])
-@test sqrtm(Atn) |> t->t*t ≈ Atn
-@test typeof(sqrtm(Atn)[1,1]) <: Complex
-@test sqrtm(Atp) |> t->t*t ≈ Atp
-@test typeof(sqrtm(Atp)[1,1]) <: Real
+@test sqrt(Atn) |> t->t*t ≈ Atn
+@test typeof(sqrt(Atn)[1,1]) <: Complex
+@test sqrt(Atp) |> t->t*t ≈ Atp
+@test typeof(sqrt(Atp)[1,1]) <: Real
 
 Areal   = randn(n, n)/2
 Aimg    = randn(n, n)/2
@@ -511,5 +511,5 @@ end
 isdefined(Main, :TestHelpers) || @eval Main include("../TestHelpers.jl")
 using TestHelpers.Furlong
 let A = UpperTriangular([Furlong(1) Furlong(4); Furlong(0) Furlong(1)])
-    @test sqrtm(A) == Furlong{1//2}.(UpperTriangular([1 2; 0 1]))
+    @test sqrt(A) == Furlong{1//2}.(UpperTriangular([1 2; 0 1]))
 end