/SketchyCGAL.jl

Primary LanguageJuliaMIT LicenseMIT

SketchyCGAL

Build Status Coverage

A high-performance implementation of SketchyCGAL with weighted average gradients.

Currently a work-in-progress

Use

This package solves trace-constrained SDPs of the form

$$min tr(CX) s.t. A(X) = b tr(X) ≦ 1 X ⪰ 0,$$

where A(X) is a linear map. To define a problem, we must specify

  • C, the cost matrix in R^{n x n}
  • b, the vector in R^d on the RHS of the equality constraints
  • A(X), the linear map from S^{n x n} to R^d.
  • A_adj(y), the adjoint of the linear map (R^d to S^{n x n}).

We specify the maps as the non-allocating functions

function A!(y, X)
        ...
end

function A_adj!(S::SparseMatrixCSC, z)
        ...
end

where the first argument is modified. Note the adjoint should add to instead of replace S (see example).

Function Calls

Scaling

  • scale_X = 1/tr(X), so that scale_X * tr(X) = 1
  • scale_C = 1 / norm(C) for numerical stability

Optional Parameters

  • R is the sketch sized used in SketchyCGAL (10 or 20 works reasonably well)
  • logging = true causes the solver to return a log of the duality gap (estimated), objective value, primal infeasibility, and time
  • logging_primal = true causes the solver to return a log with the SDP objective after reconstruction of X (instead of just the implicit objective) and the reconstruction's primal infeasibility
  • compute_cut = true causes the solver to add the MAXCUT objective (after rounding) to the primal log
  • ηt is the step size. Defaults to ηt=t->2.0/(t + 1.0)
  • print_iter controls how often the solver prints an update

We can call the CGAL method as follows

XT, yT = SketchyCGAL.cgal_full(
    C, b, A!, A_adj!; n=n, d=d, scale_X=scale_X, scale_C=scale_C,
    max_iters=250,
    print_iter=25
)

Similarly, we can call SketchyCGAL as follows

scgal_full(
    C, b, A!, A_adj!, A_uu!; n=n, d=d, scale_X=scale_X, scale_C=scale_C,
    max_iters=1_000,
    print_iter=100,
    R=R,
    logging=true,
    logging_primal=false,
    ηt=ηt,
)

Example: MAXCUT

See the MAXCUT example in the \examples folder for full code. Assume we have an adjacency matrix G. The maxcut objective is

$$min -1/4( ∑_ij G_ij*(1-xi*xj) ) = -1/4⟨diag(∑ᵢ G_ij) - G), X⟩$$

with constraints

$$diag(X) = 1 X ⪰ 0.$$

We can construct the problem as follows

C = -0.25*(Diagonal(G*ones(n)) - G)
b = ones(n)

# Scaling variables -- so trace is bounded by 1
scale_C = 1 / norm(C)
scale_X = 1 / n

# Linear map
# AX = diag(X)
function A!(y, X)
    @views y .= X[diagind(X)]
    return nothing
end

# Adjoint
# A*z = Diagonal(z)
function A_adj!(S::SparseMatrixCSC, z)
    @views S[diagind(S)] .+= z
    return nothing
 end

Additionally, we must define a function A_uu! that is the linear map applied to the rank one matrix u*u^T. (Note: This inferface may change in the future)

#u -> A(uu^T)
function A_uu!(z, u)
    z .= u.*u
end

Finally, we can call SketchyCGAL (in this case, with no logging)

soln = scgal_full(
    C, b, A!, A_adj!, A_uu!; n=n, d=d, scale_X=scale_X, scale_C=scale_C,
    max_iters=1_000,
    print_iter=100,
    R=20,
)

References

  • Yurtsever, A., Tropp, J. A., Fercoq, O., Udell, M., Cevher, V., Scalable Semidefinite Programming, SIAM Journal on Mathematics of Data Science, 3(1):171-200, 2021 (link)
  • Zhang, Y., Li, B., Giannakis, G., Accelerating Frank-Wolfe with Weighted Average Gradients, ICASSP 2021 (link)

TODOs

  • Finish implementation that only uses primitives defined in the paper
  • Add test cases
  • Add examples
  • Add example documentation