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]

• [x] ad.constraint submodule:
  • [x] 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.
  • [x] 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.
• [x] 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):

• [x] 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.
• [x] Write a base constraint class.
• [x] Compile the exported bindings.
• [x] 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.
• [x] Write a non-negative constraint class with the API as above.
• [x] 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.

• [x] 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.
• [x] 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]

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

[JamesYang007]

• [x] 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.
• [x] 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.
• [x] Change all mini-optimizers to throw exception when max_iters reached.
• [x] Change all constraint classes to take a buffer. Implement a virtual method that returns the buffer size.
• [x] 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().

[JamesYang007]

Optimizations:

• [ ] Dual is assumed to be dense currently. This is a waste of space since it is pretty sparse most of the time.
• [ ] Solver uses dual only for 2 reasons: 1) keep a contiguous storage for all duals (helps when duals are small?) 2) store a sparse-vector dual solution by converting from the dense version. Otherwise, any other interaction with the dual is through passing it to the constraint objects. This suggests the constraint object should be responsible for storing the dual however they want. The simple constraints can keep the dense storage and ConstraintLinear can use a sparse format. This will save a lot of space!
• [ ] Constraint needs an additional function like clear() to clear any cached information. The solver should do this in the beginning (do it in solve_core).
• [ ] New Constraint API?
  • [ ] solve(...): only return primal. Dual is saved internally.
  • [ ] dual_nnz(): returns number of non-zero entries of the current dual variable.
  • [ ] dual(indices, values): populate the current dual variable sparse format.
  • [ ] gradient(x, dual_indices, dual_values, out): compute gradient where primal and dual are given by user. This is useful when computing full gradient in diagnostic.
  • [ ] gradient(x, out): compute gradient where primal is given by user and dual is assumed to be the currently saved value. This is useful in the solver when the gradient is computed on screen set first then the non-screen set.
  • [ ] clear(): clears internal structure to the same state as initial construction. Solver should call this in the very beginning.