From 0cc0f48222314e9964d6e5b042b30930bf6f8455 Mon Sep 17 00:00:00 2001 From: Mladen Kolar Date: Thu, 22 Oct 2026 16:54:09 +0400 Subject: [PATCH] Allow CoordinateDescent 0.5 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Compat CoordinateDescent "0.4, 0.5"; version 0.2.1 - The second stages copied the options field by field into a new CDOptions, which would drop relativeTol, the field CoordinateDescent 0.5 adds; _warmStart copies every field and sets warmStart - CDInverseKroneckerLoss and CDInverseReducedSymKroneckerLoss define coordinateScale and responseScale (when CoordinateDescent has them), as for its quadratic loss, so the relative stopping rule of CoordinateDescent 0.5 does not depend on the units of Σx and Σy; the curvature is shared with descendCoordinate! - Tests for the copied options, the scales against those of CDQuadraticLoss with the same Hessian, and the invariance of the solution and the number of passes to the scale of Σx and Σy Co-Authored-By: Claude Opus 5.5 --- Project.toml | 4 +- src/diffEstimation.jl | 7 +-- src/invHessianEstimationASym.jl | 24 +++++++---- src/invHessianEstimationSymReduced.jl | 30 ++++++++----- src/reducedDiffEstim.jl | 14 +----- src/support.jl | 7 +-- src/util.jl | 6 +++ test/invHessian.jl | 61 +++++++++++++++++++++++++++ 8 files changed, 108 insertions(+), 45 deletions(-) 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 "] -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