diff --git a/Project.toml b/Project.toml
index d9e2993..580dfca 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "DiffPrecTest"
 uuid = "bb4bfe13-bfd7-5c6f-a2d2-8072376c6806"
 authors = ["Mladen Kolar <mkolar@chicagobooth.edu>"]
-version = "0.2.0"
+version = "0.2.1"
 
 [deps]
 CoordinateDescent = "c00ed354-b9be-54b8-8a23-c5ac771bc43a"
@@ -18,7 +18,7 @@ CovSel = {url = "https://github.com/mlakolar/CovSel.jl", rev = "master"}
 
 [compat]
 Aqua = "0.8"
-CoordinateDescent = "0.4"
+CoordinateDescent = "0.4, 0.5"
 CovSel = "0.2"
 Ipopt = "1"
 JuMP = "1"
diff --git a/src/diffEstimation.jl b/src/diffEstimation.jl
index 92f3913..3a61ab4 100644
--- a/src/diffEstimation.jl
+++ b/src/diffEstimation.jl
@@ -158,12 +158,7 @@ function diffEstimation(Sx::Symmetric, Sy::Symmetric, X, Y, λ, options=CDOption
     ##################
 
     # recompute variance
-    opt1 = CDOptions(;
-        maxIter=options.maxIter,
-        optTol=options.optTol,
-        randomize=options.randomize,
-        warmStart=true,
-        numSteps=options.numSteps)
+    opt1 = _warmStart(options)
 
     _computeVar!(ω, f, X, Y, x)
 
diff --git a/src/invHessianEstimationASym.jl b/src/invHessianEstimationASym.jl
index 8101a1b..74a9bd5 100644
--- a/src/invHessianEstimationASym.jl
+++ b/src/invHessianEstimationASym.jl
@@ -26,6 +26,19 @@ function CDInverseKroneckerLoss(Σx::Symmetric{T,S}, Σy::Symmetric{T,S}, a::Int
 end
 
 CoordinateDescent.numCoordinates(f::CDInverseKroneckerLoss) = f.p*f.p
+
+# second derivative of the loss along the coordinate Θ[r, c]
+@inline _curvature(f::CDInverseKroneckerLoss, r, c) = @inbounds f.Σx[r, r] * f.Σy[c, c]
+
+# Scales for the stopping rule of CoordinateDescent 0.5 and later (relativeTol),
+# as for its quadratic loss: the change of coordinate k is measured as
+# |Δxₖ| √aₖ, with aₖ the curvature, and compared with optTol ⋅ maxₖ |bₖ| / √aₖ,
+# with -bₖ the derivative at zero, the indicator of the entry (a, b).
+if isdefined(CoordinateDescent, :coordinateScale)
+  CoordinateDescent.coordinateScale(f::CDInverseKroneckerLoss, k::Int) =
+    sqrt(_curvature(f, Tuple(CartesianIndices((f.p, f.p))[k])...))
+  CoordinateDescent.responseScale(f::CDInverseKroneckerLoss) = 1 / sqrt(_curvature(f, f.a, f.b))
+end
 function CoordinateDescent.initialize!(f::CDInverseKroneckerLoss, x::SparseIterate)
 # compute residuals for the loss
 
@@ -70,9 +83,7 @@ function CoordinateDescent.descendCoordinate!(
   # ri, ci = ind2sub(Σx, j)
   ri, ci = Tuple(CartesianIndices(Σx)[j])
 
-  a = zero(T)
-  b = zero(T)
-  @inbounds a  = Σx[ri,ri] * Σy[ci,ci]
+  a = _curvature(f, ri, ci)
   @inbounds b = A[ri,ci]
   b = ri == f.a && ci == f.b ? b - 1. : b
 
@@ -204,12 +215,7 @@ function invAsymHessianEstimation(Sx::Symmetric, Sy::Symmetric, ri, ci, X, Y, λ
     ##################
 
     # recompute variance
-    opt1 = CDOptions(;
-        maxIter=options.maxIter,
-        optTol=options.optTol,
-        randomize=options.randomize,
-        warmStart=true,
-        numSteps=options.numSteps)
+    opt1 = _warmStart(options)
 
     _computeVar!(ω, f, X, Y, x)
     ω[(ci-1)*p + ri] = 0.
diff --git a/src/invHessianEstimationSymReduced.jl b/src/invHessianEstimationSymReduced.jl
index 0a7bb0b..aba500f 100644
--- a/src/invHessianEstimationSymReduced.jl
+++ b/src/invHessianEstimationSymReduced.jl
@@ -27,6 +27,24 @@ end
 
 CoordinateDescent.numCoordinates(f::CDInverseReducedSymKroneckerLoss) = div((f.p + 1) * f.p, 2)
 
+# second derivative of the loss along the coordinate for (r, c), r ≥ c: the
+# diagonal entry of H
+@inline function _curvature(f::CDInverseReducedSymKroneckerLoss, r, c)
+  Σx, Σy = f.Σx, f.Σy
+  @inbounds r == c ? Σx[r, r] * Σy[r, r] : (Σx[r, r] * Σy[c, c] + Σx[c, c] * Σy[r, r]) / 2 + Σx[r, c] * Σy[r, c]
+end
+
+# Scales for the stopping rule of CoordinateDescent 0.5 and later (relativeTol),
+# as for its quadratic loss: the change of coordinate k is measured as
+# |Δxₖ| √aₖ, with aₖ the curvature, and compared with optTol ⋅ maxₖ |bₖ| / √aₖ,
+# with -bₖ the derivative at zero, the indicator of ind.
+if isdefined(CoordinateDescent, :coordinateScale)
+  CoordinateDescent.coordinateScale(f::CDInverseReducedSymKroneckerLoss, k::Int) =
+    sqrt(_curvature(f, ind2subLowerTriangular(f.p, k)...))
+  CoordinateDescent.responseScale(f::CDInverseReducedSymKroneckerLoss) =
+    1 / sqrt(_curvature(f, ind2subLowerTriangular(f.p, f.ind)...))
+end
+
 
 function _A_mul_symX_mul_B_rc(
   A::Symmetric{T},
@@ -91,13 +109,10 @@ function CoordinateDescent.descendCoordinate!(
 
   ri, ci = ind2subLowerTriangular(f.p, j)
 
-  a = zero(T)
-  b = zero(T)
+  a = _curvature(f, ri, ci)
   if ri == ci
-    @inbounds a = Σx[ri,ri] * Σy[ci,ci]
     @inbounds b = A[ri,ri]
   else
-    @inbounds a = (Σx[ri,ri] * Σy[ci,ci] + Σx[ci,ci] * Σy[ri,ri]) / 2. + Σx[ri,ci] * Σy[ri,ci]
     @inbounds b = (A[ri,ci] + A[ci,ri]) / sqrt(2.)
   end
   b = j == f.ind ? b - 1. : b
@@ -258,12 +273,7 @@ function invQSymHessian(Sx::Symmetric, Sy::Symmetric, ind, X, Y, λ, options=CDO
     ##################
 
     # recompute variance
-    opt1 = CDOptions(;
-        maxIter=options.maxIter,
-        optTol=options.optTol,
-        randomize=options.randomize,
-        warmStart=true,
-        numSteps=options.numSteps)
+    opt1 = _warmStart(options)
 
     _computeReducedVar!(ω, f, X, Y, x)
     ω[ind] = 0.
diff --git a/src/reducedDiffEstim.jl b/src/reducedDiffEstim.jl
index 4bebfff..3484fe8 100644
--- a/src/reducedDiffEstim.jl
+++ b/src/reducedDiffEstim.jl
@@ -125,12 +125,7 @@ function reducedDiffEstimation(
   ##################
 
   # recompute variance
-  opt1 = CDOptions(;
-      maxIter=options.maxIter,
-      optTol=options.optTol,
-      randomize=options.randomize,
-      warmStart=true,
-      numSteps=options.numSteps)
+  opt1 = _warmStart(options)
 
   _computeVarStep1!(ω, Sx, Sy, X, Y, x)
 
@@ -243,12 +238,7 @@ function invHessianReduced(H::Symmetric,
   ##################
 
   # recompute variance
-  opt1 = CDOptions(;
-      maxIter=options.maxIter,
-      optTol=options.optTol,
-      randomize=options.randomize,
-      warmStart=true,
-      numSteps=options.numSteps)
+  opt1 = _warmStart(options)
 
   _computeVarStep2!(ω, indRow, Sx, Sy, X, Y, x)
   ω[indRow] = 0.
diff --git a/src/support.jl b/src/support.jl
index d72e4e7..b3eaeb6 100644
--- a/src/support.jl
+++ b/src/support.jl
@@ -276,12 +276,7 @@ function supportEstimate(::DTraceValidationSupport,
     loss2arr[1] = CovSel.diffLoss(SxValid, x, SyValid, 2)
     lossInfarr[1] = CovSel.diffLoss(SxValid, x, SyValid, Inf)
 
-    opt = CDOptions(;
-        maxIter=options.maxIter,
-        optTol=options.optTol,
-        randomize=options.randomize,
-        warmStart=true,
-        numSteps=options.numSteps)
+    opt = _warmStart(options)
 
     for i=2:numΛ
         g = ProxL1(Λarr[i])
diff --git a/src/util.jl b/src/util.jl
index 9924c46..ac7a41b 100644
--- a/src/util.jl
+++ b/src/util.jl
@@ -1,3 +1,9 @@
+# the options with warmStart = true, for the second stage of the estimators;
+# copies every field, including those that later versions of CoordinateDescent
+# add (relativeTol in 0.5)
+_warmStart(options::CDOptions) =
+  CDOptions((f === :warmStart ? true : getfield(options, f) for f in fieldnames(CDOptions))...)
+
 # quantile function of the standard normal distribution
 _normQuantile(α::Real) = -sqrt(2.) * erfcinv(2. * α)
 
diff --git a/test/invHessian.jl b/test/invHessian.jl
index 1127f23..8e7357f 100644
--- a/test/invHessian.jl
+++ b/test/invHessian.jl
@@ -126,4 +126,65 @@ end
   end
 end
 
+@testset "warm start options" begin
+  opt = CDOptions(; maxIter=7, optTol=1e-3, randomize=false, warmStart=false, numSteps=3)
+  w = DiffPrecTest._warmStart(opt)
+  @test w.warmStart
+  for f in fieldnames(CDOptions)
+    f === :warmStart || @test getfield(w, f) == getfield(opt, f)
+  end
+  if hasfield(CDOptions, :relativeTol)   # CoordinateDescent 0.5 and later
+    @test !DiffPrecTest._warmStart(CDOptions(; relativeTol=false)).relativeTol
+  end
+end
+
+# CoordinateDescent 0.5 and later measure convergence relative to scales of the loss
+if isdefined(CoordinateDescent, :coordinateScale)
+  @testset "stopping rule scales" begin
+    p = 5
+    X, Y = genData(p)
+    hSx, hSy = Symmetric(cov(X)), Symmetric(cov(Y))
+
+    # the scales are those of the quadratic loss with the same Hessian and linear term
+    A = kron(hSy, hSx)
+    t = (2-1)*p + 3
+    b = zeros(p*p)
+    b[t] = -1.
+    f = CDInverseKroneckerLoss(hSx, hSy, 3, 2)
+    fq = CDQuadraticLoss(A, b)
+    @test [CoordinateDescent.coordinateScale(f, k) for k=1:p*p] ≈ sqrt.(diag(A))
+    @test [CoordinateDescent.coordinateScale(f, k) for k=1:p*p] ≈ [CoordinateDescent.coordinateScale(fq, k) for k=1:p*p]
+    @test CoordinateDescent.responseScale(f) ≈ CoordinateDescent.responseScale(fq) ≈ 1 / sqrt(A[t, t])
+
+    H = DiffPrecTest.skron(hSx, hSy)
+    m = size(H, 1)
+    ind = sub2indLowerTriangular(p, 3, 2)
+    b = zeros(m)
+    b[ind] = -1.
+    f = CDInverseReducedSymKroneckerLoss(hSx, hSy, ind)
+    fq = CDQuadraticLoss(H, b)
+    @test [CoordinateDescent.coordinateScale(f, k) for k=1:m] ≈ sqrt.(diag(H))
+    @test [CoordinateDescent.coordinateScale(f, k) for k=1:m] ≈ [CoordinateDescent.coordinateScale(fq, k) for k=1:m]
+    @test CoordinateDescent.responseScale(f) ≈ CoordinateDescent.responseScale(fq) ≈ 1 / sqrt(H[ind, ind])
+
+    # so multiplying Σx and Σy by a divides the solution by a² and takes the same
+    # number of passes (powers of 2 keep the scaling exact)
+    for (loss, m) in (((Sx, Sy) -> CDInverseKroneckerLoss(Sx, Sy, 3, 2), p*p),
+                      ((Sx, Sy) -> CDInverseReducedSymKroneckerLoss(Sx, Sy, ind), m))
+      passes, sols = Int[], Vector{Float64}[]
+      for a in (0.25, 1.0, 4.0)
+        info = CoordinateDescent.CDInfo()
+        x = SparseIterate(m)
+        coordinateDescent!(x, loss(Symmetric(a * hSx), Symmetric(a * hSy)), ProxL1(0.05), CDOptions(; randomize=false), info)
+        @test info.converged
+        push!(passes, info.iterations)
+        push!(sols, a^2 * Vector(x))
+      end
+      @test passes[1] == passes[2] == passes[3]
+      @test sols[1] ≈ sols[2]
+      @test sols[3] ≈ sols[2]
+    end
+  end
+end
+
 end