Group elastic net with constraints?

[JamesYang007]

At the very least, it'd be nice to support constraints of the form Ax = b and Ax <= b.

[JamesYang007]

• ad.constraint submodule:
  • ad.constraint.linear: general linear inequality constraint Ax <= b. A should be taken as some abstract matrix type (should be fine to use MatrixNaive... since we should only require cmul and ctmul. b is obviously just a dense vector. Incorrect API.
  • All constraints must implement .update() function. It should take in c, d, alpha, Q, D, v (first 3 scalars, Q square dense matrix (small), D diagonal vector, v vector). It just solves the prox operator using nnqp solver and returns the block update.
• Need another separate solver for group elastic net with constraints. BCD just needs to check if there is constraint or not and call the respective block update (fit step). abs_grad definition needs to change (invariance step) to include the dual component. Then KKT + screening steps should remain the same. Fully integrated in existing solvers to take constraints instead.

[JamesYang007]

TODO (ordered):

• Write skeleton code for state_gaussian_pin_constraint and solve_gaussian_pin_constraint. Should generalize control flow of the unconstrained version and only the coordinate_descent function should be swapped. Code should assume a constraint class exists and the API should fall nicely.
• Write a base constraint class.
• Compile the exported bindings.
• Export the bindings to Python modules. Basically, we just need to include constraints as an argument to grpnet, state functions, and change _solve function that delegates to the C++ solver functions to properly include the case of constraints.
• Write a non-negative constraint class with the API as above.
• Write internal core routines of the constraint class in bcd.constrained for unittesting purposes. This is unnecessary since the constraint classes already form a unit for the proximal operator.

[JamesYang007]

Need better intuition for how the iterates are behaving with proximal Newton method. An issue seems to occur for ad.constraint.lower when b is positive and small (close to 0). The dual seems to move too much and enters the regime where ||v - A^T mu||_2 <= lmda. We want to confirm this behavior. Try to find a minimal problem with dimension 2 for ease of visualization.

• Write a debug method (preprocessor def) that populates primals/duals and returns them.
• Write a visualization tool that shows dual iterates with color gradient scheme.
• Reproducible example:

d = 2
eps = 1e-5
seed = 6

np.random.seed(seed)
quad = np.random.uniform(0, 1, d)
linear = np.sqrt(quad) * np.random.normal(0, 1, d)
l1 = 0.5 * np.linalg.norm(linear)
l2 = 0
Q = np.random.normal(0, 1, (d, d))
Q, _, _ = np.linalg.svd(Q)
Q = np.asfortranarray(Q)
b = np.full(d, eps)
constr = ad.constraint.lower(
    b, 
    dtype=np.float64,
)

[JamesYang007]

• Change ConstraintLowerUpper to ConstraintOneSided and make _sgn a vector of +/- 1s. The rest of the implementation should mostly be the same?
• Fix bug in multi-response GLM where constraints and dual_groups must be adjusted if intercept=True.
• Output duals as sparse matrix exactly like the betas.
• Change the diagnostic.gradient_norms implementation.
• Implement ConstraintOneSided with ADMM. Is it sensitive to rho even with warm-starts? YES very!
• Implement ConstraintBox following the template for ConstraintOneSided.
• Implement ConstraintLinear. Not sure what the right algorithm is here..
• Update write-up of what's working.

[JamesYang007]

• Change solver so that constraint buffer is "union"ed instead of having separate buffers for each constraint object (waste of space). The constraints are accessed always iteratively since CD is sequential by nature. This will save a lot of space.
• Clean-up mini-optimizers to have the same API. Just a simple class with a solve() method. Functional doesn't make much sense in this case IMO.
• Change all mini-optimizers to throw exception when max_iters reached.
• Change all constraint classes to take a buffer. Implement a virtual method that returns the buffer size.
• Change pin methods to get the largest requested buffer size from each constraint class and allocate that buffer. This must be passed into each constraint class when we invoke solve().