diff --git a/ext/PEPSKitGPUArraysExt.jl b/ext/PEPSKitGPUArraysExt.jl
index 13510fd5c..e76b0a3f2 100644
--- a/ext/PEPSKitGPUArraysExt.jl
+++ b/ext/PEPSKitGPUArraysExt.jl
@@ -12,4 +12,16 @@ function PEPSKit._column_vectors(::Type{A}, m::AbstractMatrix) where {A <: AnyGP
 end
 PEPSKit._indexable_diagview(m::AnyGPUArray) = collect(MatrixAlgebraKit.diagview(m))
 
+# Reduce into a 0-dimensional device array instead of returning a host scalar. `sdiag_pow` only
+# feeds this into a broadcast, and a 0-dim array broadcasts as a scalar, so the value never has to
+# come back to the CPU. Returning a scalar here would force a device sync on every call, and
+# `sdiag_pow` runs once per bond per weight absorption in simple update.
+function PEPSKit._maxabs(data::AnyGPUArray)
+    T = real(eltype(data))
+    acc = similar(data, T, ())
+    fill!(acc, zero(T))
+    Base.mapreducedim!(abs, max, acc, data)
+    return acc
+end
+
 end
diff --git a/src/utility/util.jl b/src/utility/util.jl
index df25d2ae5..c3bf03a25 100644
--- a/src/utility/util.jl
+++ b/src/utility/util.jl
@@ -37,7 +37,7 @@ function sdiag_pow(s::DiagonalTensorMap, pow::Real; tol::Real = eps(real(scalart
     return spow
 end
 function sdiag_pow(
-        s::AbstractTensorMap{T, S, 1, 1}, pow::Real; tol::Real = eps(real(scalartype(s)))^(3 / 4)
+       s::AbstractTensorMap{T, S, 1, 1}, pow::Real; tol::Real = eps(real(scalartype(s)))^(3 / 4)
     ) where {T, S}
     tol *= _maxabs(s)
     spow = similar(s)
@@ -53,7 +53,7 @@ function ChainRulesCore.rrule(
         ::typeof(sdiag_pow), s::AbstractTensorMap, pow::Real;
         tol::Real = eps(real(scalartype(s)))^(3 / 4),
     )
-    tol *= _maxabs(s)
+    tol *= norm(s, Inf)
     spow = sdiag_pow(s, pow; tol)
     spow_minus1_conj = scale!(sdiag_pow(s', pow - 1; tol), pow)
     function sdiag_pow_pullback(c̄_)