From 1e982b26ca3e7a0b08e14f491465faa94fe2486f Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sat, 12 Sep 2026 22:23:19 -0700 Subject: [PATCH 01/15] cleaned up the Accelerated AP algorithm --- src/NearestCorrelationMatrix.jl | 8 +- ...atingProjectionsAA.jl => AcceleratedAP.jl} | 79 +++++++++---------- 2 files changed, 43 insertions(+), 44 deletions(-) rename src/algorithms/{AlternatingProjectionsAA.jl => AcceleratedAP.jl} (56%) diff --git a/src/NearestCorrelationMatrix.jl b/src/NearestCorrelationMatrix.jl index 79dc638..fa23525 100644 --- a/src/NearestCorrelationMatrix.jl +++ b/src/NearestCorrelationMatrix.jl @@ -20,6 +20,7 @@ include("simple_interface.jl") include("algorithms/Newton.jl") include("algorithms/DirectProjection.jl") include("algorithms/AlternatingProjections.jl") +include("algorithms/AcceleratedAP.jl") include("algorithms/JuMPAlgorithm.jl") export @@ -39,9 +40,10 @@ export nearest_cor, nearest_cor!, # algorithms - Newton, - DirectProjection, + AcceleratedAP, AlternatingProjections, - JuMPAlgorithm + DirectProjection, + JuMPAlgorithm, + Newton end diff --git a/src/algorithms/AlternatingProjectionsAA.jl b/src/algorithms/AcceleratedAP.jl similarity index 56% rename from src/algorithms/AlternatingProjectionsAA.jl rename to src/algorithms/AcceleratedAP.jl index 66eeed7..022d93d 100644 --- a/src/algorithms/AlternatingProjectionsAA.jl +++ b/src/algorithms/AcceleratedAP.jl @@ -1,33 +1,40 @@ -struct AlternatingProjectionsAA{A, K} <: NCMAlgorithm +""" + AcceleratedAP(; tau=0, m=2) + +The alternating projections algorithm with Anderson acceleration applied. Should converge +in roughly half the number of steps of the standard alternating projections algorithm. +""" +struct AcceleratedAP{A, K} <: NCMAlgorithm tau::Real m::Int args::A kwargs::K end -function AlternatingProjectionsAA(args...; tau::Real = 0, m::Int = 2, kwargs...) - return AlternatingProjectionsAA(tau, m, args, kwargs) +function AcceleratedAP(args...; tau::Real = 0, m::Int = 2, kwargs...) + return AcceleratedAP(tau, m, args, kwargs) end -default_iters(::AlternatingProjectionsAA, A) = clamp(size(A, 1), 20, 200) -modifies_in_place(::AlternatingProjectionsAA) = true -supports_float16(::AlternatingProjectionsAA) = true -supports_symmetric(::AlternatingProjectionsAA) = false -supports_parameterless_construction(::Type{AlternatingProjectionsAA}) = true +default_iters(::AcceleratedAP, A) = clamp(size(A, 1), 20, 200) +modifies_in_place(::AcceleratedAP) = true +supports_float16(::AcceleratedAP) = true +supports_symmetric(::AcceleratedAP) = false +supports_parameterless_construction(::Type{AcceleratedAP}) = true +supports_mask(::AcceleratedAP) = true -function autotune(::Type{AlternatingProjectionsAA}, prob::NCMProblem) - return AlternatingProjectionsAA(; tau = eps(eltype(prob.A)), m = 2) +function autotune(::Type{AcceleratedAP}, prob::NCMProblem) + return AcceleratedAP(; tau = eps(eltype(prob.A)), m = 2) end -function CommonSolve.solve!(solver::NCMSolver, alg::AlternatingProjectionsAA; kwargs...) +function CommonSolve.solve!(solver::NCMSolver, alg::AcceleratedAP; kwargs...) A = solver.A n = size(A, 1) - size(A, 2) == n || throw(DimensionMismatch("Input matrix A must be square.")) T = eltype(A) m = alg.m - tol = solver.reltol - maxiter = solver.maxiters + tau = convert(T, alg.tau) + mask = solver.mask + A_orig = solver.A_orig # Initialize working matrices X = copy(A) @@ -35,6 +42,7 @@ function CommonSolve.solve!(solver::NCMSolver, alg::AlternatingProjectionsAA; kw S = zeros(T, n, n) R = similar(A) G = similar(A) + scratch = similar(A) # Pre-allocate memory for Anderson Acceleration history vec_dim = n * n @@ -47,36 +55,26 @@ function CommonSolve.solve!(solver::NCMSolver, alg::AlternatingProjectionsAA; kw m_eff = 0 iter = 0 - converged = false - rel_err = 0.0 + resid = Inf - while iter < maxiter + while iter < solver.maxiters iter += 1 - # R = Y - S - @. R = Y - S - - # X = P_S(R) : Project onto Positive Semidefinite Cone S+ - X .= project_s(R) - - # S = X - R : Update Dykstra correction - @. S = X - R - - # G = P_U(X) : Project onto Unit Diagonal U + R .= Y .- S + project_psd!(X, R, tau, scratch) + S .= X .- R copyto!(G, X) - for i in 1:n - G[i, i] = one(T) + + if mask !== nothing + project_fixed!(G, A_orig, mask) end - # Relative residual error check - rel_err = norm(X .- G, 2) / max(one(T), norm(X, 2)) + # project unit after projecting fixed to ensure that the unit diagonal is preserved + project_unit!(G) - if solver.verbose - println("Iter $iter: rel_err = $rel_err") - end + resid = norm(X .- G) / norm(X) - if rel_err <= tol - converged = true + if resid <= solver.reltol break end @@ -131,11 +129,10 @@ function CommonSolve.solve!(solver::NCMSolver, alg::AlternatingProjectionsAA; kw end end - # Ensure output X strictly satisfies unit diagonal and exact symmetry - for i in 1:n - X[i, i] = one(T) + if mask !== nothing + project_fixed!(Y, A_orig, mask) end - X .= (X .+ X') ./ 2 + project_unit!(X) - return build_ncm_solution(alg, X, rel_err, solver; iters = iter) + return build_ncm_solution(alg, X, resid, solver; iters = iter) end From 5c280a430224510d0a65deb74ffeb6a392fdb0ac Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sat, 12 Sep 2026 22:23:32 -0700 Subject: [PATCH 02/15] small changes to AP method --- src/algorithms/AlternatingProjections.jl | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/algorithms/AlternatingProjections.jl b/src/algorithms/AlternatingProjections.jl index 64a3d46..60002aa 100644 --- a/src/algorithms/AlternatingProjections.jl +++ b/src/algorithms/AlternatingProjections.jl @@ -54,7 +54,9 @@ function CommonSolve.solve!(solver::NCMSolver, alg::AlternatingProjections; kwar iter = 0 resid = Inf - while iter < solver.maxiters && resid ≥ solver.reltol + while iter < solver.maxiters + iter += 1 + R .= Y .- ΔS project_psd!(X, R, tau, scratch) ΔS .= X .- R @@ -68,7 +70,10 @@ function CommonSolve.solve!(solver::NCMSolver, alg::AlternatingProjections; kwar project_unit!(Y) resid = norm(Y .- X) / norm(Y) - iter += 1 + + if resid ≤ solver.reltol + break + end end return build_ncm_solution(alg, Y, resid, solver; iters = iter) From 9426516ea143c8f9760da8115853eb0eb7ba786e Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sat, 12 Sep 2026 22:28:58 -0700 Subject: [PATCH 03/15] changed default algorithm to AcceleratedAP when a mask is supplied --- src/NCMSolver.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/NCMSolver.jl b/src/NCMSolver.jl index a22bc69..85860a2 100644 --- a/src/NCMSolver.jl +++ b/src/NCMSolver.jl @@ -38,7 +38,7 @@ end Get the default algorithm type for a given input matrix. """ -default_algtype(prob::NCMProblem) = prob.mask === nothing ? Newton : AlternatingProjections +default_algtype(prob::NCMProblem) = prob.mask === nothing ? Newton : AcceleratedAP """ init(prob, alg, args...; kwargs...) From 1ffce9022dcacfe473384a4ee755744c28428443 Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sat, 12 Sep 2026 22:29:24 -0700 Subject: [PATCH 04/15] Added residual calculation to DirectProjection --- src/algorithms/DirectProjection.jl | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/src/algorithms/DirectProjection.jl b/src/algorithms/DirectProjection.jl index e8a96bf..f567a85 100644 --- a/src/algorithms/DirectProjection.jl +++ b/src/algorithms/DirectProjection.jl @@ -1,5 +1,5 @@ """ - DirectProjection(; tau=eps()) + DirectProjection(args...; tau=0, kwargs...) Single step projection of the input matrix into the set of correlation matrices. Useful when a "close" correlation matrix is needed without concern for it being the most optimal. @@ -11,10 +11,10 @@ struct DirectProjection{A, K} <: NCMAlgorithm tau::Real args::A kwargs::K -end -function DirectProjection(args...; tau::Real = 0, kwargs...) - return DirectProjection(tau, args, kwargs) + function DirectProjection(args...; tau::Real = 0, kwargs...) + return DirectProjection(tau, args, kwargs) + end end modifies_in_place(::DirectProjection) = true @@ -24,7 +24,7 @@ supports_parameterless_construction(::Type{DirectProjection}) = true autotune(::Type{DirectProjection}, prob::NCMProblem) = _autotune(DirectProjection, prob.A) -function _autotune(::Type{DirectProjection}, A::AbstractMatrix{Float64}) +function _autotune(::Type{DirectProjection}, ::AbstractMatrix{Float64}) return DirectProjection(; tau = 1.0e-12) end @@ -61,12 +61,14 @@ function _autotune(::Type{DirectProjection}, A::AbstractMatrix{Float16}) end function CommonSolve.solve!(solver::NCMSolver, alg::DirectProjection; kwargs...) - X = solver.A - T = eltype(X) - tau = max(T(alg.tau), zero(T)) + A = solver.A + X = copy(A) + tau = convert(eltype(X), alg.tau) project_psd!(X, tau) cov2cor!(X) - return build_ncm_solution(alg, X, nothing, solver; iters = 1) + resid = norm(X .- A) / norm(X) + + return build_ncm_solution(alg, X, resid, solver; iters = 1) end From f7de12bacf66d09cea14811fd647540feae25da6 Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sat, 12 Sep 2026 22:34:58 -0700 Subject: [PATCH 05/15] moved constructor outside of struct definition --- src/algorithms/DirectProjection.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/algorithms/DirectProjection.jl b/src/algorithms/DirectProjection.jl index f567a85..7a2b383 100644 --- a/src/algorithms/DirectProjection.jl +++ b/src/algorithms/DirectProjection.jl @@ -12,9 +12,10 @@ struct DirectProjection{A, K} <: NCMAlgorithm args::A kwargs::K - function DirectProjection(args...; tau::Real = 0, kwargs...) - return DirectProjection(tau, args, kwargs) - end +end + +function DirectProjection(args...; tau::Real = 0, kwargs...) + return DirectProjection(tau, args, kwargs) end modifies_in_place(::DirectProjection) = true From a6abc02ebcf4f3085516514b7573c6d42674f949 Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sat, 12 Sep 2026 23:26:01 -0700 Subject: [PATCH 06/15] Added unit tests for masking --- src/NCMAlgorithm.jl | 7 +++--- src/algorithms/AcceleratedAP.jl | 18 +++++++++++---- src/algorithms/AlternatingProjections.jl | 15 ++++++++---- test/runtests.jl | 1 + test/test_masking.jl | 29 ++++++++++++++++++++++++ 5 files changed, 58 insertions(+), 12 deletions(-) create mode 100644 test/test_masking.jl diff --git a/src/NCMAlgorithm.jl b/src/NCMAlgorithm.jl index 1ef7063..49f74e8 100644 --- a/src/NCMAlgorithm.jl +++ b/src/NCMAlgorithm.jl @@ -52,14 +52,15 @@ default_tol(::Type{<:Rational}) = 0 default_tol(::Type{<:Integer}) = 0 """ - supports_mask(alg) + supports_mask(algtype) Trait for whether an algorithm supports a fixed-element mask. When `true`, the mask is enforced during the solve (for alternating projections, this means projecting onto the fixed-element subspace each iteration). When `false`, passing a mask throws an informative error. Default is `false`. """ -supports_mask(::NCMAlgorithm) = false +supports_mask(::Type{T}) where {T<:NCMAlgorithm} = false +supports_mask(alg::T) where {T<:NCMAlgorithm} = supports_mask(T) """ default_iters(alg, A) @@ -100,7 +101,7 @@ If `false`, then a copy using the upper or lower matrix is used instead. supports_symmetric(::NCMAlgorithm) = false """ - supports_parameterless_construction(alg) + supports_parameterless_construction(algtype) Trait for if an algorithm can be constructed without any parameters (default is `false`). """ diff --git a/src/algorithms/AcceleratedAP.jl b/src/algorithms/AcceleratedAP.jl index 022d93d..487ee41 100644 --- a/src/algorithms/AcceleratedAP.jl +++ b/src/algorithms/AcceleratedAP.jl @@ -19,11 +19,19 @@ default_iters(::AcceleratedAP, A) = clamp(size(A, 1), 20, 200) modifies_in_place(::AcceleratedAP) = true supports_float16(::AcceleratedAP) = true supports_symmetric(::AcceleratedAP) = false -supports_parameterless_construction(::Type{AcceleratedAP}) = true -supports_mask(::AcceleratedAP) = true +supports_parameterless_construction(::Type{<:AcceleratedAP}) = true +supports_mask(::Type{<:AcceleratedAP}) = true -function autotune(::Type{AcceleratedAP}, prob::NCMProblem) - return AcceleratedAP(; tau = eps(eltype(prob.A)), m = 2) +function autotune(::Type{<:AcceleratedAP}, prob::NCMProblem) + T = eltype(prob.A) + tau = eps(T) + + if prob.mask !== nothing + tau = 10 * sqrt(tau) + end + + # if the problem implements masking, then err on the safe side for tau + return AcceleratedAP(; tau = tau, m = 2) end function CommonSolve.solve!(solver::NCMSolver, alg::AcceleratedAP; kwargs...) @@ -130,7 +138,7 @@ function CommonSolve.solve!(solver::NCMSolver, alg::AcceleratedAP; kwargs...) end if mask !== nothing - project_fixed!(Y, A_orig, mask) + project_fixed!(X, A_orig, mask) end project_unit!(X) diff --git a/src/algorithms/AlternatingProjections.jl b/src/algorithms/AlternatingProjections.jl index 60002aa..61b12bb 100644 --- a/src/algorithms/AlternatingProjections.jl +++ b/src/algorithms/AlternatingProjections.jl @@ -32,11 +32,18 @@ default_iters(::AlternatingProjections, A) = clamp(size(A, 1), 20, 200) modifies_in_place(::AlternatingProjections) = true supports_float16(::AlternatingProjections) = true supports_symmetric(::AlternatingProjections) = false -supports_parameterless_construction(::Type{AlternatingProjections}) = true -supports_mask(::AlternatingProjections) = true +supports_parameterless_construction(::Type{<:AlternatingProjections}) = true +supports_mask(::Type{<:AlternatingProjections}) = true -function autotune(::Type{AlternatingProjections}, prob::NCMProblem) - return AlternatingProjections(; tau = eps(eltype(prob.A))) +function autotune(::Type{<:AlternatingProjections}, prob::NCMProblem) + T = eltype(prob.A) + tau = eps(T) + + if prob.mask !== nothing + tau = 15 * sqrt(tau) + end + + return AlternatingProjections(; tau = tau) end function CommonSolve.solve!(solver::NCMSolver, alg::AlternatingProjections; kwargs...) diff --git a/test/runtests.jl b/test/runtests.jl index 3c9e63e..b38aa9b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,6 +15,7 @@ include("datadeps_registration.jl") # Algorithm Robustnes @safetestset "Constructors" include("test_constructors.jl") @safetestset "Convergence" include("test_convergence.jl") +@safetestset "Fixed Element Masking" include("test_masking.jl") @safetestset "Real World Data" include("test_real_world.jl") # Extension Packages diff --git a/test/test_masking.jl b/test/test_masking.jl new file mode 100644 index 0000000..1a0405b --- /dev/null +++ b/test/test_masking.jl @@ -0,0 +1,29 @@ +using Test +using InteractiveUtils +using NearestCorrelationMatrix +import NearestCorrelationMatrix as NCM + +include("Datasets.jl") +using .Datasets + +include("CustomTestMacros.jl") +using .CustomTestMacros + +masking_algs = filter(alg -> NCM.supports_mask(alg), subtypes(NCMAlgorithm)) +supported_types = (Float64, Float32, Float16) + +for algtype in masking_algs, T in supported_types + @testset "$(NCM.alg_name(algtype)) - $T" begin + A, m = usgs13() + A = convert(Matrix{T}, A) + X = copy(A) + prob = NCMProblem(X; mask = m) + alg = autotune(algtype, prob) + solver = init(prob, alg) + @test solver.mask !== nothing + @test Base.mightalias(solver.A_orig, A) == false + sol = solve!(solver) + @test_iscorrelation sol.X + @test sol.X[m] == A[m] + end +end From 13fad5950fb0aee9644b33b271543d960fd2b956 Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sat, 12 Sep 2026 23:27:57 -0700 Subject: [PATCH 07/15] fixed formatting --- src/NCMAlgorithm.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/NCMAlgorithm.jl b/src/NCMAlgorithm.jl index 49f74e8..38ebf26 100644 --- a/src/NCMAlgorithm.jl +++ b/src/NCMAlgorithm.jl @@ -59,8 +59,8 @@ during the solve (for alternating projections, this means projecting onto the fi subspace each iteration). When `false`, passing a mask throws an informative error. Default is `false`. """ -supports_mask(::Type{T}) where {T<:NCMAlgorithm} = false -supports_mask(alg::T) where {T<:NCMAlgorithm} = supports_mask(T) +supports_mask(::Type{T}) where {T <: NCMAlgorithm} = false +supports_mask(alg::T) where {T <: NCMAlgorithm} = supports_mask(T) """ default_iters(alg, A) From 2d6f7afbe91faeba29f9e55d5d7af36a54b91925 Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sat, 12 Sep 2026 23:47:12 -0700 Subject: [PATCH 08/15] updated autotune --- src/algorithms/AcceleratedAP.jl | 6 +++--- src/algorithms/AlternatingProjections.jl | 5 +++-- src/algorithms/Newton.jl | 18 +++++++++++++++++- 3 files changed, 23 insertions(+), 6 deletions(-) diff --git a/src/algorithms/AcceleratedAP.jl b/src/algorithms/AcceleratedAP.jl index 487ee41..d156937 100644 --- a/src/algorithms/AcceleratedAP.jl +++ b/src/algorithms/AcceleratedAP.jl @@ -24,13 +24,13 @@ supports_mask(::Type{<:AcceleratedAP}) = true function autotune(::Type{<:AcceleratedAP}, prob::NCMProblem) T = eltype(prob.A) - tau = eps(T) + tau = sqrt(eps(T)) + # if the problem implements masking, then err on the safe side for tau if prob.mask !== nothing - tau = 10 * sqrt(tau) + tau = 15 * tau end - # if the problem implements masking, then err on the safe side for tau return AcceleratedAP(; tau = tau, m = 2) end diff --git a/src/algorithms/AlternatingProjections.jl b/src/algorithms/AlternatingProjections.jl index 61b12bb..495df71 100644 --- a/src/algorithms/AlternatingProjections.jl +++ b/src/algorithms/AlternatingProjections.jl @@ -37,10 +37,11 @@ supports_mask(::Type{<:AlternatingProjections}) = true function autotune(::Type{<:AlternatingProjections}, prob::NCMProblem) T = eltype(prob.A) - tau = eps(T) + tau = sqrt(eps(T)) + # if the problem implements masking, then err on the safe side for tau if prob.mask !== nothing - tau = 15 * sqrt(tau) + tau = 15 * tau end return AlternatingProjections(; tau = tau) diff --git a/src/algorithms/Newton.jl b/src/algorithms/Newton.jl index 9b19994..cc7e31f 100644 --- a/src/algorithms/Newton.jl +++ b/src/algorithms/Newton.jl @@ -33,7 +33,23 @@ function Newton( end autotune(::Type{Newton}, prob::NCMProblem) = _autotune(Newton, prob.A) -_autotune(::Type{Newton}, A::AbstractMatrix{Float64}) = Newton(; tau = 1.0e-12) +function _autotune(::Type{Newton}, A::AbstractMatrix{Float64}) + n = size(A, 1) + + tau = if n ≤ 50 + 1.0e-12 + elseif n ≤ 100 + 1.0e-11 + elseif n ≤ 500 + 1.0e-10 + elseif n ≤ 1000 + 1.0e-8 + else + 1.0e-6 + end + + return Newton(; tau = tau) +end function _autotune(::Type{Newton}, A::AbstractMatrix{Float32}) n = size(A, 1) From f522bb6142203a7a063a71eda9cb422322b40e20 Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sun, 13 Sep 2026 00:18:29 -0700 Subject: [PATCH 09/15] changed method for ensuring PD to also call cov2cor at end --- src/NCMSolution.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/NCMSolution.jl b/src/NCMSolution.jl index 87af951..2c72573 100644 --- a/src/NCMSolution.jl +++ b/src/NCMSolution.jl @@ -83,16 +83,17 @@ function CommonSolve.solve!(solver::NCMSolver, args...; kwargs...) sol = solve!(solver, solver.alg, args...; kwargs...) if sol.solver.ensure_pd && !isposdef(sol.X) + project_psd!(sol.X, sqrt(eps(eltype(sol.X)))) + # Strict PD and exact fixed-element feasibility cannot both be guaranteed: repairing # definiteness perturbs every entry, so re-apply the mask afterwards. The fixed elements # (and unit diagonal) take precedence - the result is PD up to O(√eps). - project_psd!(sol.X, sqrt(eps(eltype(sol.X)))) - if sol.solver.mask !== nothing project_fixed!(sol.X, sol.solver.A_orig, sol.solver.mask) + project_unit!(sol.X) + else + cov2cor!(sol.X) end - - project_unit!(sol.X) end return sol From 93a7485e07cd3854cacd58fcb85736ef4209c9e3 Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sun, 13 Sep 2026 00:18:51 -0700 Subject: [PATCH 10/15] updated real world test structure --- test/test_real_world.jl | 221 +++++++++++++++++++++------------------- 1 file changed, 119 insertions(+), 102 deletions(-) diff --git a/test/test_real_world.jl b/test/test_real_world.jl index d9339af..8270128 100644 --- a/test/test_real_world.jl +++ b/test/test_real_world.jl @@ -1,5 +1,7 @@ using Test +using InteractiveUtils using NearestCorrelationMatrix +import NearestCorrelationMatrix as NCM using NearestCorrelationMatrix.Internals: iscorrelation include("Datasets.jl") @@ -8,106 +10,121 @@ using .Datasets include("CustomTestMacros.jl") using .CustomTestMacros -@testset "BCCD16" begin - A = bccd16() - @test !iscorrelation(A) - nearest_cor!(A) - @test_iscorrelation A -end - -@testset "BEYU11" begin - A = beyu11() - @test !iscorrelation(A) - nearest_cor!(A) - @test_iscorrelation A -end - -@testset "BHWI01" begin - A = bhwi01() - @test !iscorrelation(A) - nearest_cor!(A) - @test_iscorrelation A -end - -@testset "COR1399" begin - A = cor1399() - @test !iscorrelation(A) - nearest_cor!(A) - @test_iscorrelation A -end - -@testset "COR3120" begin - A = cor3120() - @test !iscorrelation(A) - nearest_cor!(A) - @test_iscorrelation A -end - -@testset "FING97" begin - # unmasked - A, _ = fing97() - @test !iscorrelation(A) - nearest_cor!(A) - @test_iscorrelation A - - A, m = fing97() - @test !iscorrelation(A) - B = nearest_cor(A, AlternatingProjections(); mask = m) - @test_iscorrelation B - @test B[m] == A[m] -end - -@testset "HIGH02" begin - A = high02() - @test !iscorrelation(A) - nearest_cor!(A) - @test_iscorrelation A -end - -@testset "MMB13" begin - A = mmb13() - @test !iscorrelation(A) - nearest_cor!(A) - @test_iscorrelation A -end - -@testset "TEC03" begin - A = tec03() - @test !iscorrelation(A) - nearest_cor!(A) - @test_iscorrelation A -end - -@testset "TYDA99R1" begin - A = tyda99r1() - @test !iscorrelation(A) - nearest_cor!(A) - @test_iscorrelation A -end - -@testset "TYDA99R2" begin - A = tyda99r2() - @test !iscorrelation(A) - nearest_cor!(A) - @test_iscorrelation A -end - -@testset "TYDA99R3" begin - A = tyda99r3() - @test !iscorrelation(A) - nearest_cor!(A) - @test_iscorrelation A -end - -@testset "USGS13" begin - A, _ = usgs13() - @test !iscorrelation(A) - nearest_cor!(A) - @test_iscorrelation A - - A, m = usgs13() - @test !iscorrelation(A) - B = nearest_cor(A, AlternatingProjections(tau = 1.0e-6); mask = m) - @test_iscorrelation B - @test B[m] == A[m] +internal_algtypes = setdiff(subtypes(NCMAlgorithm), (JuMPAlgorithm,)) +fast_algtypes = (Newton, DirectProjection) + +for alg in internal_algtypes + @info "Working on $alg" + + @testset "$(NCM.alg_name(alg))" begin + @testset "BEYU11" begin + A = beyu11() + @test !iscorrelation(A) + nearest_cor!(A, alg) + @test_iscorrelation A + end + + @testset "BHWI01" begin + A = bhwi01() + @test !iscorrelation(A) + nearest_cor!(A, alg) + @test_iscorrelation A + end + + @testset "COR1399" begin + A = cor1399() + @test !iscorrelation(A) + nearest_cor!(A, alg) + @test_iscorrelation A + end + + if alg ∈ fast_algtypes + @testset "BCCD16" begin + A = bccd16() + @test !iscorrelation(A) + nearest_cor!(A, alg) + @test_iscorrelation A + end + + @testset "COR3120" begin + A = cor3120() + @test !iscorrelation(A) + nearest_cor!(A, alg) + @test_iscorrelation A + end + end + + @testset "FING97" begin + # unmasked + A, _ = fing97() + @test !iscorrelation(A) + nearest_cor!(A, alg) + @test_iscorrelation A + + if NCM.supports_mask(alg) + A, m = fing97() + @test !iscorrelation(A) + B = nearest_cor(A, alg; mask = m) + @test_iscorrelation B + @test B[m] == A[m] + end + end + + @testset "HIGH02" begin + A = high02() + @test !iscorrelation(A) + nearest_cor!(A, alg) + @test_iscorrelation A + end + + @testset "MMB13" begin + A = mmb13() + @test !iscorrelation(A) + nearest_cor!(A, alg) + @test_iscorrelation A + end + + @testset "TEC03" begin + A = tec03() + @test !iscorrelation(A) + nearest_cor!(A, alg) + @test_iscorrelation A + end + + @testset "TYDA99R1" begin + A = tyda99r1() + @test !iscorrelation(A) + nearest_cor!(A, alg) + @test_iscorrelation A + end + + @testset "TYDA99R2" begin + A = tyda99r2() + @test !iscorrelation(A) + nearest_cor!(A, alg) + @test_iscorrelation A + end + + @testset "TYDA99R3" begin + A = tyda99r3() + @test !iscorrelation(A) + nearest_cor!(A, alg) + @test_iscorrelation A + end + + @testset "USGS13" begin + A, _ = usgs13() + @test !iscorrelation(A) + nearest_cor!(A, alg) + @test_iscorrelation A + + if NCM.supports_mask(alg) + A, m = usgs13() + @test !iscorrelation(A) + B = nearest_cor(A, alg; mask = m) + @test_iscorrelation B + @test B[m] == A[m] + end + end + end end From 691140fc9452e1587be5d33533d126029b3e5585 Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sun, 13 Sep 2026 00:19:56 -0700 Subject: [PATCH 11/15] updated AP autotune --- src/algorithms/AcceleratedAP.jl | 2 +- src/algorithms/AlternatingProjections.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/algorithms/AcceleratedAP.jl b/src/algorithms/AcceleratedAP.jl index d156937..eccbc29 100644 --- a/src/algorithms/AcceleratedAP.jl +++ b/src/algorithms/AcceleratedAP.jl @@ -28,7 +28,7 @@ function autotune(::Type{<:AcceleratedAP}, prob::NCMProblem) # if the problem implements masking, then err on the safe side for tau if prob.mask !== nothing - tau = 15 * tau + tau = 20 * tau end return AcceleratedAP(; tau = tau, m = 2) diff --git a/src/algorithms/AlternatingProjections.jl b/src/algorithms/AlternatingProjections.jl index 495df71..f9fcc36 100644 --- a/src/algorithms/AlternatingProjections.jl +++ b/src/algorithms/AlternatingProjections.jl @@ -41,7 +41,7 @@ function autotune(::Type{<:AlternatingProjections}, prob::NCMProblem) # if the problem implements masking, then err on the safe side for tau if prob.mask !== nothing - tau = 15 * tau + tau = 20 * tau end return AlternatingProjections(; tau = tau) From 6de1664538d13585e5e818a6797274dea4476607 Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sun, 13 Sep 2026 00:22:24 -0700 Subject: [PATCH 12/15] removed stray info log --- test/test_real_world.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/test_real_world.jl b/test/test_real_world.jl index 8270128..d07bd4c 100644 --- a/test/test_real_world.jl +++ b/test/test_real_world.jl @@ -14,8 +14,6 @@ internal_algtypes = setdiff(subtypes(NCMAlgorithm), (JuMPAlgorithm,)) fast_algtypes = (Newton, DirectProjection) for alg in internal_algtypes - @info "Working on $alg" - @testset "$(NCM.alg_name(alg))" begin @testset "BEYU11" begin A = beyu11() From 87816b98dd250f1dd3d611442d881f02978e27e2 Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sun, 13 Sep 2026 00:52:04 -0700 Subject: [PATCH 13/15] Added min_eigenvalue option to NCMSolver --- src/NCMSolution.jl | 2 +- src/NCMSolver.jl | 33 +++++++++++++++++++++--- src/algorithms/AcceleratedAP.jl | 10 +------ src/algorithms/AlternatingProjections.jl | 10 +------ 4 files changed, 32 insertions(+), 23 deletions(-) diff --git a/src/NCMSolution.jl b/src/NCMSolution.jl index 2c72573..0077752 100644 --- a/src/NCMSolution.jl +++ b/src/NCMSolution.jl @@ -83,7 +83,7 @@ function CommonSolve.solve!(solver::NCMSolver, args...; kwargs...) sol = solve!(solver, solver.alg, args...; kwargs...) if sol.solver.ensure_pd && !isposdef(sol.X) - project_psd!(sol.X, sqrt(eps(eltype(sol.X)))) + project_psd!(sol.X, solver.min_eigenvalue) # Strict PD and exact fixed-element feasibility cannot both be guaranteed: repairing # definiteness perturbs every entry, so re-apply the mask afterwards. The fixed elements diff --git a/src/NCMSolver.jl b/src/NCMSolver.jl index 85860a2..02a67d3 100644 --- a/src/NCMSolver.jl +++ b/src/NCMSolver.jl @@ -14,6 +14,7 @@ Common interface for solving NCM problems. Algorithm-specific cache is stored in - `maxiters`: The number of iterations allowed. Defaults to `size(A,1)` - `ensure_pd`: Checks (and corrects) that the resulting matrix is positive definite. Defaults to `false`. +- `min_eigenvalue`: The minimum eigenvalue to enforce when `ensure_pd` == true. - `verbose`: Whether to print extra information. Defaults to `false`. - `mask`: The fixed-element mask, or `nothing` if unmasked. - `A_orig`: The original values of A to be used when a mask is supplied. @@ -28,6 +29,7 @@ mutable struct NCMSolver{TA, P, Talg, Tc, Ttol, Tm} reltol::Ttol # relative tolerance for convergence maxiters::Int # maximum number of iterations ensure_pd::Bool # ensures that the resulting matrix is positive definite + min_eigenvalue::Union{Nothing, Real} # the minimum eigenvalue to enforce verbose::Bool # whether to print extra information mask::Tm # fixed-element mask, or nothing A_orig::TA # a copy of A, or an alias of A if no mask is given @@ -84,6 +86,7 @@ function CommonSolve.init( convert_f16::Bool = false, force_f16::Bool = false, ensure_pd::Bool = false, + min_eigenvalue = nothing, verbose::Bool = false, kwargs... ) @@ -113,6 +116,8 @@ function CommonSolve.init( A = prob.A p = prob.p + T = eltype(A) + A = if alias_A verbose && println("Aliasing A") A @@ -159,7 +164,7 @@ function CommonSolve.init( end end - if eltype(A) === Float16 && !supports_float16(alg) + if T === Float16 && !supports_float16(alg) if convert_f16 verbose && println( @@ -188,15 +193,35 @@ function CommonSolve.init( A_orig = mask === nothing ? A : copy(A) # Guard against type mismatch for user-specified reltol/abstol - reltol = real(eltype(A))(reltol) - abstol = real(eltype(A))(abstol) + reltol = real(T)(reltol) + reltol = max(reltol, sqrt(eps(T))) + abstol = real(T)(abstol) + abstol = max(abstol, eps(T)) + + min_eigenvalue = if min_eigenvalue === nothing + if ensure_pd + if mask === nothing + # no mask, can default to sqrt(eps(T)) + sqrt(eps(T)) + else + # be more conservative about the min eigenvalue when there is a mask + sqrt(sqrt(eps(T))) + end + else + # no checks for PD -> min_eigenvalue is not used + nothing + end + else + # user explicitly set min_eigenvalue. Just ensure that it is Real + real(T)(min_eigenvalue) + end cacheval = init_cacheval(alg, A; maxiters = maxiters, abstol = abstol, reltol = reltol, verbose = verbose) isfresh = true Tc = typeof(cacheval) solver = NCMSolver{typeof(A), typeof(p), typeof(alg), Tc, typeof(reltol), Union{Nothing, typeof(mask)}}( - A, p, alg, cacheval, isfresh, abstol, reltol, maxiters, ensure_pd, verbose, mask, A_orig + A, p, alg, cacheval, isfresh, abstol, reltol, maxiters, ensure_pd, min_eigenvalue, verbose, mask, A_orig ) return solver diff --git a/src/algorithms/AcceleratedAP.jl b/src/algorithms/AcceleratedAP.jl index eccbc29..5523dde 100644 --- a/src/algorithms/AcceleratedAP.jl +++ b/src/algorithms/AcceleratedAP.jl @@ -23,15 +23,7 @@ supports_parameterless_construction(::Type{<:AcceleratedAP}) = true supports_mask(::Type{<:AcceleratedAP}) = true function autotune(::Type{<:AcceleratedAP}, prob::NCMProblem) - T = eltype(prob.A) - tau = sqrt(eps(T)) - - # if the problem implements masking, then err on the safe side for tau - if prob.mask !== nothing - tau = 20 * tau - end - - return AcceleratedAP(; tau = tau, m = 2) + return AcceleratedAP(; tau = sqrt(eps(eltype(prob.A))), m = 2) end function CommonSolve.solve!(solver::NCMSolver, alg::AcceleratedAP; kwargs...) diff --git a/src/algorithms/AlternatingProjections.jl b/src/algorithms/AlternatingProjections.jl index f9fcc36..e7e9702 100644 --- a/src/algorithms/AlternatingProjections.jl +++ b/src/algorithms/AlternatingProjections.jl @@ -36,15 +36,7 @@ supports_parameterless_construction(::Type{<:AlternatingProjections}) = true supports_mask(::Type{<:AlternatingProjections}) = true function autotune(::Type{<:AlternatingProjections}, prob::NCMProblem) - T = eltype(prob.A) - tau = sqrt(eps(T)) - - # if the problem implements masking, then err on the safe side for tau - if prob.mask !== nothing - tau = 20 * tau - end - - return AlternatingProjections(; tau = tau) + return AlternatingProjections(; tau = sqrt(eps(eltype(prob.A)))) end function CommonSolve.solve!(solver::NCMSolver, alg::AlternatingProjections; kwargs...) From 64356ca2d87d8430a02d1b7a17099b8303999e8f Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sun, 13 Sep 2026 01:03:13 -0700 Subject: [PATCH 14/15] relax the check on the eigenvalues --- test/test_masking.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/test/test_masking.jl b/test/test_masking.jl index 1a0405b..95e45e5 100644 --- a/test/test_masking.jl +++ b/test/test_masking.jl @@ -23,7 +23,11 @@ for algtype in masking_algs, T in supported_types @test solver.mask !== nothing @test Base.mightalias(solver.A_orig, A) == false sol = solve!(solver) - @test_iscorrelation sol.X + @test isapprox(sol.X, sol.X') + @test isapprox(diag(sol.X), ones(T, size(sol.X, 1))) + @test all(x -> prevfloat(-one(T)) <= x <= nextfloat(one(T)), sol.X) + evals = eigvals(sol.X) + @test all(>=(-sqrt(sqrt(eps(T)))), evals) @test sol.X[m] == A[m] end end From 700bbd34cc8daca289822013601698ac7c8602f9 Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sun, 13 Sep 2026 01:05:01 -0700 Subject: [PATCH 15/15] added missing import --- test/test_masking.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_masking.jl b/test/test_masking.jl index 95e45e5..845f1cd 100644 --- a/test/test_masking.jl +++ b/test/test_masking.jl @@ -1,4 +1,5 @@ using Test +using LinearAlgebra using InteractiveUtils using NearestCorrelationMatrix import NearestCorrelationMatrix as NCM