Progressively fix positive definiteness

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 <alexk.706@gmail.com>"]
 
 [deps]

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

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
 

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

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!,

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