In Sec. 2.1, we outlined the general strategy to do this: the optimal values for the Lagrange multipliers are found by minimizing the Lagrange dual as **given**
by Eq. (8). For concreteness, let us go through the mathematical details for the case of a rectangular binary matrix: D ∈ {0, 1}
m×n**. There should be no**
conceptual difficulties in adapting the derivations below for other choices of D,
and for conciseness these adaptations are omitted from this paper.

For D ∈ {0, 1}
m×n**, the Lagrange dual from Eq. (8) is equal to:**

$$\begin{array}{r c l}{{L(\lambda^{r},\lambda^{c})}}&{{=}}&{{\log(Z(\lambda^{r},\lambda^{c}))-\sum_{i}\lambda_{i}^{r}d_{i}^{r}-\sum_{j}\lambda_{j}^{c}d_{j}^{c}.}}\end{array}$$

Using Z(λ r, λ c) = Qi,j Z(λ r i
, λc j
) and Z(λ r i
, λc j
) = 1 + exp(λ r i + λ c j
), this gives:

$$\begin{array}{r c l}{{}}&{{}}&{{}}\\ {{L(\lambda^{r},\lambda^{c})}}&{{=}}&{{\sum_{i,j}\log(Z(\lambda_{i}^{r},\lambda_{j}^{c}))-\sum_{i}\lambda_{i}^{r}d_{i}^{r}-\sum_{j}\lambda_{j}^{c}d_{j}^{c};}}\\ {{}}&{{}}&{{}}\\ {{}}&{{}}&{{=}}&{{\sum_{i,j}\log\left(1+\exp(\lambda_{i}^{r}+\lambda_{j}^{c})\right)-\sum_{i}\lambda_{i}^{r}d_{i}^{r}-\sum_{j}\lambda_{j}^{c}d_{j}^{c}.}}\end{array}$$

The optimal values of the parameters are easily found using standard methods for unconstrained convex optimization such as Newton's method or (conjugate) gradient descent, possibly with a preconditioner (Shewchuk, 1994; Boyd and Vandeberghe, 2004). We will report computational results for two possible choices in Sec. 7.

Gradient descent type methods rely on the gradient of L**, while Newton's method**
relies on the gradient as well as the Hessian. Both can easily be computed analytically. The gradient is determined by the first order partial derivatives:

$$\begin{array}{r c l}{{\frac{\partial L}{\partial\lambda_{i}^{r}}}}&{{=}}&{{\sum_{j}\frac{\exp(\lambda_{i}^{r}+\lambda_{j}^{c})}{1+\exp(\lambda_{i}^{r}+\lambda_{j}^{c})}-d_{i}^{r},}}\\ {{\frac{\partial L}{\partial\lambda_{j}^{c}}}}&{{=}}&{{\sum_{i}\frac{\exp(\lambda_{i}^{r}+\lambda_{j}^{c})}{1+\exp(\lambda_{i}^{r}+\lambda_{j}^{c})}-d_{j}^{c}.}}\end{array}$$

```
Note that these derivatives have a natural interpretation. Indeed, the sum
Pjexp(λ
         r
         i +λ
            c
            j
             )
   1+exp(λ
          r
          i +λ
             c
             j
              )
               is equal to the expected number of ones in the ith column for
the distribution with the current parameter values, and the partial derivative
∂L
∂λr
  i
    is equal to the difference between that expected number and the value d
                                                                         r
                                                                         i
it needs to be as required by the constraints. The Hessian is determined by the

```
