From 263e6e0ef86a0928909b99bf00eb21591c5113ee Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sun, 13 Sep 2026 19:43:32 -0700 Subject: [PATCH 1/4] Updated problem doc strings --- src/NCMProblem.jl | 24 ++++++++---------------- 1 file changed, 8 insertions(+), 16 deletions(-) diff --git a/src/NCMProblem.jl b/src/NCMProblem.jl index 1f16270..4443dbd 100644 --- a/src/NCMProblem.jl +++ b/src/NCMProblem.jl @@ -4,28 +4,19 @@ Defines the semi-definite programming problem of finding the nearest correlation matrix to a given input matrix. -To define a `NCMProblem`, you only need to provide a square `AbstractMatrix` ``A``. +To define a `NCMProblem`, you only need to provide a square matrix ``A``. Optionally, a mask of fixed element-pairs can be supplied. Only certain algorithms can make use of the mask. -## Problem Type - -### Constructors - -```julia -NCMProblem(A, mask=nothing, p=NullParameters(); kwargs...) -``` - -Parameters are optional, and if not given, then a `NullParameters()` singleton will be used, -which will throw nice errors if you try to index non-existent parameters. Any extra keyword -arguments are stored in the `kwargs` field and forwarded on to the solvers. - -### Fields +## Arguments - `A`: The input matrix. Must be square. Should be symmetric. - `p`: The parameters for the problem. Defaults to `NullParameters`. Currently unused. -- `mask`: A pattern -- `kwargs`: The keyword arguments passed on to the solvers. + +## Keyword Arguments + +- `mask`: A BitMatrix or a matrix of 1s/0s indicating which elements must remain fixed. +- `kwargs`: Additional keyword arguments passed on to the `init` function. """ struct NCMProblem{T, M, P, K} A::T @@ -36,6 +27,7 @@ struct NCMProblem{T, M, P, K} require_matrix(A) require_square(A) require_real(A) + return new{typeof(A), typeof(mask), typeof(p), typeof(kwargs)}(A, mask, p, kwargs) end end From d7274a385e9cf9fcc66b64a6676c25206a0abe0a Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sun, 13 Sep 2026 19:44:31 -0700 Subject: [PATCH 2/4] Updated solve! to progressively check/fix positive definiteness --- src/NCMSolution.jl | 91 ++++++++++++++++++++-------------------------- 1 file changed, 40 insertions(+), 51 deletions(-) diff --git a/src/NCMSolution.jl b/src/NCMSolution.jl index 9f670d6..e10cbd4 100644 --- a/src/NCMSolution.jl +++ b/src/NCMSolution.jl @@ -1,9 +1,9 @@ """ - NCMSolution(x, resid, alg, iters, solver, stats) + NCMSolution(X, resid, alg, iters, solver, stats) Representation of the solution to an NCM problem defined by a `NCMProblem` -## Fields +# Fields - `X`: The solution to the NCM problem. - `resid`: The residual of the solver. @@ -34,65 +34,54 @@ function build_ncm_solution(alg, X, resid, solver; iters = 0, stats = nothing) end """ - solve(prob, alg, args...; kwargs...) + solve!(solver, args...; kwargs...) -Solve the NCM problem with the given algorithm. +Solve the initialized NCM problem. """ -function CommonSolve.solve(prob::NCMProblem, alg::NCMAlgorithm, args...; kwargs...) - return solve!(init(prob, alg, args...; kwargs...)) -end +function CommonSolve.solve!(solver::NCMSolver, args...; kwargs...) + verbose = solver.verbose -""" - solve(prob, algtype, args...; kwargs...) + verbose && println("Beginning solve...") -Solve the NCM problem with the given algorithm type. -The algorithm will be autotuned to the problem. -""" -function CommonSolve.solve( - prob::NCMProblem, algtype::Type{<:NCMAlgorithm}, args...; kwargs... - ) - return solve!(init(prob, algtype, args...; kwargs...)) -end + # solve! must dispatch on both the solver and the algorithm + sol = solve!(solver, solver.alg; kwargs...) -""" - solve(prob, args...; kwargs...) + verbose && println("Finished solving...") -Solve the NCM problem with the default algorithm. -The algorithm will be autotuned to the problem. -""" -function CommonSolve.solve(prob::NCMProblem, args...; kwargs...) - return solve(prob, nothing, args...; kwargs...) -end + if solver.ensure_pd + verbose && println("Checking that the solution matrix is positive definite") + δ = solver.min_eigenvalue + δ = max(δ, eps(eltype(sol.X))) -""" - solve(prob, nothing, args...; kwargs...) + attempt = 0 -Solve the NCM problem with the default algorithm. -The algorithm will be autotuned to the problem. -""" -function CommonSolve.solve(prob::NCMProblem, ::Nothing, args...; kwargs...) - return solve!(init(prob, nothing, args...; kwargs...)) -end + while attempt < solver.max_pd_attempts + if isposdef(sol.X) + verbose && println("Solution matrix is positive definite") + break + end -""" - solve!(solver, args...; kwargs...) + λpre = eigmin(sol.X) -Solve the initialized NCM problem. -""" -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, 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 - # (and unit diagonal) take precedence - the result is PD up to O(√eps). - if sol.solver.mask !== nothing - project_fixed!(sol.X, solver.A_orig, solver.mask) - project_unit!(sol.X) - else - cov2cor!(sol.X) + if solver.mask === nothing + project_psd!(sol.X, δ) + cov2cor!(sol.X) + else + # 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. + project_psd!(sol.X, δ) + project_fixed!(sol.X, solver.A_orig, solver.mask) + project_unit!(sol.X) + end + + λpost = eigmin(sol.X) + attempt += 1 + + verbose && println("Attempt=$attempt, δ=$δ, λ_min_pre=$λpre, λ_min_post=$λpost") + + # if δ∈(0, 1), then repeatedly applying `√` causes δᵢ to converge to 1. + δ = sqrt(δ) end end From 63b40509b5853ec504fb83fa175639a765b5e498 Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Sun, 13 Sep 2026 19:57:04 -0700 Subject: [PATCH 3/4] updated solver with verbosity --- src/NCMSolver.jl | 99 ++++++++++++++++++++------------- src/NearestCorrelationMatrix.jl | 5 +- src/simple_interface.jl | 1 + 3 files changed, 64 insertions(+), 41 deletions(-) diff --git a/src/NCMSolver.jl b/src/NCMSolver.jl index 5f639d8..3845446 100644 --- a/src/NCMSolver.jl +++ b/src/NCMSolver.jl @@ -4,6 +4,8 @@ Common interface for solving NCM problems. Algorithm-specific cache is stored in the `cacheval` field. +# Fields + - `A`: The input matrix. Must be square. Should be symmetric. - `p`: The parameters for the problem. Defaults to `NullParameters`. Currently unused. - `alg`: The algorithm used by the solver. @@ -13,30 +15,31 @@ Common interface for solving NCM problems. Algorithm-specific cache is stored in - `reltol`: The relative tolerance. Defaults to `√(eps(eltype(A)))`. - `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`. +- `max_pd_attempts`: The maximum number of attempts to force the solution to be positive definite. - `mask`: The fixed-element mask, or `nothing` if unmasked. - `A_orig`: The original values of A to be used when a mask is supplied. +- `verbose`: Whether to print extra information. Defaults to `false`. """ -mutable struct NCMSolver{TA, P, Talg, Tc, Ttol, Tm} - A::TA # the input matrix - p::P # parameters - alg::Talg # ncm algorithm - cacheval::Tc # store algorithm cache here - isfresh::Bool # false => cacheval is set wrt A, true => update cacheval wrt A - abstol::Ttol # absolute tolerance for convergence - 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 +struct NCMSolver{TA, P, Talg, Tc, Ttol} + A::TA + p::P + alg::Talg + cacheval::Tc + isfresh::Bool + abstol::Ttol + reltol::Ttol + maxiters::Int + ensure_pd::Bool + min_eigenvalue::Ttol + max_pd_attempts::Int + mask::Union{Nothing, BitMatrix} + A_orig::TA + verbose::Bool end """ - init(prob, alg, args...; kwargs...) + init(prob, alg, args...; kwargs...)::NCMSolver Initialize the solver with the given algorithm. @@ -66,6 +69,7 @@ Initialize the solver with the given algorithm. - `min_eigenvalue`: The minimum eigenvalue to enforce when `ensure_pd` is `true`. Defaults to `nothing`, in which case it is either unused or set to a reasonable value depending on the problem parameters. +- `max_pd_attempts`: The maximum number of attempts to force the solution to be positive definite. - `verbose`: Whether to print extra information. Defaults to `false`. """ function CommonSolve.init( @@ -74,22 +78,32 @@ function CommonSolve.init( args...; mask = nothing, alias_A = default_alias_A(alg, prob.A), + # generic algorithm controls abstol = default_tol(real(eltype(prob.A))), reltol = default_tol(real(eltype(prob.A))), maxiters::Int = default_iters(alg, prob.A), + # keywords regarding symmetry fix_sym::Bool = false, uplo::Symbol = :U, + # regarding Float16 inputs convert_f16::Bool = false, force_f16::Bool = false, + # regarding positive definiteness ensure_pd::Bool = false, min_eigenvalue = nothing, + max_pd_attempts::Int = 5, + # additional keywords verbose::Bool = false, kwargs... ) + A = prob.A + p = prob.p + T = eltype(A) + # Resolve the effective mask first: an explicit `mask=` kwarg overrides any mask set on the # problem; otherwise fall back to the problem's mask, then finally default to nothing. mask = if mask !== nothing - verbose && println("Using fixed-element mask") + verbose && println("Using fixed-element mask supplied to init") mask elseif prob.mask !== nothing verbose && println("Using fixed-element mask from the problem") @@ -100,6 +114,14 @@ function CommonSolve.init( # Ensure that the mask is normalized to a BitMatrix or Nothing mask = normalize_mask(mask) + if mask !== nothing + size(A) == size(mask) || + throw( + DimensionMismatch( + lazy"The problem matrix and the mask must both be square matrices of the same size. Got $(size(A)) and $(size(mask))" + ) + ) + end # 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 @@ -112,30 +134,25 @@ function CommonSolve.init( ) end - A = prob.A - p = prob.p - - T = eltype(A) - A = if alias_A - verbose && println("Aliasing A") + verbose && println("Aliasing `A` to the matrix in the problem") A elseif A isa Symmetric if supports_symmetric(alg) - verbose && println("Creating a Symmetric copy of A") + verbose && println("Creating a Symmetric copy of `A`") copy(A) else verbose && println( "$(alg_name(alg)) does not support Symmetric types. " * - "Creating a symmetric copy of A.data" + "Creating a symmetric copy of `A.data`" ) Matrix(A) end elseif A isa Matrix - verbose && println("Creating a copy of A") + verbose && println("Creating a copy of `A`") copy(A) else - verbose && println("Creating a deep copy of A") + verbose && println("Creating a deep copy of `A`") deepcopy(A) end @@ -192,35 +209,39 @@ function CommonSolve.init( A_orig = mask === nothing ? A : copy(A) # Guard against type mismatch for user-specified reltol/abstol - reltol = real(T)(reltol) - reltol = max(reltol, sqrt(eps(T))) - abstol = real(T)(abstol) - abstol = max(abstol, eps(T)) + reltol = T(reltol) + abstol = T(abstol) min_eigenvalue = if min_eigenvalue === nothing if ensure_pd if mask === nothing - # no mask, can default to sqrt(eps(T)) - sqrt(eps(T)) + δ = eps(T) + verbose && println("Setting min eigenvalue to $δ") + # no mask, can default to eps(T) as a starting point + δ else # be more conservative about the min eigenvalue when there is a mask - sqrt(sqrt(eps(T))) + δ = sqrt(eps(T)) + verbose && println("Setting min eigenvalue to $δ") + δ end else # no checks for PD -> min_eigenvalue is not used - nothing + zero(T) end else # user explicitly set min_eigenvalue. Just ensure that it is Real - real(T)(min_eigenvalue) + δ = T(min_eigenvalue) + verbose && println("Setting min eigenvalue to $δ") + δ 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, min_eigenvalue, verbose, mask, A_orig + solver = NCMSolver{typeof(A), typeof(p), typeof(alg), Tc, T}( + A, p, alg, cacheval, isfresh, abstol, reltol, maxiters, ensure_pd, min_eigenvalue, max_pd_attempts, mask, A_orig, verbose ) return solver diff --git a/src/NearestCorrelationMatrix.jl b/src/NearestCorrelationMatrix.jl index fa23525..ecdf695 100644 --- a/src/NearestCorrelationMatrix.jl +++ b/src/NearestCorrelationMatrix.jl @@ -1,7 +1,8 @@ module NearestCorrelationMatrix using LinearAlgebra -using CommonSolve: CommonSolve, init, solve, solve! +import CommonSolve +using CommonSolve: init, solve, solve! include("internals/Internals.jl") using .Internals @@ -34,8 +35,8 @@ export autotune, # common solve interface init, - solve, solve!, + solve, # just a re-export of the default implementation # simple interface nearest_cor, nearest_cor!, diff --git a/src/simple_interface.jl b/src/simple_interface.jl index a486a57..0657b65 100644 --- a/src/simple_interface.jl +++ b/src/simple_interface.jl @@ -125,4 +125,5 @@ function nearest_cor(A, alg; kwargs...) return sol.X end + nearest_cor(A; kwargs...) = nearest_cor(A, nothing; kwargs...) From 109d7b0bcf619cf62f4ce0a4e1db21c80404d1c1 Mon Sep 17 00:00:00 2001 From: Alex Knudson Date: Wed, 16 Sep 2026 15:38:15 -0700 Subject: [PATCH 4/4] update version number --- Project.toml | 2 +- src/algorithms/DirectProjection.jl | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 3152743..ad127b9 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "NearestCorrelationMatrix" uuid = "59ddf330-608c-4938-8bc9-a4ee97bbbea6" -version = "1.2.1" +version = "1.3.0" authors = ["Alex Knudson "] [deps] diff --git a/src/algorithms/DirectProjection.jl b/src/algorithms/DirectProjection.jl index 7a2b383..31bc533 100644 --- a/src/algorithms/DirectProjection.jl +++ b/src/algorithms/DirectProjection.jl @@ -11,7 +11,6 @@ struct DirectProjection{A, K} <: NCMAlgorithm tau::Real args::A kwargs::K - end function DirectProjection(args...; tau::Real = 0, kwargs...)