Faster null fit

[ailurophilia]

Adds fit_null_repar, a couple helper functions, and a test to confirm this finds the same null fit as the augmented lagrangian approach.

Instead of using a constrained optimization method, fit_null_repar treats B_k_constr_j_constr as a function of the remaining elements of B_k_constr and solves the resulting (likelihood) optimization problem by coordinate descent. This is similar to but not the same as the coordinate descent algorithm used in full model fitting. The main differences are the following

• since B_k_constr_j_constr depends on every other element of B_k_constr, we loop through j such that j is not the convenience constraint (j_ref) or j_constr and then (conditional on z) take a fisher step in both j and j_constr
• we use optim to adjust the overall position of each row of B (relative to B_j_ref, which is forced to equal zero) in between loops through j. That is, we fix B at its current value and solve an optimization in epsilon_1, ..., epsilon_p, where epsilon_k is added to every element of B[k,] except B[k,j_ref] (and B[k_constr,j_constr] if k = k_constr). For B[k_constr,], epsilon_k_constr is added to all elements except j_constr and j_ref; B_k_constr_j_constr is then updated based on the new values of all other elements of B_[k_constr,]. This optimization is performed on the likelihood with z profiled out (since z is very informative wrt overall location of B). This helps us avoid situations where coordinate descent is slow because we have to take a bunch of small steps in B[,c(j,j_constr)] instead of one big.

In my experience so far, this outperforms the augmented lagrangian approach both in stability (I have yet to see it fail) and in scalability (I'm able to fit nulls in J = 1000 categories in a matter of 10 or 20 minutes even when B[k_constr,j_constr] has a large value under the full model fit.

However,

• as currently coded, this function assumes that the constraint function is permutation-invariant. This will be true of any reasonable measure of central tendency (e.g., mean, median, pseudo-Huber center, etc.), but it is not true when we want to constrain one element of B[k_constr,] to equal another under the null. I'm noting this for posterity, but I also don't think it's a big problem, because we can fit this type of null model via the full model coordinate descent algorithm if we just set B[k_constr,j_constr] to zero at the beginning of optimization and then skip updating it.
• also as currently coded, this fit_null_repar requires j_ref = J (because I didn't want to deal with indexing problems while I wrote it for the first time -- this requirement can be removed easily enough) and also that j_ref != j_constr. This second requirement is fundamental to the approach -- constraining our model to the null via reparametrization does not work if j_ref = j_constr. Practically, however, this should not pose a large problem since j_ref is arbitrary, so we can just choose a j_ref not equal to j_constr, but it's worth noting that we need to enforce this.
• I have written fit_null_repar as a replacement for / alternative to fit_null, and I have added a test comparing these two algorithms on a toy problem, but I have not integrated fit_null_repar into score_test, which will be the next step needed to incorporate this algorithm into radEmu.