From a049ab74a1b85f9b7de7a0e9e1b91d9caae1f412 Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sun, 13 Sep 2026 12:40:53 -0700 Subject: [PATCH 1/2] Updated default_algtype to properly select the most appropriate algorithm. Updated tests --- src/NCMProblem.jl | 24 --------------- src/NCMSolution.jl | 2 +- src/NCMSolver.jl | 68 ++++++++++++++++++++++++++++++++--------- src/simple_interface.jl | 29 ++++++++++++------ test/test_api.jl | 16 ++++++++++ test/test_simple_api.jl | 14 +++++++-- 6 files changed, 103 insertions(+), 50 deletions(-) diff --git a/src/NCMProblem.jl b/src/NCMProblem.jl index 2780af2..1f16270 100644 --- a/src/NCMProblem.jl +++ b/src/NCMProblem.jl @@ -36,30 +36,6 @@ struct NCMProblem{T, M, P, K} require_matrix(A) require_square(A) require_real(A) - - mask = normalize_mask(mask) - return new{typeof(A), typeof(mask), typeof(p), typeof(kwargs)}(A, mask, p, kwargs) end end - - -normalize_mask(::Nothing) = nothing - -function normalize_mask(B::BitMatrix) - n = require_square(B) - - for i in 1:(n - 1), j in (i + 1):n - b = B[i, j] || B[j, i] - B[i, j] = b - B[j, i] = b - end - - for ii in diagind(B) - B[ii] = false - end - - return B -end - -normalize_mask(B::AbstractMatrix{T}) where {T} = normalize_mask(BitMatrix(B)) diff --git a/src/NCMSolution.jl b/src/NCMSolution.jl index 0077752..9f670d6 100644 --- a/src/NCMSolution.jl +++ b/src/NCMSolution.jl @@ -89,7 +89,7 @@ function CommonSolve.solve!(solver::NCMSolver, args...; kwargs...) # 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). if sol.solver.mask !== nothing - project_fixed!(sol.X, sol.solver.A_orig, sol.solver.mask) + project_fixed!(sol.X, solver.A_orig, solver.mask) project_unit!(sol.X) else cov2cor!(sol.X) diff --git a/src/NCMSolver.jl b/src/NCMSolver.jl index 575d455..5f639d8 100644 --- a/src/NCMSolver.jl +++ b/src/NCMSolver.jl @@ -35,13 +35,6 @@ mutable struct NCMSolver{TA, P, Talg, Tc, Ttol, Tm} A_orig::TA # a copy of A, or an alias of A if no mask is given end -""" - default_algtype(prob) - -Get the default algorithm type for a given input matrix. -""" -default_algtype(prob::NCMProblem) = prob.mask === nothing ? Newton : AcceleratedAP - """ init(prob, alg, args...; kwargs...) @@ -94,10 +87,10 @@ function CommonSolve.init( kwargs... ) # Resolve the effective mask first: an explicit `mask=` kwarg overrides any mask set on the - # problem; otherwise fall back to the problem's (already-normalized) mask, then to nothing. + # problem; otherwise fall back to the problem's mask, then finally default to nothing. mask = if mask !== nothing verbose && println("Using fixed-element mask") - normalize_mask(mask) + mask elseif prob.mask !== nothing verbose && println("Using fixed-element mask from the problem") prob.mask @@ -105,6 +98,9 @@ function CommonSolve.init( nothing end + # Ensure that the mask is normalized to a BitMatrix or Nothing + mask = normalize_mask(mask) + # A non-empty effective mask must be enforced by an algorithm that supports it. Check the # *effective* mask (not just the kwarg) so a mask set on the problem is not silently ignored # by an unsupported algorithm. @@ -235,10 +231,9 @@ end Initialize the solver, and autotune the algorithm to the problem. """ -function CommonSolve.init( - prob::NCMProblem, algtype::Type{<:NCMAlgorithm}, args...; kwargs... - ) - return init(prob, autotune(algtype, prob), args...; kwargs...) +function CommonSolve.init(prob::NCMProblem, algtype::Type{<:NCMAlgorithm}, args...; kwargs...) + alg = autotune(algtype, prob) + return init(prob, alg, args...; kwargs...) end """ @@ -256,5 +251,50 @@ end Initialize the solver with the default algorithm autotuned to the problem. """ function CommonSolve.init(prob::NCMProblem, ::Nothing, args...; kwargs...) - return init(prob, default_algtype(prob), args...; kwargs...) + algtype = default_algtype(prob; kwargs...) + return init(prob, algtype, args...; kwargs...) +end + +""" + default_algtype(prob; kwargs...) + +Get the default algorithm type for a given input matrix. +""" +function default_algtype(prob::NCMProblem; mask = nothing, kwargs...) + mask = if mask !== nothing + mask + elseif prob.mask !== nothing + prob.mask + else + nothing + end + + return mask === nothing ? Newton : AcceleratedAP +end + +""" + normalize_mask(mask) + +Normalizes a mask by converting it to a BitMatrix, or leaves as `nothing`. If a matrix is +given, the resulting BitMatrix is forced to be symmetric. If A[i,j] or A[j,i] is true, then +the resulting BitMatrix will have a true value in both positions. The diagonal elements are +always set to false, since the algorithm should handle setting the diagonal elements to 1. +""" +function normalize_mask(B::BitMatrix) + n = require_square(B) + + for i in 1:(n - 1), j in (i + 1):n + b = B[i, j] || B[j, i] + B[i, j] = b + B[j, i] = b + end + + for ii in diagind(B) + B[ii] = false + end + + return B end + +normalize_mask(B::AbstractMatrix{T}) where {T} = normalize_mask(BitMatrix(B)) +normalize_mask(::Nothing) = nothing diff --git a/src/simple_interface.jl b/src/simple_interface.jl index d03d9dc..a486a57 100644 --- a/src/simple_interface.jl +++ b/src/simple_interface.jl @@ -12,8 +12,7 @@ positive definite, and corrected if it is not. When a fixed-element mask is passed (via the `mask` keyword or on the problem), the result has an exact unit diagonal and retains the masked elements exactly. Because repairing positive definiteness perturbs every entry, strict PD and exact fixed-element feasibility cannot both be -guaranteed: the fixed elements take precedence, and the result is positive definite up to -``O(\\sqrt{\\mathrm{eps}})``. +guaranteed: the fixed elements take precedence, and the result is positive definite up to `√ϵ`. # Examples @@ -42,8 +41,11 @@ true ``` """ function nearest_cor!(A, alg; kwargs...) - sol = solve( - NCMProblem(A), + # (#41) Pass any potential mask to the problem so that the proper algorithm can be selected + prob = NCMProblem(A; kwargs...) + + solver = init( + prob, alg; alias_A = true, fix_sym = true, @@ -52,6 +54,8 @@ function nearest_cor!(A, alg; kwargs...) kwargs... ) + sol = solve!(solver) + copyto!(A, sol.X) return A end @@ -71,8 +75,7 @@ positive definite, and corrected if it is not. When a fixed-element mask is passed (via the `mask` keyword or on the problem), the result has an exact unit diagonal and retains the masked elements exactly. Because repairing positive definiteness perturbs every entry, strict PD and exact fixed-element feasibility cannot both be -guaranteed: the fixed elements take precedence, and the result is positive definite up to -``O(\\sqrt{\\mathrm{eps}})``. +guaranteed: the fixed elements take precedence, and the result is positive definite up to `√ϵ`. # Examples @@ -101,8 +104,15 @@ true ``` """ function nearest_cor(A, alg; kwargs...) - sol = solve( - NCMProblem(A), + # (#41) Pass any potential mask to the problem so that the proper algorithm can be selected + prob = NCMProblem(A; kwargs...) + + # `nearest_cor` cannot simply call `nearest_cor!` with `alias_A=false`. For some reason + # the `alias_A` keyword is passed properly, but `A` still ends up getting aliased anyway. + # The solution is to call `init` with `alias_A=false` independently. + + solver = init( + prob, alg; alias_A = false, fix_sym = true, @@ -111,7 +121,8 @@ function nearest_cor(A, alg; kwargs...) kwargs... ) + sol = solve!(solver) + return sol.X end - nearest_cor(A; kwargs...) = nearest_cor(A, nothing; kwargs...) diff --git a/test/test_api.jl b/test/test_api.jl index 9682f21..ae451ec 100644 --- a/test/test_api.jl +++ b/test/test_api.jl @@ -46,3 +46,19 @@ cache = init(prob) @test_isdefined solve! @test_isimplemented solve!(cache) @test solve!(cache) isa NCMSolution + +# alias_A must be respected +A = default_negdef() +prob = NCMProblem(A) +@test Base.mightalias(A, prob.A) +solver = init(prob; alias_A = true) +@test Base.mightalias(prob.A, solver.A) +solver = init(prob; alias_A = false) +@test !Base.mightalias(prob.A, solver.A) +S = Symmetric(A) +prob = NCMProblem(S) +@test Base.mightalias(S, prob.A) +solver = init(prob; alias_A = true) +@test Base.mightalias(prob.A, solver.A) +solver = init(prob; alias_A = false) +@test !Base.mightalias(prob.A, solver.A) diff --git a/test/test_simple_api.jl b/test/test_simple_api.jl index af517bb..35ace54 100644 --- a/test/test_simple_api.jl +++ b/test/test_simple_api.jl @@ -40,5 +40,15 @@ A = rand(Float16, 4, 4) # (#41) uses an algorithm that supports masking when a mask is given A, m = default_negdef(; include_mask = true) -@test_broken nearest_cor(A; mask = m) -@test_broken nearest_cor!(A; mask = m) +@test_nothrow nearest_cor(A; mask = m) +@test_nothrow nearest_cor!(A; mask = m) + +# nearest_cor must not modify original matrix UNLESS user passes `alias_A=true` +A = default_negdef() +Y = nearest_cor(A, AlternatingProjections) +@test Base.mightalias(Y, A) == false +@test !isapprox(Y, A) +A = default_negdef() +Y = nearest_cor(A, AlternatingProjections; alias_A = true) +@test Base.mightalias(Y, A) == true +@test isapprox(Y, A) From b188a0d4dd19da0be287c93c4e7b1f3e24015597 Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sun, 13 Sep 2026 13:12:04 -0700 Subject: [PATCH 2/2] added missing import --- test/test_api.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_api.jl b/test/test_api.jl index ae451ec..970e4a7 100644 --- a/test/test_api.jl +++ b/test/test_api.jl @@ -1,4 +1,5 @@ using Test +using LinearAlgebra using InteractiveUtils using NearestCorrelationMatrix using NearestCorrelationMatrix.Internals: default_negdef