Venue: NeurIPS 2021 Conference

Reviews:

final_score: rating: 7: Good paper, accept
review_content: summary: This paper develops a proximal damped version of the Finito algorithm. The algorithm is proved to achieve the same convergence rate as proximal GD, for cyclic sampling, random reshuffling and shuffling-once versions of the algorithm. Further, the authors claim that this is the first* shuffling based variance reduction algorithm to achieve the convergence rate. The paper also gives a new norm that captures the optimality of sampling orders and provides a heuristic based on that for importance based reshuffling.

Besides the theoretical results, the empirical results seem to suggest that the proposed algorithm is indeed faster than other variance reduction algorithms. 

*: The authors cite a concurrent work (Malinovsky et al.) that also achieves the same convergence rates for general convex functions, but the algorithms in the two papers are different.
main_review: This paper proposes Prox-DFinito, which is a shuffling based variance reduction algorithm. 

The theoretical results show that the cyclic sampling, random reshuffling and shuffling-once versions of the algorithm achieve the same convergence rate as GD (up to logarithmic factors) on general convex and strongly convex functions. The paper also proposes a new heuristic to get good sampling orders based on a new norm.

Overall, the theoretical results look good and the empirical evaluation support the theory. 

A concern that I have is regarding the claim that the optimal cyclic sampling can achieve a sample complexity of $\widetilde{O}(L^2/\epsilon)$, which is independent of $n$. While this might indeed be correct, this might be slightly misleading. The fact is that without looking at all the $n$ functions at least once, good convergence cannot be achieved - consider the case where a fraction of the functions have minima very far away from others. Hence, even to determine the optimal cyclic order, at least $\Omega(n)$ computation must be done. The authors should clarify this.

The empirical evaluation for the optimal cycling sampling order seems to be done on an artificial dataset of quadratics that share the same minima. This does not give sufficient indication to whether the proposed optimal cyclic sampling order would work in practice. Can the authors provide evaluation on real datasets? 
limitations_and_societal_impact: 1. Please move the comparison done in Table 3, appendix K to the main paper.
2. Please also provide plots for the test error of the evaluations done on MNIST, CIFAR-10 etc.
review_title: null

final_score: rating: 6: Marginally above the acceptance threshold
review_content: summary: The paper proposes a new method called Prox-DFinito based on the proximal Finito with without-replacement sampling. The authors derive complexity bounds for the proposed method in convex (for making the squared norm of the gradient small) and strongly convex cases (for making the squared distance to the solution small) that match under some additional assumptions the rate of Gradient Descent. Moreover, under additional assumptions on the objective function, the authors show sample-size independent bound in the convex case. The proofs are non-standard but clean and easy to follow. However, the paper has several strong weaknesses.
main_review: ## Strengths

1. **Clarity.** The paper is clearly written and well-organized.

2. **Interesting proofs.** The proofs of the theoretical results of the paper are non-standard and, therefore, valuable. In addition, the proofs are easy to follow and do not contain inaccuracies.

## Weaknesses

1. **Weak guarantees in the convex case.** In the convex case (i.e., all $f_i$ are convex), the authors establish complexity bounds ensuring that $\text{dist}(\|\nabla F(x^{kn}) - \partial r(x^{kn})\|^2, 0) = \min_{g\in \partial r(x^{kn})}\|\nabla f(x^{kn}) + g\|^2 \leq \varepsilon$. In other papers and, in particular, in [18], bounds are obtained to ensure $F(x^{kn}) - F(x^*) \leq \varepsilon$ (no regularization). When the objective is convex and smooth it is possible that the gradient is small but functional suboptimality is huge. Moreover, it is more important to achieve small functional suboptimality rather than the small norm of the gradient. Although the proofs are interesting, the result in the convex case is not strong enough and cannot be directly compared with known results. The authors should clarify it in Table 1 since in the cited papers the bounds in the generally convex case are established for functional suboptimality, not for the squared gradient norm. Moreover, the authors should clarify what do they mean by such a comparison. In the non-proximal scenario, one can use classical inequality $\|F(x)\|^2 \le 2L(F(x) - F(x^{*}))$, but there are no analogs of such inequality for the composite problems. This place should be carefully clarified.

2. **Complexity bounds depend on $\frac{1}{n}\|\mathcal{z}^0 - \mathcal{z}^{*}\|_{\pi\text{ or }2}^2$.** This norm has implicit dependence on the heterogeneity of local loss functions $f_i$. Indeed, if $\mathcal{z}^0$ is close to $x^*$, then for RR $C \sim \frac{\alpha^2}{n}\sum_{i=1}^n\|\nabla f_i(x^{*})\|$ where $C$ is the factor appearing in all upper bounds.In other words, even if the starting point is extremely close to the solution the constant $C$ can be large. This observation immediately implies that all derived rates in the paper can be arbitrarily worse than known results. In contrast, the concurrent work [18] does not have this issue.

3. **Results are weaker than state-of-the-art ones.** In view of the first weakness, the complexity bounds in the convex case are weaker by default since they do not provide guarantees for the functional suboptimality. Next, in the strongly convex case, much tighter results are given in [35], where the authors use $L$ being the average of individual smoothness constants (that can be almost $n$ times smaller than the worst one) and $\mu$ being the strong convexity constant of the average and do not require individual loss functions to be strongly convex. In contrast, this paper uses the worst $L$ for all summands and relies on the assumption that all $f_i$ are $\mu$-strongly convex. Moreover, in view of weakness 2, it is even impossible to fairly compare the results without additional assumptions.

4. **Incomplete comparison with the related work.** The paper tries to create an impression that the obtained results are the current state-of-the-art. However, it is not true. As it was mentioned earlier, the derived results are weaker than the state-of-the-art ones but the authors do not write about such important details. Moreover, since [18] is the recent concurrent work, the authors should add more details in the main part of the paper. First of all, the rates from [18] should be added in Table 1 (otherwise it creates a biased impression). Next, the authors should explicitly write that the results in the convex case from [18] are stronger (see weakness 1). Moreover, [18] contains some results for the "Big Data regime" establishing $O(n \frac{L}{\mu} \log \frac{1}{\varepsilon})$ complexity in the strongly convex case and it is shown without assuming strong convexity of each summand. Finally, rates from [18] do not have the issue described in weakness 2. Therefore, the authors should write about all of these details in the paper and at least briefly mention them somewhere in the main part.

## Questions and comments

1. **Table 1, results from [33] and [5]:** The bounds are strange since their "physical dimension" is incorrect. That is, the complexity should be dimensionless quantity whereas $\frac{L^2}{\mu}$ and $\frac{L}{\mu^2}$ are not. Please, correct these bounds.

2. **line 42, "accelerated":** This word has a certain meaning (Nesterov's acceleration) in the optimization literature that differs from what the authors want to say (e.g., see d'Aspremont, Alexandre, Damien Scieur, and Adrien Taylor. "Acceleration methods." arXiv preprint arXiv:2101.09545 (2021).). The authors should rewrite the sentence.

3. **lines 46-47, "In the generally convex scenario, existing rates for without-replacement sampling with variance reduction are still far worse than GD":** this is not true because [18] exists.

4. **line 101, "Note that the damping step does not incur additional memory requirements":** In theory, it is still $O(nd)$, but in practice, there is a big difference between $nd$ and $2nd$.

5. **Algorithm 3:** The algorithm is not equivalent to Algorithm 1.

6. **inequality (54), the last line:** the identity is not correct since $\sum_{t=0}^{j-i-1} \neq \left(1 + \frac{1}{n}\right)^{j-i-1}$. However, one can stop one line earlier since it is not important for the next arguments in the proof.

## Comment after the rebuttal

I thank the authors for their detailed response. I have read it carefully.

**Bounds in the convex case.** Both authors and me agree that bounds for the norm of the gradient are weaker than bounds for functional suboptimality, so my concern is still valid. However, I acknowledge that before this work and [18], which is written independently and has different analysis, are the first results for VR methods with without-replacement sampling that match the rate of GD in the convex case under the assumption that all $f_i$ and $F$ have the same smoothness constant. Moreover, in the "highly-heterogeneous scenario," the rate derived in this paper is independent of $n$ that I find quite interesting.

**On the $\frac{1}{n}\|z^0 - z^*\|^2$.** Indeed, if $z_i^0 = x^0 - \alpha \nabla f_i(x^0)$, then there is no issue. I suggest the authors write about this initialization in the main text. I think it is important for the comparison with other results.

**Revised Table 1 and clarifications about the comparison with other works.** Now, these parts are much more transparent for the readers. I believe these corrections are very important.

**On Algorithm 3.** Thank you for the clarifications. Now I see the equivalence. Initially, I was confused because of the difference in the formulas for $z_{i_t}^t$ in Algorithms 1 and 3. However, after your clarifications, I realized that the two methods are indeed equivalent. I suggest the authors add these clarifications to the paper.

Overall, although the paper has several limitations, it does contain some valuable contributions. In particular, the idea of using order-specific norm in the analysis is of separate interest. Taking all these aspects into account, I decided to increase my score from 4 to 6. 
limitations_and_societal_impact: The authors adequately addressed the limitations and potential negative societal impact of their work.
review_title: null

final_score: rating: 6: Marginally above the acceptance threshold
review_content: summary: This paper proposes a new first-order algorithm using without-replacement strategy to solve a class of composite convex minimization problems. The main idea is to modify the well-known scheme, called Finito/MISO method by applying a without-replacement strategy and a damping step. Under the convexity and L-smoothness of f, the proposed algorithm achieves O(1/k) rate in epoch on the optimality residual. When f is additionally strongly convex, the rate is improved to linear as in standard proximal gradient methods. The analysis relies on a weight norm defined through the order of the underlying shuffling strategy. The authors also compare their method with other schemes such as a coordinate descent method with cyclic rule, and standard GD. Next, the authors also investigate the optimal cyclic rule for sampling and propose an adaptive variant. Numerical examples on standard logistic regression problem are presented to illustrate the performance of the proposed methods. 

main_review: Originality: I believe that the result of this paper is new, especially the use of a new weighted norm in this case. In fact, the algorithm can be cast into a variant of the ARock method [Peng et al, 2015], but the use of without-replacement is new. The use of damped step is not new, which has been widely used in fixed-point methods as well as in optimization. The key idea of analysis is to show that the corresponding fixed-point mapping is nonexpansive, which allows to use a damped step and achieve O(1/k) rate in epoch. However, the convergence rate is not really encouraging even it is comparable with GD. The reason is that this method requires to store n auxiliary vectors as in SAGA, making it less practical when n and d are large.

Quality/clarity: The paper is well written and well-motivated in general. It does have both algorithmic and theoretical contributions. The technicality sounds and the analysis seems to be nontrivial.
 
Significance: Overall, the paper has nice contribution in terms of new algorithm, especially when using without-replacement strategy. 

Below are some concrete comments and questions:

-- Is the damped step in Algorithm 1 really the key to achieve the O(1/k) rate? In fact, it is not needed for the strongly convex case. I have a feeling that this is only needed if the fixed-point mapping is only non-expansive. Could the authors clarify and discuss this point, because it is highlighted in the paper?

-- It seems that the nonexpansiveness of T is rather trivial since it is constituted by a proximal gradient operator which is nonexpansive for such a chosen stepsize. Does it need any other properties to guarantee the non-expansiveness of T in (11)?

-- The comparison with SVRG and SAGA using randomized reshuffling rule seems to not make sense. These variants may not work with without-replacement rules unless an appropriate variant (e.g., in [18,37]) is used. The authors may use [18] for SVRG, but for SAGA, which variant is used in the experiments?

-- The weighted norm only depends on z0 and z*, and does not depend on the landscape of the problem. This is a very interesting aspect of the new method. Can you also add more discussion on this aspect to elaborate the new contribution?

limitations_and_societal_impact: The paper does not discuss the limitation as well as social impact.
review_title: null

final_score: rating: 7: Good paper, accept
review_content: summary: In this paper, the authors study a stochastic variance reduced method under sampling rules that are without replacement for smooth and non-smooth, convex and strongly convex objectives. More precisely :
- they develop a proximal method call Prox-DFinito for which they study convergence rates for random reshuffling, cyclic and shuffling once samplings
- in the cyclic sampling, they derive an optimal fixed ordering (and a practical adaptive variant which do not require the knowledge of $z^*$
- Finally, the authors highlight the fact that, when assuming data heterogeneity, their study lead to a convergence rate independent of the number of data samples. This was up until now just known for iid sampling.    

 
main_review: First of all, the paper is very well written (especially the background detailing precious works), very very few typos, main proofs seem correct.

Minor comments which are just my opinion:
- 1) I do not find that "generally convex" is a proper way to describe a function which is just "convex". It's not standard and reading it quickly might make the reader think it's more general than convexity.
- 2) Lines 37 and 39 : i would not start a sentence with a ref, I would write "Authors of [XX] establish ..."
- 3) Prop. 1 & 3: I would not put "(proof in Appendix XX)" inside the proposition but before or after. 


Minor/medium importance comments:
- 4) Could you precise in Table 1 that you assume that the $f_i$'s are strongly convex (not $F$), it makes easier the comparison with [18]
- 5) Eq. (13) : I don't know what is the squared norm of a subdifferential. I think it's not standard. I would rewrite this to make it more mathematically correct if you meant, "the squared norm of all elements in the subdifferential".  
- 6) Remark 3: typo "One can takING" -> "take"
- 7) maybe it is worth citing again [16] in Remark 3.
- 8) line 154: why "Apparently" ? It makes the sentence look potentially false
- 9) line 232, 2 typos: $==$ and a word is missing before the parenthesis
- 10) line 242, typo : written $\nabla f (x^*)$, it should be $\nabla f_i (x^*)$ instead

Major comments:
- 11) Could you clearly explain and highlight the fact that Table 1 is comparing number of gradients evaluations (total complexity) to reach $\epsilon$-solutions, it makes easier the comparison with Table 1 in[18]
- 12) I find that the rate of [18] are missing in table 1. Even if an entire section in the appendix describe the differences, we need a more detailed comparison in the main text instead of the short 3 lines at the end of subsection 1.2. Giving the pros and cons of each method (and clearly stating the difference in memory $\mathcal{O} (nd)$ (yours) vs $\mathcal{O} (d)$ ([18])).
- 13) Lines 98-100: It's not clearly explained how you came up with this damping parameter. Why do you need it ?

Major comments on Appendix K:
- 14) You do not discuss the "Big Data regime" + strong convexity for which the authors of [18] clearly have a gain : in their terms they get an iteration complexity of $\mathcal{O} (\frac{L}{\mu} \log (\frac{1}{\epsilon}))$. This assumption lead to the same total complexity than yours for RR and worst cyclic order : $\mathcal{O} (n \frac{L}{\mu} \log (\frac{1}{\epsilon}))$ right ?
- 15) Could you precise that your memory cost is $n$ times the one for SVRG ?

Questions :
- 16) Did you come with the idea of this "order-specific norm". If not, could you please add a reference ?
- 17) I was not aware of the "data heterogeneity assumption" before reading this paper. Where does it come from ? Could you add ref about this ?  Is it met in practice ? What is the link between the assumption and the heterogeneity of the data ?
- 18) Proof of theorem 1, line 465: Why is $z_{\pi(j)}^{kn+j-1} = z_{\pi(j)}^{kn}, \ \forall j=1, \ldots, n\$ true ? For instance, if $\pi(j)=j$ and $k=0$, then it would mean that
$$z_j^{j-1} = z_j^0$$
which is not necessarily the case if the index sampled at step $j-1$ is j. Maybe I'm wrong on this one, but let me know if there is an easy way to see this.
limitations_and_societal_impact: - Authors do not clearly explain how to set $\theta$ (except in the strongly convex case) and how this affects the performance. As said before, the reader don't get why and if this parameter is important or if it is just introduce for the analysis.
- Same remark for $\gamma$ in (22)
- Numerical experiments seem to correctly recover what theory predicts for data-heterogeneity case and for recovering GD rate


**COMMENT AFTER THE REBUTTAL**

I thank the authors for there clear and exhaustive answers to my questions/concerns and those of the other reviewers. 

I read the authors' rebuttal and it made things clearer in my mind on some points : 
- the extended Table 1 clarifying the comparison with other works, especially [18] 
- I read carefully the discussion about the measure of suboptimality done in this work (pointed out by reviewer y5Mf). I also think it is a limitation of this paper. 
- After reading answers to my review and the one of Reviewer VeN3, I now understand now why the dampening parameter is required.

I this paper will be valuable to the community both for its theoretical results (sample independent rate with highly heterogeneous assumption) and its technical aspects (eg the introduction of order-specific norms). Also, fair comparison with previous work (especially [18] in/out of the Big data regime) will be clearly presented in the final version of this work. **Thus, I increase my score from 6 to 7.**



Side comment : I would have been interested to see what is the magnitude of $\rho$ for different datasets and if it is close to $\mathcal{O} (1/n)$ for  some real data sets.
review_title: null

Meta-review: This paper is favored by two of four expert reviewers and is acceptable to the other two. Reviewers agree that it is clearly written and organized. There is also general agreement that parts of the analysis are new and insightful, and possibly useful outside of this paper, in particular the "order-specific" norm that is introduced.

Some reviewers make valid points about the limitations of the end result. For instance, reviewer VeN3 points out that the convergence rate is not necessarily encouraging even if it matches plain GD since the paper's method has memory requirements that scale linearly with the dataset size. The relative advantage then happens in an optimistic case where an optimal cyclic ordering is known (in practice, this ordering is replaced by a heuristic) and when the data is sufficiently heterogeneous (a technical condition detailed in the paper).

I will add that the reviews on this paper were quite thorough, and that there was detailed discussion among reviewers and authors. Ultimately, these exchanges addressed many reviewer concerns and left the authors with suggestions for improvement that might only strengthen the writing further. They've even led to some additional empirical measurements and several references to related work.
Decision: Accept (Poster)