Description of the variational inference framework<\/p>\n
Exact computation of the truncated posterior (Eq. (4<\/a>)) becomes intractable as the number of parameters in the network grows, necessitating approximation methods. Here, we use variational inference (VI), which tackles this challenge by restricting the posterior to a parameterized distribution, known as the variational distribution. We choose a multivariate Gaussian \\({q}_{{{\\boldsymbol{\\theta }}}_{t}}({\\boldsymbol{\\omega }})\\).<\/p>\n To reduce complexity, we adopt the mean-field approximation40<\/a>, modeling each synaptic weight as an independent Gaussian:<\/p>\n $${q}_{{{\\boldsymbol{\\theta }}}_{t}}({\\boldsymbol{\\omega }})=\\mathop{\\prod }\\limits_{i=1}^{s}{q}_{{\\theta }_{t,i}}({\\omega }_{i})\\Rightarrow {q}_{{{\\boldsymbol{\\theta }}}_{t}}({\\boldsymbol{\\omega }})={\\mathcal{N}}({\\boldsymbol{\\omega }};{{\\boldsymbol{\\mu }}}_{t},{\\rm{diag}}({{\\boldsymbol{\\sigma }}}_{t}^{2})),$$<\/p>\n where s is the total number of synapses, and each \u03c9i is modeled by \\({q}_{{\\theta }_{t,i}}({\\omega }_{i})={\\mathcal{N}}({\\omega }_{i};{\\mu }_{i},{\\sigma }_{i}^{2})\\). Our goal is for \\({q}_{{{\\boldsymbol{\\theta }}}_{t}}({\\boldsymbol{\\omega }})\\) to approximate the true truncated posterior \\(p({\\boldsymbol{\\omega }}\\,| \\,{{\\mathcal{D}}}_{t-N},\\ldots,{{\\mathcal{D}}}_{t})\\).<\/p>\n Following standard VI practice9<\/a>,18<\/a>,41<\/a>,42<\/a>,43<\/a>, we minimize the Kullback\u2013Leibler (KL) divergence between \\({q}_{{{\\boldsymbol{\\theta }}}_{t}}({\\boldsymbol{\\omega }})\\) and the target posterior:<\/p>\n $${\\boldsymbol{\\theta}}_t \\, =\\, \\arg\\min_{{\\boldsymbol{\\theta}}} \\; D_{KL}\\!\\left[ q_{{\\boldsymbol{\\theta}}}({\\boldsymbol{\\omega}}) \\;\\|\\; p\\bigl({\\boldsymbol{\\omega}}\\!\\mid\\!{\\mathcal{D}}_{t-N},\\ldots,{\\mathcal{D}}_{t}\\bigr) \\right].$$<\/p>\n \n (20)\n <\/p>\n Variational free energy<\/p>\n Starting from Eq. (20<\/a>), we rewrite the target posterior \\(p({\\boldsymbol{\\omega }}\\,| \\,{{\\mathcal{D}}}_{t-N},\\ldots,{{\\mathcal{D}}}_{t})\\) using Bayes\u2019 rule:<\/p>\n $$p({\\boldsymbol{\\omega }}\\,| \\,{{\\mathcal{D}}}_{t-N},\\ldots,{{\\mathcal{D}}}_{t})=\\frac{p({{\\mathcal{D}}}_{t}\\,| \\,{\\boldsymbol{\\omega }})\\,p({\\boldsymbol{\\omega }}\\,| \\,{{\\mathcal{D}}}_{t-N-1},\\ldots,{{\\mathcal{D}}}_{t-1})}{p({{\\mathcal{D}}}_{t})}\\times \\frac{p({{\\mathcal{D}}}_{t-N-1})}{p({{\\mathcal{D}}}_{t-N-1}\\,| \\,{\\boldsymbol{\\omega }})}.$$<\/p>\n Applying the KL divergence definition \\(\\,{D}_{KL}[Q\\,| \\,| P]\\,=\\,{{\\mathbb{E}}}_{Q}[\\log (Q\/P)]\\) and ignoring terms independent of \u03b8, we obtain:<\/p>\n $${{\\boldsymbol{\\theta }}}_{t}=\t\\arg \\mathop{\\min }\\limits_{{\\boldsymbol{\\theta }}}{{\\mathbb{E}}}_{{q}_{{\\boldsymbol{\\theta }}}({\\boldsymbol{\\omega }})}\\,\\left[\\log (\\frac{{q}_{{\\boldsymbol{\\theta }}}({\\boldsymbol{\\omega }})}{p({\\boldsymbol{\\omega }}\\,| \\,{{\\mathcal{D}}}_{t-N-1},\\ldots,{{\\mathcal{D}}}_{t-1})})\\right.\\\\ \t\\left.-\\log p({{\\mathcal{D}}}_{t}\\,| \\,{\\boldsymbol{\\omega }})+ \\log p({{\\mathcal{D}}}_{t-N-1}\\,| \\,{\\boldsymbol{\\omega }})\\right].$$<\/p>\n \n (21)\n <\/p>\n As explained in the main text, we further approximate \\(\\,p({\\boldsymbol{\\omega }}\\,| \\,{{\\mathcal{D}}}_{t-N-1},\\ldots,{{\\mathcal{D}}}_{t-1})\\approx {q}_{{{\\boldsymbol{\\theta }}}_{t-1}}({\\boldsymbol{\\omega }})\\), which yields<\/p>\n $${{\\boldsymbol{\\theta }}}_{t}=\t\\arg \\mathop{\\min }\\limits_{{\\boldsymbol{\\theta }}}\\left[{D}_{KL}({q}_{{\\boldsymbol{\\theta }}}({\\boldsymbol{\\omega }})\\,\\parallel \\,{q}_{{{\\boldsymbol{\\theta }}}_{t-1}}({\\boldsymbol{\\omega }}))-{{\\mathbb{E}}}_{{q}_{{\\boldsymbol{\\theta }}}({\\boldsymbol{\\omega }})}\\,[\\log p({{\\mathcal{D}}}_{t}\\,| \\,{\\boldsymbol{\\omega }})]\\right.\\\\ \t \\left.+ {{\\mathbb{E}}}_{{q}_{{\\boldsymbol{\\theta }}}({\\boldsymbol{\\omega }})}\\,[\\log p({{\\mathcal{D}}}_{t-N-1}\\,| \\,{\\boldsymbol{\\omega }})]\\right].$$<\/p>\n \n (22)\n <\/p>\n Defining the function \\({{\\mathcal{F}}}_{t}\\) (the variational free energy42<\/a>,43<\/a>) as:<\/p>\n $${\\mathcal{F}}_t \\,=\t\\, \\underbrace{ D_{KL}\\bigl[q_{{\\boldsymbol{\\theta}}_t}({\\boldsymbol{\\omega}})\\,\\|\\, q_{{\\boldsymbol{\\theta}}_{t-1}}({\\boldsymbol{\\omega}})\\bigr] \\,-\\, {\\mathbb{E}}_{q_{{\\boldsymbol{\\theta}}}({\\boldsymbol{\\omega}})}\\!\\bigl[\\log p({\\mathcal{D}}_t\\,\\mid\\,{\\boldsymbol{\\omega}})\\bigr] }_{{\\text{Learning}}} \\,\\\\ \t+\\, \\underbrace{ {\\mathbb{E}}_{q_{{\\boldsymbol{\\theta}}}({\\boldsymbol{\\omega}})}\\!\\bigl[\\log p({\\mathcal{D}}_{t-N-1}\\,\\mid\\,{\\boldsymbol{\\omega}})\\bigr] }_{{\\text{Forgetting}}},$$<\/p>\n We recover the expression from the main text (Eq. (5<\/a>)), which explicitly separates the learning term (adapting to the current task) from the forgetting term (downweighting information from older tasks).<\/p>\n Forgetting term approximation<\/p>\n From Eq. (5<\/a>) in the main text, the main challenge lies in computing the \u201cforgetting\u201d term \\(\\,{{\\mathbb{E}}}_{{q}_{{\\boldsymbol{\\theta }}}({\\boldsymbol{\\omega }})}\\,[\\log p({{\\mathcal{D}}}_{t-N-1}| {\\boldsymbol{\\omega }})]\\). For this purpose, we assume each dataset \\({{\\mathcal{D}}}_{i}\\) has equal marginal likelihood, so<\/p>\n $$p({{\\mathcal{D}}}_{t-N-1},\\ldots,{{\\mathcal{D}}}_{t-1}| {\\boldsymbol{\\omega }})=\\mathop{\\prod }\\limits_{i=t-N-1}^{t-1}p({{\\mathcal{D}}}_{i}| {\\boldsymbol{\\omega }})={[p({{\\mathcal{D}}}_{t-N-1}| {\\boldsymbol{\\omega }})]}^{N},$$<\/p>\n thus weighting each dataset in the memory window equally. Because our variational and prior distributions are Gaussians, this assumption leads to a closed-form expression for the forgetting likelihood.<\/p>\n Lemma 1: Let \\({q}_{{\\boldsymbol{\\theta }}}({\\boldsymbol{\\omega }})\\approx p({\\boldsymbol{\\omega }}| {\\mathcal{D}})\\) be a mean-field Gaussian for a BNN, where \u03b8\u2009=\u2009(\u03bc,\u00a0\u03c3) and \u03c9\u2009=\u2009\u03bc\u2009+\u2009\u03f5 \u22c5 \u03c3, \\({\\boldsymbol{\\epsilon }} \\sim {\\mathcal{N}}(\\overrightarrow{0},{{\\bf{I}}}_{s})\\). If the prior is \\(p({\\boldsymbol{\\omega }})={\\mathcal{N}}({\\boldsymbol{\\omega }};{{\\boldsymbol{\\mu }}}_{{\\rm{prior}}},{\\rm{diag}}({{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}))\\), then \\(p({\\mathcal{D}}| {\\boldsymbol{\\omega }})\\propto {q}_{L}({\\boldsymbol{\\omega }})={\\mathcal{N}}({\\boldsymbol{\\omega }};{{\\boldsymbol{\\mu }}}_{L},{\\rm{diag}}({{\\boldsymbol{\\sigma }}}_{L}^{2})),\\) and the negative log-likelihood takes a quadratic form: \\({\\mathcal{L}}({\\boldsymbol{\\omega }})=\\frac{{({\\boldsymbol{\\omega }}-{{\\boldsymbol{\\mu }}}_{L})}^{2}}{2\\,{{\\boldsymbol{\\sigma }}}_{L}^{2}}+ \\frac{1}{2}\\,\\log (2\\pi \\,{{\\boldsymbol{\\sigma }}}_{L}^{2}),\\) where \\(\\frac{1}{{{\\boldsymbol{\\sigma }}}^{2}}=\\frac{1}{{{\\boldsymbol{\\sigma }}}_{L}^{2}}+\\frac{1}{{{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}}\\) and \\({\\boldsymbol{\\mu }}=\\frac{{{\\boldsymbol{\\mu }}}_{L}\\,{{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}+{{\\boldsymbol{\\mu }}}_{{\\rm{prior}}}\\,{{\\boldsymbol{\\sigma }}}_{L}^{2}}{{{\\boldsymbol{\\sigma }}}_{L}^{2}+{{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}}.\\)<\/p>\n A full proof of Lemma 1 appears in Supplementary Note\u00a02<\/a>. It uses Bayes\u2019 rule for Gaussian posteriors and priors to link \u03bcL and \u03c3L to \u03bc and \u03c3. In our truncated-posterior setting, \\(p({{\\mathcal{D}}}_{t-N-1}| {\\boldsymbol{\\omega }})\\propto {q}_{{L}_{t-1}}{({\\boldsymbol{\\omega }})}^{\\frac{1}{N}},\\) yielding<\/p>\n $${\\mathcal{F}}_t \\,=\\, \\underbrace{\\left[ D_{KL}\\bigl(q_{{\\boldsymbol{\\theta}}_t}\\|\\,q_{{\\boldsymbol{\\theta}}_{t-1}}\\bigr) \\,+\\, {\\mathcal{C}}_t \\right]}_{{\\text{Learning}}} \\,+\\, \\underbrace{\\frac{1}{N}\\left[ -\\,\\frac{\\bigl({\\boldsymbol{\\mu}}_t – {\\boldsymbol{\\mu}}_{L_{t-1}}\\bigr)^2}{2\\,{\\boldsymbol{\\sigma}}_{L_{t-1}}^2} \\,-\\, \\frac{{\\boldsymbol{\\sigma}}_t^2}{{\\boldsymbol{\\sigma}}_{L_{t-1}}^2} \\right]}_{{\\text{Forgetting}}},$$<\/p>\n \n (23)\n <\/p>\n where \\({{\\mathcal{C}}}_{t}=-\\,{{\\mathbb{E}}}_{{q}_{{{\\boldsymbol{\\theta }}}_{t}}({\\boldsymbol{\\omega }})}[\\log p({{\\mathcal{D}}}_{t}| {\\boldsymbol{\\omega }})].\\) the forgetting term thus \u201cde-consolidates\u201d each synapse, nudging it toward the prior to free capacity for new tasks.<\/p>\n Update rules<\/p>\n To approximate the true posterior distribution at time t, we seek to minimize the free energy \\({{\\mathcal{F}}}_{t}\\). The Kullback\u2013Leibler term in \\({{\\mathcal{F}}}_{t}\\) between two diagonal Gaussians has a known closed form:<\/p>\n $${D}_{KL}[{q}_{{{\\boldsymbol{\\theta }}}_{t}}({\\boldsymbol{\\omega }})\\parallel {q}_{{{\\boldsymbol{\\theta }}}_{t-1}}({\\boldsymbol{\\omega }})]=\\mathop{\\sum }\\limits_{i=1}^{s}\\log \\left(\\frac{{\\sigma }_{i,t-1}}{{\\sigma }_{i,t}}\\right)+ \\frac{{\\sigma }_{i,t}^{2}+{({\\mu }_{i,t-1}-{\\mu }_{i,t})}^{2}}{2\\,{\\sigma }_{i,t-1}^{2}}-\\frac{1}{2}.$$<\/p>\n \n (24)\n <\/p>\n Taking derivatives of \\({{\\mathcal{F}}}_{t}\\) with regards to \u03bct and \u03c3t yields the implicit equations:<\/p>\n $$\\frac{\\partial {{\\mathcal{F}}}_{t}}{\\partial {{\\boldsymbol{\\mu }}}_{t}} = \\frac{\\Delta {\\boldsymbol{\\mu }}}{{{\\boldsymbol{\\sigma }}}_{t-1}^{2}}+ \\frac{\\partial {{\\mathcal{C}}}_{t}}{\\partial {{\\boldsymbol{\\mu }}}_{t}}-\\frac{{{\\boldsymbol{\\mu }}}_{t}-{{\\boldsymbol{\\mu }}}_{{L}_{t-1}}}{N\\,{{\\boldsymbol{\\sigma }}}_{{L}_{t-1}}^{2}}= \\overrightarrow{0},$$<\/p>\n \n (25)\n <\/p>\n $$\\frac{\\partial {{\\mathcal{F}}}_{t}}{\\partial {{\\boldsymbol{\\sigma }}}_{t}}=-\\,\\frac{1}{{{\\boldsymbol{\\sigma }}}_{t-1}+\\Delta {\\boldsymbol{\\sigma }}}+ \\frac{{{\\boldsymbol{\\sigma }}}_{t-1}+\\Delta {\\boldsymbol{\\sigma }}}{{{\\boldsymbol{\\sigma }}}_{t-1}^{2}}+ \\frac{\\partial {{\\mathcal{C}}}_{t}}{\\partial {{\\boldsymbol{\\sigma }}}_{t}}-\\frac{2\\,({{\\boldsymbol{\\sigma }}}_{t-1}+\\Delta {\\boldsymbol{\\sigma }})}{N\\,{{\\boldsymbol{\\sigma }}}_{{L}_{t-1}}^{2}}=\\overrightarrow{0},$$<\/p>\n \n (26)\n <\/p>\n where \u0394\u03bc\u2009=\u2009\u03bct\u00a0\u2212\u00a0\u03bct\u22121 and \u0394\u03c3\u2009=\u2009\u03c3t\u00a0\u2212\u00a0\u03c3t\u22121. Solving these leads to the MESU update rule:<\/p>\n Theorem 1: Consider a stream of data \\({\\{{{\\mathcal{D}}}_{i}\\}}_{i=0}^{t}\\). Let \\({q}_{{{\\boldsymbol{\\theta }}}_{t}}({\\boldsymbol{\\omega }})\\) be a mean-field Gaussian for a BNN at time t, with \u03b8t\u2009=\u2009(\u03bct,\u00a0\u03c3t) and samples \\({\\boldsymbol{\\omega }}={{\\boldsymbol{\\mu }}}_{t}+{\\boldsymbol{\\epsilon }}\\cdot {{\\boldsymbol{\\sigma }}}_{t},{\\boldsymbol{\\epsilon }}\\, \\sim \\,{\\mathcal{N}}(\\overrightarrow{0},{\\bf{I}})\\). Suppose \\({q}_{{{\\boldsymbol{\\theta }}}_{t-1}}({\\boldsymbol{\\omega }})\\approx p({\\boldsymbol{\\omega }}\\,| \\,{{\\mathcal{D}}}_{t-N-1},\\ldots,{{\\mathcal{D}}}_{t-1})\\), and that each dataset \\({{\\mathcal{D}}}_{i}\\) has equal marginal likelihood. Under a second-order expansion of \\({{\\mathcal{C}}}_{t}\\) around \u03bct\u22121 and \u03c3t\u22121, and the small-update assumption \\(| \\frac{\\Delta {\\sigma }_{i}}{{\\sigma }_{i,t-1}}| \\ll 1\\), the parameter updates become:<\/p>\n $$\\Delta {\\boldsymbol{\\sigma }}=-\\frac{{{\\boldsymbol{\\sigma }}}_{t-1}^{2}}{2}\\,\\frac{\\partial {{\\mathcal{C}}}_{t}}{\\partial {{\\boldsymbol{\\sigma }}}_{t-1}}+ \\frac{{{\\boldsymbol{\\sigma }}}_{t-1}}{2\\,N\\,{{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}}({{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}-{{\\boldsymbol{\\sigma }}}_{t-1}^{2}),$$<\/p>\n \n (27)\n <\/p>\n $$\\Delta {\\boldsymbol{\\mu }}=-{{\\boldsymbol{\\sigma }}}_{t-1}^{2}\\,\\frac{\\partial {{\\mathcal{C}}}_{t}}{\\partial {{\\boldsymbol{\\mu }}}_{t-1}}+ \\frac{{{\\boldsymbol{\\sigma }}}_{t-1}^{2}}{N\\,{{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}}({{\\boldsymbol{\\mu }}}_{{\\rm{prior}}}-{{\\boldsymbol{\\mu }}}_{t-1}).$$<\/p>\n \n (28)\n <\/p>\n A full proof (Supplementary Note\u00a04<\/a>) uses a second-order Taylor expansion to handle the implicit dependence of \u03bct and \u03c3t. Under \\(| \\frac{\\Delta {\\boldsymbol{\\sigma }}}{{{\\boldsymbol{\\sigma }}}_{t-1}}| \\ll 1\\), these updates simplify to Eqs. (27<\/a>)\u2013(28<\/a>), balancing learning and forgetting in a single, efficient step.<\/p>\n Links with Newton\u2019s method<\/p>\n We have thus far derived MESU by minimizing the KL divergence between a variational posterior q\u03b8(\u03c9) and the (truncated) true posterior \\(p({\\boldsymbol{\\omega }}| {{\\mathcal{D}}}_{t-N},\\ldots,{{\\mathcal{D}}}_{t})\\). Here, we show how a similar update can emerge from Newton\u2019s method if we treat the free energy as the loss function to be minimized in a non-continual setting.<\/p>\n Newton\u2019s method states that for a parameter \u03c9 inducing loss \\({\\mathcal{L}}\\), the update is \\(\\Delta {\\boldsymbol{\\omega }}=-\\,\\gamma \\,{\\left(\\frac{{\\partial }^{2}{\\mathcal{L}}}{\\partial {{\\boldsymbol{\\omega }}}^{2}}\\right)}^{-1}\\,\\,\\frac{\\partial {\\mathcal{L}}}{\\partial {\\boldsymbol{\\omega }}},\\) where \u03b3 is a step size. Let \\({\\mathcal{F}}\\) be the KL divergence between two multivariate Gaussians: a variational posterior \\({q}_{{\\boldsymbol{\\theta }}}({\\boldsymbol{\\omega }})={\\mathcal{N}}({\\boldsymbol{\\omega }};{\\boldsymbol{\\mu }},{\\rm{diag}}({{\\boldsymbol{\\sigma }}}^{2}))\\) and a true posterior \\(p({\\boldsymbol{\\omega }}| {\\mathcal{D}})\\). If the true posterior is itself an (approximate) mean-field Gaussian \\(\\,p({\\boldsymbol{\\omega }}| {\\mathcal{D}})\\approx {q}_{{\\rm{post}}}({\\boldsymbol{\\omega }})={\\mathcal{N}}({\\boldsymbol{\\omega }};{{\\boldsymbol{\\mu }}}_{{\\rm{post}}},{\\rm{diag}}({{\\boldsymbol{\\sigma }}}_{\\,\\text{post}\\,}^{2})),\\), then<\/p>\n $${\\mathcal{F}}={D}_{KL}[{q}_{{\\boldsymbol{\\theta }}}({\\boldsymbol{\\omega }})\\parallel p({\\boldsymbol{\\omega }}| {\\mathcal{D}})]=\\log \\left(\\frac{{{\\boldsymbol{\\sigma }}}_{{\\rm{post}}}}{{\\boldsymbol{\\sigma }}}\\right)+ \\frac{{{\\boldsymbol{\\sigma }}}^{2}+{({\\boldsymbol{\\mu }}-{{\\boldsymbol{\\mu }}}_{{\\rm{post}}})}^{2}}{2\\,{{\\boldsymbol{\\sigma }}}_{\\,\\text{post}\\,}^{2}}-\\frac{1}{2}.$$<\/p>\n \n (29)\n <\/p>\n To find \u03bcpost and \u03c3post, we apply Bayes\u2019 rule with a Gaussian prior \\(p({\\boldsymbol{\\omega }})={\\mathcal{N}}({\\boldsymbol{\\omega }};{{\\boldsymbol{\\mu }}}_{{\\rm{prior}}},{\\rm{diag}}({{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}))\\) and a Gaussian likelihood \\(p({\\mathcal{D}}| {\\boldsymbol{\\omega }})\\propto {q}_{L}({\\boldsymbol{\\omega }})={\\mathcal{N}}({\\boldsymbol{\\omega }};{{\\boldsymbol{\\mu }}}_{L},{\\rm{diag}}({{\\boldsymbol{\\sigma }}}_{L}^{2}))\\) (cf. Lemma 1). This yields:<\/p>\n $$\\frac{1}{{{\\boldsymbol{\\sigma }}}_{\\,\\text{post}\\,}^{2}}=\\frac{1}{{{\\boldsymbol{\\sigma }}}_{L}^{2}}+ \\frac{1}{{{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}},\\qquad {{\\boldsymbol{\\mu }}}_{{\\rm{post}}}=\\frac{{{\\boldsymbol{\\mu }}}_{L}\\,{{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}+{{\\boldsymbol{\\mu }}}_{{\\rm{prior}}}\\,{{\\boldsymbol{\\sigma }}}_{L}^{2}}{{{\\boldsymbol{\\sigma }}}_{L}^{2}+{{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}}.$$<\/p>\n \n (30)\n <\/p>\n Substituting these into (29<\/a>) gives explicit forms for \\(\\partial {\\mathcal{F}}\/\\partial {\\boldsymbol{\\mu }}\\), \\(\\partial {\\mathcal{F}}\/\\partial {\\boldsymbol{\\sigma }}\\), and their second derivatives. In particular,<\/p>\n $$\\frac{{\\partial }^{2}{\\mathcal{F}}}{\\partial {{\\boldsymbol{\\mu }}}^{2}}=\\frac{1}{{{\\boldsymbol{\\sigma }}}_{L}^{2}}+\\frac{1}{{{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}}=\\frac{1}{{{\\boldsymbol{\\sigma }}}_{\\,\\text{post}\\,}^{2}},\\qquad \\frac{{\\partial }^{2}{\\mathcal{F}}}{\\partial {{\\boldsymbol{\\sigma }}}^{2}}=\\frac{1}{{{\\boldsymbol{\\sigma }}}^{2}}+ \\frac{1}{{{\\boldsymbol{\\sigma }}}_{L}^{2}}+ \\frac{1}{{{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}}.$$<\/p>\n \n (31)\n <\/p>\n Lemma 2: For a mean-field Gaussian q\u03b8(\u03c9) describing a BNN, with \u03c9\u2009=\u2009\u03bc\u2009+\u2009\u03f5 \u22c5 \u03c3, \\({\\boldsymbol{\\epsilon }} \\sim {\\mathcal{N}}(\\overrightarrow{0},{{\\bf{I}}}_{s})\\), the expected Hessian of the negative log-likelihood \\({\\mathcal{L}}\\) satisfies<\/p>\n $${H}_{D}({\\boldsymbol{\\mu }})={{\\mathbb{E}}}_{{\\boldsymbol{\\epsilon }}}\\left[\\frac{{\\partial }^{2}{\\mathcal{L}}}{\\partial {{\\boldsymbol{\\omega }}}^{2}}\\right]=\\frac{1}{{\\boldsymbol{\\sigma }}}\\frac{\\partial {\\mathcal{C}}}{\\partial {\\boldsymbol{\\sigma }}},$$<\/p>\n where \\({\\mathcal{C}}={{\\mathbb{E}}}_{{\\boldsymbol{\\epsilon }}}\\,[{\\mathcal{L}}({\\boldsymbol{\\omega }})]\\).<\/p>\n By Lemma 2, we have \\(\\,1\/{{\\boldsymbol{\\sigma }}}_{L}^{2}=(1\/{\\boldsymbol{\\sigma }})\\,(\\partial {\\mathcal{C}}\/\\partial {\\boldsymbol{\\sigma }})\\), and under N i.i.d. mini-batches, \\(\\,1\/{{\\boldsymbol{\\sigma }}}_{L}^{2}=(N\/{\\boldsymbol{\\sigma }})\\,(\\partial {\\mathcal{C}}\/\\partial {\\boldsymbol{\\sigma }})=N\\,{H}_{D}({\\boldsymbol{\\mu }}).\\) applying Newton\u2019s method then yields:<\/p>\n Theorem 2: Let \\({q}_{{\\boldsymbol{\\theta }}}({\\boldsymbol{\\omega }})={\\mathcal{N}}({\\boldsymbol{\\omega }};{\\boldsymbol{\\mu }},{\\rm{diag}}({{\\boldsymbol{\\sigma }}}^{2}))\\) be a mean-field Gaussian for a BNN, with \u03c9\u2009=\u2009\u03bc\u2009+\u2009\u03f5 \u22c5 \u03c3, \\({\\boldsymbol{\\epsilon }} \\sim {\\mathcal{N}}(\\overrightarrow{0},{{\\bf{I}}}_{s})\\), and with a prior \\(p({\\boldsymbol{\\omega }})={\\mathcal{N}}({\\boldsymbol{\\omega }};{{\\boldsymbol{\\mu }}}_{{\\rm{prior}}},{\\rm{diag}}({{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}))\\). Given \\({\\mathcal{D}}\\) split into N i.i.d. mini-batches, and defining \\({\\mathcal{C}}={{\\mathbb{E}}}_{{\\boldsymbol{\\epsilon }}}[{\\mathcal{L}}({\\boldsymbol{\\omega }})]\\), a diagonal Newton update for \u03c3 and \u03bc with learning rate \u03b3 and \u03c3\u00a0\u2248\u00a0\u03c3post becomes:<\/p>\n $$\\Delta {\\boldsymbol{\\sigma }}=\\frac{\\gamma N}{2}\\,\\left[-\\,{{\\boldsymbol{\\sigma }}}^{2}\\,\\frac{\\partial {\\mathcal{C}}}{\\partial {\\boldsymbol{\\sigma }}}+ \\frac{{\\boldsymbol{\\sigma }}}{N\\,{{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}}({{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}-{{\\boldsymbol{\\sigma }}}^{2})\\right],$$<\/p>\n \n (32)\n <\/p>\n $$\\Delta {\\boldsymbol{\\mu }}=\\gamma N\\left[-\\,{{\\boldsymbol{\\sigma }}}^{2}\\,\\frac{\\partial {\\mathcal{C}}}{\\partial {\\boldsymbol{\\mu }}}+ \\frac{{{\\boldsymbol{\\sigma }}}^{2}}{N\\,{{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}}({{\\boldsymbol{\\mu }}}_{{\\rm{prior}}}-{\\boldsymbol{\\mu }})\\right].$$<\/p>\n \n (33)\n <\/p>\n A full proof appears in Supplementary Note\u00a04<\/a>. Notice the close similarity to MESU\u2019s update rules, especially when \u03b3\u2009=\u20091\/N. This shows how Bayesian learning (with forgetting) can be viewed through the lens of second-order optimization, linking our framework both to biological synapse models and to Newton\u2019s method in variational inference.<\/p>\n Dynamics of standard deviations in the i.i.d. scenario<\/p>\n We now examine the evolution of each standard deviation \u03c3 under an i.i.d. data assumption, isolating its role as an adaptive learning rate. Starting from<\/p>\n $$\\Delta {\\boldsymbol{\\sigma }}=\\gamma \\left[-\\,{{\\boldsymbol{\\sigma }}}_{t-1}^{2}\\,\\frac{\\partial {{\\mathcal{C}}}_{t}}{\\partial {{\\boldsymbol{\\sigma }}}_{t-1}}+ \\frac{{{\\boldsymbol{\\sigma }}}_{t-1}}{N\\,{{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}}\\,({{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}-{{\\boldsymbol{\\sigma }}}_{t-1}^{2})\\right],$$<\/p>\n \n (34)\n <\/p>\n We note that \u03b3 need not be fixed at 0.5 as in our continual-learning derivation\u2014other interpretations (e.g., tempered posteriors44<\/a> or simple Newton steps) can alter its value.<\/p>\n Under Lemmas 1 and 2, we have \\(\\,\\frac{1}{{{\\boldsymbol{\\sigma }}}_{L}^{2}}=\\frac{N}{{\\boldsymbol{\\sigma }}}\\frac{\\partial {\\mathcal{C}}}{\\partial {\\boldsymbol{\\sigma }}}=N\\,{H}_{D}({\\boldsymbol{\\mu }})\\). Treating \u03c3 as a function of discrete time (iteration) t, we obtain a Bernoulli differential equation of the form<\/p>\n $${\\boldsymbol{\\sigma }}{\\prime} (t)+ a(t)\\,{\\boldsymbol{\\sigma }}(t)=b(t)\\,{\\boldsymbol{\\sigma }}{(t)}^{n},$$<\/p>\n \n (35)\n <\/p>\n with n\u2009=\u20093, \\(a(t)=-\\frac{\\gamma }{N}\\), and \\(b(t)=-\\frac{\\gamma }{N}(N{H}_{D}({{\\boldsymbol{\\mu }}}_{0})+\\frac{1}{{{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}}).\\) Solving this via Leibniz substitution yields:<\/p>\n Proposition 1: Consider the Bernoulli differential equation \\(\\sigma {\\prime} (t)+a(t)\\,\\sigma (t)=b(t)\\,\\sigma {(t)}^{n}\\) under \\({\\sigma }_{\\,\\text{prior}\\,}^{2}\\ge \\sigma (0) > 0\\), n\u2009=\u20093, \\(a(t)=-\\frac{\\gamma }{N}\\), and \\(b(t)=-\\frac{\\gamma }{N}(N\\,{H}_{D}({\\mu }_{0})+\\frac{1}{{\\sigma }_{\\,\\text{prior}\\,}^{2}})\\). Its solution is:<\/p>\n $$\\sigma (t)=\\frac{{\\sigma }_{0}\\,{e}^{\\frac{\\gamma t}{N}}}{\\sqrt{\\,1+ N\\,{\\sigma }_{0}^{2}({H}_{D}({\\mu }_{0})\\,+\\,1N\\,{\\sigma }_{\\,\\text{prior}\\,}^{2})\\,({e}^{2\\gamma tN}-1)}},$$<\/p>\n where \u03c30\u2009=\u2009\u03c3(0).<\/p>\n From this closed-form, one sees that the convergence timescale \\({t}_{c}=\\frac{N}{\\gamma }\\) does not directly depend on \u03c30. As t\u00a0\u2192\u00a0\u221e,<\/p>\n $$\\mathop{\\lim }\\limits_{t\\to \\infty }\\,{\\boldsymbol{\\sigma }}{(t)}^{2}=\\frac{1}{N}\\frac{1}{\\,{H}_{D}({{\\boldsymbol{\\mu }}}_{0})+ 1N\\,{{\\boldsymbol{\\sigma }}}_{\\,\\text{prior}\\,}^{2}}.$$<\/p>\n \n (36)\n <\/p>\n Thus, \u03c3(t)2 becomes inversely proportional to the Hessian diagonal (plus a small residual term), aligning well with the intuition that \u03c3 encodes synapse importance.<\/p>\n In practice, mini-batches only approximate i.i.d. data, adding stochasticity to the curvature estimate \\(\\,\\frac{1}{{\\boldsymbol{\\sigma }}}\\,\\frac{\\partial {\\mathcal{C}}}{\\partial {\\boldsymbol{\\sigma }}}\\).<\/p>\n Case \u2009N\u00a0\u2192\u00a0\u221e<\/p>\n When N is infinite, forgetting disappears, and the update rule reduces to \\(\\,{\\sigma }_{t+1}={\\sigma }_{t}(1-\\frac{{\\sigma }_{t}^{2}}{2\\,{\\sigma }_{L}^{2}})\\). Setting \\({\\alpha }_{t}=\\frac{{\\sigma }_{t}}{\\sqrt{2}\\,{\\sigma }_{L}}\\) yields the recurrence \\(\\,{\\alpha }_{t+1}={\\alpha }_{t}(1-{\\alpha }_{t}^{2})\\), which converges to zero but at a rate such that \\(\\,\\mathop{\\lim }\\limits_{t\\to \\infty }{\\alpha }_{t}\\sqrt{2t}=1\\) (the proof is available in Supplementary Note\u00a03<\/a>). Hence,<\/p>\n $$N=\\infty \\Rightarrow \\mathop{\\lim }\\limits_{t\\to \\infty }t\\,{\\boldsymbol{\\sigma }}{(t)}^{2}=\\frac{1}{{H}_{D}({{\\boldsymbol{\\mu }}}_{0})}.$$<\/p>\n \n (37)\n <\/p>\n In this regime, \u03c3(t)2 eventually collapses, reflecting overconfidence and vanishing plasticity\u2014mirroring phenomena seen in FOO-VB Diagonal. By contrast, any finite N preserves a stable nonzero variance, maintaining EU and ensuring that crucial weights are not overconstrained. This highlights the need for a controlled forgetting mechanism in continual learning.<\/p>\n Uncertainty in neural networks<\/p>\n In our experiments, we measure both aleatoric and epistemic uncertainties following the approach of Smith and Gal25<\/a>. They decompose the total uncertainty (TU) of a prediction into two parts: AU and EU. Formally, they use the mutual information between a model\u2019s parameters \u03c9 and a label y conditioned on input x and dataset \\({\\mathcal{D}}\\):<\/p>\n $${\\mathcal{I}}({\\boldsymbol{\\omega }},y| {\\mathcal{D}},x)=H[\\,p(\\,y| x,{\\mathcal{D}})]-{{\\mathbb{E}}}_{p({\\boldsymbol{\\omega }}| {\\mathcal{D}})}\\,H[\\,p(\\,y| x,{\\boldsymbol{\\omega }})],$$<\/p>\n \n (38)\n <\/p>\n $$\\,\\text{EU}\\,=\\,\\text{TU}\\,-\\,\\text{AU}\\,.$$<\/p>\n \n (39)\n <\/p>\n Here, H is the Shannon entropy45<\/a>, and \\(p(y| x,{\\mathcal{D}})\\) is the predictive distribution obtained by averaging over Monte Carlo samples of \u03c9 from \\(p({\\boldsymbol{\\omega }}| {\\mathcal{D}})\\). Concretely, one samples weights \\({{\\boldsymbol{\\omega }}}_{i} \\sim p({\\boldsymbol{\\omega }}| {\\mathcal{D}})\\) to compute p(y\u2223x,\u00a0\u03c9i), then averages these distributions over i.<\/p>\n In this framework, AU captures the irreducible noise in the data or measurement process, while EU reflects uncertainty about the model\u2019s parameters. EU tends to decrease as the model gathers more evidence or restricts the effective memory window (as in MESU), thereby limiting overconfidence. This separation of uncertainties has been widely adopted in the literature11<\/a>,23<\/a>, allowing a more nuanced evaluation of model predictions and their reliability.<\/p>\n MNIST and Permuted MNIST studies (Figs.\u20093<\/a>, 4<\/a>)<\/p>\n We use the standard MNIST dataset26<\/a>, consisting of 28\u2009\u00d7\u200928 grayscale images. All images are standardized by subtracting the global mean and dividing by the global standard deviation. For Permuted MNIST, each image\u2019s pixels are permuted in a fixed, unique manner to create multiple tasks.<\/p>\n Architecture<\/p>\n All networks displayed in Fig.\u20093<\/a> have 50 hidden rectified linear unit (ReLU) neurons, with 784 input neurons for the images and 10 output neurons for each class.<\/p>\n Training procedure<\/p>\n Supplementary Note\u00a05<\/a> describes the steps of the training algorithm. Training was conducted with a single image per batch and a single epoch per task across the entire dataset. For Permuted MNIST (Fig.\u20093<\/a>), accuracy on all permutations was recorded at each task trained upon. For MNIST (Fig.\u20094<\/a>), accuracy was recorded after each epoch. Each algorithm was tested over five runs with varying random seeds to account for initialization randomness.<\/p>\n Neural network initialization<\/p>\n The neural network weights were initialized differently depending on the algorithm. For SGD, SI, Online EWC Online as well as EWC Stream, weights were initialized using Kaiming initialization. Specifically, for a layer l with input dimension nl and output dimension ml, the weights \u03c9l were sampled according to:<\/p>\n $${\\omega }_{i,l} \\sim {\\mathcal{U}}\\left(-\\frac{1}{\\sqrt{{n}_{l}}},\\frac{1}{\\sqrt{{n}_{l}}}\\right).$$<\/p>\n \n (40)\n <\/p>\n For MESU and FOO-VB Diagonal algorithms, mean parameters \u03bci,l were initialized using a reweighted Kaiming initialization and \u03c3i,l were initialized as a constant value:<\/p>\n $${\\mu }_{i,l} \\sim {\\mathcal{U}}\\left(-\\frac{4}{\\sqrt{{n}_{l}}},\\frac{4}{\\sqrt{{n}_{l}}}\\right),{\\sigma }_{i,l}=\\frac{2}{\\sqrt{{n}_{l}}}.$$<\/p>\n \n (41)\n <\/p>\n For Bayesian models, the choice of initial variance strongly influences the posterior\u2019s exploration range during the first few updates, and we therefore treated the initialization scale as a hyper-parameter. The parameters in Eq. (41<\/a>) were obtained by grid search.<\/p>\n Algorithm parameters<\/p>\n Hyperparameters for SGD, EWC Online, and SI were obtained by doing a grid search to maximize accuracy on the first ten tasks of the 200-task Permuted MNIST. MESU uses Eqs. (10<\/a>) and (11<\/a>) directly, which have no learning rate: No tuning of the learning rate is performed (the same is true for FOO-VB Diagonal). For MESU, the memory window N was set to 300,000 to remember about five tasks of Permuted MNIST. The prior distribution for each synapse of MESU is set to \\({\\mathcal{N}}(0,0.06)\\) for MNIST and \\({\\mathcal{N}}(0,1)\\) Permuted MNIST. Table\u20091<\/a> lists the specific parameter values used for each algorithm considered.<\/p>\n Table 1 Hyper-parameter values for each algorithm of Figs.\u20093<\/a> and 4<\/a>Model-width ablation<\/p>\n All protocols above were repeated with four network widths (50, 256, 512, 1024 hidden ReLU units) so as to disentangle the roles of capacity and memory window N. The full results, together with the corresponding uncertainty curves and variance histograms, are provided in Supplementary Note\u00a06<\/a>.<\/p>\n CIFAR studies (Fig. 5<\/a>)<\/p>\n We use the standard CIFAR-10 dataset29<\/a>, consisting of 32\u2009\u00d7\u200932 RGB images. All images are normalized by dividing by 255. Our network comprises four convolutional layers followed by two fully connected layers with ReLU activations (see Table\u20092<\/a>).<\/p>\n Table 2 CIFAR-10\/100 model architecture and dropout parameters. m: number of tasksTraining procedure<\/p>\n In the single-split setting, we first train on CIFAR-10 (Task 1). We then sequentially train new tasks, each comprising ten classes of CIFAR-100 (which we follow consecutively). A multi-head strategy is used: each new task adds ten units (one \u201chead\u201d) to the final layer, and during training, the loss is computed only at the head corresponding to the current task. For testing, we use the head associated with the relevant task. Each task is trained for 60 epochs with a mini-batch size of 200, and dropout is applied (see Table\u20092<\/a>).<\/p>\n When splits are introduced, each task from the single-split procedure is subdivided into splits (sub-tasks). The network first learns the initial split of each task in sequence, and this process is repeated for subsequent splits.<\/p>\n Neural network initialization<\/p>\n The neural network weights were initialized differently depending on the algorithm. For Adam, SI, and EWC weights were initialized using Kaiming initialization, as described in the previous section. For MESU, mean parameters \u03bci,l were initialized using a reweighted Kaiming initialization, and \u03c3i,l were initialized as a constant value. Specifically, for a layer l with input dimension nl and output dimension ml:<\/p>\n $${\\mu }_{i,l} \\sim {\\mathcal{U}}\\left(-\\frac{\\sqrt{2}}{\\sqrt{{n}_{l}}},\\frac{\\sqrt{2}}{\\sqrt{{n}_{l}}}\\right),{\\sigma }_{i,l}=\\frac{1}{2\\sqrt{{m}_{l}}}.$$<\/p>\n \n (42)\n <\/p>\n As in the Permuted MNIST case, the parameters in Eq. (42<\/a>) were obtained by grid search.<\/p>\n Algorithm parameters<\/p>\n To determine the best \u03bb parameter for EWC, we performed the experiment with 1 split for different value, with 0.1\u2009\u2264\u2009\u03bb\u2009\u2264\u200910. Adam is the special case of EWC where \u03bb\u2009=\u20090. To determine the best c parameter for SI, we performed the experiment with one split for different value, with 0.02\u2009\u2264\u2009c\u2009\u2264\u20090.5. Furthermore, as in the original paper14<\/a>, the optimizer is reset at each new task. We set \u03b1\u03bc (for MESU) by choosing the smallest value that maximized test accuracy on the first CIFAR-10 task. Meanwhile, \u03b1\u03c3 (for MESU) was determined based on our theoretical analysis of standard-deviation dynamics in the i.i.d. scenario (see the related \u201cTheoretical Results and Methods\u201d sections). Specifically, as \u03c3 converges on a timescale \\({t}_{c}=\\frac{N}{{\\alpha }_{\\sigma }}\\), we selected \u03b1\u03c3 so that by the end of training on the first task, at least two of these timescales had elapsed, ensuring that the standard deviations had sufficient time to converge. Table\u20093<\/a> lists the specific parameter values used for each algorithm considered.<\/p>\n Table 3 Hyper-parameter values for each algorithm of Fig.\u00a05<\/a>Animals-dataset studies (Fig.\u20092<\/a>)<\/p>\n This study is based on a subset of the \u201cAnimals Detection Images Dataset\u201d from Kaggle46<\/a>. We used a ResNet18 model trained on ImageNet from PyTorch and removed its final fully connected layer to extract 512-dimensional feature vectors for each image. Training was conducted for five epochs per task with a batch size of one image. The neural network weights were initialized in the same way as in the CIFAR experiment. We then set the learning rates \u03b1 (for SGD) and \u03b1\u03bc,\u00a0\u03b1\u03c3 (for MESU) following the same procedure described for the CIFAR experiment. Table\u20094<\/a> lists the specific parameters used for each algorithm.<\/p>\n Table 4 Hyper-parameter values for each algorithm of Fig.\u20092<\/a>Architecture<\/p>\n