From aaa32c2f3472aa0a46dda20b7d5e4e747aa7a8e7 Mon Sep 17 00:00:00 2001 From: araujoms Date: Tue, 13 May 2025 15:54:09 +0200 Subject: [PATCH 1/5] optimization of exp, sqrt, ^ --- src/symmetric.jl | 23 ++++++++++++++++++++--- 1 file changed, 20 insertions(+), 3 deletions(-) diff --git a/src/symmetric.jl b/src/symmetric.jl index e07975ad..436fe27e 100644 --- a/src/symmetric.jl +++ b/src/symmetric.jl @@ -833,6 +833,14 @@ function svdvals!(A::RealHermSymComplexHerm) return sort!(vals, rev = true) end +#computes U * Diagonal(abs2.(v)) * U', destroying U +function _psd_spectral_product!(v, U) + @inbounds for j ∈ axes(U, 2), i ∈ axes(U, 1) + U[i, j] *= v[j] + end + return U * U' +end + # Matrix functions ^(A::SymSymTri{<:Complex}, p::Integer) = sympow(A, p) ^(A::SelfAdjoint, p::Integer) = sympow(A, p) @@ -848,7 +856,8 @@ function ^(A::SelfAdjoint, p::Real) isinteger(p) && return integerpow(A, p) F = eigen(A) if all(λ -> λ ≥ 0, F.values) - retmat = (F.vectors * Diagonal((F.values).^p)) * F.vectors' + map!(λ -> λ^0.5p, F.values) + retmat = _psd_spectral_product!(F.values, F.vectors) return wrappertype(A)(retmat) else retmat = (F.vectors * Diagonal(complex.(F.values).^p)) * F.vectors' @@ -860,7 +869,7 @@ function ^(A::SymSymTri{<:Complex}, p::Real) return Symmetric(schurpow(A, p)) end -for func in (:exp, :cos, :sin, :tan, :cosh, :sinh, :tanh, :atan, :asinh, :atanh, :cbrt) +for func in (:cos, :sin, :tan, :cosh, :sinh, :tanh, :atan, :asinh, :atanh, :cbrt) @eval begin function ($func)(A::SelfAdjoint) F = eigen(A) @@ -870,6 +879,13 @@ for func in (:exp, :cos, :sin, :tan, :cosh, :sinh, :tanh, :atan, :asinh, :atanh, end end +function exp(A::SelfAdjoint) + F = eigen(A) + map!(λ -> exp(0.5λ), F.values) + retmat = _psd_spectral_product!(F.values, F.vectors) + return wrappertype(A)(retmat) +end + function cis(A::SelfAdjoint) F = eigen(A) retmat = F.vectors .* cis.(F.values') * F.vectors' @@ -929,7 +945,8 @@ function sqrt(A::SelfAdjoint; rtol = eps(real(float(eltype(A)))) * size(A, 1)) F = eigen(A) λ₀ = -maximum(abs, F.values) * rtol # treat λ ≥ λ₀ as "zero" eigenvalues up to roundoff if all(λ -> λ ≥ λ₀, F.values) - retmat = (F.vectors * Diagonal(sqrt.(max.(0, F.values)))) * F.vectors' + map!(λ -> λ < 0 ? zero(λ) : fourthroot(λ), F.values) + retmat = _psd_spectral_product!(F.values, F.vectors) return wrappertype(A)(retmat) else retmat = (F.vectors * Diagonal(sqrt.(complex.(F.values)))) * F.vectors' From 786d49f49ef5ccc4ece5a22348ada2494e545e94 Mon Sep 17 00:00:00 2001 From: araujoms Date: Tue, 13 May 2025 16:25:01 +0200 Subject: [PATCH 2/5] use rmul! --- src/symmetric.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/symmetric.jl b/src/symmetric.jl index 436fe27e..24768a96 100644 --- a/src/symmetric.jl +++ b/src/symmetric.jl @@ -835,9 +835,7 @@ end #computes U * Diagonal(abs2.(v)) * U', destroying U function _psd_spectral_product!(v, U) - @inbounds for j ∈ axes(U, 2), i ∈ axes(U, 1) - U[i, j] *= v[j] - end + rmul!(U, Diagonal(v)) return U * U' end From 8958055f2fcd4e757e594d58762122631e8b264f Mon Sep 17 00:00:00 2001 From: araujoms Date: Tue, 13 May 2025 19:43:31 +0200 Subject: [PATCH 3/5] make it work for immutable types --- src/symmetric.jl | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/symmetric.jl b/src/symmetric.jl index 24768a96..f3be534f 100644 --- a/src/symmetric.jl +++ b/src/symmetric.jl @@ -833,10 +833,10 @@ function svdvals!(A::RealHermSymComplexHerm) return sort!(vals, rev = true) end -#computes U * Diagonal(abs2.(v)) * U', destroying U -function _psd_spectral_product!(v, U) - rmul!(U, Diagonal(v)) - return U * U' +#computes U * Diagonal(abs2.(v)) * U' +function _psd_spectral_product(v, U) + Uv = U * Diagonal(v) + return Uv * Uv' end # Matrix functions @@ -854,8 +854,8 @@ function ^(A::SelfAdjoint, p::Real) isinteger(p) && return integerpow(A, p) F = eigen(A) if all(λ -> λ ≥ 0, F.values) - map!(λ -> λ^0.5p, F.values) - retmat = _psd_spectral_product!(F.values, F.vectors) + rootpower = map(λ -> λ^0.5p, F.values) + retmat = _psd_spectral_product(rootpower, F.vectors) return wrappertype(A)(retmat) else retmat = (F.vectors * Diagonal(complex.(F.values).^p)) * F.vectors' @@ -879,8 +879,8 @@ end function exp(A::SelfAdjoint) F = eigen(A) - map!(λ -> exp(0.5λ), F.values) - retmat = _psd_spectral_product!(F.values, F.vectors) + rootexp = map(λ -> exp(0.5λ), F.values) + retmat = _psd_spectral_product(rootexp, F.vectors) return wrappertype(A)(retmat) end @@ -943,8 +943,8 @@ function sqrt(A::SelfAdjoint; rtol = eps(real(float(eltype(A)))) * size(A, 1)) F = eigen(A) λ₀ = -maximum(abs, F.values) * rtol # treat λ ≥ λ₀ as "zero" eigenvalues up to roundoff if all(λ -> λ ≥ λ₀, F.values) - map!(λ -> λ < 0 ? zero(λ) : fourthroot(λ), F.values) - retmat = _psd_spectral_product!(F.values, F.vectors) + rootroot = map(λ -> λ < 0 ? zero(λ) : fourthroot(λ), F.values) + retmat = _psd_spectral_product(rootroot, F.vectors) return wrappertype(A)(retmat) else retmat = (F.vectors * Diagonal(sqrt.(complex.(F.values)))) * F.vectors' From 634f7a5c5d576eff7419f36a72c673e54114460c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mateus=20Ara=C3=BAjo?= Date: Thu, 15 May 2025 15:52:02 +0200 Subject: [PATCH 4/5] Update src/symmetric.jl Co-authored-by: Steven G. Johnson --- src/symmetric.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/symmetric.jl b/src/symmetric.jl index f3be534f..60a68a15 100644 --- a/src/symmetric.jl +++ b/src/symmetric.jl @@ -836,7 +836,7 @@ end #computes U * Diagonal(abs2.(v)) * U' function _psd_spectral_product(v, U) Uv = U * Diagonal(v) - return Uv * Uv' + return Uv * Uv' # often faster than generic matmul by calling BLAS.herk end # Matrix functions From dcbf39e715cc8419c909f57e3c0845b834b58ace Mon Sep 17 00:00:00 2001 From: araujoms Date: Thu, 15 May 2025 19:22:07 +0200 Subject: [PATCH 5/5] type stability --- src/symmetric.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/symmetric.jl b/src/symmetric.jl index 60a68a15..9187caf1 100644 --- a/src/symmetric.jl +++ b/src/symmetric.jl @@ -854,7 +854,7 @@ function ^(A::SelfAdjoint, p::Real) isinteger(p) && return integerpow(A, p) F = eigen(A) if all(λ -> λ ≥ 0, F.values) - rootpower = map(λ -> λ^0.5p, F.values) + rootpower = map(λ -> λ^(p / 2), F.values) retmat = _psd_spectral_product(rootpower, F.vectors) return wrappertype(A)(retmat) else @@ -879,7 +879,7 @@ end function exp(A::SelfAdjoint) F = eigen(A) - rootexp = map(λ -> exp(0.5λ), F.values) + rootexp = map(λ -> exp(λ / 2), F.values) retmat = _psd_spectral_product(rootexp, F.vectors) return wrappertype(A)(retmat) end