Conceptio › Archive › arXiv CS

arXiv CSopen access

Normative Networks for Source Separation via Local Plasticity and Dendritic Computation

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

neuralnetworks

machine learning, deep learning, neural networks

arXiv:2605.19965v1 [cs.LG] 19 May 2026

Normative Networks for Source Separation via Local Plasticity and Dendritic Computation

Bariscan Bozkurt1,3∗ Efe Ali Gorguner2 Francesco Innocenti3,4 Rafal Bogacz3,4 1 Gatsby Computational Neuroscience Unit, University College London, UK 2 Department of Computer Science, University of Oxford, UK 3 Brain Network Dynamics Unit, University of Oxford, UK 4 MRC Centre of Research Excellence in Restorative Neural Dynamics, UK

Abstract Blind source separation (BSS) is a natural framework for studying how latent causes may be recovered from sensory mixtures, but deriving online and biologically plausible algorithms for structured (i.e., constrained to known domains) and potentially correlated sources remains challenging. Recent work has derived neural networks for BSS from maximization of an entropy measure, yet its online implementations involve complex and nonlocal recurrent dynamics. Motivated by this perspective, we propose Predictive Entropy Maximization, which achieves competitive performance in BSS, using only local weight updates. The method employs a close approximation of an entropy measure, yielding an objective function with easily interpretable components. Minimizing this objective leads to a predictive neural architecture in which feedforward synapses follow an error-driven rule (that can be realized through dendritic mechanisms), lateral inhibitory connections are learned with local Hebbian plasticity, and source-domain constraints are enforced through simple output nonlinearities. We derive explicit spectral bounds on the surrogate error, characterizing when the approximation is accurate. Empirically, Predictive Entropy Maximization remains robust under increasing source correlation and observation noise, outperforms biologically plausible algorithms that rely on stronger independence or decorrelation assumptions, and remains competitive with exact determinant- and correlativeinformation-based baselines. These results show how local plasticity and adaptive lateral inhibition can emerge from maximizing a regularized second-order entropy over structured source domains. Our implementation code is available at https: //github.com/BariscanBozkurt/Predictive-Entropy-Maximization.

Introduction and related works

Natural sensory inputs are rarely generated by a single cause. In audition, several sound sources may overlap; in vision, objects, textures, and illumination are superimposed; in olfaction, sensory responses reflect mixtures of odorants. The problem of recovering the underlying causes from mixtures alone is known as blind source separation (BSS) [1, 2]. BSS is a classical problem in signal processing and machine learning, and it has also served as a normative model of sensory representation learning in neuroscience. Sparse and independent decompositions of natural stimuli produce features reminiscent of receptive fields in early sensory systems [3, 4, 5], while the cocktail-party problem illustrates the broader challenge of source separation in the presence of strong interference [6, 7, 8]. ∗ Work completed in part while visiting the University of Oxford.

Preprint.

Because BSS is ill-posed without additional structure, identifiability must come from assumptions on the latent sources. Independent component analysis (ICA), including the informationmaximization (InfoMax) formulation, achieves this through statistical independence [3, 9, 10]. A complementary geometric line of work instead assumes that the sources lie in a structured domain— for example, a simplex (nonnegative entries summing to one) [11, 12], a sparse or bounded set [13, 14], the nonnegative orthant (all coordinates nonnegative) [15, 16], or, more generally, a polytope (a convex bounded set with finitely many vertices) [17, 18, 19]. These geometric approaches seek outputs that spread across the available dimensions of the source domain, through maximization of the determinant of the covariance matrix of output signals. An appealing feature of this approach is that it retains information about inputs as it avoids degenerate low-rank solutions. Another useful property of these methods is that it can recover dependent or correlated sources (under suitable scattering assumptions) [11, 17]. Recent work further shows that determinant-based structured source separation admits an information-theoretic interpretation through a second-order entropy measure, called correlative entropy (based on the log-determinant of the regularized output covariance matrix) [20]. A separate challenge, central to neuroscience and neuromorphic computing, is how to realize such objectives with local dynamics and plasticity. Biologically plausible BSS networks have been developed for ICA settings [21, 22, 23], for nonnegative and bounded source separation through similarity matching [24, 25], and more recently for correlated-source separation under determinant-based objectives [26, 27, 28]. These works show that source-domain structure can be implemented through neural nonlinearities, inhibition, and local learning rules. Nonetheless, the online determinant-based formulations for correlated sources still tend to induce complex recurrent dynamics, because optimizing the exact log-determinant objective naturally brings inverse-covariance into the update. Our goal is therefore to retain the robustness of determinant-based correlated-source separation while replacing these dynamics with a more directly local and biologically plausible recurrent implementation. Specifically, we strive to derive neural networks which satisfy previously proposed criteria of plausibility [29]: (i) local computation, where a neuron’s activity is only based on its inputs, and (ii) local plasticity, where weight modification is based on the variables encoded in pre- and postsynaptic neurons. To this end, we replace the exact log-determinant with an online second-order surrogate obtained by a Taylor expansion around the diagonal part of the output covariance. The resulting objective decomposes into two interpretable terms: one encourages large variance in each output dimension, while the other penalizes redundant dependencies through normalized cross-covariances. This, in turn, yields a predictive neural architecture with a direct local interpretation: feedforward synapses follow a local error-driven rule, recurrent interactions are mediated by local covariance traces and implement adaptive lateral inhibition, and source-domain constraints are enforced through simple output nonlinearities. In this way, we preserve the determinant-based geometric viewpoint while obtaining a more directly local and biologically plausible recurrent implementation. Beyond determinant-based BSS, our construction is also connected to two adjacent lines of work. First, our approach is related to mechanistic models of learning in neural circuits. These models include predictive coding which emphasizes recurrent inference and local prediction errors as a route to learning [30, 31], and supervised Correlative Information Maximization (CorInfoMax) which shows that predictive pathways and lateral interactions can arise from a log-determinant objective that encourages layer activities to spread across their available dimensions [32]. Likewise, predictive-plasticity models highlight the joint role of predictive and Hebbian mechanisms in sensory representation learning [33]. Second, more abstractly, our surrogate objective is also related to redundancy-reduction methods from self-supervised learning, such as VICReg, which combine variance preservation with penalties on cross-component correlations [34]. These latter methods are not biologically plausible models, but they provide a useful point of comparison for the variancecovariance structure of our objective. Our setting remains distinct from both lines of work, since observations arrive one at a time as linear mixtures of latent sources, and the recovered sources must satisfy explicit domain constraints. The main contributions of our Predictive Entropy Maximization model are: • A Taylor approximation of the online determinant-maximization objective, obtained by expanding the log-determinant of the regularized output covariance around its diagonal part. Unlike exact online determinant-based formulations, the resulting updates depend directly on covariance statistics rather than inverse-covariance states.

• A family of predictive neural architectures for different source domains, induced by the proposed approximation. Feedforward synapses follow local error-driven plasticity, recurrent interactions are mediated by local covariance traces, and source-domain constraints are enforced through simple output nonlinearities, as illustrated in Figure 1. • A theoretical characterization of the surrogate, including explicit error bounds (derived through an exact spectral representation of the Taylor remainder). • Empirical evidence that our method remains robust with increasing source correlation and observation noise, competes with bio-plausible baselines built on stronger decorrelation assumptions, and extends naturally to auditory source separation and sparse receptive-field learning. The remainder of the paper is organized as follows. After defining the BSS setting and the determinant-maximization objective (Section 2), we derive the online Taylor-approximated objective and the resulting predictive neural dynamics (Section 3). We then present numerical results evaluating our model (Section 4), before concluding with the limitations and future directions of this work (Section 5).

Problem setup for blind source separation

We now review the formulation of the linear BSS setting considered in this work [1, 2]. Each observation vector x(t) ∈ Rm is generated by linearly mixing an unknown source vector s(t) ∈ P ⊂ Rn through an unknown matrix A ∈ Rm×n : t ∈ {1, 2, . . . , T }.

x(t) = As(t),

Here, the set P denotes a prescribed source domain that encodes structural assumptions on the source vectors. We focus on the determined or overdetermined regime (m ≥ n) and assume that A has full column rank. Stacking the samples in the snapshot matrices S = [s(1), . . . , s(T )] and X = [x(1), . . . , x(T )], the model can be written compactly as X = AS. The goal of BSS is to recover the source vectors from the mixtures alone by learning a separator W ∈ Rn×m such that the outputs y(t) = W x(t) match the original sources. Since linear BSS is non-unique without additional structural assumptions on P, successful recovery is defined only up to permutation and sign ambiguities, as formalized next. Definition 2.1 (Ideal separation). A separator W achieves ideal separation if the recovered outputs satisfy Y = W X = ΠΛS, where Π is a permutation matrix and Λ is a diagonal sign matrix. The BSS model is identifiable over a source domain P if every exact factorization X = A′ S ′ with columns of S ′ in P agrees with the true factorization up to the ideal-separation ambiguity of Definition 2.1. For the source domains considered here, this identifiability property is established in the determinant-maximization literature for polytopic and simplex-structured source models [11, 18]. Determinant-maximization. A natural way to exploit the geometry of P is through the determinant-maximization criterion, which selects, among all feasible outputs constrained to lie in P, those with maximal second-order spread as measured by the log-determinant of their sample autocorrelation or autocovariance matrix [12, 17, 20]. The corresponding recovery guarantees additionally require the source samples to be sufficiently scattered in P. Geometrically, this means that the convex hull of the source samples must capture the shape of P more faithfully than its best ellipsoidal approximation [12, 17]. Accordingly, we consider the batch determinant-maximization problem minimize − log det(Ĉ) W ,Y

subject to

y(t) = W x(t), and y(t) ∈ P,

t = 1, . . . , T.

(1)

In the above objective, Ĉ is the sample autocovariance matrix defined as: Ĉ = where µ̂ = T1

T ⊤ 1X y(t) − µ̂ y(t) − µ̂ , T t=1

(2)

t=1 y(t) is the sample mean.

Entropy interpretation. The above optimization problem has an information-theoretic interpretation. Following recent work [20], the regularized log-determinant of a sample covariance matrix is 3

approximately proportional to a deterministic second-order entropy measure, which we refer to as the correlative entropy (CE): (ε)

ĤCE (Y ) =

1 n log det(Ĉ + εI) + log(2πe), 2 2

(3)

where ε > 0 is a small regularization constant. Under this interpretation, the optimization problem in Eq. 1 can be viewed as maximizing the output correlative entropy subject to the source-domain and mixing constraints. Further details on this second-order entropy and the related correlativeinformation definitions are given in Appendix A. Source domains. In this paper, we consider the following source domains for which the BSS model is identifiable under the scattering assumptions discussed above: • Antisparse sources: Bmax := {s ∈ Rn : maxi |si | ≤ 1}. • Nonnegative antisparse sources: Bmax,+ := Bmax ∩ Rn+ , where Rn+ is the nonnegative orthant. • Sparse sources: B1 := {s ∈ Rn : ∥s∥1 ≤ 1}, where ∥ · ∥1 denotes the ℓ1 norm. • Nonnegative sparse sources: B1,+ := B1 ∩ Rn+ . P • Simplex sources: ∆ := {s ∈ Rn : si ≥ 0 ∀i, i si = 1}.

Online surrogate objective for blind source separation

In this section, we derive an online approximation to the regularized correlative-entropy objective in Eq. 3 and show how it leads to biologically plausible learning dynamics. We also summarize the corresponding approximation error. From batch correlative entropy to an online objective. To make the correlative-entropy objective online, we replace the batch output covariance by an exponentially weighted estimate, following the standard streaming-statistics construction also used in online CorInfoMax [28]. At time t, the algorithm receives x(t), infers y(t) using pre-update weights and statistics, and then updates these state variables. Specifically, we replace Ĉ defined in Eq. 3 with: Ĉ λ (t) =

t ⊤ 1 X t−t′ λ y(t′ ) − µ̂λ (t′ ) y(t′ ) − µ̂λ (t′ ) , η(t) ′

(4)

t =1

Pt Pt ′ ′ where λ ∈ (0, 1) is a forgetting factor, η(t) = t′ =1 λ t−t , and µ̂λ (t) = η(t)−1 t′ =1 λ t−t y(t′ ). The corresponding exact finite-time recursions are given in Appendix B. Since η(t)−1 → 1 − λ as t → ∞, the estimators can be updated online based on a new sample y(t): µ̂λ (t) = λµ̂λ (t − 1) + (1 − λ)y(t),

(5)

⊤ Ĉ (t) = λĈ (t − 1) + (1 − λ) y(t) − µ̂ (t) y(t) − µ̂λ (t) . λ

(6)

Thus, the exact online counterpart of the batch problem is to minimize the negative output-entropy term, − log det(Ĉ λ (t) + εI), under the source-domain constraint y(t) ∈ P and the input-output relation y(t) = W x(t). Existing online formulations based on correlative information maximize a related mutual-information objective built from this output-entropy term together with an additional error-entropy term [28]. Direct gradients of these log-determinant objectives involve inverse covariance or correlation matrices, leading to inverse second-order states updated through the matrix inversion lemma. These inverse-state updates are less directly synapse-local than covariancetrace updates, making the local-plasticity interpretation less immediate; see Appendix A.4. In contrast, our approach introduced below avoids this step altogether. By replacing the exact regularized output-entropy term with a second-order Taylor surrogate, we obtain an objective whose gradients depend directly on covariance statistics rather than on their inverse. This preserves the entropymaximization viewpoint while leading to recurrent interactions and plasticity rules with a more direct local interpretation. Approximate entropy objective. We present here our first result - an approximate objective that can optimised for BSS. Since below we apply a Taylor expansion around the diagonal of the covariance estimate Ĉ λ (t), we decompose it as Ĉ λ (t) = D̂ λ (t) + Ô λ (t), where D̂ λ (t) is diagonal and Ô λ (t) 4

collects the off-diagonal cross-covariances. We apply a second-order Taylor expansion to the log det objective of the regularized covariance matrix, as formalized in Theorem D.1 of Appendix D, n X n n X 1X ĉij (t)2 +R2 (t), (7) − log det Ĉ λ (t) + εI = − log v̂i (t) + ε + 2 i=1 j=1 v̂i (t) + ε v̂j (t) + ε i=1 j̸=i

variance expansion

}

normalized covariance penalty

}

where R2 (t) denotes the Taylor remainder, v̂i (t) := [D̂ λ (t)]ii = [Ĉ λ (t)]ii denotes the i-th output variance, and ĉij (t) := [Ô λ (t)]ij = [Ĉ λ (t)]ij for i ̸= j denotes the corresponding crosscovariance. Up to the additive constant in the entropy definition, maximizing the output correlative entropy is therefore equivalent to minimizing the right-hand side of Eq. 7. The resulting surrogate decomposes this objective into two interpretable terms: (i) a variance-expansion term and (ii) a variance-normalized cross-covariance penalty. Importantly, it does not target exact decorrelation, as uncorrelated y(t) do not necessarily optimize this objective. Instead, it favors outputs with large component-wise variance and sufficiently small normalized cross-covariances, subject to the geometric constraint P. This is precisely the regime relevant to structured source recovery: successful separation does not require zero cross-covariances, but rather nondegenerate spread within the source domain together with controlled inter-component dependence. This also explains why the method can recover correlated sources beyond the scope of ICA-based approaches built on mutual independence, as demonstrated by our numerical results in Section 4. Bound for the surrogate approximation error. The error |R2 (t)| resulting from the above approximation can be quantified directly. Let D̂ λ,ε (t) := D̂ λ (t) + εI, and define B̂ λ,ε (t) := (D̂ λ,ε (t))−1/2 Ô λ (t)(D̂ λ,ε (t))−1/2 . Corollary D.4 in Appendix D shows that the Taylor remainder R2 (t) in Eq. 7 satisfies |R2 (t)| ≤

∥B̂ λ,ε (t)∥2F ∥B̂ λ,ε (t)∥2 . 3 1 + λmin (B̂ λ,ε (t))

(8)

Here, λmin (·) denotes the smallest eigenvalue, ∥·∥2 denotes the spectral norm, and ∥·∥F denotes the Frobenius norm. The numerator ∥B̂ λ,ε (t)∥2F ∥B̂ λ,ε (t)∥2 measures the overall size of the regularized normalized off-diagonal covariance: the Frobenius norm captures its aggregate magnitude, while the spectral norm captures its largest absolute eigenvalue. The denominator 1 + λmin (B̂ λ,ε (t)) shows that the bound worsens as the smallest eigenvalue approaches −1; hence the surrogate is accurate when the normalized off-diagonal covariance is small and the spectrum stays away from −1. In Section 4, Figure 4 reports a diagnostic of the resulting error bound as source correlation varies. Two-timescale optimization. Since biological neural networks can only learn online, we wish to generate the output y and then update the parameters W based on the currently provided input x. Therefore at each step or sample t, first the output y is optimized, and then the weights W are updated. To enable this bi-level optimization, we replace the hard relation y(t) = W x(t) with a quadratic penalty, adding to the online objective of Eq. 7: Jt = −

n X

n X n 1X ĉij (t)2 + γ∥y(t) − W x(t)∥22 , log v̂i (t) + ε + 2 v̂ (t) + ε v̂ (t) + ε i j i=1 i=1 j=1

(9)

j̸=i

and γ > 0 controls the strength of the prediction-error term. To optimize Eq. 9, we adopt a two-timescale procedure (which naturally admits a neural interpretation as we will show later). For each incoming sample, the current output y(t) is first updated on a fast timescale. During this fast relaxation, the input x(t), the weights W , and the running statistics are held fixed, and the output y(t) is updated until the value minimizing Jt over y(t) ∈ P is found. Once this fast inference stage has settled, the resulting output state is used to update the feedforward weights and the running second-order statistics on a slower timescale. Inference of output activity. We first introduce notation that will help us to describe the inference process. For each streaming sample t, we introduce a faster neural-relaxation time index τ = 0, 1, . . . , τmax . The index t labels the slow arrival of new input samples and the corresponding parameter updates, whereas τ labels the fast within-sample activity dynamics used to infer y(t). For 5

notational convenience, let µ̂k (t) := [µ̂λ (t)]k denote the k-th component of the running mean and we define the centered activity as ȳk (t, τ ) = yk (t, τ ) − µ̂k (t). During inference we change the output y to reduce the objective of Eq. 9. Using the gradient of Eq. 9, we define the fast activity-update direction as the negative truncated gradient d(t, τ ) ≈ −∇y Jt , where the truncation discards the higher-order term quadratic in the off-diagonal covariance entries for biological plausibility; see Appendix B.3. Componentwise, ! n m X X ȳk (t, τ ) ĉkj (t) ȳj (t, τ ) − γ yk (t, τ ) − dk (t, τ ) = − Wkℓ (t − 1)xℓ (t) . v̂k (t) + ε j=1 v̂k (t) + ε v̂j (t) + ε ℓ=1 | {z } j̸=k {z } | {z } | predictive correction variance drive

covariance reduction

(10) The corresponding projected activity update is yk (t, τ + 1) = σP (yk (t, τ ) + ηy (t, τ ) dk (t, τ )) ,

(11)

where ηy (t, τ ) > 0 is the inference step size. Terms in Eq. 10 have intuitive interpretation. The first term in dk (t, τ ) promotes regularized variance expansion, the second reduces covariance, and the third pulls the output activity toward its feedforward prediction. The form of σP in Eq. 11 is determined by the source domain. For box-type domains such as Bmax and Bmax,+ , it reduces to elementwise saturation, namely clipping to [−1, 1] or [0, 1]. For B1 , B1,+ , and ∆, we instead project the output on the source domain as in standard proximal algorithms [35], and the corresponding domain-specific derivations are given in Appendix C. Parameter updates. Once the fast inference dynamics have settled for sample t, the resulting output state is used to update the weights W and the running output statistics on a slower timescale. With the settled output y(t) held fixed, the only W -dependent part of Jt is the quadratic prediction loss. A gradient descent step on 12 ∥y(t) − W x(t)∥22 with learning rate αW (t) gives W (t) = W (t − 1) + αW (t)e(t)x(t)⊤ ,

e(t) := y(t) − W (t − 1)x(t),

(12)

as derived in Appendix B (Eq. B.8). The running mean is updated according to Eq. 5, and the centered activity ȳ(t) = y(t)−µ̂λ (t) is then used to update the cross-covariance and variance traces: ĉij (t) = λĉij (t − 1) + (1 − λ)ȳi (t)ȳj (t), i ̸= j,

(13)

v̂i (t) = λv̂i (t − 1) + (1 − λ)ȳi (t)2 , i = 1, . . . , n,

(14)

as derived in Appendix B (see Eq.s B.6 and B.7). We refer to the resulting procedure as Predictive Entropy Maximization (PEM) and summarize it in Algorithm 1 in Appendix C. Neural implementation. PEM can be naturally implemented in a neural network with local weight update rules. In this implementation, we assume that output signals yi are encoded in neural activity with dynamics on a fast time scale, while parameters Wij and ĉij are encoded in the weights of connections between neurons which evolve on a slower time scale. PEM naturally maps onto a two-layer network, where the input layer encodes the mixture inputs x, while the output layer encodes the outputs y. To make the network mapping explicit, we rewrite the activity-update direction from Eq. 10 as dk (t, τ ) = γ

m X

Wkℓ (t − 1)xℓ (t) −

ℓ=1

feedforward input

}

n X

ĉkj (t)ȳj (t, τ ) ȳk (t, τ ) − γyk (t, τ ) + . (15) (v̂ (t) + ε)(v̂ (t) + ε) v̂ k j k (t) + ε j=1 | {z } j̸=k leak | {z } recurrent input

Eq. 15 describes changes in the output signals, which correspond to changes in the output-layer activity in the network in Figure 1a. The first term is the feedforward input to output unit k through weights Wkℓ . The second term maps to lateral recurrent inhibition via connections with weights ĉkj . The remaining terms depend only on neuron k’s own activity (but not on any other neurons in the network), so they can be interpreted as a leak or decay of activity (dependent on parameters µ̂k and v̂k that are intrinsic to the neuron). 6

e 1 y1

e1 y1

x1 e 2 y2

e2 y2

Ĉ λ

Ĉ .. .

λL

.. .

e n yn

.. . en yn

xm (a) Antisparse network

(b) Sparse network

Figure 1: Representative Predictive Entropy Maximization architectures for two source domains. (a) Antisparse architecture. Mixture inputs are mapped through feedforward weights W to local prediction compartments, each paired with an output unit yk . Prediction errors ek are computed as the differences between somatic and dendritic activity. The output layer is coupled through adaptive recurrent inhibitory interactions driven by the running output-covariance statistics Ĉ λ , while the antisparse constraint is enforced locally by clipping. (b) Sparse architecture. The same predictive feedforward pathway and adaptive recurrent inhibitory interactions are retained, but an additional shared inhibitory unit λL sends a common suppressive signal to all outputs. This shared inhibition enforces the sparsity constraint and yields soft-thresholding dynamics.

Assuming domains with limited total activity of output neurons (B1 , B1,+ , and ∆) requires normalization of activity that could be driven by an additional neural population λL , which acts as a shared inhibitory signal enforcing the corresponding population-level constraint (as shown in Figure 1b). The parameter updates in Eqs. 12-14 are local. The feedforward rule of Eq. 12 depends only on the presynaptic activity and the postsynaptic prediction error, and corresponds to a previously proposed plasticity rule [36]. It has been postulated that biological neurons can implement such rule: If input W (t − 1)x(t) comes to dendritic compartments (denoted by grey rectangles in Figure 1a), then the neurons can compute errors e(t) as the differences between their somatic and dendritic activity, and trigger plasticity proportional to these errors [36]. The change in recurrent connections in Eq. 13 is proportional to the product of normalized activity of pre and post-synaptic neurons, hence it corresponds to Hebbian plasticity. This yields a more directly synapse-local interpretation of the recurrent updates compared to CorInfoMax [28] (where it relies on transformed inverse-covariance signals). The variance trace v̂i (t) determines the leak of neuron i and its change in Eq. 14 depends only on the centered activity of neuron i. The updates of the variance trace could correspond to homeostatic regulation [37], that seeks to maintain average neural activity in a desired range. Unnormalized Predictive Entropy Maximization. A potential challenge for the neural implementation of Eq. 15 is the dependence of recurrent inhibition between neurons k and j on v̂k and v̂j , because we assumed above that these parameters determine neurons’ leak, hence it is unclear how they could affect synaptic transmission specifically via recurrent connections. Therefore, in our subsequent experiments, we also consider a simplified and arguably even more biologically plausible variant of PEM called unnormalized PEM (u-PEM). This model is obtained by replacing the variance normalization ((v̂k (t) + ε)(v̂j (t) + ε))−1 from the ‘covariance reduction’ term (Eq. 10) with a parameter γlateral > 0 that controls the strength of the lateral coupling. Hence, this gives the following cost function and the activity-update direction:

Jt = −

n X

n X n X 1 log(v̂i (t) + ε) + γlateral ĉij (t)2 + γ∥y(t) − W (t − 1)x(t)∥22 , 2 i=1 i=1 j=1 j̸=i

(a) Correlated nonnegative antisparse (Bmax,+ )

(b) Noisy sparse (B1 )

Figure 2: Performance comparisons across source domains. (a) Mean component SNR versus correlation ρ for nonnegative antisparse sources. Predictive Entropy Maximization (PEM) and its unnormalized variant (u-PEM) remain robust and outperform online baselines relying on stronger independence or decorrelation assumptions. (b) Mean component SNR versus input SNR for sparse sources. The PEM models stay close to the batch CorInfoMax (Batch) baseline across noise levels while remaining online and biologically plausible. ! m n X X ȳk (t, τ ) Wkℓ (t − 1)xℓ (t) . − γlateral ĉkj (t) ȳj (t, τ ) − γ yk (t, τ ) − dk (t, τ ) = v̂k (t) + ε j=1 ℓ=1

j̸=k

u-PEM maps onto the same neural architecture as before, but now the lateral inhibition is directly implemented via the product of the weights and the neural activity, without needing further normalizing parameters. Further details are given in Appendix E.7.

Numerical results

We compare PEM and its unnormalized variant (u-PEM) with both batch and online baselines. The batch baselines are CorInfoMax (Batch) [20] and ICA-InfoMax [3]; among them, CorInfoMax (Batch) serves as the oracle reference, since it optimizes the exact determinant-maximization objective using the full dataset. The bio-plausible online baselines are CorInfoMax (Online) [28] and Nonnegative Similarity Matching (Online) [24]. Among these, CorInfoMax (Online) is the most directly comparable method, as it targets the exact log-determinant objective that our Taylor surrogate approximates. Our aim is therefore to match the robustness of exact determinant-based correlatedsource separation while obtaining a more directly local and biologically plausible implementation. Performance is measured by the mean component signal-to-noise ratio (mSNR) defined in Eq. E.1. Unless stated otherwise, all results are averaged over 30 independent realizations, and shaded bands denote 95% confidence intervals computed using Eq. E.2. Our Python implementation is provided with the supplementary material. Additional implementation details, extended results, and supplementary experiments, including comparisons with additional baselines, are reported in Appendix E. Synthetic source separation. We first consider synthetic linear mixtures with n = 5 sources, m = 10 mixtures, and T = 105 samples. For each realization, the mixing matrix A ∈ R10×5 has i.i.d. standard normal entries, and the observations are generated as x(t) = As(t) + ϵ(t), ϵ(t) ∼ N (0, σ 2 I10 ), t = 1, . . . , T, with σ 2 chosen to match the prescribed input SNR level SNRin . We evaluate the method in two complementary regimes. In the first, we vary the source correlation level while keeping the observation noise fixed. In the second, we vary the observation noise while keeping the source-domain geometry fixed. For the correlation experiment, we consider the nonnegative antisparse domain Bmax,+ . Sources are generated using a copula-t model, a standard way to generate correlated variables with uniform marginals; here ρ controls the dependence level and is varied over {0, 0.05, . . . , 0.5}, while the mixture SNR is fixed at 30 dB. The results are shown in Figure 2a. Both PEM and its unnormalized variant (u-PEM) remain robust as ρ increases and clearly outperform ICA-InfoMax and NSM, whose performance degrades rapidly as source correlation grows. PEM appears more resistant to increased source correlation than u-PEM, suggesting that the variance normalization inherited from the entropy objective contributes to robustness under correlated sources. 8

Relative to CorInfoMax, the PEM models preserve robust correlated-source separation while replacing inverse-covariance dynamics by the simpler surrogate-based updates introduced in Section 3. As expected, the batch CorInfoMax (Batch) baseline achieves the strongest overall performance. For the noise experiment, we consider the sparse domain B1 . Sources are generated uniformly from B1 , and the mixture SNR is varied over SNRin ∈ {30, 25, . . . , 5} dB. Figure 2b shows that both PEM and u-PEM track the batch CorInfoMax (Batch) baseline closely across the full noise range. This indicates that replacing inverse-covariance dynamics by the Taylor surrogate does not materially degrade source separation performance, while still yielding an online and bio-plausible neural implementation. Additional results for the antisparse, nonnegative sparse, and simplex domains follow the same pattern and are reported in Appendix Figure 6. Learning sparse receptive fields. We evaluate the sparse architecture illustrated in Figure 1b on prewhitened 12 × 12 natural image patches [4].2 The network is trained on vectorized patches in R144 . The learned filters, shown in Figure 3, are localized and oriented, with the familiar Gabor-like structure associated with sparse coding models of the early visual cortex [4]. This experiment shows that the same sparse neural mechanism used for source separation also captures statistically meaningful features from natural images.

Receptive Filters Learned from Natural Scenes

Auditory source separation. We also evaluate the sparse architecture on an audio separation task inspired by the cocktail-party setting [6, 8]. We use three 5second audio sources from the librosa library [38], namely fishin, pistachio, and vibeace, and transform the mixtures with a db4 discrete wavelet trans- Figure 3: Receptive fields learned by the form before running PEM. This choice is motivated sparse Predictive Entropy Maximization by the approximate sparsity of natural sounds in multi- model from natural image patches exhibit scale representations [39, 5, 40]; the transform-domain localized and oriented structure characformulation and representative temporal alignments are teristic of sparse sensory representations. reported in Appendix E. Across 30 random mixing scenarios, the three recovered sources achieve source-wise SNR values of 24.12 ± 1.98 dB, 25.07 ± 2.14 dB, and 21.54 ± 2.17 dB (95% CI), indicating that the sparse-domain dynamics remain effective on realistic audio mixtures. Diagnostic validation of the Taylor surrogate. To complement the analytical results of Section 3, we empirically assess the accuracy of the second-order Taylor surrogate used to derive the neural dynamics. We consider the antisparse synthetic setting with controlled source correlation ρ ∈ {0, . . . , 0.5}, while fixing the input SNR at 30 dB. To make the approximation nontrivial from the outset, we initialize √ the covariance statistic with substantial off-diagonal structure by setting Ĉ λ (0) = GG⊤ , G = 0.2 I +Z, Zij ∼ N (0, 1/5). At each recorded iteration, we compute both the exact Taylor remainder and the sharper spectral upper bound given by the first inequality in Eq. 8, using the current covariance estimate Ĉ λ (t). Figure 4a reports all recorded error–bound pairs across runs and correlation levels. The empirical points lie on or above the identity line, confirming that the analytical bound consistently upperbounds the observed approximation error. Figure 4b summarizes the converged error and bound as functions of ρ. As expected, the approximation error increases with source correlation, but remains small throughout the range considered and continues to be well controlled by the theoretical bound. Additional transient diagnostics, including time-evolution plots of the Taylor approximation error and its upper bound for ρ = 0 and ρ = 0.4, are reported in Appendix E.

2 The dataset and original Sparsenet implementation are available at https://www.rctn.org/bruno/sparsenet/.

103

y=x

Source Correlation ( )

Theoretical Bound

0.4 101

0.3 10 1

0.2

10 3

0.1

10 5 10 5

10 3

10 1

Actual Error

101

103

Mean Surrogate Error at Convergence vs. Correlation

0.5

0.0

Mean Converged Actual Error Mean Converged Theoretical Bound

Error Magnitude (Log Scale)

Actual Error vs. Theoretical Bound

10 1

10 2

0.0

(a) Actual error vs. theoretical bound

0.1

0.2

0.3

Source Correlation ( )

0.4

0.5

(b) Converged error vs. source correlation

Figure 4: Taylor-surrogate diagnostics. (a) Exact Taylor remainder versus the theoretical bound from (Eq. 8), over recorded runs and iterations, color-coded by source correlation ρ. The solid line denotes y = x. (b) Mean Taylor remainder and corresponding bound (Eq. 8) as functions of ρ.

Conclusion

Summary. In this work, we introduced Predictive Entropy Maximization, an online determinantbased entropy maximization framework for blind source separation over identifiable source domains. By replacing the exact log-determinant objective with a second-order Taylor surrogate, we obtained a predictive neural architecture in which feedforward plasticity is local and error-driven, while recurrent interactions are mediated by covariance traces and implement adaptive lateral inhibition. We also provided explicit spectral bounds for the surrogate error, clarifying when this approximation is accurate. Empirically, the resulting method remains robust under increasing source correlation and observation noise, compares favorably with existing biologically plausible baselines, and extends naturally to sparse receptive-field learning and auditory source separation. Taken together, these results provide a normative account of how structured source separation can give rise to simple local plasticity rules and interpretable inhibitory dynamics in recurrent networks. Limitations and future directions. The performance of the proposed model can be sensitive to hyperparameter choices and to the number of observed mixtures, as analyzed empirically in Appendix E.5-E.6. In addition, similar to other biologically plausible recurrent BSS networks, the computational cost of our method scales with the number of fast neural-dynamics iterations (see Appendix E.9), motivating future work on faster inference schemes.

Acknowledgements B.B. is supported by the Gatsby Charitable Foundation. This work was supported by the Wellcome Trust grant 313955/Z/24/Z and the Medical Research Council grant UKRI/MR/B000936/1. Part of this work was completed while B.B. was a visiting PhD student at the University of Oxford. B.B. thanks Vasco Portilheiro and Houssam Zenati for helpful discussions on an earlier version of this work.

References [1] Pierre Comon and Christian Jutten. Handbook of Blind Source Separation: Independent component analysis and applications. Academic press, 2010. [2] Andrzej Cichocki, Rafal Zdunek, Anh Huy Phan, and Shun-ichi Amari. Nonnegative matrix and tensor factorizations: applications to exploratory multi-way data analysis and blind source separation. John Wiley & Sons, 2009. [3] Anthony J. Bell and Terrence J. Sejnowski. An information-maximization approach to blind separation and blind deconvolution. Neural Comput., 7(6):1129–1159, November 1995. ISSN 0899-7667. doi: 10.1162/neco.1995.7.6.1129. URL https://doi.org/10.1162/neco. 1995.7.6.1129. [4] Bruno A. Olshausen and David J. Field. Emergence of simple-cell receptive field properties by learning a sparse code for natural images. Nature, 381(6583):607–609, 1996. doi: 10.1038/ 381607a0. [5] Michael S. Lewicki. Efficient coding of natural sounds. Nature Neuroscience, 5(4):356–363, 2002. doi: 10.1038/nn831. URL https://doi.org/10.1038/nn831. [6] E Colin Cherry. Some experiments on the recognition of speech, with one and with two ears. The Journal of the Acoustical Society of America, 25(5):975–979, 1953. [7] Adelbert W Bronkhorst. The cocktail party phenomenon: A review of research on speech intelligibility in multiple-talker conditions. Acta Acustica united with Acustica, 86(1):117– 128, 2000. [8] Simon Haykin and Zhe Chen. The cocktail party problem. Neural Computation, 17(9):1875– 1902, 2005. doi: 10.1162/0899766054322964. [9] Aapo Hyvärinen and Erkki Oja. Independent component analysis: A tutorial, 1999. URL https://api.semanticscholar.org/CorpusID:118629. [10] Aapo Hyvärinen and Erkki Oja. Independent component analysis: Algorithms and applications. Neural Networks, 13(4–5):411–430, 2000. doi: 10.1016/S0893-6080(00)00026-5. [11] Chia-Hsiang Lin, Wing-Kin Ma, Wei-Chiang Li, Chong-Yung Chi, and Arulmurugan Ambikapathi. Identifiability of the simplex volume minimization criterion for blind hyperspectral unmixing: The no pure-pixel case. IEEE Transactions on Geoscience and Remote Sensing, 53, 06 2014. doi: 10.1109/TGRS.2015.2424719. [12] Chia-Hsiang Lin, Ruiyuan Wu, Wing-Kin Ma, Chong-Yung Chi, and Yue Wang. Maximum volume inscribed ellipsoid: A new simplex-structured matrix factorization framework via facet enumeration and convex optimization. SIAM Journal on Imaging Sciences, 11(2):1651–1679, 2018. doi: 10.1137/17M114145X. URL https://doi.org/10.1137/17M114145X. [13] Eren Babatas and Alper T Erdogan. An algorithmic framework for sparse bounded component analysis. IEEE Transactions on Signal Processing, 66(19):5194–5205, August 2018. [14] Alper T Erdogan. A class of bounded component analysis algorithms for the separation of both independent and dependent sources. IEEE Transactions on Signal Processing, 61(22): 5730–5743, August 2013. [15] Pentti Paatero and Unto Tapper. Positive matrix factorization: A non-negative factor model with optimal utilization of error estimates of data values. Environmetrics, 5(2):111–126, 1994. doi: https://doi.org/10.1002/env.3170050203. URL https://onlinelibrary.wiley.com/ doi/abs/10.1002/env.3170050203. [16] Xiao Fu, Kejun Huang, Nicholas D. Sidiropoulos, and Wing-Kin Ma. Nonnegative matrix factorization for signal and data analytics: Identifiability, algorithms, and applications. IEEE Signal Processing Magazine, 36(2):59–80, 2019. doi: 10.1109/MSP.2018.2877582. 11

[17] Gokcan Tatli and Alper T. Erdogan. Polytopic matrix factorization: Determinant maximization based criterion and identifiability. IEEE Transactions on Signal Processing, 69:5431–5447, 2021. doi: 10.1109/TSP.2021.3112918. [18] Gokcan Tatli and Alper T. Erdogan. Generalized polytopic matrix factorization. In ICASSP 2021 - 2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3235–3239, 2021. doi: 10.1109/ICASSP39728.2021.9413709. [19] Bariscan Bozkurt and Alper T. Erdogan. On identifiable polytope characterization for polytopic matrix factorization. In ICASSP 2022 - 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3343–3347, 2022. doi: 10.1109/ ICASSP43922.2022.9746370. [20] Alper T. Erdogan. An information maximization based blind source separation approach for dependent and independent sources. In ICASSP 2022 - 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 4378–4382, 2022. doi: 10.1109/ ICASSP43922.2022.9746099. [21] Takuya Isomura and Taro Toyoizumi. A local learning rule for independent component analysis. Scientific Reports, 6(1):28073, Jun 2016. ISSN 2045-2322. doi: 10.1038/srep28073. URL https://doi.org/10.1038/srep28073. [22] Takuya Isomura and Taro Toyoizumi. Error-gated hebbian rule: A local learning rule for principal and independent component analysis. Scientific Reports, 8(1):1835, Jan 2018. ISSN 2045-2322. doi: 10.1038/s41598-018-20082-0. URL https://doi.org/10.1038/ s41598-018-20082-0. [23] Yanis Bahroun, Dmitri Chklovskii, and Anirvan M. Sengupta. A normative and biologically plausible algorithm for independent component analysis. In A. Beygelzimer, Y. Dauphin, P. Liang, and J. Wortman Vaughan, editors, Advances in Neural Information Processing Systems, 2021. URL https://openreview.net/forum?id=fpvUKdqcPV. [24] Cengiz Pehlevan, Sreyas Mohan, and Dmitri B. Chklovskii. Blind nonnegative source separation using biological neural networks. Neural Computation, 29(11):2925–2954, 11 2017. ISSN 0899-7667. doi: 10.1162/neco_a_01007. URL https://doi.org/10.1162/neco_a_ 01007. [25] Alper Tunga Erdogan and Cengiz Pehlevan. Blind bounded source separation using neural networks with local learning rules. ICASSP 2020 - 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pages 3812–3816, 2020. URL https: //api.semanticscholar.org/CorpusID:215745493. [26] Berfin Simsek and Alper T. Erdogan. Online bounded component analysis: A simple recurrent neural network with local update rule for unsupervised separation of dependent and independent sources. In 2019 53rd Asilomar Conference on Signals, Systems, and Computers, pages 1639–1643, 2019. doi: 10.1109/IEEECONF44664.2019.9048916. [27] Bariscan Bozkurt, Cengiz Pehlevan, and Alper Tunga Erdogan. Biologically-plausible determinant maximization neural networks for blind separation of correlated sources. In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances in Neural Information Processing Systems, 2022. URL https://openreview.net/forum?id=espX_4CLr46. [28] Bariscan Bozkurt, Ateş İsfendiyaroğlu, Cengiz Pehlevan, and Alper Tunga Erdogan. Correlative information maximization based biologically plausible neural networks for correlated source separation. In The Eleventh International Conference on Learning Representations, 2023. URL https://openreview.net/forum?id=8JsaP7j1cL0. [29] Rafal Bogacz. A tutorial on the free-energy framework for modelling perception and learning. Journal of mathematical psychology, 76:198–211, 2017. [30] Rajesh PN Rao and Dana H Ballard. Predictive coding in the visual cortex: a functional interpretation of some extra-classical receptive-field effects. Nature neuroscience, 2(1):79–87, 1999. doi: 10.1038/4580. 12

[31] James CR Whittington and Rafal Bogacz. Theories of error back-propagation in the brain. Trends in cognitive sciences, 23(3):235–250, 2019. [32] Bariscan Bozkurt, Cengiz Pehlevan, and Alper Tunga Erdogan. Correlative information maximization: A biologically plausible approach to supervised deep neural networks without weight symmetry. In Thirty-seventh Conference on Neural Information Processing Systems, 2023. URL https://openreview.net/forum?id=TUGoUNkccV. [33] Manu Srinath Halvagal and Friedemann Zenke. The combination of hebbian and predictive plasticity learns invariant object representations in deep sensory networks. Nature Neuroscience, 26(11):1906–1915, Nov 2023. ISSN 1546-1726. doi: 10.1038/s41593-023-01460-y. URL https://doi.org/10.1038/s41593-023-01460-y. [34] Adrien Bardes, Jean Ponce, and Yann LeCun. VICReg: Variance-invariance-covariance regularization for self-supervised learning. In International Conference on Learning Representations, 2022. URL https://openreview.net/forum?id=xm6YD62D1Ub. [35] Neal Parikh, Stephen Boyd, et al. Proximal algorithms. Foundations and trends® in Optimization, 1(3):127–239, 2014. [36] Robert Urbanczik and Walter Senn. Learning by the dendritic prediction of somatic spiking. Neuron, 81(3):521–528, 2014. [37] Gina Turrigiano. Homeostatic synaptic plasticity: Local and global mechanisms for stabilizing neuronal function. Cold Spring Harbor Perspectives in Biology, 4(1):a005736, 2012. doi: 10.1101/cshperspect.a005736. [38] Brian McFee, Matt McVicar, Daniel Faronbi, Iran Roman, Matan Gover, Stefan Balke, Scott Seyfarth, Ayoub Malek, Colin Raffel, Vincent Lostanlen, Benjamin van Niekirk, Dana Lee, Frank Cwitkowitz, Frank Zalkow, Oriol Nieto, Dan Ellis, Jack Mason, Kyungyun Lee, Bea Steers, Emily Halvachs, Carl Thomé, Fabian Robert-Stöter, Rachel Bittner, Ziyao Wei, Adam Weiss, Eric Battenberg, Keunwoo Choi, Ryuichi Yamamoto, CJ Carr, Alex Metsai, Stefan Sullivan, Pius Friesch, Asmitha Krishnakumar, Shunsuke Hidaka, Steve Kowalik, Fabian Keller, Dan Mazur, Alexandre Chabot-Leclerc, Curtis Hawthorne, Chandrashekhar Ramaprasad, Myungchul Keum, Juanita Gomez, Will Monroe, Viktor Andreevitch Morozov, Kian Eliasi, nullmightybofo, Paul Biberstein, N. Dorukhan Sergin, Romain Hennequin, Rimvydas Naktinis, beantowel, Taewoon Kim, Jon Petter Åsen, Joon Lim, Alex Malins, Darío Hereñú, Stef van der Struijk, Lorenz Nickel, Jackie Wu, Zhen Wang, Tim Gates, Matt Vollrath, Andy Sarroff, Xiao-Ming, Alastair Porter, Seth Kranzler, Voodoohop, Mattia Di Gangi, Helmi Jinoz, Connor Guerrero, Abduttayyeb Mazhar, toddrme2178, Zvi Baratz, Anton Kostin, Xinlu Zhuang, Cash TingHin Lo, Pavel Campr, Eric Semeniuc, Monsij Biswal, Shayenne Moura, Paul Brossier, Hojin Lee, Waldir Pimenta, Jon Petter Åsen, Shin Hyun, Iliya S, Eugene Rabinovich, Geo Lei, Jize Guo, Phillip S.M. Skelton, Matt Pitkin, Anmol Mishra, Slava Chaunin, BenedictSt, Scott VanRavenswaay, and David Südholt. librosa/librosa: 0.11.0, March 2025. URL https://doi.org/10.5281/zenodo.15006942. [39] Evan C Smith and Michael S Lewicki. Efficient auditory coding. Nature, 439(7079):978–982, 2006. [40] Ingrid Daubechies. Ten Lectures on Wavelets. SIAM, Philadelphia, PA, 1992. ISBN 978-089871-274-2. [41] Serdar Ozsoy, Shadi Hamdan, Sercan O Arik, Deniz Yuret, and Alper T Erdogan. Selfsupervised learning with an information maximization criterion. In Alice H. Oh, Alekh Agarwal, Danielle Belgrave, and Kyunghyun Cho, editors, Advances in Neural Information Processing Systems, 2022. URL https://openreview.net/forum?id=5MgZAu2NR7X. [42] Jean-Bastien Grill, Florian Strub, Florent Altché, Corentin Tallec, Pierre Richemond, Elena Buchatskaya, Carl Doersch, Bernardo Avila Pires, Zhaohan Guo, Mohammad Gheshlaghi Azar, et al. Bootstrap your own latent-a new approach to self-supervised learning. Advances in neural information processing systems, 33:21271–21284, 2020. 13

[43] Jure Zbontar, Li Jing, Ishan Misra, Yann LeCun, and Stéphane Deny. Barlow twins: Selfsupervised learning via redundancy reduction. In Marina Meila and Tong Zhang, editors, Proceedings of the 38th International Conference on Machine Learning, ICML 2021, 18-24 July 2021, Virtual Event, Proceedings of Machine Learning Research, pages 12310–12320. PMLR, 2021. URL http://proceedings.mlr.press/v139/zbontar21a.html. [44] Yann LeCun et al. A path towards autonomous machine intelligence version 0.9. 2, 2022-0627. Open Review, 62(1):1–62, 2022. [45] Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey E. Hinton. A simple framework for contrastive learning of visual representations. CoRR, abs/2002.05709, 2020. URL https://arxiv.org/abs/2002.05709. [46] Horace B. Barlow. Possible principles underlying the transformations of sensory messages. In Walter A. Rosenblith, editor, Sensory Communication, pages 217–234. MIT Press, Cambridge, MA, 1961. [47] Katrina Drozdov, Ravid Shwartz-Ziv, and Yann LeCun. Video representation learning with joint-embedding predictive architectures, 2024. URL https://arxiv.org/abs/2412. 10925.

Appendix A Review of Correlative Information Maximization

A.1 Correlative entropy and mutual information . . . . . . . . . . . . . . . . . . . . .

A.2 Batch CorInfoMax objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

A.3 Online CorInfoMax objective and neural dynamics . . . . . . . . . . . . . . . . .

A.4 Comparison with Predictive Entropy Maximization. . . . . . . . . . . . . . . . . .

B Derivation of the update rules for the online optimization objective

B.1 Recursive updates for the running mean and covariance . . . . . . . . . . . . . . .

B.2 Derivation of the feedforward weight update . . . . . . . . . . . . . . . . . . . . .

B.3 Gradient of the surrogate objective for the output dynamics . . . . . . . . . . . . .

C Domain-specific network realizations

C.1 Antisparse sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

C.2 Nonnegative antisparse sources . . . . . . . . . . . . . . . . . . . . . . . . . . . .

C.3 Sparse sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

C.4 Nonnegative sparse sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

C.5 Simplex sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

D Error analysis and analytical bounds for the second-order Taylor approximation

D.1 Pointwise remainder representation and bounds . . . . . . . . . . . . . . . . . . .

D.2 Batch surrogate optimality relative to the exact determinant objective . . . . . . . .

E Supplementary on numerical experiments

E.1 Performance evaluation metric . . . . . . . . . . . . . . . . . . . . . . . . . . . .

E.2 Transform invariance in auditory source separation . . . . . . . . . . . . . . . . .

E.3 Antisparse, nonnegative sparse, and simplex source separation examples . . . . . .

E.4 Transient diagnostics for the Taylor approximation . . . . . . . . . . . . . . . . .

E.5 Simulation hyperparameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

E.6 Ablation Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

E.7 Unnormalized Predictive Entropy Maximization . . . . . . . . . . . . . . . . . . .

E.8 Variance-covariance regularisation in self-supervised learning . . . . . . . . . . .

E.9 Computational complexity of the proposed method . . . . . . . . . . . . . . . . .

Review of Correlative Information Maximization

This section briefly reviews the correlative-information-maximization (CorInfoMax) framework underlying the batch method of Erdogan [20] and its online neural implementation in Bozkurt et al. [28]. Our purpose is to make explicit the information-theoretic quantities that motivate CorInfoMax, the exact batch and online objectives it optimizes, and how its recurrent learning rules differ from those derived in the present paper. A.1

Correlative entropy and mutual information

Following Erdogan [20], we use the correlative entropy as a second-order entropy measure. For a finite set of vectors X = [x(1), . . . , x(T )] ∈ Rm×T with sample autocorrelation matrix T

R̂x =

1X x(t)x(t)⊤ , T t=1

the corresponding deterministic correlative entropy is m 1 log det R̂x + εI + log(2πe), 2 2 where ε > 0 is a small regularization constant. This quantity was introduced in Erdogan [20] under the name log-determinant entropy. In later work, the terminology correlative entropy has been adopted for the same second-order entropy measure [28, 32, 41], and we follow that convention here. (ε)

ĤCE (X) :=

For a pair of finite sets X = [x(1), . . . , x(T )] ∈ Rm×T and Y = [y(1), . . . , y(T )] ∈ Rn×T , define T

R̂xy =

1X x(t)y(t)⊤ , T t=1

The deterministic joint correlative entropy is 1 R̂x + εI (ε) ĤCE (X, Y ) := log det 2 R̂yx

⊤ R̂yx = R̂xy .

R̂xy R̂y + εI

m+n log(2πe). 2

Using the block determinant identity, this can be written as (ε)

(ε)

ĤCE (X, Y ) = ĤCE (X) + ĤCE (Y |L X), where

n 1 log det R̂e + εI + log(2πe), 2 2 −1 R̂e := R̂y − R̂yx R̂x + εI R̂xy .

(ε)

ĤCE (Y |L X) := and

(A.1)

The notation Y |L X emphasizes that this is the correlative entropy of the residual under the best linear predictor of Y from X, not Shannon conditional entropy. The corresponding deterministic correlative mutual information is defined by (ε)

(ε)

ÎCMI (X, Y ) := ĤCE (Y ) − ĤCE (Y |L X). (ε)

For sufficiently small ε, ĤCE (Y |L X) corresponds to the correlative entropy of the error of the best linear minimum-mean-square predictor of y from x. Thus, the associated mutual information measures linear/correlative dependence rather than general statistical dependence. A.2

Batch CorInfoMax objective

Let X = [x(1), . . . , x(T )] ∈ Rm×T denote the mixture matrix and Y = [y(1), . . . , y(T )] ∈ Rn×T the separator outputs, constrained by y(t) ∈ P for all t. The batch CorInfoMax objective is the maximization of the deterministic correlative information: 1 1 (ε) maximize ÎCMI (X, Y ) = log det R̂y + εI − log det R̂e + εI , 2 2 Y ∈Rn×T 16

y(t) ∈ P,

subject to

t = 1, . . . , T.

Here R̂e is defined by Equation A.1. For sufficiently small ε, it is the sample correlation matrix of the error of the best linear minimum-mean-square predictor of y from x. In the non-zero-mean setting, the same formulation may be written with sample autocovariances in place of sample autocorrelations. Under the same identifiable-domain and sufficient-scattering assumptions discussed in Section 2, global optima of Equation A.2 recover the sources up to ideal separation. For the purposes of the present paper, the key point is that batch CorInfoMax optimizes the exact regularized log-determinant objective. Consequently, inverse second-order matrices appear directly in the gradient and become central dynamical variables in the online formulation. A.3

Online CorInfoMax objective and neural dynamics

The online CorInfoMax construction replaces the batch correlation matrices by exponentially weighted versions [28]. Using notation parallel to the main text, define t

R̂yy (t) := R̂eλe (t) :=

1 X t−t′ ′ λy y(t )y(t′ )⊤ , ηy (t) ′

(A.3)

1 ηe (t) ′

(A.4)

t =1 t X

′

λet−t e(t′ )e(t′ )⊤ ,

t =1

with forgetting factors λy , λe ∈ (0, 1) and normalization factors ηy (t) =

t X

′ λyt−t ,

ηe (t) =

t′ =1

t X

′

λet−t .

t′ =1

The error vector is e(t) = y(t) − W (t)x(t). The corresponding online CorInfoMax problem at time t is 1 1 λ maximize JtCI (y(t)) := log det R̂yy (t) + εI − log det R̂eλe (t) + εI , n 2 2 y(t)∈R subject to y(t) ∈ P. In the online formulation, the best linear predictor W (t) is treated as an adaptive parameter and updated through an online regularized least-squares criterion, which yields the usual LMS-type rule. Define the inverse-correlation states −1 −1 λ , Be (t) := R̂eλe (t) + εI . By (t) := R̂yy (t) + εI Gradient ascent on Equation A.5 yields ∇y(t) JtCI = γy (t) By (t − 1)y(t) − γe (t) Be (t − 1)e(t), where −1 γy (t) := λy ηy (t−1)+y(t)⊤ By (t−1)y(t) ,

−1 γe (t) := λe ηe (t−1)+e(t)⊤ Be (t−1)e(t) .

Depending on the domain P, CorInfoMax then combines this gradient step with either a projected ascent update (for box-type domains) or a primal-dual/proximal step (for sparse and simplex domains). The feedforward predictor is updated by the local LMS rule W (t + 1) = W (t) + µW (t)e(t)x(t)⊤ . By contrast, the lateral learning rules are based on inverse-correlation states. Applying the Matrix Inversion Lemma to Equations A.3 and A.4 yields the exact recursions By (t + 1) =

ηy (t) By (t) − γy (t) By (t)y(t)y(t)⊤ By (t) , λy ηy (t − 1) 17

(A.6)

Be (t + 1) =

ηe (t) Be (t) − γe (t) Be (t)e(t)e(t)⊤ Be (t) . λe ηe (t − 1)

(A.7)

In Bozkurt et al. [28], these exact inverse-state updates are then simplified to obtain a biologically plausible implementation. The inverse error-correlation update in Equation A.7 is not retained as a learned state; instead, the approximation R̂eλe (t) + εI ≈ εI,

hence

Be (t) ≈ ε−1 I,

(A.8)

is invoked, motivated by the expectation that the prediction error becomes small in the noiseless linear setting. In addition, for λy close to 1 and t sufficiently large, the scalar factor γy (t) is approximated by 1 − λy γy (t) ≈ . (A.9) λy Substituting Equation A.9 into Equation A.6 yields the simplified lateral update used in the CorInfoMax neural networks: 1 − λy −1 ⊤ By (t + 1) = λy By (t)y(t)y(t) By (t) . By (t) − (A.10) λy A.4

Comparison with Predictive Entropy Maximization.

The distinction from the present paper can now be stated precisely. In CorInfoMax, the exact regularized log-determinant objective produces gradients involving inverse second-order states, so recurrent plasticity is expressed through By (t) rather than through covariance traces. Even after the approximations in Equations A.8 and A.9, the lateral update in Equation A.10 remains a rank-one update built from transformed population activities. Indeed, if we define z(t) := By (t)y(t), then Equation A.10 becomes 1 − λy ⊤ z(t)z(t) . By (t + 1) = λ−1 B (t) − y y λy Hence the update of an individual lateral coefficient depends on the transformed signal z(t), not directly on the co-activity of the neuron pair it couples. In the 3 × 3 case, if # # " " y1 (t) b11 (t) b12 (t) b13 (t) y(t) = y2 (t) , By (t) = b21 (t) b22 (t) b23 (t) , y3 (t) b31 (t) b32 (t) b33 (t) then z1 (t) = b11 (t)y1 (t)+b12 (t)y2 (t)+b13 (t)y3 (t),

z2 (t) = b21 (t)y1 (t)+b22 (t)y2 (t)+b23 (t)y3 (t),

and Equation A.10 gives b21 (t + 1) = λ−1 y

1 − λy b21 (t) − z2 (t)z1 (t) . λy

Thus, the update of b21 (t) depends not only on the pair (y2 (t), y1 (t)), but also on y3 (t). In general, a single recurrent update depends on the activity of the full output population. This violates the local plasticity condition of biological plausibility, as defined in the introduction. By contrast, Predictive Entropy Maximization introduces a second-order Taylor surrogate of the regularized log-determinant and works directly with covariance traces rather than inverse-covariance states. The corresponding recurrent statistics are the exponentially weighted variance and crosscovariance traces of the centered outputs, defined by v̂i (t) := [Ĉ λ (t)]ii and ĉij (t) := [Ĉ λ (t)]ij for i ̸= j. Their recursions are v̂i (t) = λv̂i (t − 1) + (1 − λ)ȳi (t)2 ,

ĉij (t) = λĉij (t − 1) + (1 − λ)ȳi (t)ȳj (t),

so the cross-covariance update is written directly in terms of centered pairwise co-activity, while the variance update depends only on the squared centered activity of the corresponding neuron. This is the key distinction emphasized throughout the main text: the surrogate formulation removes the need to maintain inverse-covariance states and yields a more directly local interpretation of recurrent plasticity. 18

Derivation of the update rules for the online optimization objective

In this section, we derive the recursive updates for the online second-order statistics and the feedforward weights, as well as the output gradient used in the fast neural dynamics. B.1

Recursive updates for the running mean and covariance

Recall from Equation 4 that the exponentially weighted output covariance is defined by Ĉ λ (t) =

t ⊤ 1 X t−t′ λ y(t′ ) − µ̂λ (t′ ) y(t′ ) − µ̂λ (t′ ) , η(t) ′ t =1

where η(t) =

t X

′

λ t−t =

t′ =1

1 − λt . 1−λ

1 Since η(t) = λη(t − 1) + 1, we define α(t) := η(t) .

Using this identity, we isolate the most recent sample in the covariance sum and obtain " t−1 # X 1 t−1−t′ ′ λ ′ ′ λ ′ ⊤ λ λ λ y(t ) − µ̂ (t ) y(t ) − µ̂ (t ) Ĉ (t) = η(t) ′ t =1

⊤ 1 + y(t) − µ̂λ (t) y(t) − µ̂λ (t) η(t) ⊤ λη(t − 1) λ 1 = Ĉ (t − 1) + y(t) − µ̂λ (t) y(t) − µ̂λ (t) η(t) η(t) λ ⊤ = 1 − α(t) Ĉ (t − 1) + α(t) y(t) − µ̂λ (t) y(t) − µ̂λ (t) .

(B.1)

Exactly the same argument applied to the exponentially weighted mean t

µ̂λ (t) =

1 X t−t′ ′ λ y(t ) η(t) ′ t =1

yields the recursion µ̂λ (t) = 1 − α(t) µ̂λ (t − 1) + α(t)y(t).

(B.2)

For large t, the normalization converges to η(t) ≈ (1 − λ)−1 , hence α(t) ≈ 1 − λ. Under this steady-state approximation, Equations B.1 and B.2 reduce to µ̂λ (t) = λµ̂λ (t − 1) + (1 − λ)y(t),

(B.3)

⊤ Ĉ (t) = λĈ (t − 1) + (1 − λ) y(t) − µ̂ (t) y(t) − µ̂λ (t) . λ

(B.4)

For convenience, define the centered activity ȳ(t) := y(t) − µ̂λ (t). Then Equation B.4 can be written as Ĉ λ (t) = λĈ λ (t − 1) + (1 − λ)ȳ(t)ȳ(t)⊤ . We now decompose the covariance matrix as Ĉ λ (t) = D̂ λ (t) + Ô λ (t),

(B.5)

where D̂ λ (t) is diagonal and Ô λ (t) has zero diagonal. Reading off the diagonal and off-diagonal parts of Equation B.4 yields v̂i (t) = λv̂i (t − 1) + (1 − λ)ȳi (t)2 , i = 1, . . . , n, ĉij (t) = λĉij (t − 1) + (1 − λ)ȳi (t)ȳj (t), i ̸= j. 19

(B.6) (B.7)

B.2

Derivation of the feedforward weight update

Let e(t) := y(t) − W x(t) denote the instantaneous prediction error. We define the instantaneous quadratic cost 1 1 Jt (W ) ∝ ∥e(t)∥22 = ∥y(t) − W x(t)∥22 . 2 2 Differentiating with respect to W gives 1 ∂ ⊤ (y(t) − W x(t)) (y(t) − W x(t)) ∇W Jt (W ) = ∂W 2 = −(y(t) − W x(t))x(t)⊤ = −e(t)x(t)⊤ . A gradient step with learning rate αW (t) therefore yields W (t) = W (t − 1) − αW (t)∇W Jt (W (t − 1)) = W (t − 1) + αW (t)e(t)x(t)⊤ .

(B.8)

This formulation yields a biologically plausible, local error-modulated Hebbian learning rule where the synaptic update is proportional to the product of the presynaptic activity x(t) and the local postsynaptic error signal e(t). B.3

Gradient of the surrogate objective for the output dynamics

We now derive the gradient of the fast-time-scale objective with respect to the current output y(t). Recall from Equation 9 that Jt (y) = −

n X n 1X ĉij (t)2 + γ∥y(t) − W (t − 1)x(t)∥22 . log v̂i (t) + ε + 2 v̂ (t) + ε v̂ (t) + ε i j i=1 j=1 i=1

n X

j̸=i λ

At time t, the current statistics µ̂ (t), D̂ λ (t), and Ô λ (t) depend on y(t) through the steady-state updates in Equations B.3, B.6, and B.7. We first compute the necessary derivatives. Derivative of the centered activity.

For each component,

µ̂k (t) = λµ̂k (t − 1) + (1 − λ)yk (t). Hence ȳk (t) = yk (t) − µ̂k (t) = λ yk (t) − µ̂k (t − 1) , and therefore

∂ ȳi (t) = λ δik , ∂yk (t)

(B.9)

where δik is the Kronecker delta. Derivatives of the diagonal and off-diagonal statistics. Using Equation B.6 together with Equation B.9, we obtain ∂v̂i (t) = 2λ(1 − λ)ȳi (t) δik . (B.10) ∂yk (t) Similarly, from Equation B.7, ∂ĉij (t) = λ(1 − λ) δik ȳj (t) + δjk ȳi (t) , i ̸= j. (B.11) ∂yk (t) Derivative of the variance term. ∂ ∂yk (t)

Using Equation B.10, ! n n X X − log(v̂i (t) + ε) = −

1 ∂v̂i (t) v̂ (t) + ε ∂y k (t) i=1 i

i=1

=−

2λ(1 − λ)ȳk (t) . v̂k (t) + ε

(B.12)

Derivative of the cross-covariance penalty. Define n n 1 XX ĉij (t)2 . Rt := 2 i=1 j=1 v̂i (t) + ε v̂j (t) + ε j̸=i

Only terms with i = k or j = k contribute to the derivative with respect to yk (t). Using the symmetry of Ô λ (t), together with Equations B.10 and B.11, a direct calculation gives n X ∂Rt ĉkj (t)ȳj (t) = 2λ(1 − λ) ∂yk (t) v̂k (t) + ε v̂j (t) + ε j=1 j̸=k

− 2λ(1 − λ)

n X j=1 j̸=k

ĉkj (t)2 ȳk (t) 2 . v̂k (t) + ε v̂j (t) + ε

(B.13)

Derivative of the prediction term. Let ek (t) := [y(t) − W (t − 1)x(t)]k . Then ∂ γ∥y(t) − W (t − 1)x(t)∥22 = 2γek (t). ∂yk (t) Full gradient.

Combining Equations B.12, B.13, and B.14, we obtain 

(B.14)



n X   ĉkj (t)ȳj (t) ȳk (t) ∂Jt (y)  = 2λ(1 − λ)  − +  ∂yk (t) v̂k (t) + ε j=1 v̂k (t) + ε v̂j (t) + ε  j̸=k

− 2λ(1 − λ)

n X j=1 j̸=k

ĉkj (t)2 ȳk (t) 2 + 2γek (t). v̂k (t) + ε v̂j (t) + ε

(B.15)

The last term in Equation B.15 is quadratic in the off-diagonal covariance entries. Since the surrogate objective is used precisely in the regime where the normalized cross-covariances are small, this term is of higher order than the leading recurrent interaction retained in the main text. Discarding it yields the approximation   n X   ∂Jt (y) ĉkj (t)ȳj (t) ȳk (t)  ≈ 2λ(1 − λ)  − + + 2γek (t).  ∂yk (t) v̂k (t) + ε j=1 v̂k (t) + ε v̂j (t) + ε  j̸=k

Absorbing the constant prefactor λ(1 − λ) into γ, we recover the simplified gradient form used in the main text:   n X  ∂Jt (y)  ȳk (t) ĉkj (t)ȳj (t)  ∝ − + + γek (t).  ∂yk (t) v̂k (t) + ε j=1 v̂k (t) + ε v̂j (t) + ε  j̸=k

This is the descent direction used for the fast neural dynamics in Section 3. Descent property of the truncated direction. We now show that the truncated direction used in the fast neural dynamics is indeed a strict local descent direction for the exact objective under a simple pointwise condition. Define the regularized diagonal matrix D̂ λ,ε (t) := D̂ λ (t) + εI, let g(t) ∈ Rn denote the truncated direction used in the fast neural dynamics, and define −1/2 λ −1/2 B λ,ε (t) := D̂ λ,ε (t) Ô (t) D̂ λ,ε (t) . Further, let r(t) ∈ Rn be the vector with entries rk (t) :=

[B λ,ε (t)2 ]kk ȳk (t), v̂k (t) + ε 21

k = 1, . . . , n.

Proposition B.1 (Sufficient condition for descent of the truncated direction). With the notation above, the exact gradient can be written as ∇y Jt (y) = 2λ(1 − λ) g(t) − r(t) . (B.16) Consequently, −g(t) is a strict descent direction for Jt whenever ∥r(t)∥2 < ∥g(t)∥2 .

(B.17)

Moreover, ∥r(t)∥2 ≤

∥B λ,ε (t)∥22 (ε) v̂min (t)

(ε) v̂min (t) := min v̂k (t) + ε ,

∥ȳ(t)∥2 ,

(B.18)

1≤k≤n

so the simpler sufficient condition ∥B λ,ε (t)∥22 (ε)

v̂min (t)

∥ȳ(t)∥2 < ∥g(t)∥2

(B.19)

also guarantees that −g(t) is a strict descent direction. Proof. By construction, the k-th component of g(t) is n X ĉkj (t)ȳj (t) ȳk (t) + γek (t). + gk (t) = − v̂k (t) + ε j=1 v̂k (t) + ε v̂j (t) + ε j̸=k λ

Since Ô (t) has zero diagonal, [B λ,ε (t)2 ]kk =

n X

λ,ε Bkj (t)2 =

j=1

Therefore rk (t) =

n X j=1 j̸=k

ĉkj (t)2 . v̂k (t) + ε v̂j (t) + ε

n X ĉkj (t)2 ȳk (t) [B λ,ε (t)2 ]kk ȳk (t) = 2 . v̂k (t) + ε v̂j (t) + ε j=1 v̂k (t) + ε j̸=k

Comparing this identity with Equation B.15 shows that ∂Jt (y) = 2λ(1 − λ) gk (t) − rk (t) , ∂yk (t) which proves Equation B.16. Now let d(t) := −g(t). For any sufficiently small step size η > 0, Jt (y + ηd) = Jt (y) + η ⟨∇y Jt (y), d(t)⟩ + O(η 2 ) = Jt (y) − η ⟨∇y Jt (y), g(t)⟩ + O(η 2 ). Hence it is enough to show that ⟨∇y Jt (y), g(t)⟩ > 0. Using Equation equation B.16, ⟨∇y Jt (y), g(t)⟩ = 2λ(1 − λ) ∥g(t)∥22 − ⟨r(t), g(t)⟩ . By Cauchy–Schwarz, ⟨∇y Jt (y), g(t)⟩ ≥ 2λ(1 − λ) ∥g(t)∥2 ∥g(t)∥2 − ∥r(t)∥2 . Since 2λ(1 − λ) > 0, condition equation B.17 implies strict descent. Finally, ∥r(t)∥2 =

−1 D̂ λ,ε (t) diag B λ,ε (t)2 ȳ(t)

−1 D̂ λ,ε (t) diag B λ,ε (t)2

≤ 2

Moreover, for every k, [B λ,ε (t)2 ]kk ≤ ∥B λ,ε (t)2 ∥2 = ∥B λ,ε (t)∥22 , hence D̂ λ,ε (t)

−1

diag B λ,ε (t)2

≤

∥B λ,ε (t)∥22

. (ε) v̂min (t) This proves Equation B.18, and Equation equation B.19 is then immediate. 22

∥ȳ(t)∥2 .

Algorithm 1 Predictive Entropy Maximization: generic online procedure 1: Input: streaming mixtures {x(t)}T t=1 , source domain P 2: Hyperparameters: forgetting factor λ, feedforward step size schedule αW (t), neural step size

schedule ηy (t, τ ), fast-dynamics horizon τmax 3: Initialize W (0), µ̂λ (0), {v̂i (0)}n i=1 , {ĉij (0)}i̸=j , and θP (0) if needed 4: for t = 1, . . . , T do 5: Initialize y(t; 0) and θP (t; 0) if required by the domain 6: for τ = 0, . . . , τmax − 1 do 7: Compute g(t; τ ) from Equation 10 8: y(t; τ + 1) ← σP (y(t; τ ) − ηy (t, τ )g(t; τ ) ; θP (t; τ )) 9: Update θP (t; τ + 1) if required by the domain 10: end for 11: Set y(t) ← y(t; τmax ) and e(t) ← y(t) − W (t − 1)x(t) 12: W (t) ← W (t − 1) + αW (t)e(t)x(t)⊤ 13: µ̂λ (t) ← λµ̂λ (t − 1) + (1 − λ)y(t), ȳ(t) ← y(t) − µ̂λ (t) 14: v̂i (t) ← λv̂i (t − 1) + (1 − λ)ȳi (t)2 , ∀i 15: ĉij (t) ← λĉij (t − 1) + (1 − λ)ȳi (t)ȳj (t), ∀i ̸= j 16: end for

Domain-specific network realizations

The slow synaptic and statistical updates are identical across all source domains and are derived in Appendix B. The only domain-dependent part is the fast output inference step. All domains share the same descent direction, while the source domain determines the output nonlinearity and, when necessary, the dynamics of auxiliary inhibitory variables. For each sample t and neural-dynamics iteration τ , define n X ĉkj (t) ȳj (t; τ ) ȳk (t; τ ) + γek (t; τ ), gk (t; τ ) := − + v̂k (t) + ε j=1 v̂k (t) + ε v̂j (t) + ε j̸=k

where ȳk (t; τ ) = yk (t; τ ) − µ̂k (t),

ek (t; τ ) = yk (t; τ ) −

m X

Wkℓ (t − 1)xℓ (t).

ℓ=1

We further define the unconstrained pre-activation update ỹk (t; τ + 1) := yk (t; τ ) − ηy (t, τ ) gk (t; τ ). (C.1) The domain-specific architectures below are then obtained by combining Equation C.1 with the appropriate projection or proximal operator [35]. C.1

Antisparse sources

Consider the antisparse domain Bmax := {s ∈ Rn : ∥s∥∞ ≤ 1}. The output constraint is enforced directly through the clipping nonlinearity onto the box [−1, 1]n . The fast output dynamics are ỹk (t; τ + 1) = yk (t; τ ) − ηy (t, τ ) gk (t; τ ), yk (t; τ + 1) = clip[−1,1] ỹk (t; τ + 1) , k = 1, . . . , n, where  −1, u < −1, clip[−1,1] (u) := u, −1 ≤ u ≤ 1,  1, u > 1. Thus, the corresponding network consists only of the output population, driven by the feedforward estimate W (t)x(t), variance-normalized self-excitation, and covariance-dependent lateral inhibition. 23

C.2

Nonnegative antisparse sources

Consider the nonnegative antisparse domain Bmax,+ := Bmax ∩ Rn+ . The dynamics are identical to the antisparse case except that the projection is now onto the box [0, 1]n : ỹk (t; τ + 1) = yk (t; τ ) − ηy (t, τ ) gk (t; τ ), yk (t; τ + 1) = clip[0,1] ỹk (t; τ + 1) , k = 1, . . . , n, where  0, u < 0, clip[0,1] (u) := u, 0 ≤ u ≤ 1,  1, u > 1. C.3

Sparse sources

Consider the sparse domain B1 := {s ∈ Rn : ∥s∥1 ≤ 1}. We enforce the ℓ1 -constraint through a shared nonnegative variable λL , which acts as an adaptive threshold. Equivalently, one may view the fast-time update as a proximal-gradient step for the ℓ1 penalized objective associated with the Lagrangian Jt (y) + λL (∥y∥1 − 1),

λL ≥ 0.

Using the proximal mapping of the ℓ1 -norm [35], the resulting updates are ỹk (t; τ + 1) = yk (t; τ ) − ηy (t, τ ) gk (t; τ ), yk (t; τ + 1) = STλL (t;τ ) ỹk (t; τ + 1) , k = 1, . . . , n, λL (t; τ + 1) = ReLU λL (t; τ ) + ηλ (t, τ )

n X

!! |yi (t; τ + 1)| − 1

i=1

where STλ (u) := sign(u) max{|u| − λ, 0} is the soft-thresholding operator. The resulting architecture contains one additional inhibitory unit encoding λL , which delivers a common threshold to all output neurons. C.4

Nonnegative sparse sources

Consider the nonnegative sparse domain B1,+ := B1 ∩ Rn+ . Since y ≥ 0, the ℓ1 -constraint reduces to 1⊤ y ≤ 1. As in the sparse case, we introduce a shared nonnegative variable λL and interpret the fast-time step as a proximal / projected-gradient update for Jt (y) + λL (1⊤ y − 1), λL ≥ 0, followed by projection onto the nonnegative orthant [35]. The resulting updates are ỹk (t; τ + 1) = yk (t; τ ) − ηy (t, τ ) gk (t; τ ), yk (t; τ + 1) = ReLU ỹk (t; τ + 1) − λL (t; τ ) , λL (t; τ + 1) = ReLU λL (t; τ ) + ηλ (t, τ )

n X

k = 1, . . . , n, !! yi (t; τ + 1) − 1

i=1

Thus, the nonnegative sparse dynamics correspond to a thresholded ReLU nonlinearity driven by the shared inhibitory variable λL . 24

C.5

Simplex sources

Consider the simplex domain ∆ := {s ∈ Rn : s ≥ 0, 1⊤ s = 1}. The only difference with respect to the nonnegative sparse case is that the population constraint is now an equality rather than an inequality. Accordingly, the shared offset variable λL is no longer constrained to be nonnegative. The same proximal / projected interpretation then yields ỹk (t; τ + 1) = yk (t; τ ) − ηy (t, τ ) gk (t; τ ), yk (t; τ + 1) = ReLU ỹk (t; τ + 1) − λL (t; τ ) , λL (t; τ + 1) = λL (t; τ ) + ηλ (t, τ )

n X

k = 1, . . . , n, !

yi (t; τ + 1) − 1 .

i=1

Thus, the simplex network has the same output architecture as the nonnegative sparse network, except that the inhibitory unit corresponding to λL is linear rather than rectified. Summary. All five domains share the same feedforward pathway, the same covariance-driven recurrent interactions, and the same slow synaptic updates. What changes with the source domain is only the output-layer nonlinearity and, for B1 , B1,+ , and ∆, the inclusion of one shared inhibitory variable enforcing the corresponding population constraint. In particular, box-type domains require only local clipping operations, whereas sparse and simplex domains lead to piecewise-linear thresholding dynamics mediated by a single auxiliary inhibitory unit.

Error analysis and analytical bounds for the second-order Taylor approximation

In this section, we analyze the approximation error induced by the second-order surrogate in Equation 7. The analysis proceeds in two steps. We first study the Taylor remainder pointwise for a generic covariance matrix, deriving an exact spectral representation together with sharp spectral and norm-based bounds. We then specialize these results to the online covariance sequence used by the algorithm and, finally, use the same bounds to relate the batch surrogate objective to the exact regularized determinant objective. D.1

Pointwise remainder representation and bounds

We begin with a generic covariance matrix. Writing C = D + O, where D collects the diagonal entries of C and O its off-diagonal part, and introducing the regularized diagonal matrix D ε := D + εI, we define the normalized perturbation B ε := (D ε )−1/2 O(D ε )−1/2 . This yields an exact spectral representation of the Taylor remainder for log det(C + εI). Theorem D.1 (Exact spectral representation of the second-order remainder). Fix ε > 0. Let C ∈ Rn×n be a symmetric positive semidefinite matrix, and write C = D + O,

D = diag(C),

where D is the diagonal matrix formed from the diagonal of C and O is the off-diagonal part of C, so that Oii = 0 for all i. Define the regularized diagonal matrix D ε := D + εI, the normalized off-diagonal matrix B ε := (D ε )−1/2 O(D ε )−1/2 , and let λ1 , . . . , λn denote the eigenvalues of B ε . Then: 1. B ε is symmetric, Tr(B ε ) = 0, and λi > −1 for all i = 1, . . . , n; 25

2. the regularized log-determinant admits the exact decomposition n

n X

1X 2 log det(C + εI) = λ + R2 , log(Cii + ε) − 2 i=1 i i=1 where the remainder is given exactly by n X 1 2 R2 = log(1 + λi ) − λi + λi . 2 i=1

(D.1)

(D.2)

Moreover, the second-order term can be written entrywise as n X

λ2i = Tr((B ε )2 ) =

i=1

n X n X

2 Cij , (Cii + ε)(Cjj + ε) i=1 j=1

(D.3)

j̸=i

and therefore log det(C + εI) =

n X

log(Cii + ε) −

i=1

n n 2 Cij 1 XX + R2 . 2 i=1 j=1 (Cii + ε)(Cjj + ε)

(D.4)

j̸=i

Proof. Since C is symmetric positive semidefinite, each diagonal entry satisfies Cii ≥ 0. Hence D ε = D + εI is positive definite and (D ε )−1/2 is well-defined. Because O = C − D is symmetric and (D ε )−1/2 is diagonal, the matrix B ε = (D ε )−1/2 O(D ε )−1/2 is symmetric. Its diagonal entries are zero, since (B ε )ii = (Dii + ε)−1/2 Oii (Dii + ε)−1/2 = 0, Therefore Next, observe that

i = 1, . . . , n.

Tr(B ε ) = 0.

(D.5)

C + εI = D ε + O = (D ε )1/2 (I + B ε )(D ε )1/2 .

(D.6)

ε 1/2

Since C + εI is positive definite and (D ) ε

is invertible, Equation D.6 implies that

ε −1/2

I + B = (D )

(C + εI)(D ε )−1/2

is also positive definite. Hence all eigenvalues of I + B ε are strictly positive. If λ1 , . . . , λn are the eigenvalues of B ε , then the eigenvalues of I + B ε are 1 + λ1 , . . . , 1 + λn , and therefore λi > −1 for all i = 1, . . . , n. Using Equation D.6 and multiplicativity of the determinant, we obtain det(C + εI) = det(D ε ) det(I + B ε ). Taking logarithms yields log det(C + εI) = log det(D ε ) + log det(I + B ε ).

(D.7)

Since D is diagonal, log det(D ε ) =

n X

log(Cii + ε).

i=1

Now let

B ε = Q diag(λ1 , . . . , λn ) Q⊤ be an eigendecomposition of B ε , with Q orthogonal. Then I + B ε = Q diag(1 + λ1 , . . . , 1 + λn ) Q⊤ , 26

(D.8)

so log det(I + B ε ) =

n X

log(1 + λi ).

(D.9)

i=1

Substituting Equations D.8 and D.9 into Equation D.7 gives log det(C + εI) =

n X

log(Cii + ε) +

i=1

n X

log(1 + λi ).

(D.10)

i=1

P Pn P Using Equation D.5, we have i=1 λi = 0. Therefore we may add and subtract i λi − 12 i λ2i inside Equation D.10 to obtain n n n X 1X 2 X 1 2 log det(C + εI) = log(Cii + ε) − λ + log(1 + λi ) − λi + λi . 2 i=1 i i=1 2 i=1 This proves Equations D.1 and D.2. It remains to identify the quadratic term in coordinates. Since B ε is symmetric, n X

λ2i = Tr((B ε )2 ).

i=1 ε

Also, because (B )ii = 0,

n n X X (B ε )2ij . Tr((B ) ) = ε 2

i=1 j=1 j̸=i

Finally, for i ̸= j, Cij . (B ε )ij = (Dii + ε)−1/2 Oij (Djj + ε)−1/2 = p (Cii + ε)(Cjj + ε) Substituting this identity into the previous display proves Equation D.3, and Equation D.4 follows immediately. Theorem D.1 shows that the approximation error is determined by the spectrum of the regularized normalized perturbation B ε = (D ε )−1/2 O(D ε )−1/2 . This identifies B ε as the natural perturbation variable for assessing the accuracy of the regularized Taylor surrogate. The next result turns the exact representation above into a sharp two-sided bound on the remainder. Theorem D.2 (Two-sided spectral bound). Under the assumptions and notation of Theorem D.1, the remainder R2 satisfies 1 X |λi |3 1 X 3 − ≤ R2 ≤ λi . (D.11) 3 1 + λi 3 λi ≥0

λi <0

Consequently, |R2 | ≤ max

 1 X 3

λi ≥0

 3  X 1 |λi | λ3i , . 3 1 + λi 

(D.12)

λi <0

Proof. Define the scalar function 1 r(x) := log(1 + x) − x + x2 , 2

x > −1.

By Equation D.2, R2 =

n X

r(λi ).

(D.13)

i=1

Also, r(0) = 0,

r′ (x) = 27

x2 . 1+x

Hence

Z x r(x) = 0

t2 dt, 1+t

x > −1.

(D.14)

If x ≥ 0, then 1 + t ≥ 1 for all t ∈ [0, x], and Equation D.14 gives Z x 2 Z x t x3 0 ≤ r(x) = dt ≤ t2 dt = . 3 0 1+t 0 Thus 0 ≤ r(x) ≤

x3 , 3

x ≥ 0.

(D.15)

If −1 < x < 0, then Z 0 r(x) = − x

t2 dt ≤ 0. 1+t

Moreover, for t ∈ [x, 0], we have 1 + t ≥ 1 + x > 0, so 1 1 ≤ . 1+t 1+x Therefore

Z 0 |r(x)| = −r(x) = x

Hence −

1 t2 dt ≤ 1+t 1+x

|x|3 ≤ r(x) ≤ 0, 3(1 + x)

Z 0

t2 dt =

|x|3 . 3(1 + x)

−1 < x < 0.

(D.16)

Applying Equation D.15 to the nonnegative eigenvalues and Equation D.16 to the negative eigenvalues, and summing the resulting inequalities in Equation D.13, yields Equation D.11. Equation D.12 follows immediately. The preceding spectral bound immediately yields a simpler estimate expressed directly in terms of the size and conditioning of the regularized normalized perturbation B ε . Corollary D.3 (Simpler norm-based bound). Under the assumptions and notation of Theorem D.1, let λmin (B ε ) denote the smallest eigenvalue of the normalized off-diagonal matrix B ε . Then |R2 | ≤

n X ∥B ε ∥2F ∥B ε ∥2 3 . |λ | ≤ i 3 1 + λmin (B ε ) i=1 3 1 + λmin (B ε )

(D.17)

Proof. From Theorem D.2, we have |R2 | ≤ max

 1 X 3

λi ≥0

 3  X 1 |λi | λ3i , . 3 1 + λi 

(D.18)

λi <0

We bound the two terms on the right-hand side separately. First, since Tr(B ε ) = 0, either B ε = 0 or λmin (B ε ) ≤ 0. In either case, 1 + λmin (B ε ) ≤ 1. Therefore, for every eigenvalue with λi ≥ 0, λ3i ≤ Hence

λ3i 1 + λmin

(B ε )

≤

|λi |3 . 1 + λmin (B ε )

X 1 X 3 1 λi ≤ |λi |3 . ε 3 3 1 + λmin (B ) λi ≥0

λi ≥0

(D.19)

Second, for every eigenvalue with λi < 0, we have λi ≥ λmin (B ε ), so 1 + λi ≥ 1 + λmin (B ε ) > 0. Consequently, |λi |3 |λi |3 ≤ . 1 + λi 1 + λmin (B ε ) Summing over all negative eigenvalues gives X 1 X |λi |3 1 ≤ |λi |3 . 3 1 + λi 3 1 + λmin (B ε ) λi <0

(D.20)

λi <0

Combining Equations D.19 and D.20 with Equation D.18, we obtain |R2 | ≤

n X |λi |3 . 3 1 + λmin (B ε ) i=1

For the second inequality, note that n X

|λi | =

i=1

n X

|λi | |λi | ≤

i=1

X n λ2i = ∥B ε ∥2 ∥B ε ∥2F . max |λi |

1≤i≤n

i=1

Substituting this into the previous inequality proves Equation D.17. We now specialize the preceding pointwise analysis to the online covariance sequence generated by the algorithm. Corollary D.4 (Online version). Consider the online Predictive Entropy Maximization algorithm, and suppose that at time t the exponentially weighted output covariance Ĉ λ (t) is symmetric positive semidefinite. Let D̂ λ,ε (t) := D̂ λ (t) + εI,

Ĉ λ (t) = D̂ λ (t) + Ô λ (t),

and define the normalized off-diagonal matrix −1/2 λ −1/2 B̂ λ,ε (t) := D̂ λ,ε (t) Ô (t) D̂ λ,ε (t) . Let λ̂1 (t), . . . , λ̂n (t) denote the eigenvalues of B̂ λ,ε (t). Then the remainder R2 (t) of the secondorder surrogate in Equation 7 satisfies −

1 X |λ̂i (t)|3 1 X ≤ R2 (t) ≤ λ̂i (t)3 , 3 3 1 + λ̂i (t) λ̂i (t)<0

λ̂i (t)≥0

and therefore |R2 (t)| ≤ max

 1 X 3

λ̂i (t)≥0

 X |λ̂i (t)|3  1 λ̂i (t)3 , . 3 1 + λ̂i (t)  λ̂i (t)<0

Proof. Apply Theorem D.1 and Theorem D.2 with C = Ĉ λ (t),

D = D̂ λ (t),

D ε = D̂ λ,ε (t),

O = Ô λ (t).

Since Ĉ λ (t) is symmetric positive semidefinite and ε > 0, all assumptions are satisfied. Corollary D.4 is the bound used in our numerical diagnostics. It shows that the approximation error is controlled by the spectrum of the regularized normalized off-diagonal covariance matrix B̂ λ,ε (t). In particular, the surrogate is accurate when the diagonally normalized off-diagonal covariance is small and the eigenvalues of B̂ λ,ε (t) stay well away from the singular value −1. 29

D.2

Batch surrogate optimality relative to the exact determinant objective

The previous subsection controls the Taylor remainder pointwise for a fixed covariance matrix. We now use that control at the objective level. Although the surrogate in Section 3 is introduced through the online exponentially weighted covariance Ĉ λ (t), the expansion in Theorem D.1 and the remainder bounds above are pointwise statements about an arbitrary covariance matrix. It is therefore useful to consider the corresponding batch counterpart, obtained by replacing Ĉ λ (t) with the ordinary centered sample covariance. This auxiliary analysis does not address the convergence of the online dynamics; rather, it clarifies what minimizing the surrogate preserves from the exact regularized determinant objective itself. Let Y be a nonempty family of feasible output matrices Y = [y(1), . . . , y(T )] ∈ Rn×T such that y(t) ∈ P for all t, and such that the centered sample covariance Ĉy (Y ) :=

T ⊤ 1X y(t) − µ̂λ (Y ) y(t) − µ̂λ (Y ) , T t=1

µ̂λ (Y ) :=

1X y(t), T t=1

is symmetric positive semidefinite for every Y ∈ Y. For each Y ∈ Y, write D(Y ) := diag Ĉy (Y ) , D ε (Y ) := D(Y ) + εI, O(Y ) := Ĉy (Y ) − D(Y ), B ε (Y ) := D ε (Y )−1/2 O(Y )D ε (Y )−1/2 ,

v̂i (Y ) := [Ĉy (Y )]ii ,

ĉij (Y ) := [Ĉy (Y )]ij .

We consider the exact batch determinant objective batch Jdet (Y ) := − log det Ĉy (Y ) + εI ,

(D.21)

and its second-order batch surrogate batch Jsur (Y ) := −

n X

n X n 1X ĉij (Y )2 . log v̂i (Y ) + ε + 2 i=1 j=1 v̂i (Y ) + ε v̂j (Y ) + ε i=1 j̸=i

The next result shows that if the regularized normalized off-diagonal perturbation remains uniformly controlled over the feasible family Y, then every global minimizer of the surrogate is nearly optimal for the exact regularized determinant objective. Theorem D.5 (Uniform approximation and exact-objective near-optimality). Define ∥B ε (Y )∥2F ∥B ε (Y )∥2 . ε Y ∈Y 3 1 + λmin (B (Y ))

(D.22)

batch batch Jsur (Y ) − Jdet (Y ) ≤ ε̄Y ;

(D.23)

ε̄Y := sup Assume ε̄Y < ∞. Then: 1. for every Y ∈ Y,

batch 2. if Ysur is a global minimizer of Jsur over Y, then batch batch Jdet (Ysur ) ≤ inf Jdet (Y ) + 2ε̄Y . Y ∈Y

(D.24)

Proof. By Equation D.4, applied with C = Ĉy (Y ), the difference between the surrogate and exact batch objectives is exactly the Taylor remainder: batch batch Jsur (Y ) − Jdet (Y ) = R2 (Y ).

Applying Corollary D.3 with C = Ĉy (Y ) yields |R2 (Y )| ≤

∥B ε (Y )∥2F ∥B ε (Y )∥2 ≤ ε̄Y , 3 1 + λmin (B ε (Y ))

which proves Equation D.23. 30

batch Now let Ysur be a global minimizer of Jsur over Y. Then batch batch batch Jdet (Ysur ) ≤ Jsur (Ysur ) + ε̄Y = inf Jsur (Y ) + ε̄Y . Y ∈Y

Using Equation D.23 once more, batch inf Jsur (Y ) ≤ inf

Y ∈Y

batch batch Jdet (Y ) + ε̄Y = inf Jdet (Y ) + ε̄Y . Y ∈Y

Combining the last two displays proves Equation D.24. Theorem D.5 shows that the Taylor surrogate does more than approximate the exact regularized determinant objective pointwise: if the regularized normalized off-diagonal perturbation is uniformly small on the feasible family, then minimizing the surrogate yields an output matrix whose exact determinant objective value is close to the batch optimum. Corollary D.6 (Explicit spectral form). Assume that there exists ρ ∈ (0, 1) such that ∥B ε (Y )∥2 ≤ ρ, Then ε̄Y ≤

∀ Y ∈ Y.

(D.25)

nρ3 . 3(1 − ρ)

batch Consequently, every global minimizer Ysur of Jsur over Y satisfies batch batch Jdet (Ysur ) ≤ inf Jdet (Y ) + Y ∈Y

2nρ3 . 3(1 − ρ)

(D.26)

Proof. Under Equation D.25, we have λmin (B ε (Y )) ≥ −ρ,

∥B ε (Y )∥2F ≤ nρ2 ,

∀ Y ∈ Y.

Substituting these inequalities into Equation D.22 gives ε̄Y ≤

nρ2 · ρ nρ3 = . 3(1 − ρ) 3(1 − ρ)

The bound in Equation D.26 then follows directly from Theorem D.5.

Supplementary on numerical experiments

This section provides technical implementation details, extended experimental results, and additional information on the experimental protocols used throughout the paper. E.1

Performance evaluation metric

To evaluate source recovery performance, we report the mean signal-to-noise ratio (mSNR) of the estimated sources. Let si and ỹi denote the i-th ground-truth source and its corresponding permutation- and sign-corrected estimate, respectively. The mSNR is defined as n 1X ∥si ∥22 mSNR (dB) = 10 log10 . (E.1) n i=1 ∥si − ỹi ∥22 To assess statistical variability, all reported metrics are aggregated over N = 30 independent realizations with different random mixing matrices and noise seeds. We report the sample mean µ̄ together √ with the 95% confidence interval. Let σ denote the sample standard deviation and SE = σ/ N the standard error of the mean. The confidence interval is computed as σ CI95% = µ̄ ± tα/2,N −1 √ , (E.2) N where tα/2,N −1 is the critical value of the Student t-distribution with α = 0.05 and N − 1 degrees of freedom. For N = 30, this gives t0.025,29 ≈ 2.045. All shaded envelopes in the reported plots correspond to this 95% confidence interval. 31

Channel 1 Ground Truth Estimation

0.10 0.05 0.00 0.05 0.10 Channel 2

0.3

Ground Truth Estimation

0.2 0.1 0.0 0.1 0.2 Channel 3

0.08

Ground Truth Estimation

0.06 0.04 0.02 0.00 0.02 0.04 0

200

400

600 Time Samples

800

1000

Figure 5: Representative auditory source-separation result. Temporal alignment between the ground-truth and recovered sources for one trial in the cocktail-party experiment described in Section 4. E.2

Transform invariance in auditory source separation

In Section 4, we applied Predictive Entropy Maximization to audio mixtures in a sparse wavelet domain. The justification is straightforward. In the noise-free setting, let X = AS,

X ∈ Rm×T ,

S ∈ Rn×T ,

and let Φ ∈ RT ×T denote a linear transform acting on the sample axis, such as the discrete wavelet transform. The transformed mixtures satisfy X̃ = XΦ = (AS)Φ = A(SΦ) = AS̃, where S̃ := SΦ. Thus, the same mixing matrix A governs the linear relation in the transform domain. This preserves the mixing geometry while allowing us to exploit the approximate sparsity of wavelet coefficients. Accordingly, we learn a separator W such that W X̃ ≈ S̃ using the B1 architecture. After separation, the time-domain sources are recovered through the inverse transform: Ŝ = (W X̃)Φ−1 ≈ S̃Φ−1 = S. For the specific trial illustrated in Figure 5, the mixing matrix was   −1.131 0.696 −0.432  0.741 −0.478 1.386    1.149 −2.350 , A =  0.125  0.183 −0.311 −0.294 0.400 1.006 0.502 and the resulting source-wise SNR values were 25.76 dB, 24.36 dB, and 30.54 dB. E.3

Antisparse, nonnegative sparse, and simplex source separation examples

Complementing the representative results in Figure 2, Figure 6 reports three additional synthetic domains. 32

(a) Correlated antisparse (Bmax )

(b) Noisy nonneg. sparse (B1,+ )

Figure 6: Additional performance comparisons. (a) Mean component SNR (mSNR) as a function of the source correlation level ρ for antisparse sources. (b) Mean component SNR as a function of the input SNR for nonnegative sparse sources. (c) Mean component SNR as a function of the input SNR for simplex sources. These complementary experiments confirm the same pattern observed in the main text: Predictive Entropy Maximization remains robust under source correlation in antisparse domains and remains competitive with batch Det-Max baselines under increasing observation noise in nonnegative sparse- and simplex-structured domains. Shaded envelopes show 95% confidence intervals over 30 independent realizations.

For the antisparse experiment, we use the domain Bmax . As in the nonnegative antisparse case from the main text, sources are generated from a Copula-t model with four degrees of freedom, the source correlation level varies over ρ ∈ {0, 0.05, . . . , 0.5}, and the input SNR is fixed at 30 dB. Figure 6a shows that both PEM and u-PEM remain robust as correlation increases. In particular, they maintain strong performance relative to the biologically plausible baselines, while the batch LD-InfoMax method again provides the strongest overall reference performance. While u-PEM performs similarly to PEM, its performance declines slightly faster with increasing source correlation. For the nonnegative sparse experiment, we use the domain B1,+ and vary the observation noise over SNRin ∈ {30, 25, . . . , 5} dB. Sources are generated uniformly from B1,+ . Figure 6b shows that PEM and u-PEM remain competitive with the batch LD-InfoMax baseline across the full noise range, confirming that the surrogate-based online formulation extends naturally to nonnegative sparse domains as well. For the simplex experiment, we use the domain ∆ and vary the observation noise over SNRin ∈ {30, 25, . . . , 5} dB. Sources are generated uniformly from ∆. Figure 6c shows that PEM remains competitive with the batch LD-InfoMax baseline across the full noise range, confirming that the surrogate-based online formulation extends beyond box-constrained and sparse domains to simplexstructured sources as well. u-PEM also follows the other models closely, although it performs noticeably worse at high SNR values. E.4

Transient diagnostics for the Taylor approximation

Figure 7 reports the time evolution of the exact Taylor remainder and the corresponding spectral upper bound for two representative correlation levels, ρ = 0 and ρ = 0.4. In both cases, the approximation error decreases during training as the dynamics suppress normalized off-diagonal covariance structure. The correlated case exhibits a larger error scale, as expected, but remains controlled by the theoretical bound throughout training. These trajectories complement the aggregate Taylor-surrogate diagnostics reported in the main text. E.5

Simulation hyperparameters

This section summarizes the hyperparameters used in all reported simulations. We first describe the generic implementation choices for Predictive Entropy Maximization, and then list the experimentspecific settings. Generic implementation details for Predictive Entropy Maximization. All Predictive Entropy Maximization experiments were run in the online setting with a single streaming pass over the data. Unless otherwise stated, we used the regularization constant ε = 10−5 . For each incoming sample x(t), the fast inference dynamics were initialized at y(t; 0) = 0. If no custom initialization was 33

Time-Evolution of Mean Surrogate Error ( = 0.0) Error Magnitude (Log Scale)

102

Magnitude (Log Scale)

Time-Evolution of Mean Surrogate Error ( = 0.4)

Mean Actual Error Mean Theoretical Bound

101 100 10 1 10 2 10 3

5000

Mean Actual Error Mean Theoretical Bound

102 101 100 10 1 0

10000 15000 20000 25000 30000 35000 40000

Training Iterations

(a) Uncorrelated sources (ρ = 0)

5000 10000 15000 20000 25000 30000 35000 40000

Training Iterations

(b) Correlated sources (ρ = 0.4)

Figure 7: Transient Taylor-surrogate diagnostics. Exact Taylor remainder and corresponding spectral upper bound over training iterations for two representative source-correlation levels. supplied, the feedforward matrix was initialized as W (0) = I + ξW ,

(E.3)

where I denotes the rectangular identity and ξW has i.i.d. Gaussian entries with standard deviation 0.01. Likewise, if not specified explicitly, the running mean and covariance states were initialized as µ̂y (0) = 0, Ĉy (0) = 0.2 I. (E.4) The slow feedforward update uses the step size schedule selected by lr_W_rule. In the experiments below, we used either: 0 αW (t) = αW (constant), (E.5) 0 αW , 10−8 (divide_by_index), (E.6) αW (t) = max t/TW + 1 or 0 αW −8 αW (t) = max , 10 (divide_by_log_index), (E.7) 1 + log(t/TW + 2) where TW is the decay divider. The fast neural inference loop uses the step size schedule selected by neural_lr_rule. In the experiments below, we used either: ηy (τ ) = ηy0 ) ( ηy0 min ηy (τ ) = max ,η τ +1 y or

( ηy (τ ) = max

ηy0 , η min τ Ty + 1 y

(constant),

(E.8)

(divide_by_loop_index),

(E.9)

(divide_by_slow_loop_index),

(E.10)

)

where Ty is the neural learning-rate decay divider. For sparse, nonnegative sparse, and simplex domains, the scalar threshold/inhibitory variable λL was updated with step size ηλ , implemented in code as stlambda_lr. Predictive Entropy Maximization hyperparameters. The synthetic source-separation experiments use n = 5 sources, m = 10 mixtures, and T = 105 samples. The auditory experiment uses n = 3 sources and m = 5 mixtures. The sparse receptive-field experiment uses n = m = 144. 0 • Correlated antisparse sources (Bmax ). We used λlat = 0.99, γpred = 250, αW = 5 × 10−2 , 0 min −6 −7 ηy = 0.5, ηy = 10 , τmax = 250, and tolerance 10 . The feedforward learning rate used divide_by_index with divider TW = 5000, and the neural step size used divide_by_loop_index. No threshold variable was used.

• Correlated nonnegative antisparse sources (Bmax,+ ). We used λlat = 0.95, γpred = 750, 0 αW = 5 × 10−2 , ηy0 = 0.05, ηymin = 10−4 , τmax = 500, ε = 10−4 and tolerance 10−6 . The feedforward learning rate used divide_by_index with divider TW = 20000, and the neural step size used divide_by_loop_index. In this experiment we overrode the default initialization and set 1 Ĉy (0) = 2I, W (0) = 0.01 I + ξW , (E.11) 15 with ξW i.i.d. Gaussian. 0 • Noisy sparse sources (B1 ). We used λlat = 0.99, γpred = 150, αW = 5 × 10−2 , 0 min −4 −6 ηy = 0.05, ηy = 10 , ηλ = 0.5, τmax = 100, and tolerance 10 . The feedforward learning rate used divide_by_index with divider TW = 5000, and the neural step size used divide_by_loop_index. 0 • Noisy nonnegative sparse sources (B1,+ ). We used λlat = 0.99, γpred = 250, αW = 5× 10−2 , ηy0 = 0.1, ηymin = 10−4 , ηλ = 0.5, τmax = 100, and tolerance 10−7 . The feedforward learning rate used divide_by_index with divider TW = 2000, and the neural step size used divide_by_loop_index. 0 • Noisy simplex sources (∆). We used λlat = 0.99, γpred = 150, αW = 5 × 10−2 , ηy0 = 0.1, ηymin = 10−4 , ηλ = 0.05, τmax = 100, and tolerance 10−7 . The feedforward learning rate used divide_by_log_index with divider TW = 5000, and the neural step size used divide_by_loop_index. • Auditory source separation. The sparse-domain audio experiment was run after transforming the mixtures to a wavelet domain using a db4 wavelet with decomposition level 3. The Predictive 0 = 9.5 × 10−1 , Entropy Maximization hyperparameters were λlat = 0.95, γpred = 150, αW −4 −6 min 0 = 10 , ηλ = 0.5, τmax = 100, and tolerance 10 . The feedforward ηy = 0.01, ηy learning rate used divide_by_index with divider TW = 2000, and the neural step size used divide_by_loop_index. • Sparse receptive-field learning on natural images. For the 12 × 12 patch experiment we used the sparse-domain architecture with λlat = 1 −

10−3 , 7

γpred = 3,

0 αW = 10−4 ,

(E.12)

together with ηy0 = 0.05, ηymin = 10−6 , ηλ = 5 × 10−2 , τmax = 500, and tolerance 10−6 . The feedforward learning rate was kept constant, while the neural step size used divide_by_slow_loop_index with divider Ty = 100. We also used custom initializations Ĉy (0) = I +

1 ξC , 250

W (0) = I +

1 ξW , 250

(E.13)

where ξC and ξW have i.i.d. Gaussian entries. E.6

Ablation Studies

We report two additional ablations to better understand the robustness of Predictive Entropy Maximization. The first studies the sensitivity of the method to the distribution of the mixing-matrix entries. The second studies how performance varies with the number of mixtures while keeping the number of sources fixed. E.6.1

Sensitivity to the mixing-matrix distribution

We first study whether the method is sensitive to the particular law used to generate the mixing matrix. For this ablation, we consider the antisparse and nonnegative antisparse domains, fix the source correlation level at ρ = 0, and fix the input SNR at 30 dB. In both cases, we use n = 5 sources, m = 10 mixtures, and T = 105 samples. Results are aggregated over 30 random seeds. For each seed, we keep the source realization and the observation noise fixed and vary only the distribution of the entries of the mixing matrix A. We consider five centered distributions, all scaled to have unit variance: p √ √ √ N (0, 1), U(− 3, 3), Laplace(0, 1/ 2), Rad(±1), 3/5 t5 . 35

Nonnegative antisparse ablation: effect of mixing-matrix distribution 28

mSNR (dB)

Antisparse ablation: effect of mixing-matrix distribution 28

22 20

24 22 20

18 16

18 (0, 1)

(

3, 3)

Laplace(0, 1/ 2)

Rad(±1)

Mixing-matrix entry distribution

(0, 1)

3/5 t5

(a) Antisparse sources (Bmax )

(

3, 3)

Laplace(0, 1/ 2)

Rad(±1)

Mixing-matrix entry distribution

3/5 t5

(b) Nonnegative antisparse sources (Bmax,+ )

Figure 8: Ablation with respect to the distribution of the mixing-matrix entries. Box plots summarize the distribution of mSNR over 30 seeds in the uncorrelated setting (ρ = 0) at input SNR √ √ 30 dB. All candidate mixing laws are centered and scaled to unit variance: N (0, 1), U(− 3, 3), p √ Laplace(0, 1/ 2), Rad(±1), and 3/5 t5 . The performance of Predictive Entropy Maximization remains broadly stable across these choices, indicating that the method is not sensitive to the precise distributional shape of the mixing coefficients.

Here, Rad(±1) denotes a Rademacher random variable taking values ±1 with equal probability, and p t5 denotes a Student-t random variable with 5 degrees of freedom. The scaling by 3/5 ensures unit variance, so that the comparison isolates the effect of distributional shape rather than trivial scale differences. Figure 8 shows the resulting box plots for the mean component SNR. Across both domains, the central performance remains broadly stable across the five mixing laws. This indicates that the behavior of Predictive Entropy Maximization is not tied to Gaussian mixing matrices and transfers well to heavier-tailed, bounded, and discrete mixing distributions. The nonnegative antisparse case exhibits slightly larger variability across seeds, but again no systematic degradation under nonGaussian or discrete mixing is observed. E.6.2

Effect of the number of mixtures

We next study how performance varies with the number of observed mixtures. For this ablation, we consider the sparse and simplex domains, fix the number of sources at n = 5, fix the input SNR at 30 dB, and vary the number of mixtures over m ∈ {7, 8, 9, 10, 11, 12, 13}. As before, each experiment uses T = 105 samples and is repeated over 30 random seeds. For a fixed seed, we keep the source realization fixed and generate a single Gaussian mixing matrix of size 13×5. The experiment with m mixtures then uses the first m rows of this matrix. This nested construction makes the comparison more controlled, since increasing the number of mixtures corresponds to adding new observation channels rather than resampling a completely different mixing matrix. Figure 9 reports the mean mSNR together with a standard-error envelope across seeds. In both domains, performance improves as the number of mixtures increases. This is consistent with the intuition that additional observation channels make the inverse problem easier by providing richer linear views of the same latent sources. The trend is visible already in the sparse case and is even more pronounced in the simplex setting. These results confirm that the method scales favorably with observation dimension in the overdetermined regime. E.7

Unnormalized Predictive Entropy Maximization

The Predictive Entropy Maximization (PEM) model includes cross-covariance terms that are variance normalized in both the energy function and the activity gradient update (Eqs. 9-10). If we 36

Sparse ablation: effect of the number of mixtures 26.75

29.5

26.50

29.0

mSNR (dB)

Simplex ablation: effect of the number of mixtures

30.0

27.00

26.25 26.00 25.75 25.50

28.5 28.0 27.5

25.25

27.0 7

Number of mixtures

(a) Sparse sources (B1 )

Number of mixtures

(b) Simplex sources (∆)

Figure 9: Ablation with respect to the number of mixtures. Mean mSNR is plotted against the number of mixtures m ∈ {7, . . . , 13} while keeping the number of sources fixed at n = 5. The shaded envelopes show one standard error over 30 seeds. For each seed, a single Gaussian 13 × 5 mixing matrix is sampled and the experiment with m mixtures uses its first m rows, so the comparison isolates the effect of adding observation channels. In both the sparse and simplex domains, performance improves as the number of mixtures increases. remove this normalization term, the objective becomes Jt = −

n X

n X n X 1 log(v̂i (t) + ε) + γlateral ĉij (t)2 + γ∥y(t) − W (t − 1)x(t)∥22 , 2 i=1 i=1 j=1

(E.14)

j̸=i

where we also introduced an additional hyperparameter γlateral that controls the strength of the covariance term, corresponding to lateral connections. We call this unnormalized variant of the model u-PEM. Furthermore, the coordinate gradient update becomes (cf. Eq. 10) ! n m X X ȳk (t, τ ) dk (t, τ ) = −γlateral ĉkj (t) ȳj (t, τ )−γ yk (t, τ ) − Wkℓ (t − 1)xℓ (t) . (E.15) v̂k (t) + ε j=1 ℓ=1

j̸=k

The update is now directly given by the product of the lateral weight and the neural activity, with no variance scaling, making u-PEM arguably even more biologically plausible than PEM. The removal of the variance normalization is further motivated by the observation that the mean variance of the PEM model appears to stabilise during training, as shown in Figure 10. Benchmark results for u-PEM and PEM are reported in Figures 2 and 6. The hyperparameter values γlateral for the antisparse, non-negative antisparse, sparse, non-negative sparse and simplex domains were 10, 300, 50, 3200 and 100, respectively. E.8

Variance-covariance regularisation in self-supervised learning

An influential approach to self-supervised learning is to learn embeddings or latent representations of inputs (e.g. images) which are invariant to various transformations (e.g. crops). However, selfsupervised objectives are prone to learning constant or degenerate embeddings (e.g. not using all the available dimensions), a phenomenon known as “representation collapse” [42, 43, 44]. Two main methods have been developed to prevent or mitigate against such collapse: (i) contrastive methods [45], which involve training with both positive and negative examples; and (ii) non-contrastive or regularized approaches, which directly modify the loss function and which are most closely related to the Predictive Entropy Maximization model. For example, the objective of “Barlow Twins” [43] pushes the cross-correlation matrix between the embeddings of two identical networks towards the identity. This promotes both invariance to distortions (via the diagonal terms) and decorrelation (via the off-diagonal terms). The method was named after the neuroscientist H. Barlow for introducing the idea of redundancy-reduction, hypothesizing that the visual system minimizes the statistical dependency between neural firing rates [46] (which relates to our experiment on sparse receptive fields; see Figure 3). 37

(a) Antisparse sources (Bmax )

(b) Sparse sources (B1 )

Figure 10: Evolution of the variance term in PEM during training. The evolution of the mean variance (averaged across the sources) is plotted across training samples for the PEM model. Solid lines and shaded bands indicate the mean and the 95% confidence interval over 30 seeds, respectively. For each seed, a single Gaussian 10 × 5 mixing matrix is sampled. These experiments demonstrate that the variance terms stabilise to a small range of values, motivating our simplification further (as well as providing a starting point for the value of γlateral ). Building on the Barlow Twins, “VICReg” [34] introduced a variance term that encourages spread in each embedding dimension. The overall loss function therefore contains a variance, a covariance and an invariance term. This is similar to the objective of the Predictive Entropy Maximization model, which also features variance-maximization and covariance-minimization terms under structural constraints (see Eq. 9). Although the models used for self-supervised learning are not biologically plausible and address a different problem than source separation, they provide a useful point of comparison given the similarity of the objectives (and perhaps a common goal of encouraging sufficient spread in the domain). A form of variance-covariance regularization has also been used in recent large-scale models such as Video JEPA [47]. Overall, extending Predictive Entropy Maximization to self-supervised or supervised learning tasks could be an interesting future direction. E.9

Computational complexity of the proposed method

In this section, we characterize the computational cost of Predictive Entropy Maximization. Let n denote the number of sources, m the number of mixtures, and τmax the maximum number of fast neural-dynamics iterations used to compute the output for a single input sample. The dominant cost comes from the fast inference step. For one neural-dynamics iteration, the feedforward prediction W (t)x(t) requires O(nm) operations, while the covariance-driven lateral term in Equation 10 requires O(n2 ) operations. The domain-dependent output nonlinearity adds only O(n) work, and for B1 , B1,+ , and ∆, the additional update of the shared inhibitory variable λL also remains O(n). Therefore, one fast iteration costs O(nm + n2 ), and the full inference stage costs O τmax (nm + n2 ) . After the output has converged, the slow updates are cheaper. The feedforward update W (t) = W (t − 1) + αW (t)e(t)x(t)⊤ costs O(nm), the mean update costs O(n), the diagonal variance update costs O(n), and the offdiagonal covariance update costs O(n2 ). Hence the slow plasticity and statistics updates contribute O(nm + n2 ) per sample. 38

Combining both stages, the overall worst-case cost per sample is O (τmax + 1)(nm + n2 ) . In the determined or overdetermined regime m ≥ n, this simplifies to O(τmax mn), since the recurrent inference stage dominates. Thus, as in other biologically plausible recurrent BSS algorithms, the computational bottleneck is the iterative output inference rather than the synaptic updates. This is the same asymptotic order reported for the antisparse CorInfoMax network and for the determinant-maximization WSM networks, both of which are likewise dominated by recurrent output computation in digital simulations. Compute resources. The implementation is CPU-compatible and does not require GPU acceleration or high-performance computing; the reported experiments can be run on a standard personal computer with sufficient memory. The runs reported here were performed with Python 3.12.12 on a Linux x86_64 machine (Linux 6.8.0-107-generic, glibc 2.39) using a single-node, single-task CPU Slurm job with 10GB requested memory.

Record · ID 204801 · SHA-256 86a53a3b578f66e9
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.