Linear programming using diagonal linear networks

$**Response on the comments on stepsize and linear rate:**$ The upper bound on $\eta$ was explicitly written down in equation (6.1) in the appendix $\eta := \min_{ u\in B_R(0) } \Big( \min \big( \frac{1}{4\text{Res}}, \frac{1}{5L|u|_{\infty}} \big) \Big)$ where $\text{Res}:= |A^\top( A(u\circ u)-b )|_2$ and the value of $R$ can be large in the worst case (leading to a conservative constant for the convergence rate), it can be much better in cases where the matrix $A$ satisfies some regularity properties. For example, if all the entries of $A$ are i.i.d. Gaussian, then one can show that with high probability, $R$ has a polynomial dependence on the problem sizes.

The dependence of $\rho$ on $A$, $b$ and the initial point can be derived from the proof of Theorem 3.4. This dependence is hard to be written in a clean way, as we have used few constants that can be shown to exist but cannot be written in a closed form (see Appendix H.1). Furthermore, since our problem is non-convex in nature, it makes it lot more harder to delineate exact expression of the linear rate. However, we agree that if $A$ is ill-conditioned or the initial point is close to $0$, the value of $\rho$ is close to $1$.

Finally, we note that conservative constants (e.g. Hoffman constant) have been commonly used in deriving the linear convergence of first order methods for solving LP problems. Although such constants may not reflect the actual performance of the algorithm, it provides some hint on the convergence properties of the algorithm. Some of the impactful papers from the past with similar problem includes [1] and [2].

$**Response on the dependence of stepsize and running time on optimality gap:**$ This is an excellent question. We definitely like to include more details in the Appendix. To emphasize how the stepsize depends on the gap, let us define

where $x^*$ is the solution of the LP problem and $R$ is the same constant as in the response of the previous comment. Now a simple but tedious computation would show that if the regularization $\lambda(\epsilon)$ and stepsize $\eta(\epsilon)$ is chosen in way so that $2\lambda(\epsilon) C_1 + 2\eta(\epsilon) C_2 <\epsilon$, then the optimality gap after large enough time $T$ will be less than $\epsilon$.

The a number of iterations $T(\epsilon)$ to reach an $\epsilon$-optimality gap is more subtle. Since our analysis here is significantly different than usual first order methods, it is hard to write a closed-form expression of $T(\epsilon)$. This is partly because our results shows global convergence rate of GD for a non-convex problem. Unlike the analysis of a convex problem, the typical per-iteration analysis is not working here. We have to define some global constants that can be shown to exist but cannot be written in a closed form (see Appendix H.1). The dependence on $\epsilon$ is therefore encoded in a sophisticated way.

$**Response on comparison with other commercial LP solver:**$ We have run a new experiment comparing our method against the commercial LP solver Gurobi and another first-order method: an ADMM-based solver SCS. In particular, we simulate the data with the same data-generating process as in Section 4 of the paper, but with $m = 500$ and $n = 50000$. Here is a brief comparison of the methods in this example: Gurobi, SCS and our algorithm achieve respectively 0, -0.971,0.005 primal gaps and 0, 2.21e-10, 4.17e-11 feasibility gaps in 16.77, 40.49 and 9.26 seconds.

Here the primal gap is defined as $(c^\top \hat x - c^\top x^*) / c^\top x^*$, with $\hat x$ be the solution by an algorithm and $x^*$ is the optimal solution (by Gurobi). The feasibility gap is defined as $\| A \hat x - b \|_2^2 / \max\{1, \| b\|_2^2 \}$. First, it can be seen that our method is much faster than SCS, where our method computes a solution with a smaller primal gap in a much shorter time. Our method is also faster than Gurobi if a low-accuracy solution with a primal gap $0.005$ is satisfactory for the underlying application. Moreover, we note that Gurobi takes $10.11$ seconds for preprocessing, which is already similar to the runtime of our method.

Admittedly, our method is not able to obtain high-accuracy solutions as fast as Gurobi, but it is useful for quickly obtaining an approximate solution. Moreover, since Algorithm 1 only involves matrix-vector multiplications, it is easily implemented in a GPU-acceleration setting, which can potentially outperform Gurobi for larger problems.

[1] Hong, M., Luo, Z.-Q. On the, linear convergence of the alternating direction method of multipliers. (2019) Mathematical Programming.

[2] Luo, Z.-Q., Tseng, P. On the linear convergence of descent methods for convex essentially smooth minimization (1992) SIAM Journal on Control and Optimization

$**Response on the comparison with the Sinkhorn algorithm:**$ Certainly, this is an intriguing question. We acknowledge that the elucidation of the comparison between our algorithm and the Sinkhorn algorithm might not be entirely clear. Our intention was to delineate the similarities and distinctions between these two algorithms by scrutinizing the update rules in each iteration. For instance, we demonstrated in Section 2.3 that our algorithm behaves as row and column rescaling algorithms just like as the Sinkhorn algorithm.

Similar to the Sinkhorn, our neural reparametrized gradient descent initialize the coupling matrix of OT at $X^{(0)} := K$ where $K_{ij} := e^{- C_{ij}/\lambda }$. When the stepsize is small, like as in Sinkhorn, the gradient descent rescales the rows and columns of the coupling matrix. In other words, as the step size gets smaller, gradient descent iterates look like $X^{(k+1)} \approx D(1_n - g^{(k)}) X^{(k)} D(1_m - h^{(k)})$ where $g^{(k)}$ and $h^{(k)}$ are respectively the vectors of residuals of the row sums and column sums of the coupling matrix $X^{(k)}$. Hence the updates above can also be viewed as row and column rescalings, although with different rescaling rules than the Sinkhorn. Given the parallels between the Sinkhorn algorithm and the updates in our gradient descent, it's reasonable to assert that our approach introduces a Sinkhorn-like algorithm for solving linear programming problems within the context of diagonal linear networks. As far as our knowledge extends, this novel application has not been documented in prior works.

$**Resnpose on the comparison with other works on DLN:**$ Although our work bears certain similarities to the references [1] and [2] cited by the reviewer, the primary distinction lies in the diverse approaches to reparametrization. In [1] and [2], the authors have used the reparametrization $\beta= u\odot u - v\odot v$ in the basis pursuit problem while we use a slightly different reparametrization which leads to different final representations of the limit of the gradient descent in our case. However, in both cases, the limit of the gradient descent solve regularized $L_1$-norm minimization problem of $\beta$. But the nature of regularization is different in our work compared to [1] and [2]. Moreover, at the time of submitting our work, the rate of convergence result for gradient descent, as presented in [1] or [2], was not available. To the best of our knowledge, this result was initially introduced in our work. In the following, we discuss the effect of stepsize in our gradient descent in the light of [2] and [3].

$**Response on the effect of stepsize:**$ As the reviewer rightly mentioned, the stepsize play a crucial role in shaping the quality of the algorithm. Theoretical thresholds for the step size necessary for global linear convergence are detailed in Theorem 3.5. However, recent investigations, exemplified in the reference [2] pointed out by the reviewer and the recently published paper in arxiv [3], scrutinize the impact of an increasing step size. These studies, especially the latter one, reveals that, as the step size amplifies, gradient descent in quadratic regression traverses distinct phases. Beyond the theoretical bounds, the loss function oscillates around a diminishing trend, signifying the emergence of the 'catapult' phase. Notably, during the catapult phase, the application of specific ergodic averaging methods has been demonstrated to enhance generalization error. Further escalation of the step size leads to the onset of the periodic phase, marked by periodic oscillations in the loss function without a discernible downward trend. Subsequent to the periodic phase, the system transitions into the chaotic phase, characterized by a lack of convergence signals, ultimately culminating in divergence. Given that our gradient descent iterates mimic those of the least square quadratic regression, we anticipate that the insights from [3] will be applicable to our framework. We would be pleased to incorporate a comprehensive discussion of this in the revised version.

[3] Chen, X., Balasubramanian, K., Ghosal, P., Agrawalla, B. (2023). From Stability to Chaos: Analyzing Gradient Descent Dynamics in Quadratic Regression. ArXiV: 2310.01687

$**Response on the faster rate compared to mirror gradient descent:**$ That is an excellent question. We posit that our algorithm outpaces mirror descent due to the tendency of mirror descent updates to overshoot frequently, resulting in unnecessary excursions near the limit for extended periods. Nevertheless, it remains plausible that a judicious selection of the step size in mirror descent, akin to our algorithm, could potentially enhance its speed.

$**Response on the use of co-ordinate descent:**$ We are very grateful to you evaluation. This is indeed a great question. Certainly, the selection of the step-size is a subject that merits deeper exploration, and we intend to delve into this aspect, drawing insights from recent studies elucidating the impact of step size on gradient descent. Exploring unique step-sizes for coordinate-wise descent presents a fascinating avenue for investigation. We are optimistic about the potential acceleration of our algorithm through coordinate descent; nevertheless, the analysis of its convergence and rate appears to pose distinctive challenges compared to our current approach.

$**Response on the experimental comparison with the Sinkhorn algorithm and other LP solver:**$ Thank you for the question! Actually we compared our method with standard Sinkhorn, but unfortunately, we found that Sinkhorn is better for OT -- it more directly uses the structure of the problem. Nonetheless, we found that our algorithm can be useful for general LPs where Sinkhorn method is not applicable. For instance, we have performed a new experiment, comparing our method with Gurobi and another ADMM-based method. In particular, we simulate the data with the same data-generating process as in Section 4 of the paper, but with with $m = 500$ and $n = 50000$. Here is a brief comparison of the methods in this example: Gurobi, SCS and our algorithm achieve respectively 0, -0.971,0.005 primal gaps and 0, 2.21e-10, 4.17e-11 feasibility gaps in 16.77, 40.49 and 9.26 seconds.

Here the primal gap is defined as $(c^\top \hat x - c^\top x^*) / c^\top x^*$, with $\hat x$ be the solution by an algorithm and $x^*$ is the optimal solution (by Gurobi). The feasibility gap is defined as $\| A \hat x - b \|_2^2 / \max\{1, \| b\|_2^2 \}$. First, it can be seen that our method is much faster than SCS, where our method computes a solution with a smaller primal gap in a much shorter time. Our method is also faster than Gurobi if a low-accuracy solution with a primal gap $0.005$ is satisfactory for the underlying application. Moreover, we note that Gurobi takes $10.11$ seconds for preprocessing, which is already similar to the runtime of our method.

Admittedly, our method is not able to obtain high-accuracy solutions as fast as Gurobi, but it is useful for quickly obtaining an approximate solution. Moreover, since Algorithm 1 only involves matrix-vector multiplications, it is easily implemented in a GPU-acceleration setting, which can potentially outperform Gurobi for larger problems.

$**Response on the use of minibatch SGD:**$ We appreciate your inquiries. The substitution of minibatch gradient descent for gradient descent is not only feasible but can potentially offer advantages. Upon closer examination of our proofs and drawing on similarities between GD and one pass SGD on DLN, it becomes evident to us that the theoretical findings presented in our paper remain robust for minibatch GD with an appropriate choice of the batch size.

$**Response on the effect of initialization:**$ Exploring the impact of step size and initialization is indeed a crucial avenue of investigation. The initialization of our gradient descent proves to be a pivotal factor in guiding the algorithm toward the desired solution. As demonstrated in Theorem ~2.2 and 3.5, we established that initializing the gradient descent with $u^0= \exp(-c/2\lambda)$ leads to convergence toward the solution of $\min_u\{\langle c, u\odot u\rangle + \frac{\lambda}{1-\lambda}\sum_{i} u_i\odot u_i \log (u_i\odot u_i)\}.$ Indeed, selecting a small positive value for $\lambda$ results in a very small initialization. This, in turn, enhances the quality of the solution, given that the regularization term diminishes as $\lambda$ approaches $0$.

$**Response on the effects of stepsize:**$ The impact of the step size is pivotal in achieving a superior generalization error. The theoretical bounds on the step size required for global linear convergence have been established in Theorem. However, recent works, such as in [1], delve into the effect of an increasing step size. In [1], it is demonstrated that, as the step size increases, gradient descent in quadratic regression undergoes various phases. Beyond theoretical bounds, the loss function oscillates around a diminishing trend, indicating the onset of the 'catapult' phase. Notably, in the catapult phase, employing certain ergodic averaging methods has been shown to lead to improved generalization error. Further increasing the step size leads to entry into the periodic phase, characterized by periodic oscillations in the loss function with no discernible downward trend. Subsequently, the chaotic phase ensues, devoid of any signs of convergence, ultimately culminating in divergence. Given that our gradient descent iterates mirror those of least square quadratic regression, we anticipate that the observations made in [1] will hold for our framework. We would be delighted to incorporate a detailed discussion of this in the revised version.

[1] Chen, X., Balasubramanian, K., Ghosal, P., Agrawalla, B. (2023). From Stability to Chaos: Analyzing Gradient Descent Dynamics in Quadratic Regression. ArXiV: 2310.01687

$**Response on comparison with other commercial LP solver:**$ We have run a new experiment comparing our method against the commercial LP solver Gurobi and another first-order method: an ADMM-based solver SCS. In particular, we simulate the data with the same data-generating process as in Section 4 of the paper, but with $m = 500$ and $n = 50000$. Here is a brief comparison of the methods in this example: Gurobi, SCS and our algorithm achieve respectively 0, -0.971,0.005 primal gaps and 0, 2.21e-10, 4.17e-11 feasibility gaps in 16.77, 40.49 and 9.26 seconds.

Here the primal gap is defined as $(c^\top \hat x - c^\top x^*) / c^\top x^*$, with $\hat x$ be the solution by an algorithm and $x^*$ is the optimal solution (by Gurobi). The feasibility gap is defined as $\| A \hat x - b \|_2^2 / \max\{1, \| b\|_2^2 \}$. First, it can be seen that our method is much faster than SCS, where our method computes a solution with a smaller primal gap in a much shorter time. Our method is also faster than Gurobi if a low-accuracy solution with a primal gap $0.005$ is satisfactory for the underlying application. Moreover, we note that Gurobi takes $10.11$ seconds for preprocessing, which is already similar to the runtime of our method.

Admittedly, our method is not able to obtain high-accuracy solutions as fast as Gurobi, but it is useful for quickly obtaining an approximate solution. Moreover, since Algorithm 1 only involves matrix-vector multiplications, it is easily implemented in a GPU-acceleration setting, which can potentially outperform Gurobi for larger problems.

$**Response on generalization for deep/non-linear network:**$ Instead of a diagonal linear network, one might consider employing an $k$-layer deep diagonal linear network. It amounts to parametrizing the gradient descent as $x= u\odot \cdots \odot u$ over $k$ layers for $k>2$. The gradient descent initialized at $u_0 = \alpha1_n$ then $u_k\odot \cdots \odot u_k$ converges to the solution of $\min_x\{\langle 1_n, x\rangle - \frac{k \alpha^{k-2}}{2} \sum_{i=1}^d x_i^{2/k}\}$. Extending this result to intricate architectures of deep non-linear networks is open. However, the optimal choice of initialization in such scenarios remains an area of active research.