We thank the reviewer for their assessment and questions. We respond to each point inline below. We encourage the reviewer to ask any follow-up question and to comment on our responses. We hope that they will consider raising their score at the end of our exchange.

The results look promising but are limited to a single metric, accuracy, which might not accurately the performance of each method. Additional common metrics, such as, the negative log predictive density and the calibration error, are missing at the current stage.

We did not provide these metrics because they are not available for our competitors. Here we provide a selection of the results of our method using the negative log predictive density (NLPD) and the calibration error (ECE) metrics.

For the 1W1A setting on the CIFAR-10 dataset:

For the 1W32A setting on the CIFAR-10, CIFAR-100 and TinyImageNet dataset:

For the ImageNet dataset:

Could the authors elaborate on connections of their methods to natural gradient-based optimization in variational inference, as done in classical and recent literature such as [Honkela] or [Shen]?

We thank the reviewer for pointing out these references. Our approach is comparable to that of Shen in the setup of the variational learning problem and the usage of Gaussian variational distributions. However, in contrast with Shen, our method works over low-rank covariances, while Shen adopts diagonal covariance matrices. Finally, to the best of our understanding, the natural-gradient methods of Honkela and Shen are based on considering the Riemannian gradient with respect to the Fisher metric, while our method implicitly uses the Riemannian gradient over the space of low-rank Gaussian distributions defined via the embedding $\Sigma = Z Z^T.$ This has the advantage of enforcing the low-rank constraint and is the basis for the Burer-Monteiro approach to semidefinite programming. We believe that this key ingredient, together with the connection with Grothendieck’s inequality and hyperplane rounding, is completely absent from the previous work cited by the reviewer. This is to be expected as that work does not have to deal with the training of binary neural networks and is not concerned with rounding from real to binary weights. In conclusion, we cannot agree with the reviewer’s assessment regarding the limited novelty of our method and reassert its significant departure from previous work.

The authors model correlations between all units in the network (line 182) through a low-rank approximation and need to resort to a very low rank approximation of rank smaller or equal to ten. Often a block-diagonal matrix assumption is made (assuming layers are independent). Did the authors experiment with a block-diagonal structure in which blocks are low-rank? If so, how does such an approach compare? It seem that a high-rank for each block would be feasible in this setting.

We concur with the reviewer on the promise of using a block low-rank structure, a direction we have already been exploring. In our first presentation of the method we opted to model all correlations as we were interested in studying the maximum possible advantage we can obtain via a Gaussian approximation. Now that we have established the existence of this advantage, our future work will focus on how to reap it without incurring the computational costs, particularly in terms of memory, that VISPA suffers. We will include this discussion in Section 5.

After reading up about hyperplane rounding, I can see the connection to using a sign function. However, the authors only mention this without elaborating on the connection. I would recommend making those connections more explicit and elaborating on them.

While this connection is well-known in the theoretical computer science community, we agree that it deserves more elaboration for a more general audience. Unfortunately, we were strongly constrained by space limitations in the detail we could provide in Section 3.1, as the geometric view of hyperplane rounding requires introducing the vector embedding view of the underlying SDP. If the reviewer deems it satisfactory, we would be happy to include a separate section detailing the geometric derivation of hyperplane rounding in the Supplementary Material.

A critical aspect of the method seem that it relies on a low-rank approximation, however, choosing the rank might be non-trivial. Moreover, as the ablation in Fig.1 (note that this should be a table) indicates it is not clear that increasing the rank necessarily improves the results. How is the rank chosen in the experiments and how do the authors propose to chose the appropriate rank?

See FAQ 2 in the author rebuttal. Most results presented use $K=4$, as this was experimentally showed to be a good compromise between performance and computational cost. More generally, for a robust result, we would recommend averaging solutions from different choices of rank.

We thank the reviewer for their assessment and questions. We respond to each point inline below. We encourage the reviewer to ask any follow-up question and to comment on our responses.

The non-diagonal low-rank representation of the covariance adds a significant overhead. It is unclear whether the algorithm can remain scalable. From Figure 1 (Impact of K on model accuracy) it remains unclear whether setting K > 1 is really worth the additional computational overhead.

The computational overhead is mostly restricted to the usage of memory. If $n$ is the number of weights, our memory for storing relevant variables is at most $(K+1) n$ compared to $n$ for other methods. The total significance of this increase depends on the batch size and the size of the examples. For instance, in the case of ResNet18, there are roughly $9 \cdot 10^6$ weights while, for batch size=100, the total size of a data batch is $15 \cdot 10^6$, for a total memory usage of $24 \cdot 10^6$. Hence, choosing K=1 will lead to a usage of $33 \cdot 10^6$, a $37.5$ increase. This increase will typically get smaller as we consider models trained on larger images. In practice, we often observe even smaller increases, as our estimate does not include memory usage due to backpropagation, which can be very large, e.g., in the case of residual connections.

The effect on the runtime is negligible because most of the runtime is spent in forward and back-propagation anyways. In conclusion, our method is most suitable in situations where memory is not strictly constrained or can be augmented at low cost. We plan to include this discussion in the final version of the paper.

The VISPA algorithm is in some sense "heuristically" derived, and it is unclear whether for example stationary points of the method are critical points of the original training objective.

We agree, but notice that the focus of this paper is to provide a mathematically motivated algorithm that performs competitively in practice. We believe it eminently possible to obtain the suggested result for sufficiently small step size under the assumption that the loss is smooth by adapting existing results in manifold optimization. Indeed, we can think of our algorithm in the Burer-Monteiro formulation as minimizing the loss function over the manifold given by the constraints $\mu_i^2 + (ZZ^T)\_{ii} = \mu_i^2 + \sum_j Z_{ij}^2=1$ for all $i$, which is a product of spheres, a well-studied case [1]. As suggested by the reviewer, we can view step 11 of our algorithm as a projection step onto the manifold, which effectively implements a retraction from the tangent space back to the manifold. To make the convergence argument complete, it only remains to analyze the effect of minibatch sampling of the loss on the gradient estimates used by the algorithm. Again, by taking the step size sufficiently small, it should be possible to make the contribution of the variance negligible and guarantee that the algorithm obeys a descent lemma and finds a near-stationary point.

[1] Boumal, N. (2023). An introduction to optimization on smooth manifolds.

I appreciate the connection to variational inference, but in the objective function there is no entropy or prior present, so it is unclear whether this algorithm actually has any connections to Bayesian inference.

Please see FAQ 3 in the author rebuttal. We will add the entropy term in the final version and justify why we let it vanish for our final algorithm.

Can lines 11-13 in VISPA algorithm be viewed as an orthogonal projection onto the constraint? It would be interesting to view VISPA as a projected stochastic gradient method on the problem P_corr, which would help in proving convergence and making the algorithm more theoretically sound.

Yes. As explained above, Step 11 of our algorithm indeed performs a projection step over the manifold of interest. We would be happy to add this explanation in the final version.

Even for K=1, the method seems to perform better than other variational baselines like BayesBiNN. Is there any intuitive explanation as of why the algorithm works so much better in practice?

The improved behavior for $K=1$ is most likely a benefit of the randomization built into our rounding argument, even for $K=1$. While other methods deterministically round based on the current $\mu$-value, our algorithm for K=1 effectively adds noise to $\mu$ via the $Zr$ term before rounding, allowing us to better explore the space of possible assignments.

We thank the reviewer and regret to find that the connections to variational inference and semidefinite programming were not clear to them. We believe that, in the attempt to streamline the presentation of the mathematical component of our method, we may have omitted some steps that clarify such connections. A summary of our clarifications are given in FAQ 3 and FAQ 4 of the author rebuttal. Below, we aim to present more details of our clarifications, together with their proposed location in the final version of the paper. Lastly, we note that the closely related works by Ajanthan et al [1] and Meng et al [2] (see Section 3.1) also rely on the same connection to variational inference. More recently, the paper by Shen et al (pointed out by Reviewer 6N2m) [3] describes a similar framework to ours for the (non-binarized) training of LLMs under the moniker of variational learning, a variant of variational inference. We welcome further questions and suggestions that may improve the clarity of our work.

CONNECTION TO VARIATIONAL INFERENCE

ENTROPY REGULARIZATION OF PROBLEM 2 (Top of Section 3)

Motivated by stability and generalization concerns, we may consider an entropy regularized version of Problem 2 given by

$OPT(\theta) := \min_{p \in \mathcal{P}} \theta \cdot \mathbb{E}_{w \sim p, x}[L(f(x, w), y_x)] - H(p)$

where $H(p)$ is the entropy of the distribution $p.$ This can be seen to be the log-partition function of a probabilistic graphical model (PGM), albeit one with a complicated $n$-way sufficient statistic given by the loss function $L$, and the corresponding optimizer is the Gibbs distribution $p_\theta.$

Following Wainwright and Jordan [4 -Chapter 3.7 and Chapter 5.1], the task of approximating the log-partition function of a model is exactly where we may deploy the toolkit of variational inference.

MEAN FIELD APPROXIMATION (Line 145 in Section 3)

The mean field technique [4 - Chapter 5.2-5.3] is based on approximating $p_\theta$ by finding its closest approximation, in KL divergence, from a family of tractable models. Meng [2] and Ajanthan [1] consider the naive mean field approximation, which corresponds to selecting the class of graphical models with no edges, i.e., product distributions. The resulting optimization problem is:

$\min\_{p \in \mathcal{P}\_{Bernoulli}} D(p || p\_{\theta}) = \min\_{p \in \mathcal{P}\_{Bernoulli}} D(p || p\_{\theta}) = \min\_{p \in \mathcal{P}\_{Bernoulli}} \theta \cdot \mathbb{E}\_{w \sim p, x}[L(f(x, w), y\_x)] - H(p) - OPT(\theta)$ where $\mathcal{P}_{Bernoulli}$ is defined in the submission.

BNN Training via Gaussian Variational Inference (Section 3.1)

Gaussian variational inference ([5] and references therein) suggests to use as tractable models the family of multivariate Gaussian distributions. We then obtain:

$\min\_{p \in \mathcal{P}\_{\textrm{corr}}} D\_{KL}(p || p\_{\theta}) = \min\_{p \in \mathcal{P}\_{\textrm{corr}}} \theta \cdot \mathbb{E}\_{w \sim p, x}[L(f(x, w), y\_x)] - H(p) - OPT(\theta)$

At first, this formulation may seem problematic as we are fundamentally changing the reference measure from the discrete measure over the binary hypercube $\\{-1,+1\\}^n$ to the Lesbegue measure over $\mathbb{R}^n.$ However, this change is justified, when $\theta$ goes to infinity, by the approximation result of Grothendieck's inequality, which formally shows that the loss minimization problem over $\mathcal{P}\_\textrm{corr}$ approximates that over $\mathcal{P}\_{Bernoulli}$ for certain classes of objective functions. See discussion in 3.1 lines 165-178.

This reasoning leads us to consider the regime in which the entropy term vanishes ($\theta \to \infty$) to obtain the problem:

$$\min\_{p \in \mathcal{P}\_\textrm{corr}} \mathbb{E}\_{w \sim p, x}[L(f(x, w), y\_x)]$$

This is exactly the problem studied in our paper. For sake of completeness, at a preliminary stage of our investigation, we also had a version of our gradient-based SDP algorithm include the entropy term for different large values of $\theta$. In all cases, our final algorithm ($\theta \to \infty$) proved more stable, accurate and faster to run, so that we did not pursue different settings of $\theta$ further.

[1] Thalaiyasingam Ajanthan, Puneet Dokania, Richard Hartley, and Philip Torr. Proximal Mean-Field for Neural Network Quantization. ICCV, 2019.
[2] Xiangming Meng, Roman Bachmann, and Mohammad Emtiyaz Khan. Training Binary Neural Networks using the Bayesian Learning Rule. ICML, 2020.
[3] Shen et al., Variational Learning is Effective for Large Deep Networks. ICML, 2024.
[4] Wainwright, M.J. and Jordan M.I., Graphical Models, Exponential Families, and Variational Inference.
[5] Diao, M.Z., Proximal Gradient Algorithms for Gaussian Variational Inference:Optimization in the Bures–Wasserstein Space, MIT Thesis.

CONNECTION TO SEMIDEFINITE PROGRAMMING

See FAQ 4 in the author's rebuttal.

QUESTION ABOUT MEMORY USAGE

The memory consumption of this method is not reported, and I suspect that this method is way more memory-intensive than baselines.

If $n$ is the number of weights, our memory for storing relevant variables is at most $(K+1) n$ compared to $n$ for other methods. The total significance of this increase depends on the batch size and the size of the examples. For instance, in the case of ResNet18, there are roughly $9 \cdot 10^6$ weights while, for batch size=100, the total size of a data batch is $15 \cdot 10^6$, for a total memory usage of $24 \cdot 10^6$. Hence, choosing K=1 will lead to a usage of $33 \cdot 10^6$, a $37.5$ increase. This increase will typically get smaller as we consider models trained on larger images.

Hi Authors,

Thanks for the response.

Connection to Variational Inference

Could you clarify what is parameter $\theta$? You have used $\theta$ as a scaler to multiply the expectation term $\mathbb{E} [L(f(x, w), y)]$ and subsequently let $\theta$ go to infinity. On the other hand, you also use $\theta$ to index the weight distribution $\mathbf{w} \sim p_{\theta}$, which in general has to be represented by a high dimensional vector in the context of neural networks.

Also, can you clarify how you arrive at the first equation $\theta \cdot \mathbb{E}[L(f(x, w), y)] - H(p)$? I think you want to define a conditional distribution $p((x, y) \mid w) = \exp(- L(f(x, w), y))$ with a prior $p_0(w)$. Doing variational inference in $w$ gives the (negative) ELBO $-\mathbb{E}[\log p((x, y) \mid w) + \log p(w)] - H(p) = \mathbb{E} _ {p} [L(f(x, w), y)] - \mathbb{E}_p [\log p_0(w)] - H(p),$ where $p$ is the variational distribution. It's quite close, but not the same as the equation you wrote.

Connection to SDP

I have read FAQ 4. But I have to say that I still find this connection is very superficial.

Memory Usage

I am partially convinced. However, I also want to note that ResNet18 is not a huge model by today's standard. For instant, a ResNet50 model will have a higher (relative) memory overhead compared to ResNet18.

We thank the reviewer for their assessment and questions. We respond to each point inline below. We encourage the reviewer to ask any follow-up question and to comment on our responses. We hope that they will consider raising their score at the end of our exchange.

The paper can benefit from empirical evaluation over a broader class of architectures since BNNs can have applicability beyond vision as a domain. I do acknowledge that the domain is typically restricted to the class of vision based transforms, but leveraging these techniques in other major model classes of interest can have significant impact which is left unexplored.

See FAQ 1 in the author rebuttal.

The authors claim that the problem involves an SDP but the algorithm presented has very loose connections to most formalizations of semi definite programming (and variational inference). The SDP claim seems to derive only from a semi definite constraint and the problem itself fails to be convex let alone semidefinite. (See Line 163-164).

In the author rebuttal, please see FAQ 3 for the connection to variational inference and FAQ 4 for the connection to semidefinite programming.

In L242-243, the authors mention taking 40 samples and average results for predictions. How much does this impact performance (aka inference time, storage etc), since one of the primary motivations of BNNs is to minimize computational storage and runtime requirements?

This has limited impact on storage, as the algorithm can maintain the average dynamically. The impact on inference time may be as large as 40x, but may be significantly reduced by parallelization of the processing of the 40 samples. It is an active area of focus to further reduce this impact. In preliminary results, we find that the process may be sped by using a smaller number of correlated samples via Gaussian quadrature. In this case, $2K +1$ samples would suffice, which is as small as 3 for the $K=1$ version of our method. In conclusion, we note that the prime runtime concern is to reduce the cost of training, which we achieve by binarizing both weights and activations while outperforming competitors.

It is well established in quantization literature that activations (and weights) both admit low precision representations without significant performance loss. Did you try activation precision below 32 bits since low precision architectures are quite common (and perhaps even the norm)?

We conducted experiments using binary activation precision. The results of these experiments are detailed in Tables 1, 3, and 4, which show the performance under the settings of binary activation and binary weights across different datasets and network architectures. For future work, we plan to explore other low precision configurations, including 2-bit, 3-bit, and 4-bit precision.

Do you have a hypothesis for why lower K perform better for 1W1A while higher K is better for 1W32A architectures?

See FAQ 2 in the author rebuttal.

We thank the reviewer for their assessment and questions. We respond to each point inline below. We encourage the reviewer to ask any follow-up question and to comment on our responses. We hope that they will consider raising their score at the end of our exchange.

How does VISPA perform on other types of neural networks, such as transformer architectures?

See FAQ 1 in the author rebuttal. We will include this discussion in Section 5 (Limitations and Open Problems) in the final version.

We have performed preliminary experiments of VISPA on the transformer architecture Swin-T[3] using ImageNet. We binarized the weights of the linear layers except for the qkv layer before the softmax. The experiment was performed with a learning rate of 0.1, a weight decay of 1e-5, and an epoch size of 150, without fine-grained hyper-parameter tuning. Results are shown in the table within the PDF attached to the author's rebuttal.

We have identified 3 methods to which to compare VISPA's performance : BiViT[1], BiT[2] and BiBERT[4]. As summarized in the table, such methods differ in fundamental ways from our application on VISPA, making an immediate comparison somewhat tenuous. All of them use Knowledge Distillation to enhance the performance of binarized transformer architectures, putting our method at a disadvantage. Furthermore, our strongest competitor BiViT uses a special handling of the softmax operation, facilitating its task. On the other hand, our implementation of VISPA does not binarize QKV layers, gaining an advantage there. With these caveats, we find that our method outperforms all competitors in terms of top-1% test accuracy, sometimes decisively. We believe that this result strongly suggests that the performance improvements, shown in this paper for vision-based transforms, will translate successfully to transformer. A more detailed study would involve different architectural choices on the binarized transformer, which are beyond the scope of our submission.

[1] He, Y., et al. "BiViT: Extremely Compressed Binary Vision Transformers." ICCV, 2023.
[2] Liu, Z., et al. "BiT: Robustly Binarized Multi-distilled Transformer." NIPS, 2022.
[3] Liu, Z., et al. "Swin Transformer: Hierarchical Vision Transformer using Shifted Windows." ICCV, 2021.
[4] Qin, H., et al. BiBERT: Accurate Fully Binarized BERT. ICLR, 2022.

Could you provide more insights into the impact of different hyperparameter settings, especially the covariance rank 𝐾, on the performance of VISPA?

See FAQ 2 in the author rebuttal.

Can you provide more details on the computational complexity and runtime of VISPA compared to other state-of-the-art BNN training methods?

As shown in Algorithm 1, let $M$ represent the batch size, $n$ the number of weight parameters, and $K$ the embedding dimension. The computational complexity of the VISPA approach is $O(Mn+nK)$. In comparison, other state-of-the-art BNN training approaches typically require $O(Mn)$. In most cases $M >> K,$ so that the increased cost due to the nK term is a negligible fraction of the total running time, as most time is spent performing forward- and back-propagation through the neural network.

Can you elaborate on the initialization strategy for the weight deviation matrix 𝑍 and its impact on the training stability and final performance?

In Section 4 (Experiments), we discussed that the weight deviation matrix $Z$ is initialized using the Xavier normal distribution, with a mean of 0 and a standard deviation of s * \sqrt{\frac{2}{\text{fan\\_in} + \text{fan\\_out}}}. Here, fan_in represents the number of input units, fan_out the number of output units, and s is a scaling factor parameter. Xavier initialization has been empirically shown to work well in practice across a variety of neural network architectures (e.g., CNNs, ResNets, GANs, and Transformers). It is based on a theoretical analysis that ensures the variance of the input to each neuron is similar to the variance of the output. This choice balances the scale of weights and helps maintain stable gradient magnitudes [1]. For details on the impact of $Z$'s initialization on training stability and performance, see Appendix A.2. There we show that VISPA achieves consistent performance across a range of choices of scaling factor s.

[1] Xavier Glorot and Yoshua Bengio, Understanding the difficulty of training deep feedforward neural networks, AISTATS, 2010.

What are the potential limitations of the proposed method when applied to resource-constrained devices, and how can these be mitigated?

The main potential limitation of VISPA on resource-constrained device is the increase in memory usage due to having to maintain the covariance variable $Z.$ If $n$ is the number of weights, our memory for storing relevant variables is at most $(K+1) n$ compared to $n$ for other methods. The total significance of this increase depends on the batch size and the size of the examples. For instance, in the case of ResNet18, there are roughly $9 \cdot 10^6$ weights while, for batch size=100, the total size of a data batch is $15 \cdot 10^6$, for a total memory usage of $24 \cdot 10^6$. Hence, choosing K=1 will lead to a usage of $33 \cdot 10^6$, a $37.5\\%$ increase. This increase will typically get smaller as we consider models trained on larger images. In practice, we often observe even smaller increases, as our estimate does not include memory usage due to backpropagation, which can be very large, e.g., in the case of residual connections.

A possible mitigation strategy is to store our variable $Z$ at lower precision. Given that this variable is only accessed via multiplication with Gaussian noise, we believe that this will not change the behavior of our algorithm.