Refine Accelerated AP method

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`).
 """

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, 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).
-        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

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
@@ -38,7 +40,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...)
@@ -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

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

src/algorithms/AlternatingProjectionsAA.jl → src/algorithms/AcceleratedAP.jl

@@ -1,40 +1,48 @@
-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(::Type{<: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 = sqrt(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)
     Y = copy(A)
     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!(X, 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

src/algorithms/AlternatingProjections.jl

@@ -32,11 +32,11 @@ 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)
+    return AlternatingProjections(; tau = sqrt(eps(eltype(prob.A))))
 end
 
 function CommonSolve.solve!(solver::NCMSolver, alg::AlternatingProjections; kwargs...)
@@ -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)

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,6 +11,7 @@ struct DirectProjection{A, K} <: NCMAlgorithm
     tau::Real
     args::A
     kwargs::K
+
 end
 
 function DirectProjection(args...; tau::Real = 0, kwargs...)
@@ -24,7 +25,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 +62,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

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)