BLoB: Bayesian Low-Rank Adaptation by Backpropagation for Large Language Models

Authors: yibin  wang
haizhou  shi
ligong  han
dimitris  metaxas
hao  wang
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 07:34:18

Abstract: Large Language Models (LLMs) often suffer from overconfidence during inference, particularly when adapted to downstream domain-specific tasks with limited data. Previous work addresses this issue by employing approximate Bayesian estimation after the LLMs are trained, enabling them to quantify uncertainty. However, such post-training approaches' performance is severely limited by the parameters learned during training. In this paper, we go beyond post-training Bayesianization and propose Bayesian Low-Rank Adaptation by Backpropagation (BLoB), an algorithm that continuously and jointly adjusts both the mean and covariance of LLM parameters throughout the whole fine-tuning process. Our empirical results verify the effectiveness of BLoB in terms of generalization and uncertainty estimation, when evaluated on both in-distribution and out-of-distribution data.

Main paper: 
 introduction: 
 Despite the recent advancements in Large Language Models (LLMs) [9, 106, 105, 66, 15, 4, 91, 92, 77, 12, 1, 71] , addressing the challenges of reliability and responsibility remains imperative [42, 100, 99] . LLMs often produce overconfident responses detached from factual grounding, posing potential harm to users [3, 107, 44, 43, 87, 50, 7, 118, 111, 120, 33, 68, 115, 45] . Therefore, accurately estimating response confidence (or uncertainty) is crucial to preemptively intervene before harm occurs. Current research predominantly focuses on eliciting the internal capability of uncertainty estimation of LLMs. For example, studies suggest that verbalized uncertainty yields better-calibrated results compared to conditional probability [87, 45] .While effective, the aforementioned methods do not offer a universal solution for expressing LLM uncertainty across all scenarios, especially when adapted [113] to domain-specific corpora, human preferences, or downstream tasks [44] . Even a well-calibrated LLM may struggle to estimate uncertainty during fine-tuning due to catastrophic forgetting of general knowledge [84] . Moreover, when applied to limited-scale downstream tasks, excessively over-parameterized LLMs can rapidly overfit, leading to overconfidence. Thus, enabling accurate uncertainty estimation of LLMs is vital for their reliable and responsible deployment.Bayesian methods emerge as a natural solution for learning uncertainty estimation abilities among their counterparts [88, 11, 101, 98, 29, 46, 51, 61, 97, 58, 102, 20, 110] . These methods model predictive uncertainty P (y|x, D) by marginalizing the posterior parameter distribution P (θ|D) after observing the dataset D: P (y|x, D) = P (y|x, θ)P (θ|D)dθ.(1)However, adapting the Bayesian framework to LLMs poses significant challenges. LLM architectures typically incorporate complex components, including non-linear activation functions, rendering exact Bayesian inference of parameter posteriors intractable, i.e., unable to compute the integral precisely. Consequently, finding an accurate approximation algorithm for the true posterior distribution becomes a primary challenge. Additionally, modeling parameter posterior distributions demands extra memory space, imposing a prohibitive burden on systems due to the massive scale of LLMs.Contemporary methods leverage Parameter-Efficient Fine-Tuning (PEFT) to reduce the number of tunable parameters, thus alleviating computational and storage resource burdens [23, 41, 26, 121, 56, 53] . Built on this, recent research explores Bayesianizing only the PEFT module during fine-tuning to calibrate LLMs [8, 103, 116, 69] , somewhat relieving the burden of introducing more parameters for posterior approximation. However, initial investigations suggest that straightforward combinations of PEFT and basic Bayesian techniques like Monte-Carlo Dropout (MCD, [29] ) or Deep Ensemble (ENS, [51, 8, 103] ) yield only marginal improvements in generalization and uncertainty estimation. The most promising results to date involve Kronecker factorized Laplace approximation, applied after maximum a posteriori (MAP) estimation provided by any optimization algorithm [116] . Nevertheless, we argue that such post-training procedures bifurcate posterior approximation into two stages, inevitably leading to suboptimal estimation.To address this challenge, we propose Bayesian Low-Rank Adaptation by Backpropagation (BLoB), a Bayesian Deep Learning framework for fine-tuning LLMs with LoRA. BLoB jointly estimates the low-rank variational distributions' mean and covariance throughout the entire fine-tuning stage via backpropagation. Unlike methods relying on post-training approximation, BLoB enables simultaneous estimation of both the parameter mode (i.e., the mean if one assumes Gaussian distributions) and the parameter variance. Random sampling of model parameters based on variance estimation can enhance mode estimation. It thereby improves model performance in terms of accuracy and uncertainty estimation on both in-distribution and out-of-distribution datasets, as verified by our extensive experiments across multiple datasets. In summary, our contributions are:• We propose a principled Bayesianization framework for Low-Rank Adaptation (LoRA) in Large Language Models (LLMs) by assuming that full weights' approximate posterior distribution has a low-rank structure containing a linear combination of independent Gaussian distributions.• We show that, under mild conditions, optimization of the full-weight variational distribution can be done efficiently in the low-rank space of the weight update matrices.• We introduce BLoB, a variational Bayesian low-rank adaptation framework for LLMs that jointly learns the mean and covariance of the variational distribution during fine-tuning.• Extensive evaluations demonstrate the superiority of BLoB in terms of generalization and uncertainty estimation across different scenarios. 
preliminaries: 
 In this section, we describe the notation as well as some preliminaries.Notation. In this paper, scalars are denoted by lowercase letters, vectors by lowercase boldface letters, and matrices by uppercase boldface letters. Probability, expectation, and the dataset are denoted by P , E, and D, respectively. We use [m] = {1, 2, • • • , m} to denote the set of consecutive integer numbers starting from 1 and ending at m. For a matrixX = [x 1 , • • • , x n ] ∈ R m×n , we use vec(X) = [x ⊤ 1 , x ⊤ 2 , • • • , x ⊤ n ] ⊤ ∈ R (mn)×1to denote the vectorization operation; we use∥X∥ p = ij |X ij | p 1 /pto define the p-norm of a matrix. We use ⊗ and • to denote the Kronecker product and the element-wise product, respectively.Inspired by the pioneering work on identifying and leveraging the low intrinsic rank of overparameterized models during fine-tuning [55, 2] , Low-Rank Adaptation (LoRA) assumes a low rank for the network's weight updates [41] . Typically in a single linear layer, LoRA decomposes each update matrix ∆W = BA into the product of two low-rank matrices, where B ∈ R m×r and A ∈ R r×n . Here, m, n, and r denote the number of input neurons, output neurons, and the rank of the decomposition, respectively [41] . The forward pass of the linear layer with LoRA is formulated as:EQUATIONwhere h and z denote the input and output of the layer. Since the rank r ≪ min{m, n} is significantly smaller than the numbers of input and output neurons (e.g., r = 8 ≪ m = n = 4096 in the attention layer [41] ), LoRA can drastically reduce the number of trainable parameters by approximately three orders of magnitude compared to full-parameter fine-tuning, while achieving comparable performance to the full-rank fine-tuning. This also leads to a similar reduction in memory consumption for storing optimizer states, thereby reducing the hardware requirements for fine-tuning LLMs to a great extent.Bayesian Neural Networks (BNNs) estimate the posterior distributions of network parameters rather than relying on single-point estimates [10, 102] . Due to the intractability of exact inference of the true posterior, Variational Bayesian Networks (VBNs) approximate the true posterior using a variational distribution; this is done by minimizing its KL divergence from the true posterior distribution [39, 30, 11] . Specifically, if the weights W 's variational distribution q(W |θ) is parameterized by θ, minimizing the divergence KL[q(W |θ)∥P (W |D)] is equivalent to minimizing the following variational free energy with respect to θ [67, 117, 28] :EQUATIONThe final formulation of the objective function in Eqn. 3 offers another interpretation beyond minimizing the KL divergence between the variational and true posterior distributions [11] . Specifically, the first term maximizes the likelihood of the data, while the second term regularizes the variational distribution q(W |θ). We refer to the first term as the likelihood cost and the second term as the complexity cost. Optimizing these two terms involves balancing the expressiveness of the approximate posterior distribution and its simplicity.Optimizing the first term of Eqn. 3 requires integrating out the parameterized variational distribution, necessitating Monte Carlo gradient estimation [52, 81] . Using this approach, we can incorporate the re-parameterization trick to enable backpropagation of the gradient to the underlying parameter θ [72, 47, 78] . In Bayes By Backprop (BBB) [11] , the variational distribution is further simplified as a diagonal Gaussian N (µ, σ 2 ), where σ = log(1 + exp(ρ)) ensures the standard deviation is positive. Then we have the Monte-Carlo estimation of Eqn. 3 that can pass the gradient to θ:EQUATIONwhere W k = µ + log(1 + exp(ρ)) ⊙ ϵ k is the k-th sample of the weights yielded by parameterization and ϵ k ∼ N (0, I). In BBB, the authors assume the prior distributionEQUATION2 ) to be a mixture of Gaussians. Consequently, they optimize the second term based on weight sampling. In different scenarios, a simpler form of the prior, which allows for a closed-form solution, can also be considered. Although our proposed method is largely based on the existing framework of BBB, trivially combining BBB with LoRA does not yield satisfactory results. It is important to note that our specific designs are necessary to encourage the fast convergence of the variational distribution, which will be introduced later in Sec. 3. 
methodology: 
 In this section, we formally introduce our proposed method, Bayesian Low-Rank Adaptation by Backpropagation (BLoB). We begin by discussing the design choices for Bayesianizing LoRA parameters in Sec. 3.1, highlighting the assumptions BLoB makes about the approximate posterior in the full-weight space. Next, in Sec. 3.2, we explore the low-rank structure of the prior distribution in the full-weight space, which in turn motivates our choice of prior distributions in the low-rank parameter space. In Sec. 3.3, we introduce our parameterization method for the variational distributions. In Sec. 3.4, we integrate Flipout [108] into LoRA for improved sampling efficiency and faster convergence. Finally, we present the complete algorithmic description of BLoB in Sec. 3.5. Proof of the theorems and claims in this section can be found in Appendix A. 3.1 Low-Rank Variational Approximate Posterior Distribution: LoRA Bayesianization Asymmetric LoRA Bayesianization. In LoRA [41] , the weights are treated asymmetrically. A is randomly initialized, usually from the standard normal distribution or using Kaiming initialization, while B = 0 is initialized as a zero matrix to ensure that the model fully retains the capabilities of the pre-trained weights at the start of fine-tuning. The trivial solution of estimating the variational approximate posterior for the entire set of LoRA parameters can significantly hinder training convergence. For example, consider the Gaussian posteriors q(A|θ) = N (A|M A , Ω 2 A ) and q(B|θ) = N (B|0, Ω 2 B ), where Ω A and Ω B are variance estimates added to A and B, respectively. Although the expectationE A,B [(W 0 + BA)x] = W 0 x + E A,B [BAx] = W 0 xpreserves the functionality of the pre-trained model, accurate estimation requires an impractically large number of weight samples. Such variational distributions lead to significant fluctuations during the early stages of fine-tuning, unless the initial variance of B, Ω B → 0 + , is intentionally minimized towards zero. Therefore, we take an asymmetric approach to initialize Ω B = 0 and keep it fixed throughout the fine-tuning process. This, in effect, gives up Bayesian modeling of the B component and focuses only on the posterior of A in LoRA, as shown in Fig. 1 .Additional Advantages. In addition to reducing sampling noise and improving convergence speed, our Bayesianization design has two further advantages. First, compared to modeling the variational distributions of both A and B, our approach significantly reduces additional memory cost by approximately 50% per layer. Second, our design is equivalent to finding a posterior estimate for the full-weight matrix with a low-rank structure. For instance, by assuming a deterministic B and Bayesianizing P (A|θ) = N (A|M , Ω 2 ), each element of the full weight matrix W ij is calculated asEQUATIONwhereA kj ∼ N (M kj , Ω 2 kj ) is drawn independently ∀k ∈ [r].It is noteworthy that due to the low-rank structure defined in Eqn. 5, the full-weight parameters of W are no longer independent from each other. The correlation among them can be reflected by the following theorem: Theorem 3.1 (Variational Distribution of the Full-Weight Matrix in BLoB). With the pre-trained weight matrix W 0 ∈ R m×n and the low-rank weight update matrix B ∈ R m×r , suppose that the variational distribution of the other low-rank update matrix A ∈ R r×n is Gaussian withq(A|θ = {M , Ω}) = ij N (A ij |M ij , Ω 2 ij ), where M = [M ij ] ∈ R r×n and Ω = [Ω ij ]∈ R r×n are its mean and standard deviation, respectively. The equivalent variational distribution defined on the full weight matrix W as in Eqn. 3 is given by q(vec(W )|B, θ) = N (vec(W )|µ q , Σ q ), (6) where µ q = vec(W 0 + BM ),EQUATIONTheorem 3.1 shows that our asymmetric LoRA Bayesianization is equivalent to using a Gaussian variational distribution for the full weight W (i.e., Eqn. 6), with a flexible covariance matrix (i.e., Eqn. 8), to approximate the postetior distribution of the full weight W .Remark. The covariance matrix Σ q is strictly singular, which consequently inspires us to design a prior P (W ) with such low-rank structure in Sec. 3.2. Previous work on low-rank Gaussians typically considers covariance with a similar structure D 2 + Σ q , where D is diagonal [89, 70, 86, 90, 73] . However, sampling from a Gaussian with this structure requires sampling noise of the same shape as the full-weight matrix, which is not parameter-efficient; we therefore do not adopt this in our work.In Eqn. 3, optimizing the KL divergence between the variational and prior distributions in the space of full weights can be burdensome. Therefore, we assume the prior distribution of the full weights to be a low-rank Gaussian, with its mean centered at the pre-trained weights vec(W 0 ) and its covariance matrix parameterized by a rank-r ′ matrix R ∈ R (mn)×r ′ :EQUATIONΣ p = R R ⊤ .Assuming a low-rank prior distribution and designing an appropriate R allows us to optimize the KL divergence in the decomposed low-rank weight space, as suggested by the following theorem. Theorem 3.2 (Efficient Computation of Full-Weight KL Divergence). Suppose the pre-trained weights W 0 , update matrix B, and the variational distribution q(A|θ) are defined as in Theorem 3.1, and the prior distribution of the full-weight matrix P (vec(W )) is defined as Eqn. 9. Consider the Gaussian prior distributionEQUATIONp ); we then have:EQUATIONif R = [σ p I n ⊗ R], where R satisfies RR ⊤ = BB ⊤ .Theorem 3.2 shows that with a proper R, one can compute the KL divergence for the high-dimensional full weight vec(W ) simply by computing the KL divergence for A, which is much lower-dimension, more parameter-efficient, more memory-efficient, and faster. Note that the Gaussian distributions we define for both the prior and the posterior are degenerate. However, they are valid for probabilistic inference [83] , as (i) their probability density is well-defined, and (ii) their KL divergence is computable under the assumptions of Theorem 3.2. See Appendix A.1 for a detailed discussion.Concretely, we assume that the prior distribution in BLoB follows the low-rank structure described in Theorem 3.2 and minimize the KL divergence term for the low-rank component A using its analytical solution in Eqn. 3:EQUATIONThe parameterization of the Gaussian variational distribution q(A|θ) significantly affects the convergence speed of the KL term in Eqn. 11. The mean matrix M of q(A|θ) has no additional constraints, we therefore parameterize it directly as the output of a neural network. Each entry of q(A|θ)'s diagonal covariance matrix Ω (i.e., standard deviation) is non-negative; we therefore use element-wise parameterizationΩ ij = G 2 ij , where G = [G ij ] ∈ R r×nis the real parameter matrix that determines the standard deviation Ω. Since Ω is usually initialized with small positive values close to zero, our parameterization method provides large gradients initially, contributing to the rapid decrease of the KL term. We further show, both theoretically and empirically, that our parameterization method, unlike BBB's softplus function log(1 + exp(•)), is crucial for the fast convergence of Ω when q(A|θ) is close to the prior distribution P (A) (see more analysis in Appendix A.2).Improving Sample Efficiency with Flipout. One main challenge in estimating the variational distribution (i.e., the approximate posterior) during fine-tuning lies in the sample efficiency of the Algorithm 1 Bayesian Low-Rank Adaptation by Backpropagation (BLoB) Require: dataset D, pre-trained weight W 0 , low-rank component B, θ = {M , G} for parameterizing the mean and variance of A; Require: prior standard deviation σ p , initialization hyperparameter ϵ, number of input features n; Require: number of samples during training K, number of iterations T , learning rate η;1: G ∼ U( ϵ √ 2 , ϵ), M ∼ U -6 n , 6 n ▷ Initialization of A's parameters. 2: B ← 0 ▷ Initialization of B. 3: for t = 1, • • • , T do 4:Sample a mini-batch of data D t ∼ D containing b samples. 5:for k = 1, • • • , K do 6: Sample batched noise E k ∼ N (0, I).▷ Sample the noise.Let{ E kj } b j=1 ← BLoBFlipout(E k ).▷ Eqn. 12 8:Let{A kj = M + G 2 • E kj } b j=1 . 9:end for 10:LetF t = -1 Kb K k=1 b j=1 log P (D t |A kj , B) + 1 2σ 2 p (∥M ∥ 2 2 + ∥G∥ 4 2 ) -2 ij log G ij .11:▷ Eqn. 13 and 11.12:Calculate the gradient w.r.t. the parameters:∆ M = ∂ Ft /∂M, ∆ G = ∂ Ft /∂G, ∆ B = ∂ Ft /∂B. 13:Update the parameters:M ← M -η∆ M ; G ← G -η∆ G ; B ← B -η∆ B . 14: end forweights [109, 25, 58] . During mini-batch stochastic gradient descent, a batch of examples typically share the same weights drawn from the variational distribution. This can lead to slow convergence of the likelihood cost in Eqn. 3. Drawing inspiration from [108] , we introduce the technique of flipout to speed up the sampling procedure of our low-rank variational distributions q(A|θ).LoRA Flipout. Unlike the original approach, which applies rank-1 random flipping to the full weights, we apply flipout exclusively to the low-rank component A. Specifically, suppose we have a mini-batch of input vectors H ∈ R n×b , where b represents the batch size. We randomly sample two low-rank flipping matrices S ∈ {-1, +1} n×b and T ∈ {-1, +1} b×r . Denoting as E ∈ R r×n the weight noise sampled for this mini-batch, the batched output Z after applying flipout is thenEQUATIONIt is crucial that the independent noises added to the low-rank weight noise ∆A ≜ E • Ω ensure sampling independence across examples within a mini-batch, thereby enhancing the sampling efficiency of the algorithm. This is done without violating the assumptions outlined in Theorem 3.1 and 3.2. As illustrated in algorithm 1, we use E kj to represent the equivalent noise applied to parameter A for the j-th example in the k-th batch after BLoBFlipout. Due to the low-rank structure of our Bayesianization method, the computational overhead of employing flipout in BLoB is also minimal.We are now ready to present our full BLoB algorithm.During training, under the assumptions outlined in Theorem 3.1 and 3.2, optimizing the evidence lower bound on the full weight W can be efficiently done in the low-rank space, using the following final objective function:EQUATIONwhere θ = {M , Ω} denotes the set of the parameters underlying the variational distribution of the low-rank matrix A. Additionally, to trade off between data fitting and posterior approximation, During inference, for an input x, we approximate the expected output distribution P (y|x) of BLoB by drawing N samples from the variational distribution q(W |θ). Empirically, N = 10 provides a good balance between estimation quality and computational efficiency:EQUATIONQuestion: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact.• Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
experiments: 
 In this section, we compare our BLoB with existing methods on real-world datasets. Sec. 4.1 introduces the experimental settings, including baselines, fine-tuning, and evaluation protocols.We then evaluate BLoB's generalization and uncertainty estimation abilities in both in-distribution (Sec. 4.2) and out-of-distribution scenarios (Sec. 4.3).Fine-tuning and Evaluation. We implement BLoB in the PEFT library [63] and fine-tune the LlaMA2-7B [92] model on common-sense reasoning tasks. Following Laplace-LoRA [116] , we apply LoRA to the output layer as well as the queries and values of all the attention layers. For hyperparameters, we strictly adhere to the default settings in the PEFT library and the original LoRA paper [63, 41] to ensure maximal reproducibility. This includes the number of training steps, learning rate, and LoRA rank r (see Appendix B.1 for details). For common-sense reasoning tasks, we select the next token logits corresponding to possible answers from each dataset and fine-tune the LLM to maximize the likelihood of the correct token. For evaluation, in addition to Accuracy (ACC), we use Expected Calibration Error (ECE [31] ) and Negative Log-Likelihood (NLL) to assess the models' uncertainty estimation ability (see Appendix B.2 for details).Baselines and Implementation Details. We compare BLoB with state-of-the-art uncertainty estimation methods applied to the LoRA adapters of LLMs, including Monte-Carlo Dropout (MCD) [29] , Bayes By Backprop (BBB) [11] , Deep Ensemble (ENS) [51, 8, 103] , and the latest Laplace-LoRA (LAP) [116] . We also report the performance of two standard PEFT baseline methods for reference: Maximum Likelihood Estimation (MLE) [41] and Maximum A Posteriori (MAP).For MLE, we use the LoRA implementation. For MAP, we use a weight decay rate of 1e -5. For MCD, we use an ensemble of 10 LoRAs with a dropout rate of p = 0.1. For ENS, we independently fine-tune 3 LoRAs and average their logits during evaluation. For BBB, we adopt the default settings from the Bayesian-Torch library [48] and only Bayesianize the A matrix, similar to BLoB. We sample N = 10 times for BBB during test. We re-implement LAP and apply it to the MAP checkpoints. We keep all BLoB-specific hyperparameters consistent across all datasets. Typically, we set the number of samples K = 1 during training for all our BLoB experiments, which highlights BLoB's sampling efficiency. As shown in Table 1 , we also report BLoB's performance with different numbers of samples during Bayesian inference, where N = 0 indicates directly using the mean of the weight distribution for prediction.We fine-tune Llama2-7B on six common-sense reasoning tasks: Winogrande-small (WG-S), Winogrande-medium (WG-M) [82] , ARC-Challenge (ARC-C) [18] , ARC-Easy (ARC-E) [18] , Open-BookQA (OBQA) [65] , and BoolQ [17] . For all baseline methods, using the same pre-trained LLM backbone, we maintain consistent hyperparameters across all datasets and do not use additional validation sets to achieve higher performance (See Appendix B.3 for detailed settings). showcasing the advantage of jointly optimizing the mean and covariance of the variational weight distributions during fine-tuning. However, there remains considerable room for improvement.BLoB consistently achieves better or comparable performance across all datasets. With the number of samples during inference set to N = 10, the same as MCD, BLoB provides the best uncertainty estimation performance, significantly reducing NLL and ECE, and greatly mitigating overconfidence while maintaining comparable or better ACC than MLE. Even with half the number of samples, N = 5, BLoB still delivers performance comparable to that of N = 10 and outperforms other baselines on most datasets. By abandoning the modeling of the posterior distribution, prediction using the mean of the weight distribution, i.e., BLoB (N=0) sacrifices some degree of calibration in exchange for improved accuracy. Appendix C.5 presents the trade-off between accuracy and calibration, which is controlled by the standard deviation of the prior Gaussian distribution).Besides Llama2-7B, we also include additional results for RoBERTa-base [60] on text classification tasks in Appendix C.1. Our method consistently achieved either the best or runner-up performance across nearly all datasets, demonstrating its versatility across different architectures.1 LAP encounters training failures on the ARC-C dataset under a unified setting, where the same hyperparameter configuration is applied across all datasets. To achieve competitive performance with LAP on the ARC-C dataset, we deviate from this unified setting, allowing dataset-specific hyperparameters for LAP. Note that this introduces an unfair advantage for LAP. Specifically, for the LAP on ARC-C dataset, we set the dropout rate to 0.1, the learning rate to 5e-5, and applied early stopping at the 5,000-th iteration (out of 10,000). We use models fine-tuned on OBQA [65] to evaluate the generalization ability of different methods under distributional shifts. OBQA consists of multiple-choice elementary-level science questions. We categorize the distributional shifts into two types: smaller and larger shifts. The ARC [18] dataset, which also consists of multiple-choice science questions, represents a smaller distributional shift. The college-level chemistry and physics subsets of MMLU [38] represent larger distributional shifts.  
related work: 
 Parameter-Efficient Fine-Tuning (PEFT) for LLMs. Due to the prohibitively large size of LLMs, parameter-efficient fine-tuning has become a trending topic. Computational paradigms in this area include adapter-based fine-tuning [40, 36, 80, 75, 62] , prompt-based fine-tuning [34, 54, 59, 57, 94, 6] , and partial fine-tuning [119, 124, 5, 112, 32] . Among these, LoRA [41] has gained significant attention due to its simplicity and effectiveness. Building on LoRA, numerous studies have aimed to further optimize parameter efficiency when fine-tuning large models [26, 35, 22, 21] . For instance, KronA models weight updates as the Kronecker product of two smaller matrices without decreasing the update rank [26] , and SVDiff performs Singular Value Decomposition (SVD) on the original weight matrices, fine-tuning only the singular values [35] . However, in this paper, we focus solely on Bayesianizing LoRA due to its widespread application in existing works. We also note that BLoB can be naturally adapted to handle different LoRA variants.Uncertainty Estimation in Large Language Models. Large-scale pre-trained models are wellcalibrated during pre-training [44] , but fail to accurately express predictive uncertainty during inference [3, 107, 44, 43, 87] , especially after fine-tuning [8, 103, 116, 69] . This indicates that measures effective during pre-training [93, 114, 16, 122, 14] may lose their power of uncertainty estimation after fine-tuning for domain-specific knowledge. To address this issue, [27, 123] define priors and approximate posteriors on the full attention weights during fine-tuning, achieving better uncertainty estimation but at a significant cost in time and space. Consequently, recent work integrates Bayesian methods and PEFT for efficient uncertainty estimation. For instance, [8, 103] train and store multiple copies of different LoRAs, ensembling their outputs during inference to achieve somewhat better results. [116] applies Kronecker factorized Laplace approximation on fine-tuned LoRA. However, such post-training procedures bifurcate posterior approximation into two stages, leading to suboptimal estimation. In contrast, our BLoB enables simultaneous estimation of both the mean and covariance of LLM parameters in a single fine-tuning stage, substantially improving performance. 
conclusion: 
 In this work, we propose a principled Bayesianization framework for parameter-efficiently fine-tuning LLMs. Our theoretical analysis shows that a full-weight variational distribution can be efficiently optimized by approximately using a low-rank space of the weight update matrices. Our empirical evaluations corroborate this theoretical insight, demonstrating superior generalization and uncertainty estimation capabilities across diverse scenarios compared to various baseline methods. Building on LoRA, our approach seamlessly integrates with existing LLM architectures while imposing minimal additional memory overhead and training time. Our method highlights that jointly learning the mean and covariance of the variational distribution during fine-tuning can mutually improve both, underscoring the powerful potential of Bayesian methods in enhancing the reliability and generalization of LLMs. 
limitations: 
 The main limitations of our proposed BLoB method are: (i) BLoB is confined to fine-tuning scenarios and is not applicable to training-free tasks, such as direct uncertainty estimation during inference [64] .(ii) As a typical mean-field variational inference method, BLoB requires multiple sampling iterations during inference, which challenges stable and efficient deployment. (iii) While BLoB's effectiveness has been empirically demonstrated for downstream classification tasks, its application to generation tasks requires further investigation.r i=1 b i • a ⊤ i .We then rewrite vec(W ) in the form of Kronecker product ⊗:EQUATIONWe then calculate the covariance matrix Σ q asEQUATIONEQUATIONEQUATIONcompleting the proof.It is crucial to note here, the final covariance matrix of q(vec(W )) follows a block-diagonal structure, which will be further utilized for the proof of Theorem 3.2. Defining Σ i = diag(Ω 2 i: ), we have:EQUATIONAnother important fact about Σ q is its singularity. It can be seen directly as we consider the rank of any one of the block matrixEQUATION)where r is the rank of LoRA, strictly smaller than the output dimension of m. Theorem 3.2 (Efficient Computation of Full-Weight KL Divergence). Suppose the pre-trained weights W 0 , update matrix B, and the variational distribution q(A|θ) are defined as in Theorem 3.1, and the prior distribution of the full-weight matrix P (vec(W )) is defined as Eqn. 9. Consider the Gaussian prior distributionEQUATIONp ); we then have:KL[q(vec(W )|B, θ)∥P (vec(W ))] = KL[q(A|θ)∥P (A)], if R = [σ p I n ⊗ R], where R satisfies RR ⊤ = BB ⊤ .Proof. We start by assuming the low-rank structure of the prior P (vec(W )), and then reveal the conditions reaching to our final conclusion step by step.Typically, for two Gaussian distributions q and p whose covariance matrices Σ q ∈ R d×d and Σ p ∈ R d×d are both full-rank, and their means as µ q ∈ R d and µ p ∈ R d , we have their KLdivergence asEQUATIONThe singularity of the covariance matrices of P (vec(W )) and q(vec(W )), i.e., |Σ q | = |Σ p | = 0, can cause issues when computing the KL-divergence as it includes the log-determinant term. Therefore in this proof, we consider the alternative of the covariance matrices, where an extremely small diagonal elements are added.For the prior distribution, following the alternative form of a degenerate Gaussian [83] , as suggested in Eqn. 9, we assumeEQUATIONwhereΣ p = λI + R R ⊤ , (λ → 0 + ).By default, we assume that the low-rank tall matrix R ∈ R (mn)×r ′ has the full column rank r ′ . Otherwise if r( R) = r ′′ < r ′ , then we can in effect consider a new matrix component R ′ ∈ R (mn)×r ′′ that has the same rank as r ′′ , which satisfies our assumption of full column rank. Therefore, we have the SVD decomposition of R is given byEQUATIONwhereU R ∈ R (mn)×(mn) and V R ∈ R r ′ ×r ′ are orthonormal, i.e., U R U ⊤ R = U ⊤ R U R = I (mn) and V R V ⊤ R = V ⊤ R V R = I r ′ .D R is a tall matrix where its upper part is diagonal and the lower part is a zero matrix, denoted asD R = [D * R , O] ⊤ = [diag([d R1 > 0, d R2 > 0, • • • , d R r ′ > 0]), O] ⊤ .For the approximate posterior q(vec(W )|B, θ), we considerEQUATIONwhereΣ q = λI + BΣ B ⊤ , (λ → 0 + ),where we simplify the notation for the covariance matrix Σ q by defining B = [I n ⊗B] ∈ R (mn)×(mr) and Σ = diag(vec(Ω) 2 ). Likewise, we have the SVD-decomposed matrices for B where they are defined in the same way as Eqn. 29:EQUATIONwhere U B and V B are orthogonal matrices, andD B = [D * B , O] ⊤ = [diag([d B1 > 0, d B2 > 0, • • • , d Bmr > 0]), O] ⊤ .First, we calculate the log-determinant part of the KL-divergence. For the log-determinant of the covariance matrix of the prior distribution Σ p , by applying SVD decomposition in Eqn. 29, we haveEQUATIONFollowing (almost) the same idea, we now have the log-determinant of the variational distribution's covariance Σ q asEQUATIONEQUATION)We make two observations when λ → 0 + : (i) compare the terms that contain log λ on both sides, to make sure the log-determinant in the divergence term bounded, we have to set r ′ = mr; (ii) the last term in Eqn. 39,EQUATIONNext we calculate tr(Σ -1 p Σ q ) in Eqn. 27. Following the same assumptions and notations above, we have the inverse of the covariance of the prior distribution asEQUATIONHence we haveEQUATIONBy using the conditionEQUATIONand there exists an orthogonal matrix P ∈ R (mr)×(mr) , such thatEQUATION)The SVD decomposition on R can then be formulated as:EQUATIONEQUATIONSubstituting U R , D R , V R back to Eqn. 42 and applying λ → 0 + , we haveEQUATIONEQUATION= (mn -nr) + σ -2 p tr(Σ).For the quadratic term in Eqn. 27, the pre-trained weights W 0 cancel out, and we can calculate it asvec(BM ) ⊤ Σ -1 p vec(BM ) (50) =[M ⊤ :1 B ⊤ , • • • , M ⊤ :1 B ⊤ ]    (BΣ 1 B ⊤ ) -1 . . . (BΣ n B ⊤ ) -1       BM :1 . . . BM :n    (51) = n i=1 M ⊤ :i B ⊤ (BΣ i B ⊤ ) -1 BM :i (52) = n i=1 M ⊤ :i (V B      d 2 B 1 σ 2 p (d 2 B 1 +λ). . .EQUATIONFinally, proof is completed by combining Eqn. 40, Eqn. 49, and Eqn. 55. 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA]Justification: The paper does not pose any risks since there is no release of models.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes]Justification: we perform all the experiments under the license, and have cited work properly in the paper.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.• If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA] Justification: No crowdsourcing experiments are involved in this paper.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA]Justification: No IRB is involved in this paper.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.• We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: General Analysis on Parameterization. Consider a path of parameterization for a single variable:EQUATIONwhere ρ is the real parameter we perform update on, f is our parameterization choice for the variable σ, and l represents the loss function we aim to minimize. When comparing two different parameterization methods, we consider the same initial conditions of σ = σ 0 , and we assume the same step size η on the real parameter ρ. To show the influence of the choice of parameterization, we calculate the decrease of the loss value by performing one step of gradient descent. First, by the chain rule, the gradient w.r.t. ρ 0 is calculated asEQUATIONAfter one step of the gradient descent, we have ρ 1 asEQUATIONand the loss value decreased at ρ 1 can be approximated by the first-order Taylor expansion,EQUATIONEQUATIONEQUATION)Since the initialization of the different parameterization variable ρ 0 is set to ensure the same initial condition of σ 0 for different f s, we can see that the amount the loss decreases by after one step of update is proportional to the squared gradient ∆L ∝ (l′ (σ 0 )) 2 • (f ′ (ρ)) 2 = ( dL /dρ) 2 .Parameterization with log(1 + exp(•)) or (•) 2 ? Previous VBNs [11, 49] typically use a softplus function σ q = log(1 + exp(ρ)) to parameterize the standard deviation. For a single element ρ, the derivative of the closed-form solution of the KL divergence in Eqn. 11 is calculated asEQUATIONDue to the fact that σ q is typically initialized to a small value close to 0 to ensure stable optimization of the likelihood cost term (e.g., 1e -3), and in order to ensure that the model obtains good uncertainty, σ p is usually set to a larger value close to 1 (e.g., 1e -1). In this case, the derivative of ρ in Eqn. 63 is almost always a constant -1, which, based on our previous analysis, leads to slow convergence for large σ p values when σ q is small. Therefore, we parameterize σ q with quadratic function: σ q = ρ 2 . In this case, the derivative of the KL divergence with respect to ρ in Eqn. 11 becomes:EQUATIONUnder the same initialization conditions, the derivative in Eqn. 64 is approximately of the order of ρ -1 , leading to rapid convergence towards larger σ p values when σ q is small. Building on this, we use SGD without momentum to optimize the complexity loss term, thereby achieving the natural convergence of σ q .Visualization. To visually demonstrate the differences in the convergence of KL divergence during training with these two parameterizations, we set σ p = 1 and employ gradient descent to optimize the KL divergence. As introduced in B.1, in practical mini-batch gradient descent, the KL divergence is weighted by 1 /#mini-batches. Therefore, assuming there are 100 mini-batches, the learning rate is set to 0.01, which translates to an actual learning rate of 1e -4 for ρ. We initialize σ q = 0.01 for both parameterizations. The growth of σ q during the KL training, for σ q = log(1 + e ρ ) and σ q = ρ 2 is shown in Fig. 2 . In the same setting, the softplus parameterization takes nearly 100,000 gradient steps to converge, while the square parameterization takes only about 5,000 gradient steps. This modification makes it suitable for our fine-tuning setup, where the number of gradient steps is relatively small. The growth curve of σ q = log(1 + e ρ ) and σ q = ρ 2 during the optimization of KL divergence (without data likelihood). The number of gradient steps (5000) is marked with the red line.For the i-th input vector h i in a mini-batch, we randomly sample two flipping vector s ∈ {-1, +1} n and t ∈ {-1, +1} r . Denoting A as the weight matrix sampled from posterior distribution, and ∆A as the batched noise for sampling A, the output vector z i after applying flipout is:EQUATIONSimilarly, for a mini-batch input matrix H ∈ R n×b with batch size b, we randomly sample two low-rank flipping matrices S ∈ {-1, +1} n×b and T ∈ {-1, +1} b×r . The batched output matrix Z after applying flipout is then: EQUATION)In this section, we first introduce the implementation details of BLoB in Appendix B.1, including the KL Re-weighting scheme, initialization of the parameters, and learning scheduling, etc. Next, we introduce the two evaluation metrics for uncertainty estimation in Appendix B.2. Finally, we present some statistics of the adopted datasets in Appendix B.3.KL Re-weighting. In mini-batch SGD, the training data D is randomly divided into M equally sized subsets: D 1 , D 2 , . . . , D M . For mini-batch i = 1, 2, . . . , M , the cost function is:EQUATIONwhere λ i ∈ [0, 1] andM i=1 λ i = 1.There are various approaches for controlling the weight of KL divergence. [30] utilizes λ i = 1 /M, while [11] adopts λ i = 2 M -i /2 M -1.In fine-tuning tasks, we found that using a scheduler with λ i = 2 i /2 M -1 performs well. This allows the model to find good fits to the data points within the early stages and then optimize the complexity cost in later stages.In multiple epochs of mini-batch SGD, larger datasets require more iterations to complete one epoch, resulting in delayed convergence of the complexity cost. To enhance the stability of BLoB's performance across datasets with varying sizes, we pseudo-rescaled the size of the training dataset to make smaller datasets slightly larger and larger datasets slightly smaller. For the portions of the dataset that required expansion, we incorporated additional mini-batches from subsequent epochs. Conversely, for the datasets needing reduction, we deferred the excess mini-batches to subsequent epochs. We denote the size of original dataset as L 0 , The rescaled dataset size is:EQUATION)where γ is a coefficient used to control the scaling magnitude, and we set it to 8 in all experiments.The pseudo-rescaling does not affect the likelihood cost in practical mini-batch gradient descent. In fact, it only changes the warm-up period in KL reweighting from M to L * /batch size, thereby facilitating more consistent optimization of the complexity cost across datasets of different sizes.Additional Details. We initialize standard deviation parameterization matrix G by element-wise sampling from a uniform distribution with a range of [ ϵ √ 2 , ϵ], while keeping the remaining initialization settings consistent with LoRA. To maintain consistency, we use the same learning rate scheduler and warmup ratio for the optimizer of the KL term as we do for the likelihood term. We sample only once during the training process. During inference, we sample N times, then take the average of the logits obtained after passing through the softmax function. Detailed hyperparameter settings are provided in the Table 4 . Table 3 provides the hyperparameters for fine-tuning with LoRA shared with other baselines. Our experiments on Llama2-7B were conducted using 2 NVIDIA RTX A5000 GPUs for parallel training, while experiments on RoBERTa-base were conducted using 4 NVIDIA RTX A5000 GPUs for parallel training.Negative Log-Likelihood (NLL) and Expected Calibration Error (ECE [31] ) are two prevalent metrics for assessing uncertainty estimation. NLL calculates the sum of the negative expected log probability of predicting the actual label. Suppose this predicted probability is given by the model P θ , and we have a test dataset {x n , y n } N n=1 of size N . Then the NLL measured on this dataset isEQUATION)This metric prefers models that assign higher probabilities to correct labels. If the model exhibits overconfident in an incorrect prediction, the probability assigned to the correct label will be diminished, thereby increasing the NLL.On the other hand, ECE measures how well the model's confidence matches its accuracy. This is done by binning the predictions based on their confidence levels and then computing a weighted average of the absolute difference between accuracy and confidence within each bin:EQUATIONwhere acc(B m ) and conf(B m ) denote the average accuracy and confidence within bin B m , respectively. These are given by:EQUATIONwhere |B m | is the number of samples in bin m. We set |B m | = 15 across all experiments.Table 5 summarizes the size of the training set and the number of labels for each dataset. Table 6 summarizes the prompt templates used for common sense reasoning tasks. This section provides additional experimental results omitted from the main body of the paper due to space limitations. First, we present the results of BLoB when applied to RoBERTa, another pre-trained language model, in Appendix C.1. Next, in Appendix C.2, we conduct the ablation Table 7 : Performance of different methods applied to LoRA on RoBERTa-base pre-trained weights. The evaluation is undertaken on five GLUE [96] and SuperGLUE [95] tasks, with a shared hyper-parameter setting without using individual validation dataset. "↑" and "↓" represent that higher and lower values are preferred, respectively. The boldface and underline are used to denote the best and runner-up performance, respectively. The asterisk " * " denotes training failure.RTE [19] MRPC [24] WiC [76] CoLA [104] BoolQ [17] ACC We also evaluate different methods on RoBERTa-base, which has approximately 1 /50 the parameter count of Llama2-7B. Table 7 shows the results. Compared to MLE, MAP shows minor improvements in NLL and ECE, while MCD, ENS, and LAP enjoy significant improvements. The convergence difficulty observed with the BBB algorithm is further exacerbated on the smaller model, resulting in significant decreases in ACC across all datasets, and even training failures on RTE and WiC. In contrast, our method demonstrates the best or runner-up performance in uncertainty estimation on almost all datasets. Only a slight decrease in ACC is observed on BoolQ and CoLA. We suspect that such decrease is caused by RoBERTa-base's small model size compared to the large size of these datasets BoolQ and CoLA (i.e., underfitting). Using a larger pretrained model, e.g., Llama2-7B, would potentially address this issue.We perform an ablation study on the Llama2-7B model to showcase the effects of a range of techniques we designed: KL Re-Weighting (RW, Appendix B.1), Re-Parameterization (RP, Sec. 3.3), and Asymmetric Bayesianization (AB, Sec. 3.1). In the scenarios w/o AB, we Bayesianize both matrices, A and B. In practice, using identical initialization and prior for the standard deviation 1 . We set the number of samples during training K = 1 and the number of samples during inference N = 10 across the variants (denoted by "-") of the BBB [11] and BLoB for fair comparison. We denote by " * " experiments with the scaled standard deviation matrix. The hyphen "-" in the table denotes training failure caused by "NaN" loss. "↑" and "↓" indicate that higher and lower values are preferred, respectively. Boldface and underlining denote the best and the second-best performance, respectively. We further evaluate the inference time and maximum memory usage for standard LoRA, Laplace-LoRA (LAP), and our BLoB, as shown in Table 10 . The experiments are conducted on two NVIDIA A100 GPUs. These results show that compared to LAP, our BLoB can achieve comparable or better ECE and NLL with less inference time and less memory usage. Notably, our BLoB's memory overhead compared to standard LoRA is minimal, while LAP introduces significant memory overhead.Here we provide a more detailed empirical study on the sample size of BLoB during inference. Specifically, we report the results for different number of samples from N = 1 to N = 160 on the WG-S dataset, demonstrating improved uncertainty estimation with increased number of samples, as shown in Fig. 3 . The empirical trade-off between accuracy and calibration caused by different model architectures is observed in [85] . By controlling the standard deviation of the prior Gaussian distribution, we observed a similar trade-off between accuracy and calibration. Specifically, we report results for different prior Gaussian standard deviations, ranging from 0.05 to 0.25, while proportionally scaling the learning rate of KL divergence from its original value of 0.01 to values between 0.0025 and 0.0125. This highlights the trade-off between accuracy and calibration, as shown in Fig. 4 . We fine-tune the Llama2-7B model on the WG-S dataset for 5,000 gradient steps, evaluating the model's performance with different prior Gaussian standard deviations and learning rates of KL divergence.Estimating the uncertainty of LLMs in the embedding space has recently garnered significant attention in the community [13] . Expressing models' uncertainty via their generated embeddings can benefit both discriminative (the focus of this paper) and generative models. In this section, we present a preliminary study on uncertainty estimation in the embedding spaces of different models, as illustrated in Fig. 5 . We compare BLoB with two baseline models, BBB and MCD, which can generate embedding samples and effectively estimate uncertainty. We exclude LAP from this section due to its excessive memory consumption, which consistently results in Out-Of-Memory (OOM) errors during inference. The experiment is conducted on the OBQA dataset [65] , which consists of four categories.For each input sequence s, we use the last token's embedding generated by the final transformer block in Llama2-7B as the final embedding. Given the weights W , we denote the embedding as ϕ(s; W ). Generally, three types of embeddings can be generated using the Bayesian approach:(a) Embeddings generated by the mean of the weights (these embeddings are shown as "★" in Fig. 5 ):EQUATION(b) Embedding samples generated by sampling different weights from the approximate posterior, whose distribution is plotted by the solid line (-).(c) The expectation of the embedding, which is approximated by averaging the sampled embeddings:EQUATIONwhere N denotes the number of samples during inference, and E n denotes the n-th sampled noise for the weight matrix. We show this expectation as "▼" in Fig. 5 .To visually demonstrate the confidence calibration effect of the Bayesian treatment, we adopt the following pipeline of visualization, which we believe can be further applied in visualizing other frameworks' embedding uncertainty quality.(1) Acquire high-dimensional embeddings produced by the weight mean for the given test dataset, as decribed in (a) and Eqn. 76 above.NeurIPS Paper Checklist Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.Justification: There is no new assets introduced in this paper.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 


Decision: The paper proposes an algorithm to learn a (low-rank) variational posterior distribution over (a subset of) LLM weights during fine-tuning by combining Low-Rank Adaptation (LoRA) and Bayes by Backprop (BBB). To make the algorithm work in practice, several nontrivial modifications are introduced, such as certain parameterisations and optimisation strategies. Finally, the proposed method is empirically evaluated on a variety of in-distribution and out-of-distribution benchmarks and compared to relevant baselines.

All reviewers provided accept rating after the author effective rebuttal, but one is “borderline” and two are “weak”.
The most quoted strength of this work is the thorough experimental comparison that shows strong results compared to other LORA-based fine-tuning approach. The approach is also simpler that other Bayesian methods.

The most quoted weakness is about the contribution. 
- Reviewer Adn3: “Blob is  mean-field Variational inference on LoRA parameters”
- Reviewer Lm4v:  “Novelty is limited compared to Yang et al. [1], which implements a very similar idea via Laplace approximations instead of the Bayes-by-Backprop variational posterior.”

But the author detailed rebuttal show that their method is a substantial progress over the Laplace LoRA, being much more stable and practical. This is confirmed by reviewer AGd6: 
“While the proposed algorithm seems like a straightforward combination of LoRA and BBB, the paper mentions that certain modifications are necessary for the algorithm to work in practice. These modifications are a core part of the paper's contribution, and they are discussed thoroughly and justified with a mix of theoretical and empirical arguments.”

Reviewers strongly disagree about the presentation, with ratings ranging from 1 (poor) to 4 (excellent). The AC, who is not an expert in LoRA or Bayesian methods, found the paper easy to read and understand.

Reviews: review: paper_summary: paper_summary: This work introduces a Bayesian Deep Learning framework for finetuning LLMs with probabilistic LoRA. Unlike existing work by Yang et al. [1] which uses a Laplace approximation, the variational distribution is parameterized using diagonal Gaussians as done in [2] (”Bayes by Backprop”). 

Specifically, the authors design the variational distribution as follows:
Given the standard LoRA updates $W = W_0 + BA$, only $A$ is modeled probabilistically while $B$ is deterministic. In combination with an appropriate choice of prior, this allows for a a low-rank parameterization of the variational posterior and fast KL-computation in the lower dimensional parameter space.
Experiments compare the methods to existing baselines in terms of accuracy of the mean predictions, calibration (ECE) and NLL on a range of datasets. 

[1] https://arxiv.org/pdf/2308.13111
[2] https://proceedings.mlr.press/v37/blundell15.pdf

____

I have updated my overall score +1 after the rebuttal.
strength_weakness: strength_weakness: Novelty is limited compared to Yang et al. [1], which implements a very similar idea via Laplace approximations instead of the Bayes-by-Backprop variational posterior. Additionally, the method by Yang et al. ("LAP") is not replicated successfully in Table 1. For example, for the LAP method on ARC-C datasets the results in this work fall behind the results in [1] by 50%: 
- This work: Accuracy: 29.73±12.02, ECE: 14.24±1.65 and NLL: 1.53±0.01
- [1]: Accuracy: 66.91±1.1, ECE:  7.5±1.2 and NLL: 0.86±0.02, respectively. 

This makes it difficult to gauge how much the approach improves compared to [1].
questions: questions: Questions:
- Table 1: As $N \rightarrow \infty$, would we expect the accuracy to converge to those obtained with $N=0$ (mean prediction)? Is BLoB an unbiased estimator of the mean prediction?
- l. 240: “For MCD, we use an ensemble of 10 LoRAs with a dropout rate of p = 0.1”. Why is an ensemble used here? Shouldn’t the uncertainties come from the dropout variational posterior? Why are multiple models necessary? 
- Will you open source the code on publication?

Suggestions & Typos: 
- l. 33: I.e., unable to … (sentence is not grammatical, delete) 
- l. 110: yielded by parameterization —> yielded by reparameterization  
- Algorithm 1, line 11: Updater —> Update 
- l. 237: Maximize A Posteriori —> Maximum A Posteriori  
- l. 239: “For MAP, we use a weight decay rate of 1e-5.” I suggest to formulate this in terms of prior variance instead. 
- l. 638 covaraince —> covariance
limitations: limitations: No, I recommend a more explicit treatment of the limitations, perhaps in a dedicated section.
paper_summary: paper_summary: Proposes to have a Bayesian version of LoRA by placing priors over the low rank matrices. This evolves into placing prior over one of the low rank matrices. Inference is via variational inference. Results are mixed, though it shows the method is competitive with others.
strength_weakness: strength_weakness: The flow of the paper is rather poor, and the presentation of the method overly complicated. I prefer the authors to simply start from (5) and then say the deltas to the pre-trained weights are Gaussians because they are sum of Gaussians. The consideration to also put prior on $B$ can be deferred to discussions, and the trade-offs discussed there. This make the model presentation cleaner, and avoids the problem of having to involve improper priors in Theorem 3.1 (since B is low rank, $\Sigma_{q}$ is low rank, the Gaussian is degenerate) and Theorem 3.2, and the subsequent problem of reader having to think about how the KL with generate Gaussians would work out in theory. In short, the theorems are unnecessary and need fixing because of improper distributions.

1. Line 137: What is it that requires accurate estimation?
2. Section 4: How is more than one sample (for N>0) used in generating a single output for evaluation and then influence the evaluation metric?

*Minor*
- Eq (3) to remove the $\min_{\theta}$ prefix, and to introduce the prior and likelihood.
- Section 4: Shouldn't ECE and NLL be the metrics in focus because this is about "Bayesian"?
- Section 4: The numbers in the table does not justify saying BLoB is the best.
- "Bayesianization" is a mouthful --- can we think of a better derived word?
questions: questions: Q1 and Q2 above (in weaknesses)
limitations: limitations: No. Paper need to include some discussion on the limitations of their approach/paper. I do not see any potential negative societal impact to be addressed.
paper_summary: paper_summary: This paper proposes a method called Blob, which uses mean-field variational inference on LoRA parameters to obtain Bayesian estimation on LLM. This will result in a richer posterior structure on the weight than diagonal, while maintaining the same computational cost as diagonal covariance. During training, a new parameterization is used for posterior covariance to accelerate convergence and Flipout is used to make training more stable. Empirical results show Blob obtains competitive performance compared with other Bayesian methods.
strength_weakness: strength_weakness: - Not clear about the method's limitations. Since it uses VI, common disadvantages of VI also apply: (1) it needs to do sampling when making predictions, which can be computationally expensive given the size of LLM. (2) Given the complex posterior structure, there will be a very high chance of getting noise samples, which leads to suboptimal performance. This is already demonstrated in Table 1 that using MC estimation decreases the ACC compared with only using the mean. In additional, it seems in practice many extra efforts are needed to get the method to work, e.g., the complex KL re-weighting and rescale size of training data mentioned in Appendix B.1.

- The whole paper feels a bit overselling: (1) I don't see the need to make Theorem 3.1 and Theorem 3.2 theorems, nothing surprising is stated there. (2) In Line 318, the authors claim their method "demonstrates superior generalization and uncertainty estimation capabilities", I don't think the results justify the "superior" claim. (3) In essence, Blob is mean-field VI on LoRA parameters trained with better strategies to make it suffer less from VI's high variance. I have no problem at explaining methods in great detail, but the whole writing on page 4 and page 5 feels like the authors are trying to make it sound as fancy as possible. I like the idea of the paper, but I think it can be presented in a more honest way.
questions: questions: Major:
- I'm not sure what it mean for the mean matrix of q(A) to be an output of a neural network (Line 186-187). Does this mean you are using another neural network to fit the mean of posterior? 

- In Eq. (70), it seems different mini-batch has different KL re-weighting hyperparamers lambda? How sensitive is the training to the value of lambda?

- What is the difference between BBB and Blob? From my understanding, the only difference is how you do the training. For BBB you use original parameterization and for Blob you use new parameterization plus Flipout? Also how many samples did you use to evaluate BBB in Table 1?

- I would be interested to see a figure plot of model performance versus the number of MC samples. Specifically, the y-axis is accuracy/ECE/NLL, the x-axis is the number of samples used during inference. This helps the reader to understand the sufficient sampling size for inference.

- I would be interested to see Laplace added to Table 8 to show a fair comparison of the computation cost compared with Blob.

- Why is Laplace performing so poorly on ARC-C dataset? The original Bayesian LoRA paper it seems to perform much better on that dataset. 

Minor:
- Line 236, MLE and MAP are not uncertainty estimation methods.
- I don't entirely agree with the claim made in Lines 46-47. The possible suboptimal estimation of Laplace approximation comes from (1) LA, which assumes we reach MAP estimation while in practice you never know if the complex NN has really converged to MAP; (2) setting prior precision. It doesn't come from the post-training procedures.
- the notation "n" refers to two different things, at Eq. (4) it's the number of parameters, and at Line 153 it's the dimension of the weight matrix.
- Bolding is missing on ACC in Table 7.
limitations: limitations: See weaknesses and questions above.
paper_summary: paper_summary: The paper proposes an algorithm to learn a (low-rank) variational posterior distribution over (a subset of) LLM weights during fine-tuning by combining Low-Rank Adaptation (LoRA) and Bayes by Backprop (BBB). To make the algorithm work in practice, several nontrivial modifications are introduced, such as certain parameterisations and optimisation strategies. Finally, the proposed method is empirically evaluated on a variety of in-distribution and out-of-distribution benchmarks and compared to relevant baselines.
strength_weakness: strength_weakness: - The proposed method requires additional compute time and memory compared to standard LoRA. While this information is provided in Appendix C.2 (and seems to be small), I encourage the authors to be transparent about this by briefly mentioning it in the main paper (if that is not already the case, perhaps I missed it).
questions: questions: - Appendix B.1 mentions weights for the KL divergence term. How important is this? Please correct me if I am wrong, but I believe this was not explicitly discussed in the main paper. If this is important for the overall performance of the algorithm (and typically these kind of adjustments tend to be...), I suggest to at least mention it in the main paper (even if details remain in the appendix).
- Is there any particular reason why $\mathbf{A}$ is modelled in a probabilistic way instead of $\mathbf{B}$? In other words, would it also be possible to do the reverse, i.e. assume a prior over $\mathbf{B}$ and treat $\mathbf{A}$ deterministically? The choice seems arbitrary to me, but perhaps there is an underlying reason to it.
limitations: limitations: Theoretical assumptions are stated as part of the provided theorems. As far as I am aware, general limitations of the proposed method (beyond the ones which are being addresses by the proposed method itself) are not explicitly discussed and could at least be briefly mentioned in the conclusion.

Comments: title: Thank You for Your Further Feedback
Thank You for Your Further Feedback
Response to rebuttal
Improved Outline of the Method
Thank you
Thank You for Your Further Feedback
Degenerate Gaussians are Proper for Probabilistic Inference
Remaining Questions for Reviewer Adn3
comment: Thank you for your further feedback and for maintaining an open line of communication. We are glad that you found our response helpful and that it addressed your concerns. Below we address your follow-up comment on Q1.

**[Higher Accuracy of Mean Prediction Compared to Bayesian Prediction]**

This is a good question. By abandoning the modeling of the posterior distribution, mean prediction sacrifices some degree of calibration in exchange for improved accuracy. This empirical trade-off between accuracy and calibration has been noted in [3]. We will include this discussion in our revision.

Lastly, we would like to express our gratitude once again for your insightful and constructive comments. They have significantly contributed to the improvement of our paper. 

[3] Stengel-Eskin, Elias, and Benjamin Van Durme. "Calibrated interpretation: Confidence estimation in semantic parsing." Transactions of the Association for Computational Linguistics 2023
Thank you for your continued feedback and for keeping the communication channel open. We are glad that you found our response helpful and that it addressed most of your concerns. Below we address your remaining comments in detail. 

**W3.1 [the Term "Theorem"]**

We appreciate your insightful feedback. In light of your suggestion, we will revise Theorem 3.1 and Theorem 3.2 to Proposition 3.1 and Proposition 3.2, respectively. We believe this change can more accurately reflect the nature of our contributions. Nonetheless, we maintain that the discussion presented here offers valuable insights to the broader community. The theoretical foundation of Asymmetric Bayesianization laid out in this work can serve as a springboard for the development of new methods in the field.

In addition, inspired by your comments, we plan to reorganize Section 3.1 to enhance clarity. We will start with the calculation of the full weight matrix $W\_{ij}=W\_{0,ij} + \sum\_{k=1}^r B\_{ik}A\_{kj}$ (i.e., Eqn. 5 of the paper) and subsequently introduce the advantages of Asymmetric Bayesianization over bayesianizing both $A$ and $B$.

**[Related Work]**

Thank you for bringing the ICML 2024 paper "Variational Learning is Effective for Large Deep Networks" to our attention. It was not available at the time of our submission. We agree that it is highly relevant to our work, and we will ensure it is cited and discussed in the revised version.

Finally, we would like to express our sincere gratitude once again for your insightful and constructive comments, which have greatly helped improve our paper.
Thanks for the additional details, especially the derivation and experiment in Q1. I see that 

$$E_{q(W|\theta)}[P(Y | W, X)]  \neq P(Y | E_{q(W|\theta)} [W], X).$$

Follow up: What's the intuition behind the mean prediction (RHS) achieving higher accuracy than the (approximate) Bayesian prediction (LHS), also for large $N$? 

Overall, the rebuttal partly addresses my concerns around novelty (Q2 and Q3 in "General response"), hence I increase my overall score by one point.
Thank you for the rebuttal, it addresses most of my concerns. Still, I am not fully convinced that a Theorem is needed to state results in Theorem 3.1 and Theorem 3.2, maybe this is because I come from a more Bayesian background and thus the results do not look surprising.

On a side note, I recently came across an ICML 2024 paper "Variational Learning is Effective for Large Deep Networks". In essence, they do natural gradient descent in natural parameter space for mean-field VI, and they make it work for large-scale neural networks including GPT2. This is highly relevant to the paper so I strongly suggest the author to discuss it in the paper.
Thank you for your valuable suggestions on enhancing the presentation quality of our work. Following your suggestions, we have refined Sections 3.1 and 3.2. Below is the outline of our updated version: 

**Section 3.1: Low-Rank Variational Approximate Posterior Distribution**

1. Main Method: Introduction of Asymmetric Bayesianization Scheme
   - Calculating each entry of the full weight matrix: $W\_{ij}=W\_{0,ij} + \sum\_{k=1}^r B\_{ik}A\_{kj}$ (i.e., Eqn. 5 of the paper)
   - $A$ modeled by independent Gaussians: $P(A|\theta)=\mathcal{N}(A|M,\Omega^2)$
   - $B$ modeled as deterministic values

2. Presentation of Theorem 3.1
   - Statement: Asymmetric Bayesianization corresponds to assuming a low-rank Gaussian variational distribution for the full weight $W$

3. Discussion on the Choice of Asymmetric Bayesianization
   - Advantages over Bayesianizing both $A$ and $B$:
     + Stable training at the early stage
     + Faster convergence of parameters during training
     + Lower memory cost

**Section 3.2: Low-Rank Prior Distribution**

1. Presentation of the Prior Distribution on the Low-Rank Component $A$

2. Corresponding Prior Distribution on the High-Dimensional Space of the Full Weight Matrix $W$

3. Presentation of Theorem 3.2
   - Statement: Full-weight KL divergence can be efficiently computed on the low-rank component

4. Dedicated Remark on the Legitimacy of using Degenerate Gaussians for Probabilistic Inference (as in Our Previous Discussion)

Thank you once again for your suggestions. We believe this revision will make the paper more accessible for our audience and help address any potential areas of confusion.
Thank you very much for the encouraging comments and for acknowledging our contribution. We are glad that our response has been helpful and convincing. We will include the discussion and ablation study in the revision based on your suggestion.
Thank you for your further feedback and for keeping the communication channel open. We are glad that you found our response helpful and that it addressed your concerns. Below we address your remaining comments in detail. 

**W1+W2 [Presentation]**: Thank you again for your suggestion regarding the presentation. We will re-organize the presentation to start with Eqn. 5 in the revision as you suggested.

**W5 ["Generally Better"]**: Thank you for your suggestion. We will adjust the claim from "the best" to "generally better" to more accurately reflect our contribution.

**Q2 [ACC and N=10 Samples]**: Thank you for your recognition and insightful question. For ACC, we first draw N=10 samples, i.e., N=10 output probabilities, from our BLoB, and then compute the average of these N samples as *one* single final prediction. 

This approach aligns with the method we previously described: "approximate this expectation through sampling: $E_{q(W|\theta)}[P(Y|W,X)] \approx \frac{1}{N} \sum_{n=1}^{N} P(Y|W_n, X)$." Based on this averaged output probability (i.e., the approximation of the expectation), we then calculate the ACC. The ACC is *not* calculated based on whether *any* individual sample is correct. Instead, N=10 samples are aggregated to make *only one prediction*. Therefore there is *no discount* needed. 

We will follow your suggestion to provide further clarification on the output of variational inference and clarify that "experimental setup strictly adheres to the original LoRA [2] framework" in our revised version.

Last but not least, we would like to thank you again for your insightful and constructive comments. They have greatly help improve our paper.
Thank you for your response and for providing concrete numbers for inference time and memory consumption. I encourage you to include those in your manuscript. I am also glad to hear that you will present an ablation study about your KL reweighing schedule, as this is "crucial to BLoB's final performance". I am maintaining my score as I believe it is appropriate.
W1 + W2. I understand that the method is sound because your degenerate Gaussians are on the same subspace. It is now clear to me that you also understand this constraint. I'll increase the soundness score, and the overall score.

However, I still believe that the presentation can be vastly simplified (and hence the paper more accessible) by simply starting from (5). The relation to full Bayesianization can be deferred to discussion.

W5.  It is "generally better" than the other methods. You can also say it is "the best on average" if you do some averaging either over the score or by counting wins. You cannot just say it is "the best".

Q2. I see that the "experimental setup strictly adheres to the original LoRA [2] framework". It is good to say this explicitly in the paper. For ACC, do you generate, say N=10, and then count if *any* is correct? Is there some discount because you have N chances?
In this separate Official Comment, we will elaborate on two key reasons why degenerate Gaussian distributions are valid for probabilistic inference:

1. **Their probability density is well-defined.**

According to Schoeman et al. [1], the *probability density* of a degenerate Gaussian distribution can be factorized into the product of the density of a non-degenerate Gaussian distribution in a lower-dimensional linear subspace and a Dirac delta function (Eq. 12 in [1]). This density function has an alternative expression in the form of limit (Eq. 16-18 in [1]), which aligns with our main theorem 3.1 and 3.2:
$$
p(x) = \lim_{a\rightarrow 0}\mathcal{N}(x | \mu, Q(\Lambda^{-1} - aI)Q^T + a I) = \mathcal{N}(x | \mu, \Sigma=\lim_{a\rightarrow 0} Q(\Lambda^{-1} - aI)Q^T + a I). 
$$
Here $\Lambda^{-1}$ is the precision matrix of the lower-dimensional non-degenerate Gaussian distribution, and $Q$ is the linear transformation matrix that expands this to high dimensions. In our main paper, we use the simplified notation $\mathcal{N}(x | \mu, \Sigma)$ for readability, omitting the limit. The full limit notation and proof appear in Appendix A.1. We will clarify this simplification in our revision.

2. **The KL divergence between two degenerate Gaussians is computable under certain conditions.**

Secondly, the *KL divergence* between two degenerate Gaussians can be computed under the condition that they "have support on the same lower-dimensional manifold" (Eq. 44 in Schoeman et al. [1]). This condition has also been clearly stated as satisfied in our Theorem 3.2 ($RR^{T} = BB^{T}$, line 177), and we also include a detailed derivation of how we reach this condition in Appendix A.1 (line 676-678).

[1] Schoeman, Johannes Cornelius, Corné E. van Daalen, and Johan A. du Preez. "Degenerate Gaussian factors for probabilistic inference." IJAR 2022.
**W1. Not clear about the method's limitations.**

Please refer to **Q1. [Discussion of Limitations]** in **General Response**.


**Q1. What it mean for the mean matrix of q(A) to be an output of a neural network (Line 186-187)?**

We are sorry for the typo. $q(A)$ is not the output of a neural network. It is directly modeled as a Gaussian distribution $\mathcal{N}(M, \Omega)$, with $M$ and $\Omega$ as learnable parameters. 

**Q2. In Eq. (70), it seems different mini-batch has different KL re-weighting hyperparamers lambda? How sensitive is the training to the value of lambda?**

Yes, different mini-batches has different $\lambda$ values. When we use the typical settings of $\lambda_i = \frac{2^{M-i}}{2^M-1}$ [1] or $\lambda_i = \frac{1}{M}$ [2], the model struggles to fit the data distribution while converging to the prior distribution, resulting in an "NaN" NLL loss. Our setting of $\lambda_i = \frac{2^i}{2^M-1}$ allows the model to fit to the data in the early stages of training and achieve a reasonable level of complexity cost in the later stages. This re-balance scheme of data likelihood and complexity cost has proven effective across multiple datasets and can be potentially applied to other VI methods. 

We will include the discussion above in the revision as suggested. 

[1] Blundell, Charles, et al. "Weight uncertainty in neural network." International conference on machine learning. PMLR 2015.

[2] Graves, Alex. "Practical variational inference for neural networks." NeurIPS 2011.


**Q4. I would be interested to see a figure plot of model performance versus the number of MC samples. ...**

This is a good suggestion. Following your suggestion, we have run additional experiments for model performance versus the number of MC samples on WG-S dataset. The results are presented in the table below. The corresponding figure is provided as **Figure 1** in the PDF attached to the **General Response**.

|                        | **N=0** | **N=1** | **N=2** | **N=3** | **N=4** | **N=5** | **N=10** | **N=20** | **N=40** | **N=80** | **N=160**  |
| ---------------------- | ------- | ------- | ------- | ------- | ------- | ------- | -------- | -------- | -------- | -------- | ------ |
| **ACC ($\uparrow$)**   | 71.44   | 65.20   | 66.28   | 66.91   | 67.63   | 68.20   | 68.31    | 68.04    | 68.15    | 68.31    | 68.15 |
| **ECE ($\downarrow$)** | 19.71   | 19.72   | 14.60   | 12.82   | 11.27   | 10.58    | 9.51    | 9.47    | 9.25     | 8.75     | 8.75   |
| **NLL ($\downarrow$)** | 0.84    | 0.8617  | 0.7355  | 0.6971  | 0.6766  | 0.6620  | 0.6395   | 0.6313   | 0.6233   | 0.6213   | 0.6215 |

We can see that in general, the performance improves with more MC samples, but saturates after around $N=10$. 

**Q6. Why is Laplace performing so poorly on ARC-C dataset? The original Bayesian LoRA paper it seems to perform much better on that dataset.**

Please refer to **Q2. [Reproducing LAP's Results on the ARC-C Dataset]** in **General Response**. 

**Q7. Line 236, MLE and MAP are not uncertainty estimation methods.**

We apologize for the confusion. We will remove "uncertainty estimation" when describing MLE and MAP in the revision as suggested. 

**Q8. I don't entirely agree with the claim made in Lines 46-47. The possible suboptimal estimation of Laplace... never know if the complex NN has really converged to MAP; (2) setting prior precision...**

This is a good point. We agree that these are also major reasons that lead to Laplace approximation's poor performance. We will include these two additional reasons in our revision as suggested. 

**Q9. Notation issue and missing Bolding**

We are sorry for these typos and will correct them in the revision as suggested.