Foundations of Independent Component Analysis

Foundations of Independent Component Analysis Patrick Forré AI4Science Lab Korteweg-de Vries Institute for Mathematics University of Amsterdam p.d.forre@uva.nl

arXiv:2608.13229v1 [math.ST] 13 Aug 2026

Abstract We present the mathematical foundations of linear independent component analysis (ICA) models based on standard literature in a self-contained note. It is aimed at readers with a background in measure-theoretic probability theory. We first develop the theory of the characteristic functions of probability measures on Rd , including their analyticity and the way in which they determine and characterise the distributions. We then focus on several identifiability results of ICA models with successively strengthened assumptions on the sources: from merely non-constant, to non-Gaussian, to Gaussian-free independent sources. Under the strictest assumptions, we show that the independent sources are identifiable up to translation, permutation, scales and signs, and this even in the presence of additive Gaussian noise. Furthermore, we present the online equivariant gradient descent ICA algorithm for recovering the independent sources from data, in the standard complete noiseless non-Gaussian ICA setting.

Contents 1 Introduction

2

2 Notation and Conventions

4

3 Characteristic Functions of Probability Measures on Rd 3.1 Definition and elementary properties . . . . . . . . . . . . . . . . . . . . . . . 3.2 How characteristic functions determine distributions . . . . . . . . . . . . . . 3.3 Moments, the distinguished logarithm, and cumulants . . . . . . . . . . . . . 3.4 Analytic characteristic functions . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Gaussian, sub-Gaussian and super-Gaussian distributions . . . . . . . . . . . 3.6 Why non-Gaussianity is what makes ICA work . . . . . . . . . . . . . . . . .

5 5 6 8 10 12 15

4 Identifiability for Non-Constant Independent Sources

17

5 Identifiability for Non-Gaussian Independent Sources

20

6 Identifiability for Gaussian-Free Independent Sources 6.1 Gaussian-free random variables . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Identifiability in the presence of additive Gaussian noise . . . . . . . . . . . .

27 27 31

7 Complete Noiseless ICA: Estimation by Equivariant Gradient Descent 7.1 The complete noiseless model and its identifiability . . . . . . . . . . . . . . . 7.2 The maximum likelihood objective . . . . . . . . . . . . . . . . . . . . . . . . 7.3 The relative gradient: what the “preconditioner” really is . . . . . . . . . . . 7.4 Stationary points, stability, and the choice of nonlinearity . . . . . . . . . . .

35 35 38 40 42

1

7.5 The algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.6 LiNGAM: a causal order removes the permutation . . . . . . . . . . . . . . . 7.7 Generalisations: a brief overview . . . . . . . . . . . . . . . . . . . . . . . . .

49 49 52

A Proofs of the Classical Results on Characteristic Functions A.1 Uniqueness and inversion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Lévy’s continuity theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3 Entire characteristic functions . . . . . . . . . . . . . . . . . . . . . . . . . . . A.4 Marcinkiewicz’ theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.5 Cramér’s decomposition theorem . . . . . . . . . . . . . . . . . . . . . . . . .

53 54 55 56 57 58

B A Proof of the Kagan–Linnik–Rao Theorem B.1 Finite differences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2 Proof of the column dichotomy . . . . . . . . . . . . . . . . . . . . . . . . . .

59 60 62

References

65

1

Introduction

The linear independent component analysis (ICA) model postulates that an observed pdimensional random vector X is an affine image of a k-dimensional random vector Z with mutually independent components: X = AZ + µ,

A ∈ Rp×k ,

µ ∈ Rp .

(1)

The components Z1 , . . . , Zk are called the sources and A is called the mixing matrix; both are unobserved. The identifiability question asks how much of the pair (A, L(Z)) is determined by the law of X alone (Comon, 1994; Hyvärinen and Oja, 2000; Eriksson and Koivunen, 2004). Some ambiguity is unavoidable and built into the model: we may rescale a column of A and inversely rescale the corresponding source, we may permute the columns of A together with the sources, and we may shift a source into the offset µ. Written out, for any permutation matrix P , any invertible diagonal matrix Λ and any c ∈ Rk , AZ + µ = (AP Λ) Λ−1 P −1 (Z − c) + (Ac + µ),

(2)

and the right hand side is again a representation of the form Eq. (1) with independent sources. The content of the classical identifiability theory is that, under suitable assumptions, this is the only ambiguity – with one important exception: Gaussian sources. That exception is genuine and not an artefact of the proofs. If Z ∼ N (0, Ik ) then QZ ∼ N (0, Ik ) for every orthogonal Q, so a purely Gaussian source vector can be rotated arbitrarily without changing the law of X. The results collected here go back to the characterisation theory of the Gaussian law that grew out of the theorems of Cramér (1936), Marcinkiewicz (1939), Darmois (1953) and Skitovich (1954), and that is developed systematically in the monographs of Kagan et al. (1973) and Linnik and Ostrovskii (1977). The ICA community rediscovered and popularised these results, notably through Comon (1994) and Eriksson and Koivunen (2004), and they underpin causal discovery methods such as LiNGAM (Shimizu et al., 2006). Prerequisites. These notes are written for readers who have taken a measure-theoretic probability course. We assume laws and weak convergence, product measures and Fubini’s theorem, and basic complex analysis up to Morera’s theorem and the identity theorem. Everything else – in particular the entire theory of characteristic functions that the arguments rest on – is developed in Section 3, with precise references for the standard results and full proofs for everything specific to our setting. 2

Outline and contributions. Section 2 fixes notation. Section 3 develops characteristic functions of probability measures on Rd : their elementary properties, how they determine and characterise distributions (uniqueness, inversion, the Cramér–Wold device, Lévy’s continuity theorem), why local agreement of two characteristic functions is not enough, moments and cumulants via the distinguished logarithm, analyticity in a strip and the resulting rigidity, and finally the characterisation of the Gaussian law together with the sub- and super-Gaussian (platykurtic and leptokurtic) distributions that ICA is built on. It closes with the precise sense in which kurtosis acts as a contrast function, and with the observation that the identifiability theory needs strictly less than that: non-normality, not non-zero kurtosis. Section 4 treats sources that are merely assumed to be non-constant and states the fundamental result of Kagan et al. (1973) in the form we need; we take the opportunity to reinstate a proportionality constant that is lost in the statement given there (see Footnote 1). Section 5 specialises to non-Gaussian sources. The main technical device is Lemma 5.2, which shows how a Gaussian noise vector can be traded against extra columns of the mixing matrix, and back; this is what allows us to reduce the non-Gaussian case to the non-constant case. Section 6 introduces the stronger requirement that the sources be Gaussian-free, i.e. that no independent Gaussian noise can be split off them at all. We show in Theorem 6.7 that every real-valued random variable decomposes, essentially uniquely, into a Gaussian-free part and independent Gaussian noise, and we then prove the strongest statement of these notes, Theorem 6.12: for full column rank mixing matrices the sources are then identifiable up to permutation, scale and translation even in the presence of additive Gaussian noise with an arbitrary, possibly non-diagonal and degenerate covariance matrix. Section 7 turns to estimation in the case that dominates applications – a square invertible mixing matrix and no noise. We derive the classical identifiability statement as Corollary 7.4, set up the maximum likelihood objective, show that the customary “preconditioner” in the gradient step is not an approximate inverse Hessian but the exact gradient for a natural right-invariant metric, and prove in Theorem 7.20 exactly when a separating solution is a locally asymptotically stable equilibrium of the mean dynamics underlying the online algorithm – which, for correctly specified sources, is precisely when Corollary 7.4 says the model is identifiable. The same section then shows in Corollary 7.32 how a causal ordering removes the last remaining ambiguities, giving the linear non-Gaussian acyclic model LiNGAM, and closes with a short survey of how the model has been generalised. Section A collects the proofs of the four classical results about characteristic functions that Section 3 quotes: the uniqueness and inversion theorem, Lévy’s continuity theorem, Marcinkiewicz’ theorem and Cramér’s decomposition theorem. Section B proves Theorem 4.2, the identifiability theorem of Kagan et al. (1973) on which Section 4 onwards rest. The proof is a self-contained finite-difference argument: it needs only the distinguished logarithm, Marcinkiewicz’ theorem, and Fréchet’s functional equation, which we also prove. With these appendices in place, every result in these notes is proved here, with one exception: Bochner’s Theorem 3.3, which we quote and which is used only once, in Remark 3.8, where it could be replaced by a direct computation. Standard background is of course used throughout and is listed at the start of Section A; Section 7 additionally uses standard matrix calculus and, in Theorem 7.20, the unstable-manifold theorem for maps. How to read this note. The logical skeleton is short. Everything downstream of Section 4 rests on the single result Theorem 4.2, proved in Section B, and branches from there into the two lines of development that the rest of the note pursues:

3

Theorem 5.5

Theorem 6.12

non-Gaussian sources

Gaussian-free sources

Theorem 4.2 non-constant sources

Theorem 7.20 algorithmic stability

Corollary 7.4 complete noiseless model

Corollary 7.32 LiNGAM

The upper branch strengthens the hypotheses on the sources and asks how much of the model this buys back; the lower branch fixes the hypotheses at their most convenient and asks how to estimate. A reader interested only in the algorithm may go directly from Section 3 to Section 7, taking Corollary 7.4 on faith.

2

Notation and Conventions

Notation 2.1 (Random objects versus deterministic objects). All random objects are defined on one underlying probability space (Ω, F, P) and are denoted by capital Latin letters: X for the observed random vector, Z for the source vector, E, V , G, H, N , U for the various (mostly Gaussian) noise and remainder terms. Deterministic matrices are capital Latin or Greek letters (A, B, P , Q, R, W , Θ, Λ, Σ; in particular W always denotes the deterministic unmixing matrix of Section 7, never a random vector), deterministic vectors and scalars are lower case (aj , bl , µ, ν, ξ, λj , t). Components of a random vector carry a subscript, Z = [Z1 , . . . , Zk ]T , and ej denotes the j-th standard basis vector. Superscripts in brackets, A(1) , Z (2) , index the two competing representations of the same observed vector and are never exponents. Notation 2.2 (Multi-indices and weak convergence). For a multi-index α = (α1 , . . . , αd ) ∈ Nd0 we write |α| := α1 + · · · + αd for its order, and, for x ∈ Rd and a sufficiently differentiable f : Rd → C, x := α

d Y αj

∂ α f :=

xj ,

j=1

∂ |α| f

, ∂xα1 1 · · · ∂xαd d

(3)

with the convention x0 := 1 and ∂ 0 f := f . The symbol ⇒ denotes weak convergence of laws: for random vectors Xn , X in Rd we write Xn ⇒ X if E[f (Xn )] → E[f (X)] for every bounded continuous f : Rd → R. Notation 2.3 (Distributional relations). For a random vector X in Rp we write L(X) for its law and φX (t) := E exp(i tT X) , 

t ∈ Rp ,



(4)

for its characteristic function. We reserve • X ∼ N (µ, Σ) for “X is distributed according to the law N (µ, Σ)”, and d

• X = Y for “X and Y have the same law”, i.e. L(X) = L(Y ), equivalently φX = φY . Stochastic independence is written X ⊥⊥ Y . Mutual independence of a family {X1 , . . . , Xk } always means independence of the whole family, not merely pairwise independence. Statements such as “Zj is non-constant” are always understood P-almost surely, i.e. L(Zj ) is not a Dirac measure δc . 4

Notation 2.4 (Degenerate Gaussian distributions). We call E a (possibly degenerate) Gaussian random vector, E ∼ N (µ, Σ) with µ ∈ Rp and a positive semi-definite Σ ∈ Rp×p , if 



φE (t) = exp i tT µ − 21 tT Σt ,

(5)

t ∈ Rp .

This includes the degenerate cases: Σ is allowed to be singular, and for p = 1, σ 2 = 0, we have N (µ, 0) = δµ , so that constants count as Gaussian. Whenever we want to exclude this case we say non-degenerate.

Characteristic Functions of Probability Measures on Rd

3

Every argument in these notes is carried out at the level of characteristic functions. This section collects the facts we need, in the generality in which we need them. We assume familiarity with measure-theoretic probability – laws, weak convergence, product measures, Fubini’s theorem – and we cite the standard results rather than reproving them; everything that is specific to the independent component analysis (ICA) setting is proved in full. Standard references are Klenke (2020, Chapter 15) and Kallenberg (2021, Chapters 5–6) for the general theory, Lukacs (1970) for characteristic functions specifically, and Kagan et al. (1973, Chapter 1) for the characterisation-theoretic material.

3.1

Definition and elementary properties

Definition 3.1 (Characteristic function). Let µ be a probability measure on Rd , B(Rd ) . Its characteristic function (c.f.) is the map

φµ : Rd → C,

φµ (t) :=

Z

exp i tT x µ(dx).

Rd

(6)

For a random vector X in Rd with law L(X) = µ we write φX := φµ , so that φX (t) =   T E exp(i t X) . The integral is well defined because |exp(i tT x)| = 1, so the integrand is bounded and measurable. Proposition 3.2 (Elementary properties). Let X be a random vector in Rd with characteristic function φX . Then: (i) φX (0) = 1 and |φX (t)| ≤ 1 for all t ∈ Rd ; (ii) φX is uniformly continuous on Rd ; d

(iii) φX (−t) = φX (t), and φX is real valued if and only if X = −X; (iv) φX is positive semi-definite: for all n ∈ N, all t1 , . . . , tn ∈ Rd and all z1 , . . . , zn ∈ C, n X

zj zl φX (tj − tl ) ≥ 0;

(7)

j,l=1

(v) (affine images) for A ∈ Rm×d and b ∈ Rm , φAX+b (t) = exp i tT b · φX AT t ,

t ∈ Rm ;

(8)

(vi) (independent sums) if X ⊥ ⊥ Y are random vectors in Rd , then φX+Y (t) = φX (t) · φY (t), 5

t ∈ Rd ;

(9)

(vii) (one-dimensional projections) for u ∈ Rd and s ∈ R, φuT X (s) = φX (su).

(10)

Proof. (i) is immediate. For (ii) note that |φX (t + h) − φX (t)| ≤ E |exp(i hT X) − 1| , whose right hand side does not depend on t and tends to 0 as h → 0 by dominated convergence. (iii): the identity exp(i tT X) = exp(−i tT X) gives φX (−t) = φX (t), and hence φX is real as soon 



d

as X = −X, since then φX (−t) = φ−X (t) = φX (t). Conversely, if φX is real then φ−X (t) = d

φX (−t) = φX (t) = φX (t) for every t, so X = −X by Theorem 3.4 (i) below, whose proof does
P P 2 T not use this proposition. For (iv) we argue j,l zj zl exp i(tj −tl )T X = j zj exp(i tj X) ≥ 0 after taking expectations. For (v), E[exp(i tT (AX + b))] = exp(i tT b) E[exp(i (AT t)T X)]. (vi) is the factorisation of the expectation of a product of independent bounded random variables, and (vii) is the special case A = uT , b = 0 of (v). Properties (8), (9) and (10) are the three identities that carry all the ICA arguments: a linear ICA model is by definition an independent sum of affine images, so its characteristic function factorises into the characteristic functions of the sources evaluated along the columns of the mixing matrix. Theorem 3.3 (Bochner’s theorem; see Klenke, 2020, Chapter 15; Lukacs, 1970, Section 4.2). A function ϕ : Rd → C is the characteristic function of a probability measure on Rd if and only if ϕ is continuous, ϕ(0) = 1, and ϕ is positive semi-definite in the sense of Proposition 3.2 (iv). Beyond the standard background listed at the start of Section A, Theorem 3.3 is the only result from the theory of characteristic functions that any proof in these notes uses without our having proved it: the difficult direction is a genuinely deep theorem, and reproducing it here would take us too far afield. (Two further classical facts are quoted in remarks – Pólya’s criterion in Remark 3.8 and the Lukacs characterisation in Remark 3.16 – but no proof depends on either.) It is what makes statements of the form “ϕ is a characteristic function” verifiable without exhibiting the underlying law. We shall repeatedly need such statements – for instance in Theorem 6.7, where a supremum over admissible Gaussian factors is taken – although in practice we usually certify them through the more convenient Theorem 3.7 rather than by verifying positive semi-definiteness. The one place where we do verify positive semi-definiteness directly is Remark 3.8.

3.2

How characteristic functions determine distributions

Theorem 3.4 (Uniqueness and inversion). Let µ, ν be probability measures on Rd . (i) If φµ = φν on all of Rd , then µ = ν. Equivalently, for random vectors X, Y in Rd , d

X = Y ⇐⇒ φX = φY .

(11)

(ii) If in addition φµ ∈ L1 (Rd ), then µ has a bounded continuous Lebesgue density given by the inversion formula f (x) = (2π)−d

Z

exp −i tT x φµ (t) dt.

Rd

(12)

The proof is in Section A.1; see also Klenke, 2020, Chapter 15 and Kallenberg, 2021, Chapter 5. 6

Theorem 3.5 (Cramér–Wold device; see Klenke, 2020, Chapter 15). Let X, Y be random d

d

vectors in Rd . Then X = Y if and only if uT X = uT Y for every u ∈ Rd . Likewise, Xn ⇒ X weakly if and only if uT Xn ⇒ uT X for every u ∈ Rd . Proof. By Eq. (10) we have φuT X (1) = φX (u), so knowing the laws of all one-dimensional projections is the same as knowing φX pointwise, and Theorem 3.4 (i) applies. For the second statement, suppose first Xn ⇒ X. By Theorem 3.7 (i) in Rd , φXn → φX pointwise, so φuT Xn (s) = φXn (su) → φX (su) = φuT X (s) for every s; the limit is continuous at s = 0, so Theorem 3.7 (ii) in d = 1 gives uT Xn ⇒ uT X. Conversely, suppose uT Xn ⇒ uT X for every u. By Theorem 3.7 (i) in d = 1 and Eq. (10) at s = 1 we get φXn (u) → φX (u) for every u ∈ Rd , and φX is continuous at the origin, so Theorem 3.7 (ii) yields a random vector d

X ′ with φX ′ = φX and Xn ⇒ X ′ ; by Theorem 3.4 (i), X ′ = X and hence Xn ⇒ X. Remark 3.6 (Two different theorems of Cramér). Theorem 3.5 is the Cramér–Wold device, a statement about reducing Rd to R. It should not be confused with Cramér’s decomposition theorem, Theorem 3.20 below, which says that a Gaussian law has only Gaussian factors. Both are used in these notes, for entirely different purposes. Theorem 3.7 (Lévy’s continuity theorem). Let (Xn )n∈N be random vectors in Rd . (i) If Xn ⇒ X weakly, then φXn → φX pointwise. (ii) Conversely, if φXn (t) → ϕ(t) for every t ∈ Rd and ϕ is continuous at t = 0, then ϕ is the characteristic function of a random vector X and Xn ⇒ X weakly. The proof is in Section A.2; see also Klenke, 2020, Section 15.3. The continuity requirement in (ii) is not a technicality: for Xn ∼ N (0, n) we have φXn (t) = exp(−nt2 /2) → 1{t=0} pointwise, and the limit – discontinuous at the origin – is not a characteristic function; the mass has escaped to infinity. Remark 3.8 (Local equality of characteristic functions is not enough). Theorem 3.4 requires φµ = φν on all of Rd . Agreement on a neighbourhood of the origin does not suffice. The classical counterexample goes back to Pólya (1949); in the form given by Feller (1971, Chapter XV) and Lukacs (1970, p. 85) it reads as follows. Let ϕ1 (t) := 1 − |t| + ,

(13)

the “tent function”. It is a characteristic function. Indeed, with h := 1[−1/2,1/2] one has, for every t ∈ R, (h ∗ h)(t) =

Z R

h(u) h(t − u) du = 1 − |t| + = ϕ1 (t),

and since h is real-valued and even, ϕ1 (tj − tl ) = t1 , . . . , tn ∈ R and z1 , . . . , zn ∈ C n X j,l=1

zj zl ϕ1 (tj − tl ) =

(14)

R h(tj − u)h(tl − u) du, whence for all n ∈ N,

R

Z X n

2

zj h(tj − u) du ≥ 0.

(15)

R j=1

So ϕ1 is continuous, ϕ1 (0) = 1 and positive semi-definite, and Theorem 3.3 certifies that ϕ1 = φµ1 for a probability measure µ1 on R. (Alternatively one may invoke Pólya’s criterion, Feller (1971, Section XV.2): every continuous even ϕ : R → R that is convex on [0, ∞) and satisfies ϕ(0) = 1 and ϕ(t) → 0 as t → ∞ is a characteristic function; the tent function is the

7

standard example, being the maximum of two affine functions on [0, ∞).) Now Theorem 3.4 (ii) may be applied, since ϕ1 ∈ L1 (R): the law µ1 is absolutely continuous with density f (x) =

1 2π

Z 1 −1

1 − |t| e−itx dt =

1 π

Z 1 0

(1 − t) cos(tx) dt =

1 − cos x . πx2

(16)

Let ϕ2 be the 2-periodic extension of 1 − |t| from [−1, 1] to all of R. Expanding the resulting triangular wave into its Fourier series gives ϕ2 (t) =

X 1 4 + cos(nπt), 2 2 n odd π n2

(17)

which is the characteristic function of the discrete law with mass 12 at 0 and mass 2/(π 2 n2 ) at P 2 each of ±nπ for odd n ≥ 1; the masses do sum to 12 + π42 n odd n−2 = 12 + π42 · π8 = 1. Then ϕ1 = ϕ2 on [−1, 1], while the two laws are as different as they could be: one is absolutely continuous, the other purely atomic. This is why every local statement in these notes is stated with the explicit qualifier “in a neighbourhood of the origin”, and why such statements are never silently upgraded to statements about laws. The upgrade is legitimate under an analyticity hypothesis; see Corollary 3.14.

3.3

Moments, the distinguished logarithm, and cumulants

Proposition 3.9 (Moments and derivatives). Let X be a random vector in Rd and n ∈ N   n with E ∥X∥ < ∞. Then φX ∈ C n (Rd ; C) and, for every multi-index α with |α| ≤ n and every t ∈ Rd , T

∂ α φX (t) = i|α| E X α ei t X ; 



in particular

∂ α φX (0) = i|α| E X α . 



(18)

In the scalar case d = 1 this gives the Taylor expansion φX (t) =

n X E[X m ] m=0

m!

(it)m + o |t|n ,

t → 0.

(19)

Proof. We prove Eq. (18) by induction on |α|, the case α = 0 being the definition. Suppose it holds for α with |α| < n and let j ∈ {1, . . . , d}. For h ̸= 0, ∂ α φX (t + hej ) − ∂ α φX (t) eihXj − 1 T = i|α| E X α ei t X , h h "

#

(20)

and the integrand is dominated in modulus by |X α ||Xj | ≤ ∥X∥|α|+1 , which is integrable because |α| + 1 ≤ n; here we used |eihx − 1| ≤ |hx|. Since (eihXj − 1)/h → iXj pointwise as T h → 0, dominated convergence gives ∂j ∂ α φX (t) = i|α|+1 E[X α Xj ei t X ], which is the claim for the multi-index α + ej . The same domination shows that each such derivative is continuous in t, so φX ∈ C n . Eq. (19) is Taylor’s theorem with Peano remainder applied to the C n function φX , together with Eq. (18) at t = 0. A partial converse holds, and it needs no moment hypothesis at all: for a real-valued random variable Y and m ∈ N, if φY is 2m times differentiable at the single point 0, then E[Y 2m ] < ∞; the odd moments up to order 2m are then finite by Lyapunov’s inequality. This is Lemma A.3, proved in Section A. All our arguments manipulate logarithms of characteristic functions. Since a characteristic function is complex valued and may vanish, this requires a small amount of care, which the following lemma settles once and for all. 8

Lemma 3.10 (Distinguished logarithm). Let X be a random vector in Rd . (Uniqueness.) If U ∋ 0 is a connected open subset of Rd and ψ, ψ̃ : U → C are continuous with ψ(0) = ψ̃(0) = 0 and exp ψ = exp ψ̃ = φX on U , then ψ = ψ̃; so the germ at the origin is well defined and we may compare distinguished logarithms on any such U . (Existence.) There exist δ > 0 and a continuous function ψX : Bδ (0) → C with ψX (0) = 0 and φX (t) = exp ψX (t) ,

t ∈ Bδ (0),

(21)

where Bδ (0) := {t ∈ Rd : ∥t∥ < δ}. We call ψX the cumulant generating function of X (near  the origin). If E ∥X∥n < ∞, then ψX ∈ C n on Bδ (0). Proof. By Proposition 3.2 (i)–(ii) the function φX is continuous with φX (0) = 1, so there is a δ > 0 with |φX (t) − 1| < 1 for all t ∈ Bδ (0). Hence φX maps Bδ (0) into the open disc of radius 1 around 1, which is contained in the slit plane C \ (−∞, 0] on which the principal branch Log of the complex logarithm is defined and continuous. Put ψX := Log ◦ φX ; then ψX is continuous, ψX (0) = Log 1 = 0 and exp ◦ ψX = φX . For uniqueness, let U ∋ 0 be any connected open subset of Rd and let ψ, ψ̃ : U → C be continuous with ψ(0) = ψ̃(0) = 0 and exp ψ = exp ψ̃ = φX on U . Then exp(ψ̃ − ψ) ≡ 1 on U , so (ψ̃ − ψ)/(2πi) is a continuous integer-valued function on the connected set U vanishing at the origin, hence identically 0. This covers in particular the case U = Bδ (0), ψ = ψX . The regularity statement follows from Proposition 3.9 because Log is holomorphic on the slit plane. Definition 3.11 (Cumulants). Let X be a real-valued random variable with E |X|n < ∞ and let ψX be as in Lemma 3.10. The cumulants κ1 (X), . . . , κn (X) of X are the coefficients in the expansion 

ψX (t) =

n X κm (X) m=1

m!

(it)m + o(tn ),

t → 0,



(22)

(m)

equivalently κm (X) = i−m ψX (0). The first four cumulants are κ1 (X) = E[X],

κ2 (X) = Var(X),

κ3 (X) = E (X − EX) , 

3

κ4 (X) = E (X − EX) 

(23) 4

− 3 Var(X) . 2

(24)

Proposition 3.12 (Cumulant calculus). Let X, Y be real-valued random variables with finite n-th absolute moments and let a, b ∈ R. Then, for m ≤ n, (i) (additivity) if X ⊥ ⊥ Y then κm (X + Y ) = κm (X) + κm (Y ); (ii) (homogeneity) κm (aX) = am κm (X); (iii) (translation invariance) κm (X + b) = κm (X) for m ≥ 2, and κ1 (X + b) = κ1 (X) + b. Proof. By Eq. (9) and Lemma 3.10 we have ψX+Y = ψX + ψY on a neighbourhood of the origin, since both sides are continuous, vanish at 0 and exponentiate to φX+Y ; comparing Taylor coefficients gives (i). For (ii), φaX (t) = φX (at) by Eq. (8), so ψaX (t) = ψX (at) and the m-th coefficient picks up am . For (iii), ψX+b (t) = ψX (t) + ibt, which changes only the coefficient of (it)1 . Additivity is the reason cumulants, rather than moments, are the natural bookkeeping device for ICA: a mixture is an independent sum, so its cumulants are sums of the source cumulants. 9

3.4

Analytic characteristic functions

Theorem 3.13 (Analyticity in a strip). Let X be a random vector in Rd and suppose that there is an a > 0 with E exp a∥X∥ 



(25)

< ∞.

Then the map z 7−→ E exp i zT X 

(26)



is well defined and holomorphic on the tube Sa := {z ∈ Cd : ∥Im z∥ < a}, and it restricts to φX on Rd . If Eq. (25) holds for every a > 0, then φX extends to an entire function on Cd . Proof. We give the argument for d = 1, which is the only case used below: every appeal to this theorem in these notes is to a real-valued random variable, and the one place where a multivariate statement is wanted, Corollary 3.14, is reduced to one variable via Theorem 3.5. For general d the same argument applied in each coordinate separately, with |·| replaced by ∥·∥, yields separate holomorphy together with local boundedness on Sa ; joint holomorphy then follows from Osgood’s lemma (or, dispensing with boundedness, from Hartogs’ theorem). For z = s + iu with |u| < a we have |exp(izX)| = exp(−uX) ≤ exp(a|X|), which is integrable by Eq. (25); so F (z) := E[exp(izX)] is well defined on the strip Sa , and it is continuous there by dominated convergence. For any closed triangle ∆ ⊆ Sa , Fubini’s theorem applies because (z, ω) 7→ exp(izX(ω)) is jointly measurable and dominated by the integrable exp(a|X|) uniformly on the compact set ∂∆; hence I ∂∆

F (z) dz = E

I ∂∆



exp(izX) dz = 0,

(27)

the inner integral vanishing because z 7→ exp(izX) is entire for each fixed value of X. By Morera’s theorem F is holomorphic on Sa , and F |R = φX by construction. Corollary 3.14 (Rigidity of analytic characteristic functions). Let X, Y be random vectors in Rd both satisfying Eq. (25) for some a > 0. If φX = φY on some neighbourhood of the d

origin in Rd , then X = Y . Proof. Replacing a by the smaller of the two constants we may assume that X and Y satisfy Eq. (25) with the same a > 0, and let δ > 0 be such that φX = φY on Bδ (0). Fix u ∈ Rd \ {0}. By Eq. (10), φuT X (s) = φX (su) = φY (su) = φuT Y (s),

|s| < δ/∥u∥.

(28)

Moreover |uT X| ≤ ∥u∥∥X∥, so the real-valued random variables uT X and uT Y satisfy Eq. (25) with the constant a/∥u∥, and by Theorem 3.13 their characteristic functions extend holomorphically to a strip in C. Two functions holomorphic on a connected open subset of C which agree on a real interval agree on all of it, by the identity theorem in one complex variable; d hence φuT X = φuT Y on all of R and therefore uT X = uT Y by Theorem 3.4 (i). Since u was d

arbitrary, Theorem 3.5 gives X = Y . We deliberately reduced to one complex variable here. In several variables the naive form of the identity theorem fails: agreement of two holomorphic functions on a set with accumulation points is not sufficient, as f (z1 , z2 ) = z1 , vanishing on the uncountable set {0} × C, shows. A d-dimensional argument is nevertheless available – at a real point of a 10

real neighbourhood on which h vanishes, every complex partial derivative of h agrees with the corresponding real one and hence vanishes, so h ≡ 0 on a polydisc and then on all of the connected Sa – but the reduction to one variable is shorter and is all we need. Corollary 3.14 is exactly what Remark 3.8 rules out in general: the tent function of that remark is the characteristic function of a law whose density is of order x−2 , which has no exponential moments at all. Proposition 3.15 (Ridge property). Let d = 1 and let X satisfy Eq. (25) for some a > 0, and denote again by φX the holomorphic extension of Theorem 3.13. Then for all s ∈ R and |u| < a, |φX (s + iu)| ≤ φX (iu) = E exp(−uX) ∈ R>0 , 

(29)



that is, on each horizontal line the modulus is maximal on the imaginary axis. Proof. On the strip we have φX (z) = E[exp(izX)], so |φX (s + iu)| = E eisX e−uX 



≤ E eisX e−uX = E e−uX = φX (iu), 







(30)

which is real and strictly positive because e−uX > 0. Remark 3.16 (Ridge property and the moment generating function). Evaluating the holomorphic extension of Theorem 3.13 at z = −iu gives φX (−iu) = E[exp(uX)] = MX (u), the moment generating function of X; Eq. (29) then says that |φX | is maximal on each horizontal line at its purely imaginary point. If φX is analytic in some neighbourhood of the origin – equivalently, if 0 lies in the interior of {u : MX (u) < ∞} – then its maximal strip of analyticity is exactly the one cut out by that interior (Lukacs, 1970, Theorem 7.1.1); the moment generating function may well be finite at a boundary abscissa while φX is analytic only on the open strip. Without the hypothesis at the origin the statement is false: a law with an exponential right tail and a |x|−3 left tail has {MX < ∞} = [0, 1) with non-empty interior, yet E[X 2 ] = ∞, so φX is not even twice differentiable at the origin, let alone analytic there. The following theorem is the single most important analytic input for these notes. It says that the exponential of a polynomial is a characteristic function only in the Gaussian case. Theorem 3.17 (Marcinkiewicz’ theorem, Marcinkiewicz, 1939). Let Y be a real-valued random variable whose characteristic function is of the form φY (t) = exp g(t)

(31)

in a neighbourhood of the origin, with a (complex) polynomial g. Then deg(g) ≤ 2; more precisely, after subtracting from g the constant g(0) ∈ 2πi Z, which changes exp(g) not at all and the degree only in the trivial case of a constant g, g(t) = − 12 σ 2 t2 + iµt,

σ 2 ∈ R≥0 ,

µ ∈ R,

(32)

and consequently Y has a (possibly degenerate) Gaussian distribution, Y ∼ N (µ, σ 2 ), where the degenerate case σ 2 = 0 means Y ∼ N (µ, 0) = δµ . The proof is in Section A.4; see also Kagan et al. (1973, Lemma 1.4.2) and Lukacs (1970, Section 7.3). In the language of Definition 3.11: there is no probability distribution whose cumulant generating function equals a polynomial of degree 3 or higher near the origin. The superficially similar statement that the vanishing of all cumulants of order > n already forces the vanishing 11

of all cumulants of order > 2 is also true, but it is not a restatement of Theorem 3.17: vanishing Taylor coefficients do not by themselves make ψY a polynomial, and an additional analytic argument is required, namely a growth bound on the moments strong enough to make φY entire. For n = 2 – the only case we shall need – such a bound is immediate, because then the moments are Gaussian moments; we give the argument inside the proof of Theorem 3.19, in the implication (v)⇒(ii). For general n one argues instead as follows. The number of set partitions of {1, . . . , m} all of whose blocks have size at most n is bounded by C(n)m (m!)1−1/n , so the moment–cumulant formula gives |E[Y m ]| ≤ C ′m (m!)1−1/n with a constant C ′ ≥ 1 depending on n and on the finitely many non-zero cumulants. Lyapunov’s inequality turns this into a bound of the same shape for the absolute moments, with a
m/(m+1)
1−1/n larger constant: E|Y |m ≤ E[Y m+1 ] ≤ (C ′ )m+1 (m + 1)! for odd m, and since
1−1/n

(m + 1)! /m! = (m + 1)1−1/n (m!)−1/n one gets m≥0 am E|Y |m /m! < ∞ for every a > 0. By Tonelli’s theorem that sum equals E[exp(a|Y |)], so Eq. (25) holds for every a > 0 and Theorem 3.13 makes φY entire; then ψY is holomorphic near 0 and equal to its Taylor polynomial there, and Theorem 3.17 applies. We shall not use this generalisation. P

Remark 3.18 (How Theorem 3.17 will be used). We shall use Theorem 3.17 in the following two guises. (i) If two real-valued random variables Y1 , Y2 satisfy φY2 (λt) = φY1 (t) exp g(t)

(33)

near the origin with a polynomial g and a scalar λ ∈ R \ {0}, then Y2 is Gaussian if and only if Y1 is Gaussian. Indeed, by Eq. (8) the left hand side is φλY2 (t), so it suffices to treat λ = 1 and then to note that λY2 is Gaussian if and only if Y2 is, again by Eq. (8) and λ = ̸ 0. For λ = 1: if Y1 is Gaussian then φY1 is itself the exponential of a polynomial of degree ≤ 2, hence so is φY2 near the origin and Theorem 3.17 applies; the converse follows by exchanging the roles and replacing g by −g. (ii) The sum of an independent non-Gaussian random variable and a Gaussian random variable is again non-Gaussian: if Y = Y1 + Y2 with Y1 ⊥⊥ Y2 and Y2 Gaussian, then φY1 (t) = φY (t)φY2 (t)−1 , and φ−1 Y2 is the exponential of a polynomial of degree ≤ 2; so if Y were Gaussian, Y1 would be Gaussian as well. Note that φY2 has no zeros for Gaussian Y2 , so the above quotients are well defined on all of R.

3.5

Gaussian, sub-Gaussian and super-Gaussian distributions

Theorem 3.19 (Characterisation of the Gaussian law). Let X be a random vector in Rd . The following are equivalent. (i) X ∼ N (µ, Σ) is (possibly degenerate) Gaussian in the sense of Notation 2.4; (ii) φX (t) = exp i tT µ − 12 tT Σt for all t ∈ Rd , for some µ ∈ Rd and some positive semidefinite Σ;

(iii) every one-dimensional projection uT X, u ∈ Rd , is a (possibly degenerate) Gaussian real-valued random variable. For d = 1 these are further equivalent to each of: (iv) ψX coincides with a polynomial on a neighbourhood of the origin; (v) X has moments of all orders and κm (X) = 0 for all m ≥ 3. 12

Proof. (i)⇔(ii) is Notation 2.4 combined with Theorem 3.4 (i). (ii)⇒(iii): Eq. (10) gives 



φuT X (s) = exp is uT µ − 12 s2 uT Σu ,

(34)

which is the characteristic function of N (uT µ, uT Σu). (iii)⇒(ii): put µu := E[uT X] and σu2 := Var(uT X), which are finite by (iii); then φX (u) = φuT X (1) = exp(iµu − 12 σu2 ), and u 7→ µu is linear while u 7→ σu2 is a positive semi-definite quadratic form, so µu = uT µ and σu2 = uT Σu for µ := E[X] and Σ := Cov(X). (ii)⇒(iv) is clear, and (iv)⇒(ii) is Marcinkiewicz’ Theorem 3.17. (ii)⇒(v) is immediate from Eq. (22). For (v)⇒(ii) some care is needed, because Eq. (22) is only an asymptotic expansion: knowing all Taylor coefficients of ψX does not by itself make ψX a polynomial, so Theorem 3.17 is not directly applicable. Instead, put µ := κ1 (X) and σ 2 := κ2 (X). The moment–cumulant relations express E[X m ] as a universal polynomial in κ1 (X), . . . , κm (X); with κm (X) = 0 for all m ≥ 3 these are exactly the moments of N (µ, σ 2 ). Since exp(a|x|) ≤ exp(ax) + exp(−ax), monotone convergence gives, for every a > 0, E exp(a|X|) ≤ 2 



X a2k k≥0

(2k)!

E X 2k = 2 



X a2k k≥0

(2k)!

E N 2k < ∞, 



(35)

where N ∼ N (µ, σ 2 ). The last series converges because it is dominated by 2 E[exp(a|N |)], which is finite. So X satisfies Eq. (25) for every a > 0, and by Theorem 3.13 its characteristic function is entire, hence equal to its everywhere convergent Taylor series about the origin. By Proposition 3.9 that series is determined by the moments of X, which are those of N , so φX = φN , which is (ii). Theorem 3.20 (Cramér’s decomposition theorem, Cramér, 1936). If Y1 ⊥⊥ Y2 are real-valued random variables such that Y1 + Y2 is Gaussian, then both Y1 and Y2 are Gaussian (possibly degenerate). The proof is in Section A.5; see also Kagan et al., 1973, Theorem 1.1.1. Theorems 3.17 and 3.20 are dual to each other and together explain the special role of the Gaussian law: it cannot be built out of non-Gaussian independent pieces, and it cannot be approached by “almost polynomial” cumulant generating functions. The next definition quantifies the deviation from normality that ICA exploits. Definition 3.21 (Excess kurtosis; sub- and super-Gaussian). Let X be a real-valued random variable with E[X 4 ] < ∞ and σ 2 := Var(X) ∈ (0, ∞). Its excess kurtosis is E (X − EX)4 κ4 (X) kurt(X) := = − 3. σ4 σ4 



(36)

We call X • super-Gaussian (or leptokurtic) if kurt(X) > 0, • sub-Gaussian (or platykurtic) if kurt(X) < 0, • mesokurtic if kurt(X) = 0. This is the convention of the ICA literature (Comon, 1994; Hyvärinen et al., 2001); see Remark 3.25 for the clash with the concentration-theoretic use of the same words, and Remark 3.24 for what the sign of kurt does and does not say about tails. For the laws customarily used as ICA source models the picture is the familiar one: the super-Gaussian 13

ones are peaked at the mode and heavy tailed – the Laplace and Student laws used to model the “sparse” sources of natural images and audio – while the sub-Gaussian ones are flat-topped, the uniform law being the prototype. Remark 3.24 shows that this correspondence is a property of those particular families rather than a general implication. Proposition 3.22 (Kurtosis through the characteristic function). Let X be a real-valued random variable with E[X 4 ] < ∞. Then (4)

κ4 (X) = ψX (0) =

d4 , Log φX (t) dt4 t=0

(37)

and for independent real-valued random variables X1 , . . . , Xk with finite fourth moments and weights w ∈ Rk , κ4

k X



wj Xj =

j=1

k X

wj4 κ4 (Xj ).

(38)

j=1

Proof. The first identity is Definition 3.11 with m = 4, using i−4 = 1 and Lemma 3.10. The second follows from Proposition 3.12 (i)–(ii). Remark 3.23 (kurt = 0 does not mean Gaussian). By Theorem 3.19 (v) a non-degenerate Gaussian random variable is mesokurtic, but the converse fails badly. Consider the symmetric three-point law P[X = a] = P[X = −a] = p,

P[X = 0] = 1 − 2p,

(39)

with a > 0 and p ∈ (0, 21 ). Then φX (t) = 1 − 2p + 2p cos(at), σ 2 = 2pa2 , E[X 4 ] = 2pa4 , and hence κ4 (X) = 2pa4 (1 − 6p),

kurt(X) =

1 − 6p . 2p

(40)

For p = 16 this vanishes, yet X is a three-point law and certainly not Gaussian. Consequently, non-Gaussianity is a strictly weaker requirement than kurt ̸= 0 – which is precisely why the identifiability theorems of Sections 5 and 6 are formulated in terms of non-normality rather than in terms of kurtosis. Remark 3.24 (Kurtosis and tail weight). Excess kurtosis is routinely glossed as “tail weight”. The gloss is convenient, but it is not a theorem, and it is worth being precise about what kurt does and does not measure. What it measures. Let Z := (X − EX)/σ be the standardisation of X, so that E[Z] = 0 and E[Z 2 ] = 1. Since Var(Z 2 ) = E[Z 4 ] − (E[Z 2 ])2 = E[Z 4 ] − 1, we have the identity kurt(X) = E Z 4 − 3 = Var Z 2 − 2 = E (Z 2 − 1)2 − 2. 







(41)

The quantity Z 2 − 1 vanishes exactly at Z = ±1, that is at X = EX ± σ, so excess kurtosis is the dispersion of X about those two points: probability mass moved away from EX ± σ increases it, whether it moves outwards into the tails or inwards towards the centre. This is Moors’ reading of kurtosis as a measure of dispersion around µ ± σ (Moors, 1986); Balanda and MacGillivray (1988) phrase it as the location- and scale-free movement of probability mass from the “shoulders” of a distribution into its centre and its tails, and argue that the notion is irreducibly vague, admitting several inequivalent formalisations. Eq. (41) also yields the sharp bound kurt ≥ −2, with equality if and only if Z 2 = 1 almost surely, i.e. exactly when X takes the two values EX ± σ with probability 12 each – any affine image of the Rademacher 14

law, whose entry in Table 1 records this value. For the same reason kurtosis is not a measure of peakedness either (Kaplansky, 1945; Darlington, 1970). What it does not measure. Read “tail weight” as the asymptotic decay of the density, or of P[|X| > x], relative to a Gaussian. Under that reading the sign of kurt and tail weight are logically independent, in both directions. (a) Light tails, arbitrarily large kurtosis. The three-point law Eq. (39) is supported on {−a, 0, a}, so its tails are as light as tails can be; yet kurt(X) = (1 − 6p)/(2p) → ∞ as p ↓ 0, by Eq. (40). (b) Heavy tails, negative kurtosis. Let X have the mixture law L(X) = (1 − ε) Unif[−1, 1] + ε Laplace(1),

ε = 10−3 .

(42)

Then σ 2 = 67/200 and kurt(X) = −4515/4489 ≈ −1.006, so X is sub-Gaussian in the sense of Definition 3.21. Its density equals 12 εe−|x| for |x| > 1, and therefore exceeds the density of the Gaussian law of the same variance for |x| ≥ x0 with x0 ≈ 2.56, by a factor of about 3.1 · 105 at |x| = 4. (c) Heavy tails, no kurtosis at all. The Cauchy law, and the Student tν laws with ν ≤ 4, have E[X 4 ] = ∞, so kurt is not even defined; cf. Table 1. Under a different reading of “tail” the verdict changes. If tail weight is taken to mean the propensity to produce observations far from the mean measured in units of σ, rather than asymptotic decay, then kurtosis is a tail functional, and Westfall (2014) argues that this is its only unambiguous interpretation. Example (a) shows what separates the two readings: the √ three-point law has compact support, but its atoms sit at ±a = ±σ/ 2p, that is arbitrarily many standard deviations from the mean. Convention. In these notes “sub-Gaussian” and “super-Gaussian” always mean Definition 3.21, the sign of kurt, and never a statement about tails. Where tails are meant – in Table 1, in Remark 3.25 and in Corollary 6.9 – they are named explicitly and always refer to asymptotic decay. Remark 3.25 (Two incompatible meanings of “sub-Gaussian”). In concentration of measure and high-dimensional statistics, a real-valued random variable X is called sub-Gaussian when its tails are dominated by those of a Gaussian law, equivalently – up to the value of the constant σ – when E exp λ(X − EX) 







≤ exp 12 λ2 σ 2 ,

λ ∈ R,

(43)

for some σ > 0. By Theorem 3.13 this forces φX to be entire, and the displayed bound  
gives |φX (z)| ≤ E e|Im z||X| ≤ 2 exp |Im z||EX| + 12 σ 2 |Im z|2 , so φX is of order at most 2. This is not the notion of Definition 3.21, and neither notion implies the other. Indeed, the three-point law Eq. (39) with p < 16 is bounded, hence sub-Gaussian in the concentration sense, while Eq. (40) gives kurt(X) = (1 − 6p)/(2p) > 0, so it is super-Gaussian in the sense of Definition 3.21. Conversely, the mixture of Remark 3.24 (b) has exponential tails and so is not sub-Gaussian in the concentration sense, while its excess kurtosis is negative, so it is sub-Gaussian in the sense of Definition 3.21. Throughout these notes, “sub-Gaussian” and “super-Gaussian” always refer to Definition 3.21.

3.6

Why non-Gaussianity is what makes ICA work

The classical heuristic behind ICA is the central limit theorem: a normalised mixture wT Z = P j wj Zj of many independent sources looks “more Gaussian” than any individual source, so 15

Law of X

φX (t)

N (0, σ 2 )

e−σ t /2 sin(at) at cos t

kurt(X)

type

analyticity of φX

0 6 − 5 −2 1 − 6p 2p 3

Gaussian

entire, order 2

sub-Gaussian

entire, order 1

sub-Gaussian

entire, order 1

sub- iff p > 61

entire, order 1

super-Gaussian

strip |Im z| < b−1

6 ν−4

super-Gaussian

not analytic at 0

undefined

no moments

not differentiable at 0

2 2

Uniform on [−a, a] Rademacher on {±1} Three-point Eq. (39)

1 − 2p + 2p cos(at) 1 1 + b2 t2 √ ( ν|t|)ν/2 Kν/2 ( ν|t|) Γ(ν/2) 2ν/2−1 e−γ|t|

Laplace, scale b

√

Student tν , ν > 4 Cauchy, scale γ

Table 1: Characteristic functions, excess kurtosis and domain of analyticity for the standard examples. Kν/2 denotes the modified Bessel function of the second kind. Power-law tails destroy analyticity at the origin; exponential tails give a strip; bounded or Gaussian tails give an entire function. one may hope to recover the sources by searching for the directions in which the mixture is least Gaussian. Eq. (38) turns this heuristic into a precise statement. Proposition 3.26 (Kurtosis as a contrast function). Let Z1 , . . . , Zk be independent real-valued random variables with E[Zj ] = 0, Var(Zj ) = 1 and E[Zj4 ] < ∞, and let w ∈ Rk with ∥w∥2 = 1. Then Var(wT Z) = 1 and kurt wT Z

=

k X

wj4 kurt(Zj ) ≤ max kurt(Zj ) . 1≤j≤k

j=1

(44)

If maxj |kurt(Zj )| > 0, then equality holds if and only if w = ±ej for an index j attaining the maximum; so the maximisers of w 7→ |kurt(wT Z)| on the unit sphere are exactly the signed coordinate directions of maximal |kurt|. If instead kurt(Zj ) = 0 for every j – in particular if all sources are Gaussian – then both sides vanish for every w and the criterion carries no information at all. Proof. Since Var(Zj ) = 1 and ∥w∥2 = 1 we get κ2 (wT Z) = 1 by Proposition 3.12, so kurt(wT Z) = κ4 (wT Z), and Eq. (38) gives the stated identity. Then X j

wj4 kurt(Zj ) ≤

X

wj4 kurt(Zj ) ≤ max kurt(Zj ) · j

j

X j

wj4 ≤ max kurt(Zj ) , j

(45)

2

2 where the last step uses j wj4 ≤ = 1. Write M := maxj |kurt(Zj )| and assume j wj P M > 0. Then equality in the last inequality, which reads M j wj4 ≤ M , holds if and only P if j wj4 = 1, i.e. if and only if exactly one coordinate of w is non-zero and thus w = ±ej ; equality in the middle inequality then forces |kurt(Zj )| = M . If M = 0 all three quantities vanish for every w, which is the last assertion.

P

P

Proposition 3.26 is the theoretical justification of the kurtosis-based contrast functions of Comon (1994) and Hyvärinen et al. (2001), and it already displays the two phenomena that the rest of these notes make exact: (i) non-Gaussian sources are recoverable, and the recovery is only ever determined up to sign, scale and the ordering of the coordinates – the ambiguities of Eq. (1); 16

(ii) Gaussian sources are not recoverable at all, since the contrast degenerates. Compare Remark 3.27. What the following sections add is that neither fourth moments nor any moments at all are actually needed: the correct hypothesis is non-normality, in the sense of Theorem 3.19, and the correct tool is Theorem 3.17 rather than a contrast function. Remark 3.27 (Gaussian sources are genuinely unidentifiable). Let Z ∼ N (0, Ik ) and let A ∈ Rp×k . For every orthogonal Q ∈ O(k) we have QZ ∼ N (0, Ik ) by Eq. (8), so QZ again has mutually independent components, while d

AZ = (AQT )(QZ).

(46)

So A and AQT are indistinguishable from the law of X = AZ, and the ambiguity is a whole orthogonal group rather than only the discrete permutations and the diagonal rescalings of Section 1. This is why every identifiability statement below either excludes Gaussian sources, or explicitly quantifies the remaining Gaussian ambiguity.

4

Identifiability for Non-Constant Independent Sources

We begin with the weakest set of assumptions on the sources: they are mutually independent and non-constant, and nothing else. Before any uniqueness statement can be made we have to normalise the representation, because two trivial mechanisms produce genuinely different representations of the same random vector: zero columns of the mixing matrix, and columns that are proportional to each other. Remark 4.1 removes both. Remark 4.1 (Normalising a representation). Consider a representation of a p-dimensional random vector X ∈ Rp : X = AZ + µ =

k X

aj Zj + µ,

(47)

j=1

where Z = [Z1 , . . . , Zk ]T ∈ Rk is a random vector with mutually independent components {Z1 , . . . , Zk }, where µ ∈ Rp is a (non-random) column vector, and where A = [a1 , . . . , ak ] ∈ Rp×k is a (non-random) matrix with column vectors aj ∈ Rp . Step 1: removing zero columns and constant sources. By deleting the zero columns aj = 0 and by absorbing every almost surely constant summand aj Zj into µ, we arrive at a representation of X in which A has no zero column and in which every component of Z is almost surely non-constant. Step 2: removing proportional columns. The columns of A may still be proportional to each other, say aj = λal with a proportionality constant λ ∈ R \ {0} and j ̸= l. This is a true ambiguity of the representation, since we can always merge the two corresponding sources, aj Zj + al Zl = al λZj + Zl =: al Zel .

(48)

For a first identifiability result we therefore have to aggregate these variables. Proportionality is an equivalence relation on Rp \ {0}, so we can choose columns ã1 , . . . , ãk′ of A forming a
system of representatives of {a1 , . . . , ak } \ {0} /∝; that is, a set of non-zero columns of A that are pairwise non-proportional and such that every non-zero aj is proportional to a (then necessarily unique) ãl , written aj ∝ ãl . Denoting the corresponding proportionality constants

17

by λj , i.e. aj = λj ãl , we obtain X = AZ + µ =

k X

(49)

aj Zj + µ

(50)

j=1 ′

=

k X

!

λj Zj +µ

X

ãl

(51)

aj ∝ãl

l=1

|

{z

}

el =: Z

=

k̃ X

′

k X

ãl Zel +

l=1

(52)

ãl Zel + µ

l=k̃+1

|

{z

=: µ̃

}

= AeZe + µ̃,

(53)

where – after possibly re-indexing – k̃ ≤ k ′ denotes the number of indices l for which Zel is not almost surely constant, so that the sum inside µ̃ collects exactly the almost surely constant Zel and µ̃ is (a.s. equal to) a deterministic vector, and where Ae := [ã1 , . . . , ãk̃ ] and Ze := [Ze1 , . . . , Zek̃ ]T . If Step 1 has already been carried out, so that every Zj is non-constant, then in fact k̃ = k ′ and the second sum is empty: an independent sum of non-constant random variables is never a.s. constant, because |φY1 φY2 | ≡ 1 forces |φY1 | ≡ 1 and hence Y1 degenerate. We nevertheless state Step 2 in the general form, so that it can be applied on its own. We have thus found a representation of X in which the columns of Ae are non-zero and pairwise non-proportional and in which every component of Ze is almost surely non-constant. e since each Moreover, if the components of Z are mutually independent then so are those of Z, Zj enters exactly one Zel . Now that the existence of such normalised representations is established, we restrict attention to them and investigate their uniqueness – up to the transformations that are unavoidable. Theorem 4.2 (Identifiability for independent non-constant sources; Kagan et al., 1973, Chapter 10, Lemma 10.2.3 and Theorem 10.3.1). Let X ∈ Rp be a p-dimensional random vector with two representations A(1) Z (1) + µ(1) = X = A(2) Z (2) + µ(2) ,

(54)

with the following properties for i = 1, 2: (i)

(i) A(i) ∈ Rp×k is a (non-random) matrix whose columns are non-zero and pairwise non-proportional; (ii) µ(i) ∈ Rp is a (non-random) column vector; (iii) Z (i) ∈ Rk

(i)

is a random vector such that (i)

(i)

(a) its k (i) components {Z1 , . . . , Zk(i) } are mutually independent, and (i)

(i)

(b) each component Zj is almost surely non-constant, i.e. L(Zj ) is not a Dirac measure, j = 1, . . . , k (i) .

18

Then µ(2) − µ(1) ∈ im A(1) = im A(2) ,

rank A(1) = rank A(2) .

(1)

In particular there exist c(1) ∈ Rk and c(2) ∈ Rk Furthermore, the following statements hold.

(2)

(55)

with µ(2) − µ(1) = A(1) c(1) = A(2) c(2) . (2)

(1) If the l-th column of A(2) is not proportional to any column of A(1) , then Zl

is Gaussian.

(2) Assume that the l-th column of A(2) is proportional to the j-th column of A(1) with (2) (1) proportionality constant1 λ ∈ R \ {0}, i.e. al = λ · aj . Then there exists a (com(2)

plex) polynomial g such that the characteristic functions of Zl neighbourhood of the origin,

(1)

and Zj

satisfy, in a

φZ (2) (λt) = φZ (1) (t) · exp g(t) .

(56)

j

l

(2)

In particular, by Marcinkiewicz’ Theorem 3.17 together with Remark 3.18, Zl (1) Gaussian if and only if Zj is Gaussian.

is

The proof is in Section B. It is a self-contained finite-difference argument that uses nothing from the sections in between, so the reader may turn to it at any point; we have deferred it only because it is considerably longer than the statement and would interrupt the development here. The first part of the theorem is elementary, and we record it right away as Proposition 4.4 below. Remark 4.3 (Reading Theorem 4.2). Theorem 4.2 is best read as a dichotomy at the level of columns. Every column of A(2) either (i) is proportional to a column of A(1) , in which case, by Eq. (56), the two associated sources agree up to that proportionality factor and a factor exp(g) with g a polynomial; or (ii) is not, in which case its source is forced to be Gaussian, and by Remark 3.27 we should not have expected to recover it in the first place. In case (i), Theorem 4.2 by itself bounds neither deg(g) nor the direction in which the perturbation acts; under the additional hypotheses of Theorem 5.5 the factor becomes an independent Gaussian perturbation of one of the two sources. All the work in the following sections consists of ruling out the second alternative by strengthening the assumptions on the sources, and of bookkeeping the Gaussian perturbation in the first. Theorem 4.2 is the engine of everything that follows. Its first part – the statement about images, ranks and offsets – needs none of the machinery of Section B, so we prove it here; only the column dichotomy (1)–(2) is postponed. Proposition 4.4 (The first part of Theorem 4.2). Under the hypotheses of Theorem 4.2, the affine hull of the support of L(X) equals µ(i) + im A(i) for i = 1, 2. Consequently µ(2) − µ(1) ∈ im A(1) = im A(2) ,

rank A(1) = rank A(2) .

(57)

Proof. Fix i and abbreviate A := A(i) , Z := Z (i) , µ := µ(i) , k := k (i) . Since the components of Z are mutually independent, L(Z) is the product of the laws of its components, and the 1

Note that this proportionality constant was not reintroduced in (Kagan et al., 1973, Theorem 10.3.1) after it had been removed “without loss of generality” in (Kagan et al., 1973, Lemmata 10.2.4 and 10.2.5). It matters here, since we do not normalise the columns of A(i) .

19

support of a finite product of Borel probability measures on second countable spaces – in Q particular on R – is the product of the supports: supp L(Z) = kj=1 supp L(Zj ). Each factor contains at least two points, because Zj is almost surely non-constant. Hence, for each j, choosing two distinct points of supp L(Zj ) and fixing an arbitrary point of supp L(Zl ) in every other coordinate l ̸= j produces two elements of supp L(Z) whose difference is a non-zero multiple of ej . The direction space of the affine hull of supp L(Z) therefore contains every ej , so aff supp L(Z) = Rk .

(58)

Now let T (z) := Az + µ. For a continuous map T one has supp L(T (Z)) = T supp L(Z) : the inclusion “⊇” holds because T −1 (U ) is an open neighbourhood of any preimage point and therefore has positive mass, and “⊆” because the closed set T (supp L(Z)) has L(T (Z))measure 1. Affine hulls are unchanged by taking closures, since affine subspaces are closed, and they commute with affine maps. Hence

aff supp L(X) = T aff(supp L(Z)) = A Rk + µ = µ + im A.

(59)

Applying this to i = 1 and i = 2 and equating the two descriptions of the same affine subspace gives µ(1) + im A(1) = µ(2) + im A(2) , whence im A(1) = im A(2) (the two direction spaces coincide) and µ(2) − µ(1) ∈ im A(1) . Equality of ranks is equality of the dimensions of these images.

5

Identifiability for Non-Gaussian Independent Sources

Theorem 4.2 exhibits a clear distinction between Gaussian and non-Gaussian components in a representation X = AZ + µ: Gaussian components tend to be non-identifiable. It therefore makes sense to separate them off from the beginning and to work with a representation in which the Gaussian part appears as one additive noise vector. Remark 5.1 (Separated representations). Start from a normalised representation as in Remark 4.1 and write it as X = AZ + µ =

k X

aj Zj + µ

(60)

j=1

=

k̂ X

k X

al Zl +

l=1

al Zl + µ

(61)

l=k̂+1

|

{z

}

e) e ∼ N (µ̃,Σ =: E e = AeZe + E,

(62)

where – after possibly re-indexing – k̂ is the number of non-Gaussian components and the second sum collects all Gaussian components Zl into a single p-variate (possibly degenerate) e where Z e ∼ N (µ̃, Σ), e := [Z1 , . . . , Z ]T is a vector of mutually Gaussian random vector E k̂ e and where the mixing independent, non-constant, non-Gaussian components with Ze ⊥⊥ E, e matrix A := [a1 , . . . , ak̂ ] has non-zero, pairwise non-proportional columns. Before stating the general identifiability result – which will be based on Theorem 4.2 – we investigate how the above outsourcing of the Gaussian components can be reverted. This is the technical heart of this section: we need to be able to move Gaussian mass back and forth between the noise vector and additional columns of the mixing matrix. 20

Lemma 5.2 (Trading Gaussian noise against columns of the mixing matrix). Let X ∈ Rp be a p-variate random vector with a representation X = AZ + E,

(63)

where • A ∈ Rp×k is a (non-random) matrix with non-zero, pairwise non-proportional columns, • Z ∈ Rk is a random vector with mutually independent, non-constant, non-Gaussian components, and • E ∈ Rp is a (possibly degenerate) Gaussian random vector, E ∼ N (µ, Σ), with mean vector µ ∈ Rp and positive semi-definite covariance matrix Σ ∈ Rp×p of rank r ≤ p, such that Z ⊥ ⊥ E. Then there exist a matrix B ∈ Rp×r of rank r and a random vector V ∼ N (0, Ir ) such that Σ = BBT ,

E = BV + µ a.s.,

Z ⊥⊥ V,

(64)

which yields the representation X = AZ + BV + µ.

(65)

If r = 0 then B is the empty p × 0 matrix, V the empty vector and E = µ a.s. If r ≥ 2 then B can moreover be chosen such that no column of B is proportional to any column of A. Furthermore, for either choice of B and for every r ≥ 0, we obtain a normalised representation in the sense of Remark 4.1: " # i Z e e + µ, B e V

h

X= A

(66)

where e ∈ Rp×r̃ with r̃ ≤ r consists of those columns of B that are not proportional to any • B column of A, and Ve ∼ N (0, Ir̃ ) collects the corresponding components of V ;

• Ze = Z + Vb with Vb ∼ N (0, Γ) for a diagonal matrix Γ (possibly with zeros on the diagonal), where {Z, Vb , Ve } is mutually independent; • every component of Ze is non-Gaussian, the random vector [ZeT , Ve T ]T has mutually e has non-zero, pairwise independent non-constant components, and the matrix [A, B] non-proportional columns. Proof. The first statement follows directly from the spectral decomposition of the symmetric positive semi-definite matrix Σ (equivalently, from its singular value decomposition): "

#

∆r 0 T Σ=Q Q , 0 0

(67)

with an orthogonal matrix Q ∈ O(p) ⊆ Rp×p and a diagonal matrix ∆r ∈ Rr×r with strictly positive diagonal entries. We may then put "

1/2

∆r B := Q 0

#

h

i

T V := ∆−1/2 0 Q (E − µ), r

,

21

(68)

which yields exactly the desired properties: rank(B) = r,

Σ = BBT ,

E = BV + µ a.s.,

V ∼ N (0, Ir ),

Z ⊥⊥ V.

(69)

Note that if G ∈ O(r) is orthogonal, then replacing B by BG and V by GT V preserves all of these properties. For r ≥ 2 this freedom can be used to make every column of BG non-proportional to every column of A; the technical details are the content of Lemma 5.4. e be the matrix consisting of those columns of B = [b1 , . . . , br ] that are In any case, let B not proportional to any column of A, and let Ve be the random vector of the corresponding components of V . After re-indexing we may assume these to be the first r̃ columns of B, e = [b1 , . . . , br̃ ]. For l > r̃ we then have bl = λl aj for a unique index jl and some so that B l λl ∈ R \ {0}. With this we obtain X = AZ + E

(70)

= AZ + BV + µ =

k X

aj Zj +

j=1

=

k X

k X

aj Zj +

=

r̃ X

bl Vl +

l=1

aj Zj +

j=1 k X

bl Vl + µ

(72)

l=1

j=1

=

r X

(71)

r̃ X

j=1

|

bl Vl + µ

(73)

ajl λl Vl + µ

(74)

l=r̃+1

bl Vl +

l=1

aj Zj +

r X

r X l=r̃+1

bj =: V zX}| {!

r̃ X

bl ∝aj

l=1

λl Vl +

{z

bl Vl + µ

(75)

}

ej =: Z

=

k X

aj Zej +

j=1

h

= A

r̃ X

bl Vl + µ

(76)

l=1

e B

" # i Z e + µ. Ve

(77)

Each Zej is non-Gaussian, being the sum of the non-Gaussian Zj and independent Gaussian noise, by Marcinkiewicz’ Theorem 3.17, cf. Remark 3.18 (ii). The components of [ZeT , Ve T ]T are mutually independent because the components of [ZT , V T ]T are mutually independent e are and each component Vl enters exactly one component of [ZeT , Ve T ]T . The columns of B pairwise non-proportional because rank(B) = r makes the columns of B linearly independent, e So the matrix and they are non-proportional to the columns of A by the very choice of B. e has non-zero, pairwise non-proportional columns. [A, B] The following lemma provides the rotation used in the proof above. It is the only place where we need a genuinely geometric argument. Lemma 5.3 (Avoiding finitely many hyperplanes). Let k ≥ 1 and let A ⊆ Rk \ {0} be a finite set. Then there exists q ∈ Rk with qT a = ̸ 0 for every a ∈ A. In fact the set of such q is open and its complement is a Lebesgue null set. Proof. For a ̸= 0 the orthogonal complement a⊥ = {v ∈ Rk : vT a = 0} is a linear subspace of dimension k − 1, hence closed and a Lebesgue null set. A finite union of null sets is a null set, 22

so Rk \ a∈A a⊥ has full Lebesgue measure and is in particular non-empty; it is open as the complement of a finite union of closed sets. S

Lemma 5.4 (Rotating a frame away from finitely many directions). Let p ≥ k ≥ 2, let A ⊆ Rp \ {0} be a finite set of non-trivial (column) vectors in Rp , and let B ∈ Rp×k be a matrix with rank(B) = k. Then there exists an orthogonal matrix Q ∈ O(k) ⊆ Rk×k such that for every a ∈ A and every column vector c of C := BQ ∈ Rp×k the set {a, c} is linearly independent, i.e. a is not proportional to c (and vice versa). Proof. Put B † := (BT B)−1 BT ∈ Rk×p , which is a left inverse of B by the assumption rank(B) = k. We first discard from A those a with BT a = 0. Such an a is automatically linearly independent of every column c of BQ, whatever Q ∈ O(k): indeed c ∈ im B, so a = λc with λ= ̸ 0 would give a ∈ im B, whereas BT a = 0 says a ⊥ im B, and together these force a = 0, which is excluded. Note also that B † a = 0 if and only if BT a = 0, since BT B is invertible. Replacing A by {a ∈ A : BT a = ̸ 0} – if this is empty any Q ∈ O(k) will do – we may therefore † assume B a ̸= 0 for every a ∈ A, and consider the sets Ae := B † a | a ∈ A ⊆ Rk \ {0}, 

A′ := Ae ∪ tã | ã ∈ Ae ⊆ Rk , 

where for ã ∈ Ae the column vector tã ∈ Rk is any fixed non-trivial vector with tTã ã = 0; such a vector exists because k ≥ 2. By Lemma 5.3 there exists a vector q1 ∈ Rk that is not orthogonal to any vector in A′ , and without loss of generality we may scale it to unit norm, ∥q1 ∥2 = 1. e otherwise q1 = λã with λ ∈ R, and It follows that q1 is not proportional to any ã ∈ A: T T hence tã q1 = λ tã ã = 0, contradicting the choice of q1 . e the set A e is disjoint from the orthogonal Since q1 is also not orthogonal to any element of A, complement q1⊥ := v ∈ Rk | vT q1 = 0 , 

which is a (k − 1)-dimensional subspace of Rk and thus admits an orthonormal basis q2 , . . . , qk . We can now build the orthogonal matrix Q := [q1 , . . . , qk ] ∈ O(k) ⊆ Rk×k . Assume, by way of contradiction, that there exist a ∈ A and a column vector c of C := BQ, say c = Cej = Bqj , such that {a, c} is linearly dependent. Since a, c = ̸ 0 there is then a † λ ∈ R \ {0} with a = λc, and multiplying by B gives ã := B † a = λ · B † c = λ · B † BQej = λ · qj . e For j ≥ 2 it contradicts ã ∈ A e together For j = 1 this contradicts q1 ̸= λ−1 ã for all ã ∈ A. with Ae ∩ span{q2 , . . . , qk } = ∅. So {a, c} is linearly independent for every a ∈ A and every column vector c of BQ, which proves the claim.

Theorem 5.5 (Identifiability for independent non-Gaussian sources). Let X ∈ Rp be a random vector and assume that we have two representations A(1) Z (1) + E (1) = X = A(2) Z (2) + E (2) , with the following properties for i = 1, 2: 23

(78)

(i)

(i) A(i) ∈ Rp×k is a (non-random) matrix whose columns are non-zero and pairwise non-proportional; (ii) E (i) ∈ Rp is a p-variate Gaussian random vector, E (i) ∼ N (µ(i) , Σ(i) ), possibly degenerate, with mean vector µ(i) ∈ Rp and positive semi-definite covariance matrix Σ(i) ; (iii) Z (i) ∈ Rk

(i)

is a random vector such that (i)

(i)

(a) its k (i) components {Z1 , . . . , Zk(i) } are mutually independent, and (i)

(b) each component Zj 1, . . . , k (i) ;

is a (non-constant) non-Gaussian random variable, for j =

(iv) E (i) is independent of Z (i) : E (i) ⊥⊥ Z (i) . Then k (1) = k (2) =: k, and there exist a permutation matrix P = P (ρ) ∈ Rk×k , given by P ej = eρ(j) for a permutation ρ, and an invertible diagonal matrix Λ = diag(λ1 , . . . , λk ) ∈ Rk×k such that A(2) = A(1) P Λ,

rank A(2) = rank A(1) ,

(79)

and such that for every j = 1, . . . , k there exists a (complex) polynomial gj with φZ (2) (λj t) = φZ (1) (t) · exp gj (t)

(80)

j

ρ(j)

in a neighbourhood of the origin. Furthermore, and in particular, A(2) has a left inverse (rank(A(2) ) = k) if and only if A(1) has a left inverse (rank(A(1) ) = k). If this is the case, then for each j = 1, . . . , k separately we have deg(gj ) ≤ 2 and there exist a constant νj ∈ R and a (possibly degenerate) random variable Gj ∼ N (0, σj2 ) such that at least one of the following two cases holds, according to the sign of the degree-two coefficient of gj : (2)

λj Zj

d

(1)

+ νj = Zρ(j) + Gj , (1)

d

(2)

Zρ(j) = λj Zj

(1)

if gj′′ ≤ 0,

(81)

(2)

if gj′′ ≥ 0.

(82)

Gj ⊥⊥ Zρ(j) ,

+ νj + G j ,

Gj ⊥⊥ Zj ,

If gj′′ ̸= 0 exactly one of the two cases applies; if gj′′ = 0 then Gj = 0 and both hold and coincide. The polynomial gj is determined near the origin only up to an additive constant in 2πi Z, because φZ (1) does not vanish there; in particular deg(gj ) and gj′′ are well defined. ρ(j)

Proof. Apply Lemma 5.2 to both sides to obtain the representations h

A(1)

" # i Z h e (1) (1) e (1) + µ = X = B A(2) Ve (1)

" # i Z e (2) e (2) + µ(2) , B Ve (2)

(83)

which satisfy the requirements of Theorem 4.2. (2) (2) Since Zej is non-Gaussian, Theorem 4.2 implies that the j-th column aj of A(2) is e (1) ]. All components of Ve (1) are Gaussian whereas proportional to some column of [A(1) , B (2) those of Ze (1) are non-Gaussian, so Theorem 4.2 forces aj to be proportional to a column (1)

of A(1) , say to aρ(j) . The index ρ(j) is unique because the columns of A(1) are pairwise (2)

non-proportional, and j 7→ ρ(j) is injective: ρ(j) = ρ(j ′ ) would give aj

(1)

(2)

∝ aρ(j) ∝ aj ′ ,

contradicting pairwise non-proportionality of the columns of A(2) ; hence k (2) ≤ k (1) . Since the same argument applies to all columns of A(2) , and also with the roles of A(1) and A(2) 24

interchanged, we obtain a bijection between the columns of A(1) and those of A(2) , together with their proportionality constants. This already proves the first part of the claim, namely k (1) = k (2) =: k and the representation A(2) = A(1) P Λ for some permutation matrix P = P (ρ) with permutation ρ and some invertible diagonal matrix Λ = diag(λ1 , . . . , λk ). Indeed, multiplying by the unit vector ej from the right gives (2)

aj

(1)

= A(2) ej = A(1) P Λej = A(1) P λj ej = A(1) λj eρ(j) = λj aρ(j) .

Again by Theorem 4.2 there exists a polynomial g̃j with φZe(2) (λj t) = φλ Ze(2) (t) = φZe(1) (t) · exp g̃j (t)

(84)

j

j

j

ρ(j)

in a neighbourhood of the origin. Lemma 5.2 gives us the representations (2)

Zej

(2)

+ Vbj ,

(2)

Zj

(2)

⊥⊥ Vbj ,

(2)

(85)

(1)

(1)

Zρ(j) ⊥⊥ Vbρ(j) ,

(1)

(1)

(86)

= Zj

(1)

Zeρ(j) = Zρ(j) + Vbρ(j) , (2)

with Gaussian Vbj

(1)

and Vbρ(j) . Taking characteristic functions yields φλ Ze(2) (t) = φλ Z (2) (t) · φλ Vb (2) (t),

(87)

φZe(1) (t) = φZ (1) (t) · φVb (1) (t).

(88)

j

j

j

ρ(j)

j j

j

ρ(j)

ρ(j)

Combining this with the previous display we get φZ (2) (λj t) = φλ Z (2) (t) = φλ Ze(2) (t) · φλ Vb (2) (t)−1 j

j

j

j

(89)

j j

j

= φZe(1) (t) · exp g̃j (t) · φλ Vb (2) (t)−1

(90)

j j

ρ(j)

= φZ (1) (t) · φVb (1) (t) · exp g̃j (t) · φλ Vb (2) (t)−1

ρ(j)

j j

ρ(j)

|

{z

=: exp(gj (t))

(91)

}

= φZ (1) (t) · exp gj (t) ,

(92)

ρ(j)

where we used that the logarithm of the characteristic function of a (possibly degenerate) Gaussian distribution is a polynomial of degree ≤ 2 and that such characteristic functions have no zeros. Hence gj is a well-defined polynomial. This proves the main claim. Now assume in addition that rank(A(1) ) = k ≤ p. Then A(1) has the left inverse A(1)† := A(1)T A(1)

−1 (1)T A .

(93)

Multiplying the original representation Eq. (78) by eTρ(j) A(1)† from the left gives (1)

Zρ(j) + eTρ(j) A(1)† E (1) = eTρ(j) A(1)† A(1) Z (1) + E (1)

(94)

= eTρ(j) A(1)† A(2) Z (2) + E (2)
= eTρ(j) A(1)† A(1) P ΛZ (2) + E (2)

(95)

= eTρ(j) P ΛZ (2) + eTρ(j) A(1)† E (2)

(97)

= eTj ΛZ (2) + eTρ(j) A(1)† E (2)

(98)

(2)

= λj Zj 25

+ eTρ(j) A(1)† E (2) ,

(96)

(99)

with Gaussian (i)

(i)

(i)

Nj := eTρ(j) A(1)† E (i) ∼ N νj , (τj )2 ,

(100)

which are independent of the corresponding Z (i) . Taking characteristic functions in (1)

(1)

Zρ(j) + Nj

(2)

(2)

= λj Zj

+ Nj

(101)

gives φZ (1) (t) · φN (1) (t) = φλ Z (2) (t) · φN (2) (t). j

j

ρ(j)

j

(102)

j

Now put σj :=

r

(2)

(1)

(2)

(τj )2 − (τj )2 , (1)

Treating the two cases τj

(2)

≥ τj

(1)

νj := νj

− νj ,

(1)

separately we obtain

and τj

(2)

≤ τj

Gj ∼ N (0, σj2 ).

(103)

φλ Z (2) (t) = φZ (1) (t) · φN (1) (t) · φN (2) (t)−1 j

j

j

ρ(j)

(104)

j

(1)

= φZ (1) (t) · φGj −νj (t),

τj

ρ(j)

(2)

(105)

≥ τj ,

φZ (1) (t) = φλ Z (2) (t) · φN (2) (t) · φN (1) (t)−1 j

ρ(j)

(106)

j

j

j

(1)

= φλ Z (2) (t) · φGj +νj (t), j

(2)

≤ τj .

(107)

(1)

≥ τj ,

(2)

(108)

(1)

≤ τj ,

(2)

(109)

τj

j

This shows that we are in one of the two cases (2) d

λ j Zj

(1)

(1)

(1)

τj

(2)

τj

= Zρ(j) + Gj − νj ,

Gj ⊥⊥ Zρ(j) ,

d

Gj ⊥⊥ Zj ,

(2)

Zρ(j) = λj Zj

+ Gj + νj ,

which is the claim, once the case distinction is expressed through gj . For that, compare the two displays above with Eq. (56). Near the origin φZ (1) does not vanish, so dividing gives ρ(j)

exp gj (t) =

φλ Z (2) (t) j

j

φZ (1) (t) ρ(j)

=

 φG −ν (t), φ

(1)

≥ τj ,

(1)

≤ τj ,

j

τj

−1 Gj +νj (t) ,

τj

j

(2)

(110)

(2)

near the origin. By Notation 2.4 the right hand side equals exp ∓ 12 σj2 t2 − iνj t in the two cases, so, applying the uniqueness part of Lemma 3.10 to φGj ∓νj and its exponent ±gj ,

gj (t) ≡ ∓ 12 σj2 t2 − iνj t

mod 2πi Z ,

gj′′ = ∓σj2 .

(111)

In particular deg(gj ) ≤ 2, the coefficient gj′′ is real, σj2 = |gj′′ |, and gj′′ ≤ 0 holds exactly in the first case and gj′′ ≥ 0 exactly in the second. Remark 5.6 (The four remaining ambiguities). Theorem 5.5 states that, under the assumption of independent non-Gaussian sources and full column rank of the mixing matrix, the distributions of the sources can be recovered up to translation, scale, permutation and componentwise additive Gaussian noise. These are the four ambiguities the theorem leaves. Each of the first three ambiguities can be removed by a normalisation convention: 26

(1) the translation ambiguity by centring all random variables, e.g. by subtracting their means (assuming these exist); (2) the scale ambiguity by rescaling, e.g. by dividing by the standard deviations (assuming these exist and are finite). Beware that this pins the scale down to a sign only when the fourth ambiguity is absent: in the setting of Theorem 5.5 the first case gives (2) (1) λ2j Var(Zj ) = Var(Zρ(j) ) + σj2 , so standardising both source vectors leaves λj = q

± 1 + σj2 . In the noiseless model σj = 0 and only the sign survives, as Corollary 7.4 records with Λ = diag(±1, . . . , ±1); (3) the permutation ambiguity by enforcing a recognisable ordering of the sources, as is done for causal models such as the linear non-Gaussian acyclic model (LiNGAM) of Shimizu et al. (2006); there the ordering is not merely a convention but is itself identified, see Corollary 7.32. The fourth ambiguity – componentwise additive Gaussian noise – is of a different nature, since it changes the law of the sources rather than just their parametrisation. Removing it is the subject of Section 6.

6

Identifiability for Gaussian-Free Independent Sources

Theorem 5.5 left us with one irreducible ambiguity: the sources are determined only up to an additive Gaussian perturbation. That ambiguity is caused by sources that still contain some Gaussian noise which could equally well be attributed to the noise vector. In this section we remove it by strengthening the hypothesis on the sources from “non-Gaussian” to “no Gaussian noise can be split off at all”.

6.1

Gaussian-free random variables

Definition 6.1 (Gaussian-free). A real-valued random variable Z is called Gaussian-free if for every decomposition d

Z = U1 + U2 ,

U1 ⊥⊥ U2 ,

(112)

neither U1 nor U2 is a non-degenerate Gaussian random variable. Equivalently: there is no d σ > 0 and no real-valued random variable Y with Y ⊥⊥ G, G ∼ N (0, σ 2 ) and Z = Y + G; that is, L(Z) has no non-degenerate Gaussian convolution factor. Remark 6.2 (Gaussian-freeness is affine invariant). Whether Z is Gaussian-free depends only on L(Z), and it is invariant under invertible affine maps: if Z is Gaussian-free, λ = ̸ 0 and c ∈ R, d

then λZ + c is Gaussian-free as well. Indeed, a decomposition λZ + c = U1 + U2 with U1 ⊥⊥ U2 d and U2 ∼ N (m, τ 2 ), τ > 0, would give Z = λ−1 (U1 −c)+λ−1 U2 with λ−1 U2 ∼ N (m/λ, τ 2 /λ2 ) again non-degenerate, contradicting Gaussian-freeness of Z. We use this silently whenever a rescaled source is called Gaussian-free, for instance in the proof of Theorem 6.12. Remark 6.3 (Gaussian-free versus non-Gaussian). A non-constant Gaussian-free random variable is automatically non-Gaussian: a non-degenerate Gaussian Z fails to be Gaussian-free, d because the trivial decomposition Z = Z +0 already exhibits a non-degenerate Gaussian factor. The converse fails, by Example 6.10 (d), and constants are Gaussian-free but degenerate – which is why non-constancy has to be assumed separately. Consequently, in Theorem 6.12 the Gaussian-free hypothesis re-derives the non-Gaussianity demanded by Theorem 5.5 for those components that are assumed Gaussian-free, and for them only the non-constancy has 27

to be carried over. In part (1) of that theorem, where only Z (1) is assumed Gaussian-free, the non-Gaussianity of Z (2) remains a genuine hypothesis and cannot be dropped: with p = 2, A(1) = [e1 ], Z (1) Rademacher and E (1) = e2 N for N ∼ N (0, 1), the same X also equals A(2) Z (2) with A(2) = [e1 , e2 ], Z (2) = [Z (1) , N ]T and E (2) = 0, and the conclusion fails because (2) Z2 is Gaussian. The interpretation is that from a Gaussian-free random variable one cannot shave off any further Gaussian noise: it carries, in this sense, the cleanest possible signal. Remark 6.4 (Terminology). The classical name for the property in Definition 6.1 is that L(Z) has no Gaussian component, or equivalently no normal component, where “component” means convolution factor. This is the vocabulary of the arithmetic of probability distributions of Linnik and Ostrovskii (1977), whose chapters are devoted to distributions “with a Gaussian component”, and it is used in exactly this sense in the ICA literature (Eriksson and Koivunen, 2006). We prefer the adjective Gaussian-free for two reasons: it does not collide with the “components” Z1 , . . . , Zk of a random vector, which we would otherwise have to call components having no components, and it does not collide with the “Gaussian components” of a Gaussian mixture model, which are mixture summands rather than convolution factors. A related but genuinely different classical notion is the Khinchin–Linnik class I0 of distributions having no indecomposable factors. Here a law is a factor of µ if µ is its convolution with some law – so every law is a factor of itself – and µ is indecomposable if it is non-degenerate and in every factorisation µ = µ1 ∗ µ2 one of µ1 , µ2 is degenerate. The two notions are logically independent: N (0, 1) lies in I0 by Cramér’s Theorem 3.20, since all its factors are Gaussian and every non-degenerate Gaussian is decomposable, but it is not Gaussian-free; whereas a non-degenerate Bernoulli law is Gaussian-free by Corollary 6.9 below yet is itself indecomposable – a convolution of two laws with |A| and |B| support points has at least |A| + |B| − 1 of them, so two non-degenerate factors would force at least three – and hence is not in I0 . To make Definition 6.1 quantitative we measure how much Gaussian noise can be split off. Definition 6.5 (Splittable Gaussian scales). For a real-valued random variable Z and s ∈ R≥0 set 



gs (t) := φZ (t) exp 12 s2 t2 ,

(113)

t ∈ R,

and put S(Z) := s ∈ R≥0 : gs is a characteristic function ,

(114)

σmax (Z) := sup S(Z) ∈ [0, ∞].

(115)



as well as

We call σmax (Z) the maximal Gaussian scale of Z. The name is not standard. It is justified by the following reading: s ∈ S(Z) if and only if a N (0, s2 ) factor can be split off Z, so σmax (Z) is the supremum of the standard deviations of the Gaussians that Z contains as convolution factors – and by Lemma 6.6 (iii) below that supremum is attained, so it really is the largest such standard deviation. Lemma 6.6 (Structure of S(Z)). Let Z be a real-valued random variable. Then: (i) for every s ∈ S(Z) we have the tail bound 



|φZ (t)| ≤ exp − 12 s2 t2 , 28

t ∈ R;

(116)

(ii) σmax (Z) < ∞, always; (iii) S(Z) = 0, σmax (Z) ; in particular the supremum is attained; 



(iv) Z is Gaussian-free if and only if σmax (Z) = 0. Proof. (i) If s ∈ S(Z) then gs is a characteristic function, so |gs (t)| ≤ 1 by Proposition 3.2 (i), which is Eq. (116). (ii) Suppose S(Z) were unbounded and pick sn ∈ S(Z) with sn → ∞. Fix t = ̸ 0. By 1 2 2 Eq. (116), |φZ (t)| ≤ exp(− 2 sn t ) → 0, so φZ (t) = 0. Thus φZ vanishes on R \ {0} while φZ (0) = 1, contradicting the continuity of φZ (Proposition 3.2 (ii)). (iii) First, 0 ∈ S(Z), and S(Z) is downward closed: if s ∈ S(Z) and 0 ≤ s′ ≤ s, then 



gs′ (t) = gs (t) · exp − 12 s2 − s′2 t2

(117)

is a product of the characteristic function gs with that of a N (0, s2 − s′2 ) variable, hence a characteristic function by Eq. (9). So S(Z) is an interval containing 0, and it is bounded by (ii). It remains to see that σ := σmax (Z) itself lies in S(Z). Choose sn ∈ S(Z) with sn ↑ σ. For every t ∈ R we have the pointwise convergence gsn (t) → gσ (t), and gσ is continuous at t = 0 because φZ and exp are. By Lévy’s continuity Theorem 3.7 (ii), gσ is a characteristic function, i.e. σ ∈ S(Z). (iv) If σmax (Z) = σ > 0 then, writing φZ (t) = gσ (t) exp(− 21 σ 2 t2 ) with the characteristic function gσ from (iii), we exhibit a non-degenerate Gaussian convolution factor of L(Z), so Z d

is not Gaussian-free. Conversely, if Z is not Gaussian-free, say Z = U + V with U ⊥⊥ V and V ∼ N (µ, τ 2 ), τ > 0, then 



φZ (t) exp 21 τ 2 t2 = φU (t) exp iµt = φU +µ (t)

(118)

is a characteristic function, so τ ∈ S(Z) and σmax (Z) ≥ τ > 0. Lemma 6.6 puts us in a position to prove the following, which turns the informal wish of the introduction to this section into a theorem: every real-valued random variable splits into a Gaussian-free signal and independent Gaussian noise, and it does so in essentially one way. Theorem 6.7 (Gaussian splitting). Let Z be a real-valued random variable and put σ := σmax (Z) ∈ [0, ∞). Then there exist a Gaussian-free random variable Y and a random variable G ∼ N (0, σ 2 ) with Y ⊥ ⊥ G such that d

Z = Y + G.

(119) d

The decomposition is maximal and unique up to a translation: if Z = Y ′ + G′ with Y ′ ⊥⊥ G′ , d Y ′ Gaussian-free and G′ ∼ N (m, τ 2 ), then τ = σ and Y ′ = Y − m. Proof. By Lemma 6.6 (iii) we have σ ∈ S(Z), so that t 7→ gσ (t) = φZ (t) exp( 12 σ 2 t2 ) is the characteristic function of some random variable Y , and 



φZ (t) = gσ (t) · exp − 12 σ 2 t2 = φY (t) · φG (t)

(120)

for G ∼ N (0, σ 2 ) chosen independent of Y ; this is Eq. (119) by Eq. (9) and Theorem 3.4 (i). d

We allow the degenerate case σ = 0, where G = 0 and Y = Z.

29

d

To see that Y is Gaussian-free, suppose Y = U + V with U ⊥⊥ V and V ∼ N (µ, τ 2 ), τ > 0. Then, as in the proof of Lemma 6.6 (iv), φY (t) exp( 12 τ 2 t2 ) is a characteristic function, hence so is 







φZ (t) exp 12 σ 2 + τ 2 t2 = φY (t) exp 21 τ 2 t2 , so that s̄ :=

√

(121)

σ 2 + τ 2 ∈ S(Z). But then σmax (Z) = σ <

p

σ 2 + τ 2 = s̄ ∈ S(Z),

(122)

contradicting the definition of σmax . d For the uniqueness statement, let Z = Y ′ + G′ be as stated. Then φZ (t) exp( 12 τ 2 t2 ) = φY ′ (t) exp(imt) is a characteristic function, so τ ∈ S(Z) and hence τ ≤ σ by Lemma 6.6 (iii). If we had τ < σ, then 







φY ′ (t) exp 12 σ 2 − τ 2 t2 = φZ (t) exp 12 σ 2 t2 exp(−imt) = φY −m (t)

(123)

would be a characteristic function, exhibiting a N (0, σ 2 − τ 2 ) factor of L(Y ′ ) and contradicting that Y ′ is Gaussian-free. So τ = σ, and then φY ′ (t) = φZ (t) exp( 12 σ 2 t2 ) exp(−imt) = φY −m (t). Remark 6.8 (A dichotomy, not a trichotomy). For every real-valued random variable Z exactly one of the following holds. (1) σmax (Z) = 0. Equivalently, Z is Gaussian-free. This is the degenerate case of Theod

rem 6.7, with G = 0 and Y = Z. d

(2) 0 < σmax (Z) < ∞. Here Theorem 6.7 gives a proper decomposition Z = Y + G with a non-degenerate Gaussian G, a Gaussian-free Y , and Y ⊥⊥ G. One might expect a third case, in which arbitrarily large Gaussian factors can be shaved off; Lemma 6.6 (ii) shows that this cannot happen. The reason is Eq. (116): splitting off a N (0, s2 ) factor forces the characteristic function to decay at least as fast as exp(− 21 s2 t2 ), and no characteristic function can decay faster than every Gaussian without vanishing identically off the origin. The tail bound Eq. (116) is also the most convenient practical criterion, alongside a support argument. Corollary 6.9 (Two sufficient criteria). Let Z be a real-valued random variable. Each of the following implies that Z is Gaussian-free. (i) supp L(Z) ̸= R; this holds in particular if Z is bounded, or bounded from one side. (ii) exp( 12 s2 t2 )|φZ (t)| is unbounded in t for every s > 0; in particular, this holds whenever
 2 lim inf |t|→∞ − log|φZ (t)| t = 0. Proof. (i) The support of a convolution of two laws is the closure of the sum of their supports. A non-degenerate Gaussian law has support R, so any law with a non-degenerate Gaussian convolution factor has support R. (ii) is the contrapositive of Eq. (116), together with Lemma 6.6 (iv). Example 6.10 (Which laws are Gaussian-free?). (a) By Corollary 6.9 (i), every law whose support is not all of R is Gaussian-free: the uniform, Rademacher and three-point laws of Table 1, every Bernoulli law, and also the exponential, Gamma, Poisson and χ2 laws, which are supported on a half line. 30

(b) By Corollary 6.9 (ii), every law whose characteristic function decays more slowly than any Gaussian is Gaussian-free. This covers the Laplace law, with φZ (t) = (1 + b2 t2 )−1 decaying polynomially; the Cauchy law, with φZ (t) = exp(−γ|t|); and the Student tν laws, whose characteristic functions decay exponentially. So the standard super-Gaussian sources of ICA are all Gaussian-free. (c) A non-degenerate Gaussian law is of course not Gaussian-free, with σmax N (µ, σ 2 ) = σ.

(d) The instructive example is Z := R + G with R Rademacher, G ∼ N (0, 1) and R ⊥⊥ G, i.e. the Gaussian mixture 12 N (−1, 1) + 12 N (1, 1). Here φZ (t) = cos(t) exp(− 21 t2 ), so φZ (t) exp( 12 s2 t2 ) = cos(t) exp( 12 (s2 −1)t2 ) is unbounded for s > 1 and hence σmax (Z) = 1, with Gaussian-free part R – determined, as Theorem 6.7 says, only up to a translation. By Remark 3.18 (ii) the variable Z is non-Gaussian, so it is an admissible source for Theorem 5.5, but it is not Gaussian-free and therefore not admissible for Theorem 6.12 below. This shows that the hypothesis of this section is strictly stronger than that of Section 5. Remark 6.11 (The infinitely divisible case). For readers familiar with the Lévy–Khintchine representation (Kallenberg, 2021, Chapter 7), there is a clean characterisation. Let Z be infinitely divisible with triplet (a, ν, γ), i.e. φZ (t) = exp



iγt − 12 at2 +

Z 

e

itx





− 1 − itx1{|x|≤1} ν(dx) ,

(124)

R

with a ≥ 0 and a Lévy measure ν satisfying min(1, x2 ) ν(dx) < ∞. Then Z is Gaussian-free if and only if a = 0, and in fact σmax (Z)2 = a. Indeed, if a > 0 then splitting the factor exp(− 12 at2 ) off φZ leaves the characteristic √ function of the infinitely divisible law with triplet (0, ν, γ), so a ∈ S(Z). Conversely, assume a = 0 and let s ∈ S(Z). Taking moduli, R

− log|φZ (t)| =

Z

1 − cos(tx) ν(dx),

(125)

R

and this is o(t2 ) as |t| → ∞: on {|x| ≤ 1} we have (1−cos(tx))/t2 ≤ x2 /2, which is ν-integrable there and tends to 0 pointwise, so dominated convergence applies; on {|x| > 1} the integral is bounded by 2ν(|x| > 1) < ∞. On the other hand Eq. (116) gives − log|φZ (t)| ≥ 12 s2 t2 . Comparing the two forces s = 0, so σmax (Z) = 0. The same comparison gives the identity σmax (Z)2 = a for general a ≥ 0, not just the two implications. Indeed, without assuming a = 0, − log|φZ (t)| = 12 at2 +

Z R

1 − cos(tx) ν(dx) = 21 at2 + o t2 ,

|t| → ∞,

(126)

by the estimate just given, while Eq. (116) bounds the left hand side below by 12 s2 t2 for every s ∈ S(Z). Dividing by t2 and letting |t| → ∞ yields s2 ≤ a, hence σmax (Z)2 ≤ a; combined √ with a ∈ S(Z) this is equality. Note that this argument never assumes that the cofactor is itself infinitely divisible – which matters, since factors of infinitely divisible laws need not be infinitely divisible.

6.2

Identifiability in the presence of additive Gaussian noise

In the following we explicitly assume that the sources of our ICA model have been separated according to Theorem 6.7, so that they are mutually independent and Gaussian-free. Note that Gaussian random variables are considered noise and are usually not properly identifiable 31

in ICA models anyway, cf. Remark 3.27. Assuming that the proper signals in an ICA model are Gaussian-free is therefore a reasonable assumption, even though it is more restrictive than the usual assumption of mere non-normality – and by Example 6.10 (b) it is satisfied by all the standard super-Gaussian source models. The reward for this stronger assumption is a stronger identifiability result: it identifies the sources inside a mixture even in the presence of additive and dependent Gaussian noise. Theorem 6.12 (Identifiability for independent Gaussian-free sources with additive Gaussian noise). Let X ∈ Rp be a random vector and assume that we have two representations of X satisfying all the assumptions of Theorem 5.5: A(1) Z (1) + E (1) = X = A(2) Z (2) + E (2) .

(127)

Assume furthermore that A(1) (equivalently, by Theorem 5.5, A(2) ) has full column rank, rank(A(1) ) = k ≤ p. (1) If all components of Z (1) are Gaussian-free, then there exist a permutation matrix P ∈ Rk×k , an invertible diagonal matrix Λ ∈ Rk×k and a (possibly degenerate) Gaussian random vector H ∼ N (ξ, ∆) with mean vector ξ ∈ Rk and diagonal, possibly degenerate, covariance matrix ∆ ∈ Rk×k , independent of Z (1) and of E (2) , such that A(2) = A(1) P Λ, Z

(128)

(2) d

= Λ−1 P −1 Z (1) + H ,

H ⊥⊥ Z (1) , E

(2)

,

d

E (1) = A(1) H + E (2) ,

(130)

µ(1) = A(1) ξ + µ(2) , Σ

(1)

=A

(1)

∆A

(1)T

+Σ

(129) (131)

(2)

(132)

.

(2) If, in addition, all components of Z (2) are Gaussian-free as well, then H = ξ and ∆ = 0, and we further get d

d

Z (2) = Λ−1 P −1 Z (1) + ξ ,

E (2) = E (1) − A(1) ξ,

Σ(2) = Σ(1) .

(133)

Proof. Part (1). By Theorem 5.5 we already have A(2) = A(1) P Λ, and for each j = 1, . . . , k separately one of the two cases (2)

λ j Zj

d

(1)

+ νj = Zρ(j) + Gj , (1)

d

(2)

Zρ(j) = λj Zj

(1)

(134)

(2)

(135)

Gj ⊥⊥ Zρ(j) ,

+ νj + Gj ,

Gj ⊥⊥ Zj ,

(1)

with some νj ∈ R and Gj ∼ N (0, σj2 ). Since Zρ(j) is Gaussian-free, the second case exhibits a (1)

Gaussian convolution factor of L(Zρ(j) ) and can therefore only occur with σj2 = 0, i.e. Gj = 0; but then it is subsumed by the first case. We may therefore assume the first case for all j = 1, . . . , k, which we rearrange into (2) d

eTj Z (2) = Zj

(1)

(1)

−1 −1 = λ−1 j Zρ(j) − λj νj + λj Gj ,

Gj ⊥⊥ Zρ(j) ,

(136)

= eTj Λ−1 P −1 Z (1) − eTj Λ−1 ν + eTj Λ−1 G,

(137)

= eTj Λ−1 P −1 Z (1) − P ν + P G ,

(138)

= eTj Λ−1 P −1 Z (1) + H ,

(139)

32

where we used the abbreviations ν := [ν1 , . . . , νk ]T ,

ξ := −P ν,

(140)

G := [G1 , . . . , Gk ] ∼ N (0, Γ),

Γ := diag(σ12 , . . . , σk2 ), T

(141)

T

H := ξ + P G ∼ N (ξ, ∆),

∆ := P ΓP .

(142)

Note that the Gj were constructed componentwise, so we may and do take G to have independent components; this implies that ∆ = P ΓP T is again a diagonal matrix, and hence that H has independent components as well. Together with the fact that the components of Z (1) and of Z (2) are each mutually independent, we may gather all components into the single distributional equation d

Z (2) = Λ−1 P −1 Z (1) + H ,

H ⊥⊥ Z (1) .

(143)

By the continuous mapping theorem we may multiply Eq. (143) by A(2) = A(1) P Λ and obtain d

A(2) Z (2) = A(1) Z (1) + A(1) H,

H ⊥⊥ Z (1) ,

(144)

which in terms of characteristic functions reads φA(2) Z (2) (t) = φA(1) Z (1) (t) · φA(1) H (t).

(145)

Plugging this into the characteristic functions of the original model equation, φA(2) Z (2) (t) · φE (2) (t) = φA(1) Z (1) (t) · φE (1) (t),

(146)

we obtain an identity that we would like to divide by φA(1) Z (1) . Here we may not simply cancel the factor φA(1) Z (1) , since a characteristic function may well have zeros: for p = k = 1, A(1) = 1 and Z (1) ∼ Unif[−1, 1] – an admissible Gaussian-free, non-Gaussian source – one has φZ (1) (t) = sin(t)/t, which vanishes at t = π. We argue locally instead. By continuity and φA(1) Z (1) (0) = 1 there is a δ > 0 with φA(1) Z (1) = ̸ 0 on the ball Bδ (0) ⊆ Rp , and cancelling there gives φA(1) H (t) · φE (2) (t) = φE (1) (t),

t ∈ Bδ (0).

(147)

Both sides of Eq. (147) are characteristic functions of (possibly degenerate) multivariate Gaussian distributions, hence of the form exp(q(t)) with polynomials q of degree ≤ 2 and q(0) = 0; here we may in addition take H independent of E (2) : all assertions here are statements about laws, so we may enlarge the underlying probability space (replacing Ω by Ω × [0, 1]) and realise an H with law N (ξ, ∆) independent of the pair (Z (1) , E (2) ). Writing the two sides as exp(q1 ) and exp(q2 ), we get exp(q1 − q2 ) ≡ 1 on Bδ (0), so the continuous function (q1 − q2 )/(2πi) takes values in Z on the connected set Bδ (0) and vanishes at the origin; hence q1 = q2 on Bδ (0). Two polynomials that agree on a non-empty open subset of Rp agree everywhere, so q1 = q2 identically. Hence Eq. (147) in fact holds for every t ∈ Rp , and Theorem 3.4 gives d

A(1) H + E (2) = E (1) ,

H ⊥⊥ Z (1) , E (2) ,

(148)

and comparing mean vectors and covariance matrices yields A(1) ξ + µ(2) = µ(1) ,

(149)

A(1) ∆A(1)T + Σ(2) = Σ(1) .

(150)

33

This proves all claims of part (1). (2) Part (2). If the components Zj are Gaussian-free as well, then so are the affinely (2)

transformed variables λj Zj (2)

λ j Zj

+ νj by Remark 6.2, and the only possibility left is
d

(1)

(1)

+ νj = Zρ(j) + Gj ,

Gj ⊥⊥ Zρ(j) ,

(151)

with σj2 = 0 and thus Gj = 0. This shows H = ξ and ∆ = 0, and everything else follows from part (1). Remark 6.13 (What Theorem 6.12 buys us). Compared with Theorem 5.5, the stronger assumption on the sources removes the fourth ambiguity of Remark 5.6 – but only part (2) removes it entirely. Part (1) makes it one-sided: any residual Gaussian must sit in Z (2) and is exactly compensated in E (2) , which is what the vector H records. It is not vacuous. Take p = k = 1, A(1) = A(2) = [1], Z (1) Rademacher (Gaussian-free) with E (1) ∼ N (0, 1), against Z (2) = Z (1) + E (1) (non-constant and, by Remark 3.18 (ii), non-Gaussian) with E (2) = 0: both representations satisfy every hypothesis of part (1), yet L(Z (2) ) is absolutely continuous while L(Z (1) ) has two atoms, so no affine map carries one to the other. Three aspects deserve emphasis. (i) The Gaussian noise vectors E (i) are not assumed to have independent components, nor a non-degenerate covariance matrix. They may be arbitrarily dependent across the p observed coordinates, and they may be partially deterministic. Nevertheless, under the hypothesis of part (2), the sources are pinned down up to permutation, scale and translation. (ii) Part (2) is a genuine two-sided statement: as soon as both candidate source vectors are Gaussian-free, not only the sources but also the noise laws agree, Σ(1) = Σ(2) , so the whole model (A, L(Z), L(E)) is identified up to the unavoidable permutation, scaling and shift. (iii) The Gaussian-free assumption is not a restriction on which observations can be modelled, only on how the model is parametrised. Indeed, by Theorem 6.7 every source splits as d Zj = Yj + Gj with Yj Gaussian-free and Gj ∼ N (0, σj2 ) independent of Yj , and moving the Gaussian part into the noise is the trivial direction of Lemma 5.2: writing aj for the j-th column of A, d

aj Zj = aj Yj + aj Gj ,

aj Gj ∼ N 0, σj2 aj aTj ⊥⊥ Yj ,

(152)

so that aj Gj may simply be added to E – which is why we allow E to have an arbitrary, possibly degenerate covariance matrix. Since the Zj are mutually independent, the pairs (Yj , Gj ) may be taken jointly independent, so all k Gaussian parts can be moved at once. (Lemma 5.2 proves the harder converse direction, moving Gaussian noise out of E and into extra columns of A.) One caveat: if Zj was itself Gaussian, then its Gaussian-free part Yj is a constant, and the corresponding column must then be deleted and absorbed into the offset µ, exactly as in Step 1 of Remark 4.1, before Theorem 6.12 can be applied to the reduced model. What the assumption buys is that this splitting and pruning has been carried out, so that no Gaussian mass is left sitting ambiguously between Z and E. In practice, part (1) is the statement one applies when a candidate solution is compared against the ground truth, and part (2) is the statement one applies when both are outputs of a procedure that is guaranteed to return Gaussian-free sources.

34

7

Complete Noiseless ICA: Estimation by Equivariant Gradient Descent

The previous sections answered the question of what is identifiable. This section answers the complementary question of how one estimates it, in the special case that is by far the most common in applications: a square invertible mixing matrix and no noise. We first specialise the identifiability theory to that case (Section 7.1), then set up the maximum likelihood objective (Section 7.2), then explain what the customary “preconditioner” in the gradient step really is (Section 7.3), and finally analyse when the resulting algorithm converges to a separating solution (Section 7.4). Example 7.1 (The cocktail party problem). Two loudspeakers in a room play different pieces of music, modelled as real-valued signals Z1 (t) and Z2 (t) indexed by time t. Two microphones at different positions each record a mixture of the two, X1 (t) = a11 Z1 (t) + a12 Z2 (t),

X2 (t) = a21 Z1 (t) + a22 Z2 (t).

(153)

The task is to recover the two pieces of music Zj (t) from the two recorded mixtures Xi (t). What makes this more than a matrix inversion is that both the mixing matrix A = (aij ) and the sources Zj are unknown; all that is assumed is that the two sources are statistically independent. Corollary 7.4 below says that this is enough – provided at most one source is Gaussian.

7.1

The complete noiseless model and its identifiability

Assumption 7.2 (Complete noiseless ICA model). We observe X(1), . . . , X(T ) in Rp , all with the same law, and assume: (i) Independent sources. There are random vectors Z(t) = [Z1 (t), . . . , Zk (t)]T in Rk whose components {Z1 (t), . . . , Zk (t)} are mutually independent for each fixed t, and almost surely non-constant. Nothing is assumed about the dependence across t. (ii) Linear noiseless mixing. X(t) = A Z(t) for all t, with a matrix A ∈ Rp×k that does not depend on t. (iii) Completeness. p = k and A is invertible. We write W := A−1 and denote the rows of W by wT1 , . . . , wTk . (iv) Non-Gaussianity. At most one of the components Z1 (t), . . . , Zk (t) is Gaussian. For the estimation part we shall in addition normalise: (v) Centring and scale. E[Zj (t)] = 0 for all j and t, and the scale of each Zj is fixed either by Var(Zj ) = 1 or by the fixed-point normalisation Eq. (182); see Remark 7.3. Since the law of X(t) does not depend on t, we drop the argument and write X = AZ for a generic sample whenever only the common law matters. Remark 7.3 (The two scale conventions). Rescaling the sources is without loss of generality: it
is the reparametrisation X = (AD) D−1 Z with D diagonal and invertible, which preserves Assumption 7.2 (i)–(iv). Centring them, however, requires centring the data as well, since Assumption 7.2 (ii) carries no offset; that is step (1) of Remark 7.9. Two scale conventions appear below, and they are not the same. The first is Var(Zj ) = 1, used wherever only secondorder quantities matter (Remark 7.9); it requires the extra assumption E[Zj2 ] < ∞, which Assumption 7.2 does not otherwise impose. The second is the fixed-point scale of Lemma 7.18, which is the one Theorem 7.20 needs and which is generically not the unit-variance scale. 35

The identifiability statement usually quoted in the data-science literature is the following. It is exactly the specialisation of Theorem 4.2 to a square invertible mixing matrix, with the “at most one Gaussian” clause handled by an extra argument at the end of the proof. Corollary 7.4 (Identifiability of the complete noiseless model; Comon, 1994, Theorem 11). Let X be a random vector in Rk with two representations A(1) Z (1) + µ(1) = X = A(2) Z (2) + µ(2) ,

(154)

where for i = 1, 2: (i) A(i) ∈ Rk×k is invertible and µ(i) ∈ Rk ; (ii) the components of Z (i) are mutually independent and almost surely non-constant; and, for i = 1 only, (iii) at most one component of Z (1) is Gaussian. Then there exist a permutation matrix P ∈ Rk×k , an invertible diagonal matrix Λ ∈ Rk×k and a vector c ∈ Rk such that A(2) = A(1) P Λ,

Z (1) = P Λ Z (2) + c

(155)

almost surely.

If moreover both source vectors are centred with unit variances, then Λ = diag(±1, . . . , ±1) and c = 0: the sources are determined up to permutation and sign, and the columns of the mixing matrix up to permutation and sign. Proof. An invertible matrix has non-zero, pairwise non-proportional columns, so both representations satisfy the hypotheses of Theorem 4.2. (i) Step 1: matching the non-Gaussian columns. Let Ni := {j : Zj is non-Gaussian}. Let (2)

l ∈ N2 . By Theorem 4.2 (1) the column al

(1)

must be proportional to some column aρ(l) of A(1) ,

(1)

and by Theorem 4.2 (2) the source Zρ(l) is then non-Gaussian, i.e. ρ(l) ∈ N1 . The index ρ(l) is unique because the columns of A(1) are pairwise non-proportional, and the map ρ : N2 → N1 is (2) (1) (2) injective: if ρ(l) = ρ(l′ ) then al ∝ aρ(l) ∝ al′ , contradicting pairwise non-proportionality of the columns of A(2) . Hence |N2 | ≤ |N1 |, and by symmetry |N1 | = |N2 | =: n, with ρ : N2 → N1 a bijection. In particular the two representations have Gaussian sources in equal number, k − n each. This is why hypothesis (iii) needs to be imposed on one of the two representations only: assuming it for i = 1 gives n ≥ k − 1, and hence at most one Gaussian component of Z (2) as well. In either case n ∈ {k − 1, k}. Step 2: reduction to a single matrix. Put C := (A(1) )−1 A(2) , which is invertible, so that (2) A = A(1) C; it therefore suffices to prove C = P Λ. Multiplying Eq. (154) by (A(1) )−1 gives, almost surely, Z (1) = C Z (2) + d,

d := A(1)

−1

µ(2) − µ(1) ,

(156)

so a factorisation C = P Λ gives the second assertion as well, with c := d. For l ∈ N2 we have (2) (1) al = λl aρ(l) with λl = ̸ 0, hence Cel = λl eρ(l) : every column of C indexed by N2 is a non-zero multiple of a standard basis vector. Step 3: the case n = k. Here N2 = {1, . . . , k} and ρ is a permutation, so every column of C is of that form and C = P (ρ)Λ with Λ = diag(λ1 , . . . , λk ) invertible. Step 4: the case n = k − 1. Let l0 ∈ / N2 and k0 ∈ / N1 be the two exceptional indices, so (2) (1) that Zl0 and Zk0 are Gaussian, and non-degenerate because they are non-constant. By Step 36

2 every column of C except the l0 -th is a non-zero multiple of a standard basis vector, the corresponding indices ρ(l) exhausting {1, . . . , k} \ {k0 }. Reading off rows, Eq. (156) becomes (1)

(2)

Zk 0 = β Z l 0 + d k 0 , (1)

Zj

(2)

β := Ck0 l0 , (2)

= λρ−1 (j) Zρ−1 (j) + βj Zl0 + dj ,

βj := Cjl0 ,

(157) j ̸= k0 .

(158)

Invertibility of C forces β = ̸ 0. Suppose, for contradiction, that βj ̸= 0 for some j ̸= k0 . (2) (2) Abbreviate U := Zρ−1 (j) , λ := λρ−1 (j) ̸= 0 and G := Zl0 ∼ N (m, τ 2 ) with τ > 0; note U ⊥⊥ G. (1)

The components Zj

(1)

and Zk0 are independent, so for all s, t ∈ R, φU (sλ) φG (sβj + tβ) = φU (sλ) φG (sβj ) φG (tβ),

(159)

where the constants dj , dk0 have cancelled. Since φU is continuous with φU (0) = 1, it is non-zero in a neighbourhood of the origin, so we may cancel it there and take logarithms of the nowhere vanishing Gaussian characteristic functions. The two exponents are the quadratics of Notation 2.4; their difference is continuous, takes values in 2πi Z on the connected set where the identity holds, and vanishes at (0, 0), hence is identically zero. This gives, for all small s, t, im sβj + tβ − 12 τ 2 sβj + tβ

2

= im sβj + tβ − 12 τ 2



sβj

2

+ tβ

2 

,

(160)

i.e. τ 2 βj β st = 0 for all small s, t. As τ > 0 and β ̸= 0 this forces βj = 0, a contradiction. Hence βj = 0 for all j = ̸ k0 , so the l0 -th column of C equals βek0 and C = P Λ for a permutation matrix P and an invertible diagonal Λ. With Eq. (156) this proves the first display, and A(2) = A(1) C = A(1) P Λ. The normalised form. If both source vectors are centred, taking expectations in Z (1) = P ΛZ (2) + c gives c = 0. If in addition all variances equal 1, then each component of Z (1) equals λ times a component of Z (2) , where λ is the corresponding diagonal entry of Λ; taking variances gives 1 = λ2 , so λ = ±1. Remark 7.5 (Two Gaussian sources already destroy identifiability). Hypothesis (iii) of Corollary 7.4 cannot be weakened. Let k = 2 and Z ∼ N (0, I2 ), A(1) = I2 , so X = Z. For any orthogonal Q ∈ O(2) put A(2) := QT and Z (2) := QZ. Then Z (2) ∼ N (0, I2 ) again has independent components and A(2) Z (2) = QT QZ = X, so the whole rotation group is compatible with the observed law; cf. Remark 3.27. The same construction applies inside any two-dimensional Gaussian block, which is why at most one Gaussian source is allowed. This degeneracy reappears analytically in Proposition 7.22 below: two Gaussian sources force the stability quantities to satisfy ζj ζl = 1, so the linearised dynamics acquires a neutral direction. With the true score the implication reverses as well, by Corollary 7.23 (ζj ≥ 1 with equality only in the Gaussian case). When in addition the two source variances agree, the neutral direction is the antisymmetric one, i.e. the infinitesimal version of the rotation Q above; for unequal variances it is (σj2 , −σl2 ) instead. Remark 7.6 (Relation to the general theory). Corollary 7.4 is the noiseless, square, full-rank corner of the theory of Sections 5 and 6. Comparing: Theorem 5.5 allows a rectangular A and additive Gaussian noise but must then leave the sources undetermined up to an extra additive Gaussian term; Theorem 6.12 removes that term at the price of the Gaussian-free hypothesis. In the complete noiseless model there is no noise to hide anything in, so the weakest hypothesis – at most one Gaussian source – already suffices, and the conclusion holds almost surely rather than merely in distribution.

37

7.2

The maximum likelihood objective

From here on we work under Assumption 7.2 and additionally assume that each source Zj has a strictly positive Lebesgue density pZj . The model is parametrised by the unmixing matrix W = A−1 together with k model source densities p1 , . . . , pk , which are a modelling choice and need not equal the true densities pZ1 , . . . , pZk : θ = W, p1 , . . . , pk ,

pZ (z) = p1 (z1 ) · · · pk (zk ).

(161)

Proposition 7.7 (Model density and log-likelihood). Under Assumption 7.2, the density of X = W −1 Z induced by the parameter θ is qθ (x) = |det W |

k Y

pj wTj x ,

(162)

x ∈ Rk .

j=1

Consequently the population log-likelihood, with X distributed according to the true law q, is L(θ) := EX∼q log qθ (X) = log|det W | + 



k X

EX∼q log pj wTj X , 



(163)

j=1

and its empirical counterpart is LT (θ) :=

T T X k

1X 1X log pj wTj X(t) . log qθ X(t) = log|det W | + T t=1 T t=1 j=1

(164)

Moreover, if q has a density and finite differential entropy, L(θ) = − KL q ∥ qθ − H(q),

H(q) := −EX∼q [log q(X)].

(165)

Proof. Eq. (162) is the change-of-variables formula for the diffeomorphismRz = W x with Jacobian determinant det W : for a Borel set S, P[X ∈ S] = P[Z ∈ W S] = W S pZ (z)dz = R Q S pZ (W x)|det W | dx. Eqs. (163) and (164) follow by taking logarithms and using log j = P j log. For Eq. (165), KL(q∥qθ ) = Eq [log q] − Eq [log qθ ] = −H(q) − L(θ). Since H(q) does not depend on θ, Eq. (165) says that maximising L is the same as minimising the Kullback–Leibler divergence from the true law to the model. Remark 7.8 (The empirical objective does not need independence over time). Eq. (164) is the exact log-likelihood only if the samples X(1), . . . , X(T ) are independent. Assumption 7.2 (i) deliberately does not assume this: in the cocktail-party application consecutive samples of an audio signal are strongly dependent. Nothing below breaks, because we never use Eq. (164) as a likelihood. We use it as an M-estimation criterion whose population version Eq. (163) depends only on the common law of X(t), and all our statements – the stationarity conditions of Proposition 7.17 and the stability analysis of Theorem 7.20 – are statements about that population criterion. Provided (Z(t))t is stationary and ergodic – enough for a law of large numbers, so that LT → L almost surely – dependence across t affects only the rate of that convergence, not the location of the optima. Some such hypothesis is needed: if Z(t) ≡ Z(1) for every t, then all five conditions of Assumption 7.2 hold, yet LT = log qθ (x(1)) for every T , and this is unbounded above (send a row of W towards a direction orthogonal to x(1) while log|det W | → ∞). Remark 7.9 (Centring and whitening). Two preprocessing steps are customary.

38

(1) Centring: replace X(t) by X(t) − E[X], so that E[X] = 0; since A is invertible this is equivalent to E[Z] = 0, and being a deterministic shift of each Zj separately it preserves Assumption 7.2 (i). In practice one subtracts the sample mean X̄, which is an approximation to this whose validity again needs the law of large numbers of Remark 7.8: subtracting X̄ mixes the time points and can destroy the independence of the components at a fixed t. (2) Whitening: assume E∥Z∥2 < ∞ (so that ΣX := Cov(X) exists and is invertible) and −1/2 −1/2 replace X by X ′ := ΣX X, whose mixing matrix is A′ := ΣX A. With centred unit-variance sources this gives Ik = Cov(X ′ ) = A′ Cov(Z) (A′ )T = A′ (A′ )T ,

(166)

−1/2

with the symmetric square root ΣX , which exists and is invertible because the Zj are non-constant with finite variance; that is, the mixing matrix of the whitened data is orthogonal. For algorithms that constrain W to O(k) this reduces the number of free parameters from k 2 to the k(k − 1)/2 dimensions of O(k), and it is what most ICA algorithms exploit. The relative-gradient update Eq. (172) does not preserve orthogonality, so it gains no such reduction. The algorithm developed below needs neither step in order to be well defined. By Proposition 7.15 its trajectory depends on A only through the initial global system matrix R0 = W0 A, so whitening cannot change the shape of the dynamics – but it does change R0 , replacing it by −1/2 an orthogonal matrix when W0 = ΣX , and a better conditioned starting point is a genuine practical benefit. Centring is advisable in any case, because Theorem 7.20 assumes E[Z] = 0. Proposition 7.10 (Gradient of the log-likelihood). Fix source densities p1 , . . . , pk that are differentiable and strictly positive, and write ηj := log pj

′

p′j = pj

(167)

for the model score of the j-th source (Remark 7.11). Put C := W X, and write η(C) for the vector with components ηj (Cj ). Then, under integrability conditions permitting differentiation under the expectation,
  ∂L = E ηi (Ci ) Xj + W −1 ji , ∂Wij

that is

∇W L = W −T + E η(C)XT . 



(168)

Proof. Differentiating Eq. (163) term by term, ∂Wij log pl (wTl X) = δli ηl (wTl X) Xj by the chain rule, and summing over l leaves the single term ηi (Ci )Xj . For the determinant, Jacobi’s formula gives ∂Wij log|det W | = (W −1 )ji , i.e. ∇W log|det W | = W −T . Remark 7.11 (Two scores: the model score and the true score). From here to the end of the section, two different objects compete for the name “score”, and almost every subtlety in what follows is a statement about their relationship. We keep them notationally apart. • ηj := (log pj )′ of Eq. (167), the score of the model density pj . This is a modelling choice, fixed by the user before the algorithm is run, and it is the function the algorithm actually evaluates – in Algorithm 7.28 it is applied to the current estimate Cj of the j-th source. We also call it the nonlinearity. Nothing forces it to be the score of the true source density. The three used in practice are in fact scores of perfectly good densities – just not, in general, of the right one: − tanh is the score of the hyperbolic 2 secant density π1 sech, with kurt = +2; tanh −id is the score of ∝ cosh(u)e−u /2 , i.e. of 39

4

the mixture 12 N (−1, 1) + N (1, 1) , with kurt = − 12 ; and −u3 is the score of ∝ e−u /4 , with kurt ≈ −0.81. Choosing ηj thus amounts to positing a source law, and the three choices posit a leptokurtic, a mildly platykurtic and a platykurtic one respectively. Not 4 every nonlinearity is a score, though: +u3 would require the non-integrable eu /4 .

′

• ηj⋆ := log pZj , the score of the true density pZj of the source Zj . This is a property of the data generating process, and it is unknown. The maximum likelihood derivation of Section 7.2 would use ηj⋆ if it could; the algorithm uses ηj because it must. Theorem 7.20 is a statement about an arbitrary ηj , and this is exactly what makes it useful: it says how far ηj may deviate from ηj⋆ before the separating solution stops attracting. Corollary 7.23 is the special case ηj = ηj⋆ . Note that ηj is evaluated along the trajectory at Cj = (RZ)j , a mixture of the sources, whereas ηj⋆ belongs to Zj alone; the two are brought into contact only through their values on the sources, i.e. at the fixed point R = Ik , where Cj = Zj . A third score makes a single appearance, after Proposition 7.17: (log pCj )′ , the true score of the current estimate, which coincides with ηj⋆ at R = Ik .

7.3

The relative gradient: what the “preconditioner” really is

A plain gradient ascent step W ← W + α∇W L requires the matrix inverse W −T at every step, which is both expensive and numerically delicate. The standard remedy is to multiply the gradient by W T W on the right, giving the update W ← W + α ∇W L W T W.

(169)

This is often presented as a Newton step with an approximate inverse Hessian. It is not, and Remark 7.14 explains what it is instead: an exact gradient, taken with respect to a different – and for this problem far more natural – Riemannian metric on the group of invertible matrices. Definition 7.12 (Relative gradient). Let F be a differentiable real-valued function on the group GL(k) of invertible k × k matrices. The relative gradient of F at W is the matrix ∇rel F (W ) determined by the first-order expansion under multiplicative perturbations, F (Ik + ε)W = F (W ) + ∇rel F (W ), ε + o ∥ε∥ ,

(170)

ε ∈ Rk×k ,

where ⟨M, N ⟩ := tr(M T N ) is the Frobenius inner product. Proposition 7.13 (Relative gradient of the ICA likelihood). For any differentiable F on
T rel GL(k) we have ∇ F (W ) = ∇F (W ) W , so that a relative gradient ascent step W ← (Ik + α∇rel F (W ))W is exactly Eq. (169). For the ICA log-likelihood Eq. (163), ∇rel L(W ) = Ik + E η(C) CT , 



C = W X,

(171)




(172)

and the update rule reads 

W ← Ik + α Ik + E η(C)CT 

W.

Proof. By the ordinary chain rule, F ((Ik + ε)W ) = F (W + εW ) = F (W ) + ⟨∇F (W ), εW ⟩ +

o(∥ε∥), and ⟨∇F (W ), εW ⟩ = tr ∇F (W )T εW = tr W ∇F (W )T ε = ∇F (W )W T , ε , which is Eq. (170) with ∇rel F (W ) = ∇F (W )W T . Substituting Eq. (168) and using W −T W T = Ik and XT W T = CT , 



∇rel L(W ) = W −T + E η(C)XT W T = Ik + E η(C) XT W T = Ik + E η(C)CT . 



40







(173)

Remark 7.14 (It is a metric, not a Hessian). Differentiating Eq. (168) once more gives the true Hessian of L,  

∂2L = δil E ηi′ (Ci )Xj Xm − W −1 jl W −1 mi , ∂Wlm ∂Wij

(174)

using ∂(W −1 )ji /∂Wlm = −(W −1 )jl (W −1 )mi . The operator implicitly used in Eq. (169) is, in coordinates, K(lm),(ij) = δli W T W jm = δli

k X

Wrj Wrm ,

(175)

r=1

since (∇L)W T W lm = ij K(lm),(ij) (∇L)ij . Two observations settle the matter. First, K does not depend on the data at all, whereas Eq. (174) does, through E[ηi′ (Ci )Xj Xm ]; so K cannot be the inverse of Eq. (174), not even approximately, uniformly in the source distribution. Second, K has an exact interpretation that has nothing to do with curvature: it is the inverse metric tensor of the right-invariant Riemannian metric on GL(k) in which distances are measured multiplicatively, and Eq. (172) is the corresponding steepest-ascent step. This is Amari’s natural gradient (Amari et al., 1996; Amari, 1998) and Cardoso and Laheld’s relative gradient (Cardoso and Laheld, 1996); MacKay’s derivation of exactly the update Eq. (172) calls it the covariant algorithm (MacKay, 1996). The essential point for what follows is that Eq. (172) involves no approximation whatsoever: it is the exact gradient of the exact objective for the right-invariant metric just described. (Not for a different chart: a gradient in a chart is the gradient for the flat metric that chart induces, and no flat left- or right-invariant metric exists on GL(k) for k ≥ 2.) The only approximation in the whole algorithm is the choice of the model scores ηj – that is, the extent to which they fail to be the true scores ηj⋆ – and that is the subject of Section 7.4. 



P

The reason to prefer this particular gradient is the following invariance, which is why the algorithm is used at all. Proposition 7.15 (Equivariance). Write R := W A for the global system matrix, so that C = W X = RZ. Then the update Eq. (172) induces on R the update 

R ← Ik + α Ik + E η(RZ)(RZ)T 




R,

(176)

which does not involve A. Consequently, if (Wn )n≥0 and (Wn′ )n≥0 are the iterates produced from data X = AZ and X ′ = A′ Z with the same sources and with initialisations satisfying W0 A = W0′ A′ , then Wn A = Wn′ A′ and Wn X = Wn′ X ′ for all n. Proof. Multiplying Eq. (172) on the right by A and using C = W X = W AZ = RZ turns the update for W into Eq. (176), in which A no longer appears. The second statement follows by induction: the hypothesis Wn A = Wn′ A′ makes both iterates apply the same map Eq. (176) ′ to the same matrix Rn , so Wn+1 A = Wn+1 A′ ; and Wn X = Rn Z = Wn′ X ′ . The trajectory of the recovered signal C therefore depends on the mixing matrix only through R0 = W0 A: two mixings A and M A give identical C-trajectories provided the initialisations are co-transformed, however badly conditioned M is. This is not the statement that the conditioning of A is irrelevant; see Remark 7.16. Remark 7.16 (Equivariant, relative, natural, covariant). Four names circulate for the same update, and they name four different things about it. Equivariant (Cardoso and Laheld,

41

1996) names the property proved in Proposition 7.15: writing Φn (x, W0 ) for the n-th iterate produced from data x and initialisation W0 , that proposition says Φn M x, W0 M −1 = Φn (x, W0 ) M −1 ,

M ∈ GL(k),

(177)

which is equivariance of the map (x, W0 ) 7→ Wn under the action x 7→ M x, W0 7→ W0 M −1 of GL(k). Note that the initialisation is co-transformed: for a fixed W0 , such as the default W0 = Ik of Algorithm 7.28, the iterates are not equivariant, which is the point already made in Remark 7.9. Relative gradient (Cardoso and Laheld, 1996) names the construction: the gradient with respect to multiplicative perturbations, Definition 7.12. Natural gradient (Amari, 1998) names the same object as steepest ascent for a Riemannian metric on GL(k). Covariant (MacKay, 1996) borrows the physicists’ sense of the word – form-invariant under a change of parametrisation. We prefer equivariant. It names a property of the algorithm that we have actually proved, rather than a property of the metric used to derive it; it is the established term in the sourceseparation literature; and “covariant” is doubly unfortunate here, since it invites confusion with covariance matrices, and since in the tensor-calculus sense of the word the object K of Remark 7.14 raises an index and so produces a contravariant quantity. When the emphasis is on the derivation rather than on the algorithm, “relative gradient” and “natural gradient” are both accurate.

7.4

Stationary points, stability, and the choice of nonlinearity

We now address the question that the derivation so far has left open: the true scores ηj⋆ are unknown, so the algorithm is run with a fixed guess ηj in their place (Remark 7.11). Why should it still work? Proposition 7.17 (Stationary points are nonlinear decorrelations). W is a stationary point of the relative gradient, i.e. ∇rel L(W ) = 0, if and only if E ηj (Cj ) Cl = −δjl , 



j, l = 1, . . . , k.

(178)

Proof. Immediate from Eq. (171): the (j, l) entry of Ik + E[η(C)CT ] is δjl + E[ηj (Cj )Cl ]. The k diagonal equations of Eq. (178) fix the scale of each recovered component; the k(k − 1) off-diagonal equations say that each transformed component ηj (Cj ) is orthogonal to every untransformed Cl , l ̸= j, and hence – the data being centred – uncorrelated with it. This is a nonlinear decorrelation condition. If the model score ηj happens to coincide with the true score (log pCj )′ of the current estimate Cj , and c pCj (c) → 0 as |c| → ∞, then the j-th diagonal R equation holds automatically: integration by parts gives E[ηj (Cj )Cj ] = c p′Cj (c) dc = −1. By Proposition 7.15 the dynamics lives on the global system matrix R, and by Corollary 7.4 the separating solutions are exactly the matrices R = P Λ. Relabelling and rescaling the sources, it suffices to analyse the fixed point R = Ik . When the ηj differ across components this relabelling matters: the conditions below are to be read for the pairing of ηj with the source that P sends to position j, and a different pairing has different ζ’s and its own verdict. That is exactly what makes a wrong assignment of model scores to components unstable. We first record that the scale can always be arranged. Lemma 7.18 (The fixed-point scale exists and is unique). Let Y be a real-valued random variable that is not almost surely 0 and let the model score η be one of η(u) = − tanh(u),

η(u) = tanh(u) − u, 42

η(u) = −u3 ,

(179)

with E[|Y |] < ∞ in the first case, E[Y 2 ] < ∞ in the second and E[Y 4 ] < ∞ in the third. Then there is exactly one c > 0 with E η(cY ) cY = −1. 

(180)



Proof. In all three cases write h(c) := −E[η(cY )cY ] = E[φ(cY )] with φ(u) := u tanh(u),
φ(u) := u u − tanh(u) and φ(u) := u4 respectively. Each φ is non-negative, even, vanishes only at 0 and is strictly increasing in |u|; for the second this is because u − tanh u has the sign of u and |u − tanh u| increases in |u|. Hence c 7→ φ(cY (ω)) is non-decreasing in c > 0 for every ω, and strictly increasing whenever Y (ω) ̸= 0. So h is continuous by dominated convergence on compacts, strictly increasing, h(0+ ) = 0, and h(c) → ∞ as c → ∞ by monotone convergence. The intermediate value theorem gives a c > 0 with h(c) = 1, and strict monotonicity makes it unique. Remark 7.19 (Why the normalisation makes ζj well defined). The quantities βj , γj and ζj of Eq. (183) below are not invariant under rescaling the source: replacing Zj by cZj changes βj = E[ηj′ (Zj )] and σj2 = E[Zj2 ] separately, and in general changes their product. This is not a defect but a consequence of the fact that the fixed point under analysis is R = Ik , which presupposes a choice of scale for the sources. Lemma 7.18 says that for the three nonlinearities in use there is exactly one such choice compatible with stationarity up to sign, namely the one making Eq. (182) hold; the sign is invisible to βj , γj , ζj , because ηj′ is even for all three, and it is exactly the sign ambiguity of Corollary 7.4. Once that normalisation is imposed, ζj is a well-defined functional of L(Zj ) and of the model score ηj alone, and the stability conditions Eq. (186) are statements about the source law rather than about an arbitrary scaling. This is what makes the entries of Table 2 meaningful. Theorem 7.20 (Local stability of the separating solution). Let k ≥ 2, let Z have mutually independent components with E[Zj ] = 0 and σj2 := E[Zj2 ] ∈ (0, ∞), let the model scores η1 , . . . , ηk of Remark 7.11 be arbitrary twice continuously differentiable functions, and assume sup ηj′′ (u) < ∞

(181)

E ∥Z∥3 < ∞ 

for all j,

u∈R



(see Remark 7.21). Assume finally the normalisation E ηj (Zj )Zj = −1, 

j = 1, . . . , k,



(182)

holds. Then R = Ik is a stationary point of Eq. (176). Define βj := E ηj′ (Zj ) , 

γj := E ηj′ (Zj )Zj2 ,





ζj := −βj σj2 .



(183)

Writing R = Ik + ε, the update Eq. (176) linearised at ε = 0 decouples into the k scalar recursions εjj 7−→ 1 + α(γj − 1) εjj

(184)

and the

k
2 two-dimensional recursions

"

#

"

#

"

#

β σ 2 −1 Mjl := j l , −1 βl σj2

εjl εjl 7−→ I2 + αMjl , εlj εlj

j ̸= l.

(185)

Consequently the conditions γj < 1

for all j,

ζj > 0

for all j,

ζ j ζl > 1

for all j ̸= l

(186)

are sufficient for R = Ik to be locally asymptotically stable for all sufficiently small step sizes α > 0; if any one of them is reversed strictly, then R = Ik is unstable for every α > 0. 43

Proof. Stationarity at R = Ik : the (j, l) entry of Ik + E[η(Z)ZT ] is δjl + E[ηj (Zj )Zl ], which vanishes for j = l by Eq. (182) and for j = ̸ l because independence and E[Zl ] = 0 give E[ηj (Zj )Zl ] = E[ηj (Zj )]E[Zl ] = 0. P Now let R = Ik + ε and C = RZ, so Cj = Zj + m εjm Zm . Write G(R) := Ik + E[η(C)CT ]. Since Eq. (176) reads R 7→ (Ik + αG(R))R and G(Ik ) = 0, the induced map on ε is ε 7→ ε + αG(Ik + ε) + αG(Ik + ε)ε = ε + αG1 (ε) + O(∥ε∥2 ), where G1 is the derivative of G at Ik . Expanding to first order, ηj (Cj ) = ηj (Zj ) + ηj′ (Zj )

X

εjm Zm + O ∥ε∥2 ,

(187)

m

so that E ηj (Cj )Cl = E ηj (Zj )Zl + 





|

{z

=−δjl



X

}

m

εlm E ηj (Zj )Zm + 



X

εjm E ηj′ (Zj )Zm Zl + O ∥ε∥2 . (188) 



m

In the second sum E[ηj (Zj )Zm ] equals −1 for m = j and 0 otherwise, contributing −εlj . In the third sum, independence and E[Zm ] = 0 leave only m = l when l ̸= j, contributing εjl E[ηj′ (Zj )]E[Zl2 ] = βj σl2 εjl ; and only m = j when l = j, contributing γj εjj . Hence Gjl (Ik + ε) =

 (γj − 1) εjj ,

l = j,

β σ 2 ε − ε , j l jl lj

l ̸= j,

(189)

up to O(∥ε∥2 ), which is exactly Eqs. (184) and (185). For the stability criterion, Ik + αM is a contraction for all sufficiently small α > 0 as soon as every eigenvalue of M has strictly negative real part, and it fails to be one, for every α > 0, as soon as some eigenvalue has strictly positive real part; eigenvalues on the imaginary axis are undecided. For the scalar blocks the relevant condition is γj − 1 < 0. Each Mjl is real symmetric, so its eigenvalues are real and both are negative if and only if tr Mjl < 0 and det Mjl > 0, that is βj σl2 + βl σj2 < 0,

βj βl σj2 σl2 > 1.

(190)

In terms of ζj = −βj σj2 the second condition reads ζj ζl > 1, and the first reads ζj σl2 /σj2 + ζl σj2 /σl2 > 0. If ζj ζl > 1 then ζj and ζl have the same sign; if both were negative the first condition would fail, so both are positive, and conversely ζj , ζl > 0 implies the first condition. This is Eq. (186). It remains to check that each strict reversal produces an eigenvalue with strictly positive real part. For γj > 1 the diagonal multiplier 1 + α(γj − 1) exceeds 1 for every α > 0. For ζj ζl < 1 we get det Mjl = ζj ζl − 1 < 0, so Mjl has one positive and one negative eigenvalue. For ζj < 0, pick any l = ̸ j (possible since k ≥ 2): if ζj ζl < 1 we are in the previous case; otherwise ζj ζl ≥ 1 forces ζl < 0, and then tr Mjl = −ζj σl2 /σj2 − ζl σj2 /σl2 > 0, so at least one eigenvalue of Mjl is positive. Finally, if Ik + αM has an eigenvalue 1 + αλ with Re λ > 0, then |1 + αλ|2 = 1 + 2α Re λ + α2 |λ|2 > 1 for every α > 0, so the linearised map is expanding in the corresponding direction. Since the remainder in Eqs. (184) and (185) is O(∥ε∥2 ) and the linear part is symmetric, the unstable-manifold theorem for maps then makes R = Ik unstable for the nonlinear iteration as well. Remark 7.21 (Reading Theorem 7.20). Four comments on the statement. (i) What the regularity condition is for. The perturbation ∆j := m εjm Zm that appears in the proof is not bounded almost surely – only |∆j | ≤ ∥ε∥∥Z∥ – so the mean value theorem has to be applied at an intermediate point that can be far from Zj . A bound on ηj′′ that is uniform on all of R removes the difficulty, and E∥Z∥3 < ∞ then makes the Taylor remainder O(∥ε∥2 ). P

44

(ii) What it covers. Eq. (181) holds for ηj = − tanh and ηj = tanh −id, but not for the cubic ηj (u) = −u3 , whose second derivative ηj′′ (u) = −6u is unbounded. For the cubic one uses instead the polynomial variant of the same estimate, valid whenever E∥Z∥4 < ∞: there E[ηj (Cj )Cl ] is an explicit polynomial in ε of degree at most four whose coefficients are moments of Z of order at most four. (iii) The borderline. The only cases left undecided by the linearisation are γj = 1 and ζj ζl = 1 with ζj , ζl > 0. The apparent third borderline case ζj = 0 is decided: it gives ζj ζl = 0 < 1, hence instability. (iv) The case k = 1. There is then nothing to separate, and only the condition γ1 < 1 survives. Theorem 7.20 is the stability analysis of Amari et al. (1997), recast in the notation of these notes; see also Cardoso (1998b, Section VI). The quantity ζj therefore decides everything. The next three results compute it in the cases of interest. The first explains, in one line, why Gaussian sources are the exact borderline – for every choice of model score, correctly specified or not. Proposition 7.22 (Gaussian sources lie exactly on the stability boundary). In the setting of Theorem 7.20, suppose Zj ∼ N (0, σj2 ) with σj > 0, and suppose in addition that ηj (z) p(z) → 0 as |z| → ∞, where p denotes the N (0, σj2 ) density (see Remark 7.24). Then ζj = 1, whatever the model score ηj – correctly specified or not. Consequently, if two sources are Gaussian, then ζj ζl = 1 and Eq. (186) fails. Proof. Let p be the N (0, σj2 ) density, so that z p(z) = −σj2 p′ (z). Integration by parts, whose boundary term [ηj p]∞ −∞ vanishes by hypothesis, gives Stein’s identity E ηj (Zj )Zj = 



Z

ηj (z) z p(z) dz = −σj2

Z

ηj (z)p′ (z) dz = σj2

Z

ηj′ (z)p(z) dz = σj2 βj . (191)

The normalisation Eq. (182) says that the left hand side is −1, so ζj = −βj σj2 = 1. Corollary 7.23 (The true score: stability matches identifiability exactly). In the setting of Theorem 7.20, but requiring of each ηj only one continuous derivative and assuming neither Eq. (181) nor Eq. (182), suppose the model is correctly specified, that is
′

ηj = ηj⋆ := log pZj ,

j = 1, . . . , k,

(192)

where pZj is the density of Zj . Assume that each p is strictly positive and twice differentiable
′′ Zj 1 2 with continuous second derivative, that 1 + z pZj ∈ L (R), that the expectations in Eqs. (182) and (183) converge absolutely, and that the boundary conditions z pZj (z) −→ 0,

z 2 p′Zj (z) −→ 0,

|z| → ∞

(193)

hold. Then Eq. (182) holds automatically, ζj = I(Zj ) σj2 ≥ 1,

I(Zj ) := E ηj⋆ (Zj )2 , 



(194)

with equality if and only if Zj is Gaussian, and γj < 1. Hence the inequalities Eq. (186) hold if and only if at most one of the sources is Gaussian. Proof. Write p := pZj and η := ηj = ηj⋆ , so that ηp = p′ and η ′ p = p′′ − (p′ )2 /p; it is Eq. (192) that makes the first of these identities available, and it is what turns each of the expectations below into an integral of p and its derivatives. The second boundary condition in Eq. (193) 45

forces p′ (z) → 0, and p′′ = limR→∞ p′ (R) − p′ (−R) = 0; the integral exists because p′′ ∈ L1 . Normalisation: integration by Rparts, with vanishing boundary term by the first condition in R ′ Eq. (193), gives E[η(Zj )Zj ] = zp (z) dz = − p = −1. For βj ,

R

βj =

(p′ )2  p − 2 p= p p

Z  ′′

Z

p′′ − E η(Zj )2 = −I(Zj ). 

(195)



The inequality is Cauchy–Schwarz applied to the normalisation:
2

1 = E[η(Zj )Zj ]

≤ E η(Zj )2 E Zj2 = I(Zj )σj2 = ζj , 

 

(196)



with equality if and only if η(Zj ) and Zj are almost surely proportional, i.e. (log p)′ (z) = −z/σj2 , i.e. p is the N (0, σj2 ) density. For γj , integrating by parts twice – the first boundary term vanishing by Rthe second condition in Eq. (193) and the second by the first condition – we get R 2 ′′ z p = −2 zp′ = 2, so γj =

Z

z2

 p′′

p



− η 2 p = 2 − E η(Zj )2 Zj2 ≤ 2 − E[η(Zj )Zj ] 



2

= 1,

(197)

again by Cauchy–Schwarz, with equality only if η(Zj )Zj is almost surely constant, i.e. p(z) ∝ |z|−1 , which is not integrable; so γj < 1 strictly. Finally, ζj ≥ 1 for all j with equality exactly at the Gaussian sources, so ζj ζl > 1 for all j ̸= l if and only if at most one ζj equals 1. Remark 7.24 (On the hypotheses of Proposition 7.22 and Corollary 7.23). The decay condition ηj (z)p(z) → 0 of Proposition 7.22 is what makes the boundary term in Stein’s identity vanish; it holds whenever ηj has at most polynomial growth, and in particular for all three nonlinearities of Lemma 7.18. In Corollary 7.23 the second derivative of pZj is needed because ηj⋆ ′ = p′′ /p − (p′ /p)2 , so Theorem 7.20 presupposes it anyway, and its continuity is what makes ηj continuously differentiable; Theorem 7.20 asks for one derivative more. Note also what Corollary 7.23 does and does not deliver: it gives the inequalities Eq. (186), and these coincide with the hypothesis of Corollary 7.4; but to turn them into local stability one still needs the regularity Eq. (181), which a true score need not satisfy. Corollary 7.23 is the conceptual answer to “why does gradient descent find the right answer”. With correctly specified scores, ηj = ηj⋆ , the inequalities Eq. (186) hold exactly when the identifiability theory declares the model identifiable; subject to the regularity caveat of Remark 7.24, a separating solution is then a locally stable equilibrium of the mean dynamics. The converse is false, and it is worth being explicit about it: stable non-separating equilibria can exist even under a perfectly specified model. For k = 2 with Z1 , Z2 i.i.d. from the smooth, 2 2 strictly positive, non-Gaussian density 0.4 N (−1.5, 0.252 ) + 0.2 " N (0,#0.25 ) + 0.4 N (1.5, 0.25 ) 1 1 and ηj = ηj⋆ its exact score, the matrix R⋆ = 0.80993 · √12 is a stationary point of 1 −1 Eq. (176) whose linearisation has all four eigenvalues strictly negative, so it is locally attracting – and it is not a generalised permutation matrix, so C = R⋆ Z has dependent components. Theorem 7.20 says which separating solutions attract; it does not say that nothing else does. The remaining question is what happens when ηj ̸= ηj⋆ , which in practice is always. Corollary 7.25 (The cubic nonlinearity: kurtosis is exactly the criterion). Let ηj (u) := −u3 and let Zj have E[Zj ] = 0, E[Zj4 ] < ∞ and be normalised as in Lemma 7.18, so that E[Zj4 ] = 1. Then γj = −3 and ζj = 3σj4 =

3 . 3 + kurt(Zj )

(198)

Hence ζj > 1 if and only if kurt(Zj ) < 0, i.e. if and only if Zj is sub-Gaussian in the sense of Definition 3.21. 46

Proof. ηj′ (u) = −3u2 , so βj = −3σj2 and ζj = 3σj4 , while γj = −3E[Zj4 ] = −3. By Definition 3.21, kurt(Zj ) = E[Zj4 ]/σj4 − 3 = σj−4 − 3, so σj4 = (3 + kurt(Zj ))−1 and ζj = 3/(3 + kurt(Zj )), which exceeds 1 exactly when kurt(Zj ) < 0. Remark 7.26 (The hyperbolic tangent, and what the kurtosis rule really is). For the two tanhbased model scores the criterion ζj > 1 becomes an explicit inequality. With the normalisation Eq. (182): • ηj (u) = − tanh(u), so that E[tanh(Zj )Zj ] = 1 and βj = −E[sech2 (Zj )]; the condition ζj > 1 reads E Zj2 E sech2 (Zj ) > E tanh(Zj ) Zj . 

 







(199)

• ηj (u) = tanh(u) − u, so that βj = −E tanh2 (Zj ) ; here the condition ζj > 1 reads 



E Zj2 E tanh2 (Zj ) > 1. 

 



(200)

Eq. (199) is a Stein discrepancy: by the identity in the proof of Proposition 7.22, its two sides are equal for every Gaussian Zj and for every test function, so the inequality measures departure from Gaussianity in the direction picked out by tanh. It is exactly the switching statistic of Lee et al. (1999, eq. (2.28)). It is tempting – and it is what the ICA literature usually says – to summarise this as “use − tanh for super-Gaussian sources and tanh −id for sub-Gaussian ones”. That rule is a heuristic, not a theorem: the sign of kurt(Zj ) is exactly the criterion only for the cubic nonlinearity (Corollary 7.25), and for tanh the correct criterion is Eq. (199), which must be checked. It does hold for the standard source models; Table 2 lists the values. In general it fails in both directions, and the failures are not exotic. For the symmetric three-point family P[Y = ±a] = p/2, P[Y = 0] = 1 − p, normalised as in Lemma 7.18, write u := ca for the rescaled atom; the normalisation reads p u tanh u = 1 and one computes ζ=

u − 1, tanh u

kurt = u tanh u − 3.

(201)

So ζ = 1 happens at u = 2 tanh u, i.e. at kurt ≈ −1.1664, and every member of the family with −1.1664 < kurt < 0 is sub-Gaussian yet has ζ > 1, so that − tanh separates it. In the other direction, mass placed far out in units of σ inflates kurt while saturating tanh: the four-atom law with P[Y = ±M ] = q/2 and P[Y = ±1] = (1 − q)/2 has kurt ≈ 451 but ζ ≈ 0.46 for q = 5 · 10−6 , M = 100, so two such super-Gaussian sources make the separating solution unstable for − tanh. Note finally that the sign of the discrepancy in Eq. (199) is scale dependent – it is proportional to kurt only in the small-scale limit – so the criterion is meaningful only at the fixed-point scale of Lemma 7.18; cf. Remark 7.19. Remark 7.27 (What is approximate, and why it does not spoil convergence). It is worth separating three things that are easy to confuse. (i) The step Eq. (172) is exact: it is the gradient of the exact objective for the metric (Proposition 7.13). Nothing is approximated there, and in particular no Hessian is being approximated (Remark 7.14). (ii) The model scores ηj are an approximation to the true scores ηj⋆ , since the true source densities are unknown. Their entire effect on whether the algorithm converges to a separating solution is captured by Theorem 7.20: a model score with ηj ̸= ηj⋆ is harmless as long as γj < 1, ζj > 0 and ζj ζl > 1 – and for the three nonlinearities of Lemma 7.18 the first two are automatic, since ηj′ < 0 off a Lebesgue null set and Zj is non-degenerate, 47

Law of Zj

kurt(Zj )

ζj for η = − tanh

ζj for η = tanh −id

ζj for η(u) = −u3

Gaussian

0

1.000

1.000

1.000

Laplace

+3

1.414

0.726

0.500

Student t5

+6

1.254

0.808

0.333

Student t8

+1.5

1.126

0.895

0.667

Logistic

+1.2

1.134

0.889

0.714

Uniform

−1.2

0.762

1.273

1.667

Table 2: The stability quantity ζj of Eq. (183) for the three standard nonlinearities. Each law is taken in its standard normalisation – N (0, 1); Laplace and logistic with unit scale parameter; tν standard; uniform on [−1, 1] – and is then rescaled by the unique factor c > 0 of Lemma 7.18, which is what makes ζj well defined (Remark 7.19); the tabulated values do not depend on the starting normalisation. Stability of a separating solution requires ζj ζl > 1 for all j ̸= l, for which ζj > 1 for all j is sufficient; the remaining condition γj < 1 is automatic here, since ηj′ ≤ 0 for all three nonlinearities and hence γj ≤ 0. Every Gaussian entry equals 1 exactly, as Proposition 7.22 predicts, and the cubic column is 3/(3 + kurt) exactly, as Corollary 7.25 predicts. The values were obtained by numerical quadrature. so γj < 0 and ζj > 0. This is why one only needs to know each source qualitatively: in practice only the sign of ζj − 1, rather than the density. For the source families of Table 2 that sign is determined by the sign of the excess kurtosis, for each of the three nonlinearities separately, but the correspondence does not hold in general – see the warning and the counterexamples in Remark 7.26, and note that ζj is a statement about neither kurtosis nor tails. (iii) What ηj ̸= ηj⋆ does cost is statistical efficiency. The maximum likelihood estimator with the true scores attains the Cramér–Rao bound asymptotically; with misspecified scores the separating solutions remain stationary points, at the scale fixed by Lemma 7.18, and remain locally attracting whenever Theorem 7.20 applies, but the asymptotic variance of the resulting estimator is larger. This is the semiparametric picture of Amari and Cardoso (1997); see also Cardoso (1998b). Two caveats should be stated plainly. First, Theorem 7.20 is a local statement. The objective Eq. (163) is not concave in W , and Eq. (176) has stationary points that are not separating, some of which can themselves be attracting. The theorem says when a separating solution attracts, not that the algorithm finds one from an arbitrary start. Second, the online form Algorithm 7.28 below replaces the expectation in Eq. (172) by a single-sample estimate. It is thus a Robbins–Monro stochastic approximation scheme (Robbins and Monro, 1951) for the mean dynamics Eq. (176) – provided the samples are i.i.d. across t, or at least mixing enough that the single-sample estimate is a martingale difference plus a vanishing bias; with step P P sizes satisfying n αn = ∞ and n αn2 < ∞, and under the usual regularity and boundedness conditions, its trajectories track the associated ordinary differential equation 



Ṙ = Ik + E η(RZ)(RZ)T R 

(202)

and converge almost surely to one of its locally stable equilibria – which Theorem 7.20 identifies. A constant step size α, as used in practice, gives convergence only to a neighbourhood whose size is O(α).

48

7.5

The algorithm

Algorithm 7.28 (Equivariant online ICA for the complete noiseless model). Input: data x(1), . . . , x(T ) in Rk , centred; a step size α > 0; model scores η1 , . . . , ηk chosen according to Theorem 7.20, in practice as in Remark 7.29; an invertible initial W ∈ Rk×k , e.g. W = Ik . Repeat over the data, until convergence: for each sample x(t), (1) c := W x(t)

(current estimate of the sources) T

(2) y := η1 (c1 ), . . . , ηk (ck ) 

(nonlinearity applied componentwise)

(3) W ← W + α W + y cT W = W + α Ik + y c W .
T

Output: the unmixing matrix W , and the reconstructed sources obtained by running step (1) once more with the final W , that is c(t) = W x(t); they are determined up to permutation and scale. Remark 7.29 (Reading the algorithm). In practice one takes ηj = − tanh for super-Gaussian and ηj = tanh −id for sub-Gaussian sources, with the caveats of Remark 7.26. Step (3) is Eq. (172) with the expectation replaced by the current sample, since E[η(C)CT ] is estimated by y cT . Written out as W + α(W + η(W x) xT W T W ) it is exactly the update of Amari et al. (1996), Cardoso and Laheld (1996) and MacKay (1996) – see also MacKay (2003, Chapter 34) for the same derivation in an archival source; the infomax algorithm of Bell and Sejnowski (1995) is the same update in its plain-gradient form, and the two objectives coincide, as observed by Cardoso (1997). Each step costs O(k 2 ) operations and no matrix inversion or decomposition: y cT W is computed as y (cT W ), a vector–matrix product followed by an outer product, each O(k 2 ). By Proposition 7.15 the whole dynamics is a dynamics on the global system matrix R = W A, so the mixing matrix enters only through R0 = W0 A: once started, the algorithm behaves identically for every A giving the same R0 , however badly conditioned A itself is. This is what distinguishes the update from plain gradient ascent on L, whose trajectory depends on A throughout.

7.6

LiNGAM: a causal order removes the permutation

Corollary 7.4 leaves the sources determined only up to permutation, sign and scale. In causal discovery one adds structural assumptions that remove exactly those remaining ambiguities. The basic such model is the linear non-Gaussian acyclic model, LiNGAM, of Shimizu et al. (2006). We write Θ for the matrix of structural coefficients, which Shimizu et al. (2006) call B; the letter B is already taken in these notes, and R is the global system matrix of Proposition 7.15. Assumption 7.30 (LiNGAM). The observed random vector X ∈ Rk satisfies the structural equations X = Θ X + Z,

(203)

where Θ ∈ Rk×k and the disturbance vector Z ∈ Rk satisfy: (i) Acyclicity. There is a bijection π : {1, . . . , k} → {1, . . . , k}, the causal order, such that Θab = 0 whenever π(a) ≤ π(b). Equivalently, Θ becomes strictly lower triangular after simultaneously permuting rows and columns by π; the associated directed graph, with an edge b → a whenever Θab ̸= 0, is acyclic. (ii) Independent disturbances, causal sufficiency. The components Z1 , . . . , Zk of the disturbance vector are mutually independent and almost surely non-constant. 49

(iii) Non-Gaussianity. At most one Zj is Gaussian. Lemma 7.31 (LiNGAM is a complete noiseless ICA model). Under Assumption 7.30 (i) the matrix Ik − Θ is invertible with det(Ik − Θ) = 1, and Eq. (203) is equivalent to X = A Z,

A := Ik − Θ

−1

W := A−1 = Ik − Θ.

,

(204)

This is a complete noiseless ICA model in the sense of Assumption 7.2 (i)–(iv), and the unmixing matrix W has all diagonal entries equal to 1. Proof. Let Pπ be the permutation matrix of the causal order, so that N := Pπ ΘPπT is strictly lower triangular and hence N k = 0; then Θk = PπT N k Pπ = 0, so Θ is nilpotent, Ik − Θ is P j invertible with inverse k−1 j=0 Θ , and det(Ik − Θ) = det(Ik − N ) = 1 because Ik − N is lower triangular with unit diagonal. Rearranging Eq. (203) gives (Ik − Θ)X = Z, which is Eq. (204). The diagonal of W = Ik − Θ is 1 because Θaa = 0 by (i) with b = a. The normalisation Assumption 7.2 (v) is deliberately not imposed here: LiNGAM fixes the scale of the sources through diag(W ) = Ik rather than through either of the conventions of Remark 7.3. Corollary 7.32 (Identifiability of LiNGAM). Suppose X satisfies Assumption 7.30 with two parameter sets, that is Θ(1) X + Z (1) = X = Θ(2) X + Z (2) ,

(205)

where each Θ(i) satisfies (i) and each Z (i) satisfies (ii) and (iii). Then Θ(1) = Θ(2) ,

Z (1) = Z (2)

almost surely.

(206)

Proof. By Lemma 7.31 we have two representations X = A(i) Z (i) with A(i) = (Ik − Θ(i) )−1 invertible, so Corollary 7.4 applies (with µ(i) = 0) and yields a permutation matrix P and an invertible diagonal Λ with A(2) = A(1) P Λ, together with Z (1) = P ΛZ (2) + c. Inverting, and writing W (i) := Ik − Θ(i) , W (2) = Λ−1 P −1 W (1) .

(207)

The matrix M := Λ−1 P −1 is a generalised permutation matrix: there are a permutation σ and non-zero scalars m1 , . . . , mk with Mac = ma δc,σ(a) , so that Eq. (207) reads (2)

(1)

Wab = ma Wσ(a),b ,

a, b = 1, . . . , k.

(208)

Both W (i) have unit diagonal by Lemma 7.31, so taking b = a in Eq. (208) gives (1)

(2) 1 = Waa = ma Wσ(a),a ,

(1)

Wσ(a),a ̸= 0 for all a.

hence (1)

(209)

(1)

Now use acyclicity of Θ(1) : for c ̸= b we have Wcb = −Θcb , which is non-zero only if π1 (c) > π1 (b), where π1 is the causal order of Θ(1) . So Eq. (209) forces, for every a, either

σ(a) = a

or

π1 σ(a) > π1 (a).

(210)

Suppose σ = ̸ id and let a lie on a non-trivial cycle of σ, of length L ≥ 2. Then σ j (a) ̸= σ j+1 (a) for every j, so π1 (a) < π1 σ(a) < · · · < π1 σ L (a) = π1 (a),

(211)

(1)

a contradiction. Hence σ = id, and then Eq. (209) gives ma = 1/Waa = 1, so M = Ik . Therefore W (2) = W (1) , i.e. Θ(2) = Θ(1) , and consequently Z (2) = W (2) X = W (1) X = Z (1) . 50

So the entire structural model – graph, coefficients and disturbance distributions – is determined by the law of X; the causal order itself is identified up to the ordering of variables that are not ancestrally related, which is exactly the freedom a topological sort leaves. Remark 7.33 (What each assumption does). It is worth seeing which ambiguity of Corollary 7.4 is removed by which structural assumption. The scale ambiguity Λ is removed by the form of Eq. (203): the coefficient of Xa in its own equation is fixed to 1, which is what makes W have unit diagonal. The permutation ambiguity P is removed by acyclicity, through the cycle argument above: a non-trivial relabelling would have to increase the causal order all the way around a cycle. Assumption 7.30 (ii) is causal sufficiency: mutual independence of the disturbances says exactly that there are no hidden common causes. Non-Gaussianity is what makes Corollary 7.4 available in the first place; without it one recovers the graph at best up to its Markov equivalence class, and then only under a faithfulness assumption. Note also that Shimizu et al. (2006) assume all disturbances non-Gaussian, whereas Corollary 7.32 needs only “at most one”, because that is all Corollary 7.4 needs. Algorithm 7.34 (LiNGAM by equivariant ICA, after Shimizu et al., 2006). Input: data x(1), . . . , x(T ) in Rk . (1) Centre the data. c of the unmixing matrix. By Corollary 7.4 it (2) Run Algorithm 7.28 to obtain an estimate W −1 T estimates Λ P W for the true W = Ik − Θ and some unknown generalised permutation.

(3) Undo the permutation. Find the row permutation P1 that makes all diagonal entries of c non-zero; in practice one minimises P (P1 W c )jj −1 over permutations, which is a P1 W j c. linear assignment problem. Put W ′ := P1 W (4) Undo the scaling. Divide each row of W ′ by its own diagonal entry, that is, put ′ , . . . , W ′ ) and W ′′ := D −1 W ′ , so that W ′′ has unit diagonal. D := diag(W11 kk b := I − W ′′ . (5) Read off the coefficients. Θ k b as close to strictly lower (6) Read off the causal order. Find the permutation π making Θ triangular as possible, e.g. by repeatedly setting the smallest-magnitude entries to zero until a causal order exists. Optionally prune the remaining small coefficients. b the causal order π, and the estimated disturbances Output: the coefficient matrix Θ, ′′ zb(t) := W x(t). c = Remark 7.35 (Why steps (3) and (6) are well posed). In the population limit, where W −1 T Λ P W exactly, the proof of Corollary 7.32 is precisely the statement that step (3) has a solution and that it is unique: the cycle argument shows that P1 = P is the only row permutation producing a nowhere-zero diagonal, and step (4) then recovers W exactly. Likewise b is permutable to strictly lower triangular form. Note step (6) succeeds because the true Θ also that step (2) invokes Algorithm 7.28, which was derived under the standing assumption of Section 7.2 that each source has a strictly positive differentiable density; Assumption 7.30 does not require this, so the assumption has to be added if the ICA-based route is taken. With finite samples neither search is exact, which is why both become combinatorial problems with tolerances; see Shimizu et al. (2006, Sections 4–5) for the details and Shimizu et al. (2011) for DirectLiNGAM, a later algorithm that estimates the causal order directly by iterated regression and independence testing, and so avoids the ICA step – and with it both the local convergence caveats of Remark 7.27 and the two combinatorial searches above.

51

7.7

Generalisations: a brief overview

The model of these notes – one linear mixture, instantaneous, with independent components – has been extended in many directions. This subsection is a short guide to the landscape, with references rather than derivations. Relaxing linearity. The obvious generalisation replaces the mixing matrix by a diffeomorphism, X = f (Z) with Z having independent components. This model is never identifiable: Hyvärinen and Pajunen (1999) show that any random vector X with a positive density admits a representation X = f (Z) with mutually independent Zj , obtained from the conditionalquantile construction of Darmois (1953), which returns uniform sources; composing that solution with the measure-preserving automorphisms of the unit cube then produces infinitely many further ones. Identifiability has to be bought back with extra structure. The dominant device is an observed auxiliary variable U – a time index, a segment label, an experimental condition – with respect to which the sources are conditionally independent and sufficiently variable, pZ|U (z | u) =

k Y

pZj |U (zj | u).

(212)

j=1

Under such assumptions the nonlinear model becomes identifiable up to permutation and componentwise transformations; see Hyvärinen and Morioka (2016) for time-contrastive learning, where U is a time segment label; Hyvärinen and Morioka (2017) for temporally dependent sources, where the role of U is played by the previous time point; and Khemakhem et al. (2020) for the identifiable variational autoencoder, which unifies these results. A different route constrains the geometry of f rather than the statistics of Z: independent mechanism analysis (Gresele et al., 2021) requires the columns of the Jacobian of f to be orthogonal, which by Gresele et al. (2021, Theorem 4.7) excludes the Darmois solutions of Hyvärinen and Pajunen (1999) without an auxiliary variable. It is not, however, known to deliver full identifiability on its own. Relaxing independence. If the sources fall into groups that are independent across groups but dependent within, one obtains multidimensional ICA, or independent subspace analysis (Cardoso, 1998a): the mixing matrix is identified only up to a permutation of equal-dimensional blocks and an invertible transformation within each block – the analogue, for a general block, of the orthogonal ambiguity inside a Gaussian block in Remark 7.5. Several datasets at once. Given M related datasets X (m) = A(m) Z (m) , m = 1, . . . , M , one may, in addition to independence within each dataset, tie corresponding sources across datasets. Independent vector analysis (Kim et al., 2006) keeps the components independent within a dataset but couples the j-th components across datasets, j = 1, . . . , k, into a single source component vector; this aligns the permutations across datasets, leaving only one global permutation of the k component vectors. Multi-view ICA goes further and splits the sources into a part shared by all views and parts specific to each, "

X

(m)

=A

(m)

S

Z (m)

#

+ E (m) ,

m = 1, . . . , M,

(213)

with a shared source vector S and view-specific sources Z (m) . In Pandeva and Forré (2023a) each A(m) is square and invertible and the noise is isotropic Gaussian in the latent space, 52

so that E (m) = A(m) ϵ(m) with ϵ(m) ∼ N (0, σ 2 Ikm ) and the observed noise covariance is T σ 2 A(m) A(m ) . See that paper for identifiability results and an estimation procedure, Pandeva and Forré (2023b) for its use in integrating omics modalities, and Pandeva et al. (2025) for a related multi-view latent-variable model that replaces source independence by a shared sparse precision matrix and targets co-expression network inference. Noise and overcompleteness. The noisy and overcomplete models – X = AZ + E with A ∈ Rp×k and possibly k > p – are exactly the setting of Sections 5 and 6; note that identifiability of A there does not entail that the sources can be recovered pointwise. Additive Gaussian noise leaves them determined only up to a Gaussian summand, by Theorem 5.5, even when A has a left inverse; and in the overcomplete case k > p the matrix A has no left inverse at all. Causal discovery. Beyond LiNGAM, dropping linearity while keeping additive noise gives the additive noise models of Hoyer et al. (2008), in which Xa = fa (Xpa(a) ) + Za and identifiability again comes from a non-Gaussianity- or nonlinearity-induced asymmetry between cause and effect. Further reading. Hyvärinen et al. (2001) remains the standard textbook; Comon and Jutten (2010) is a comprehensive edited handbook covering algorithms, identifiability and applications.

A

Proofs of the Classical Results on Characteristic Functions

This appendix proves the four results of Section 3 that were stated there without proof: the uniqueness and inversion theorem (Theorem 3.4), Lévy’s continuity theorem (Theorem 3.7), Marcinkiewicz’ theorem (Theorem 3.17) and Cramér’s decomposition theorem (Theorem 3.20). Together with Section B, which proves Theorem 4.2, this makes every proof in these notes refer only to other results of these notes, together with the following standard background, which we use without proof: (a) measure theory: monotone and dominated convergence (hence differentiation under the integral sign), Fatou’s lemma, Fubini’s theorem, and the fact that a Borel probability measure on Rd is determined by its integrals against continuous functions of compact support; (b) Prokhorov’s theorem: a tight sequence of Borel probability measures on Rd has a weakly convergent subsequence (Klenke, 2020, Chapter 13), (Kallenberg, 2021, Chapter 5); (c) complex analysis of one variable: Morera’s theorem, the identity theorem, the mean value property of harmonic functions, and the fact that a nowhere vanishing entire function is exp of an entire function; (d) elementary real analysis: Taylor’s theorem with Peano remainder, Jensen’s and Lyapunov’s inequalities, existence of a median, the orthogonality relations for the trigonometric system together with termwise integration of uniformly convergent series, and the sub-subsequence criterion for convergence of real sequences. Lemma A.3 below supplies the converse half of Proposition 3.9. Its proof does use the forward half; but the forward half is proved where it is stated and does not use the converse, so there is no circularity. Two further standard facts are used outside Section 3 and its appendix: Osgood’s lemma in the several-variable half of Theorem 3.13 (which is never applied), and, in 53

Section 7, matrix calculus – Jacobi’s formula for ∂ log|det W | – together with the unstablemanifold theorem for maps in Theorem 7.20. Apart from these, the only result that a proof in these notes uses without our having proved it is Bochner’s Theorem 3.3, and it is used exactly once, in Remark 3.8, where it could be replaced by an appeal to Pólya’s criterion or by a direct Fourier computation. (Two more classical facts are quoted inside remarks – Pólya’s criterion itself and the Lukacs characterisation in Remark 3.16 – but no proof depends on them.)

A.1

Uniqueness and inversion

Lemma A.1 (Gaussian Fourier identity). For every ε > 0 and every x ∈ R, √ Z    x2  2π −itx 1 2 2 e exp − 2 ε t dt = exp − 2 . ε 2ε R

(214)

Consequently, for x ∈ Rd , (2π)−d

Z

T

Rd





e−i t x exp − 12 ε2 ∥t∥2 dt = γε (x) := 2πε2

−d/2

 ∥x∥2 

exp −

2ε2

(215)

,

the Lebesgue density of N (0, ε2 Id ). Proof. Write I(x) for the left hand side of Eq. (214). Differentiation under the integral sign is legitimate, since |t| exp(− 12 ε2 t2 ) is integrable, and an integration by parts gives I ′ (x) =

Z

i ε2

2 2

(−it)e−itx e−ε t /2 dt =

R

Z

e−itx

R

d  −ε2 t2 /2  x e dt = − 2 I(x), dt ε

(216) 2 2

−ε t /2 dt = 2 2 the √ boundary terms vanishing. Hence I(x) = I(0) exp(−x /(2ε )), and I(0) = R e 2π/ε. This is Eq. (214); Eq. (215) follows by taking the product over the d coordinates.

R

Proof of Theorem 3.4. Let X have law µ and let G ∼ N (0, Id ) be independent of X, both defined on one probability space. Fix ε > 0 and put Xε := X + εG. Conditioning on X R shows that Xε has the Lebesgue density fε (x) = Rd γε (x − y) µ(dy). Inserting Eq. (215) and applying Fubini’s theorem – legitimate because the integrand has modulus exp(− 12 ε2 ∥t∥2 ), which is integrable on Rd × Rd against the probability measure µ – we obtain fε (x) = (2π)−d = (2π)−d

Z

Z d

ZR

Rd

T

2

2

e−i t (x−y) e−ε ∥t∥ /2 dt µ(dy)

(217)

T

(218)

Rd 2

2

e−i t x e−ε ∥t∥ /2 φµ (t) dt.

So fε , and hence the law of Xε , is determined by φµ alone. Now let g : Rd → R be bounded and continuous. Since Xε → X pointwise on the underlying probability space as ε ↓ 0, dominated convergence gives Z Rd

g(x)fε (x) dx = E g(Xε ) −−→ E g(X) = 







ε↓0

Z

g dµ. Rd

(219)

The left hand side depends on µ only through φµ , so g dµ is determined by φµ for every bounded continuous g, and therefore µ is determined by φµ . This proves part (i). For part (ii) assume in addition φµ ∈ L1 (Rd ) and define R

f (x) := (2π)

−d

Z Rd

54

T

e−i t x φµ (t) dt.

(220)

Then f is bounded by (2π)−d ∥φµ ∥L1 and continuous, by dominated convergence. Comparing with Eq. (217), Z

sup fε (x) − f (x) ≤ (2π)−d

2

Rd

x∈Rd

2

e−ε ∥t∥ /2 − 1 |φµ (t)| dt −−→ 0

(221)

ε↓0

by dominated convergence, with dominating function 2|φµ |. Since every fε is a probability density and fε → f pointwise, the limit f is real valued and f ≥ 0, so ν(dx) := f (x) dx defines a Borel measure on Rd . Let now Rg be continuous with compact support. Uniform convergence R R R on that compact support gives gfε dx → gf dx, while Eq. (219) gives gfε dx → g dµ, so that Z Rd

g dµ =

Z

g ∈ Cc (Rd ).

g dν, Rd

(222)

Choosing gk ∈ Cc (Rd ) with 0 ≤ gk ↑ 1 pointwise and applying monotone convergence on both sides of Eq. (222) gives ν(Rd ) = µ(Rd ) = 1, so ν is a Borel probability measure. By Eq. (222) and background fact (a), µ = ν; that is, µ has the density f .

A.2

Lévy’s continuity theorem

Lemma A.2 (Tail bound from the characteristic function near the origin). Let Y be a real-valued random variable and δ > 0. Then h

P |Y | ≥

2i 1 ≤ δ δ

Z δ −δ

1 − Re φY (t) dt.

(223)

Proof. By Fubini’s theorem, using Re φY (t) = E[cos(tY )], Z δ −δ

1 − Re φY (t) dt = 2

Z δ 0



E 1 − cos(tY ) dt = 2δ E 1 − 



sin(δY ) , δY 

(224)

with the convention sin(u)/u := 1 at u = 0. The integrand 1 − sin(u)/u is non-negative for all u ∈ R, and for |u| ≥ 2 it is at least 1 − 1/|u| ≥ 12 . Hence 

2δ E 1 −

h   2i sin(δY ) ≥ 2δ · 12 P |δY | ≥ 2 = δ P |Y | ≥ , δY δ 

(225)

which is Eq. (223). Proof of Theorem 3.7. (i) If Xn ⇒ X then, for each fixed t, the functions x 7→ cos(tT x) and x 7→ sin(tT x) are bounded and continuous, so φXn (t) = E[cos(tT Xn )] + i E[sin(tT Xn )] converges to φX (t). (ii) Assume φXn → ϕ pointwise with ϕ continuous at 0. Since φXn (0) = 1 for all n we have ϕ(0) = 1. Tightness. Let ε > 0. By continuity of ϕ at 0 and ϕ(0) = 1 we may choose δ > 0 so small that 1 δ

Z δ −δ

1 − Re ϕ(sej ) ds <

ε , 2d

j = 1, . . . , d.

(226)

By Eq. (10) the characteristic function of the j-th coordinate eTj Xn is s 7→ φXn (sej ), and |1 − Re φXn (sej )| ≤ 2, so dominated convergence on the bounded interval [−δ, δ] gives 1 δ

Z δ −δ

1 − Re φXn (sej ) ds −−−→

n→∞

55

1 δ

Z δ −δ

1 − Re ϕ(sej ) ds <

ε . 2d

(227)

By Lemma A.2 there is therefore an N with P[|eTj Xn | ≥ 2/δ] < ε/d for all n ≥ N and all j = 1, . . . , d, whence P[∥Xn ∥∞ ≥ 2/δ] < ε for n ≥ N . Enlarging the compact set [−2/δ, 2/δ]d so as to accommodate the finitely many laws of X1 , . . . , XN −1 as well, we conclude that (L(Xn ))n∈N is tight. Identification of the limit. Let (Xnl )l be an arbitrary subsequence. By Prokhorov’s theorem it has a further subsequence converging weakly to the law of some random vector X. By part (i) the characteristic function of that limit is liml φXnl = ϕ; in particular ϕ is a characteristic function. By Theorem 3.4 (i) every subsequential weak limit has the same law, namely the one with characteristic function ϕ. Now fix a bounded continuous
g : Rd → R. By the previous paragraph, every subsequence of the real sequence E[g(Xn )] n has a further subsequence converging to E[g(X)], and a real sequence with that property converges to E[g(X)]. As g was an arbitrary bounded continuous function, Xn ⇒ X.

A.3

Entire characteristic functions

Both remaining proofs rest on the following observation: a characteristic function cannot be analytic “by accident”. As soon as it agrees with an entire function near the origin, the underlying law has Gaussian-type integrability. Lemma A.3 (Moments from differentiability at the origin). Let Y be a real-valued random variable, n ∈ N, and suppose that φY is 2n times differentiable at t = 0. Then E[Y 2n ] < ∞. (2n+2)

(2n+1)

Proof. We induct on n, the case n = 0 being trivial. Assume φY (0) exists. Then φY exists on a neighbourhood of the origin, so φY is in particular 2n times differentiable at 0 and the induction hypothesis gives E[Y 2n ] < ∞. By Proposition 3.9 the function φY is 2n times (2n) continuously differentiable on R with φY (t) = (−1)n E[Y 2n eitY ]. Put c := E[Y 2n ]. If c = 0 then Y = 0 almost surely and there is nothing to prove, so assume c > 0 and let ϱ be the (2n) probability measure with ϱ(dy) := c−1 y 2n L(Y )(dy). Then (−1)n c−1 φY is the characteristic function φϱ of ϱ, and by assumption it is twice differentiable at 0. Writing h := φϱ and using R h(t) + h(−t) − 2h(0) = −2 (1 − cos(ty))ϱ(dy) we get, by Fatou’s lemma, 2 1 − cos(ty) h(t) + h(−t) − 2h(0) y ϱ(dy) ≤ lim inf ϱ(dy) = − lim = −h′′ (0) < ∞, 2 2 t→0 t→0 t t R R (228)

Z

2

Z

because 2(1 − cos(ty))/t2 → y 2 pointwise as t → 0, and because the second symmetric difference quotient converges to h′′ (0) Rwhenever the latter exists, by Taylor’s theorem with Peano remainder. Hence E[Y 2n+2 ] = c y 2 ϱ(dy) < ∞. Lemma A.4 (Entire characteristic functions have all exponential moments). Let Y be a real-valued random variable and suppose that there are δ > 0 and an entire function F with φY (t) = F (t) for all t ∈ (−δ, δ). Then E exp(a|Y |) < ∞ 



for every a > 0,

(229)

the characteristic function φY extends to the entire function z 7→ E[eizY ], and this extension equals F on all of C. Proof. Since F is entire, φY is infinitely differentiable on (−δ, δ), so by Lemma A.3 all even moments of Y are finite, hence by Lyapunov’s inequality all moments are, and by

56

(m)

Proposition 3.9 φY (0) = im E[Y m ] for every m. These are also the derivatives of F at 0, so the Taylor series of F at the origin is F (z) =

X im E[Y m ]

m!

m≥0

(230)

zm,

and, F being entire, this series converges absolutely for every z ∈ C. In particular, for every a > 0, X a2k k≥0

(2k)!

(231)

E Y 2k < ∞. 



Using exp(a|y|) ≤ exp(ay) + exp(−ay) and monotone convergence, E exp(a|Y |) ≤ E eaY + E e−aY = 2 











X a2k k≥0

(2k)!

(232)

E Y 2k < ∞. 



By Theorem 3.13 the map z 7→ E[eizY ] is therefore entire and restricts to φY on R. It agrees with F on the interval (−δ, δ), hence on all of C by the identity theorem.

A.4

Marcinkiewicz’ theorem

Proof of Theorem 3.17. By assumption φY = eg on a neighbourhood of the origin, with a polynomial g; replacing g by g − g(0) – which changes eg not at all, and the degree only in the trivial case of a constant g, since eg(0) = φY (0) = 1 forces g(0) ∈ 2πi Z – we may and do assume g(0) = 0. The function eg is entire, so Lemma A.4 applies: Y has all exponential moments, and φY (z) = E eizY = exp g(z) , 



(233)

z ∈ C.

Write g(z) = nm=0 cm z m with cn ̸= 0 and suppose, aiming at a contradiction, that n ≥ 3. Since Y has all exponential moments, Proposition 3.15 applies on every horizontal strip and gives, for z = x + iy with x, y ∈ R, P

exp Re g(x + iy) = φY (x + iy) ≤ φY (iy) = exp Re g(iy) ,

(234)

where we used that φY (iy) = E[e−yY ] > 0 is real, so that φY (iy) = exp(Re g(iy)). Taking logarithms, u(x, y) := Re g(x + iy),

u(x, y) ≤ u(0, y),

(235)

x, y ∈ R.

Write cn = |cn |eiψ and, for θ ∈ [0, 2π) and r > 0, put z = reiθ . Since g is a polynomial of degree n, u r cos θ, r sin θ = Re cn rn einθ + O rn−1 = |cn | cos(ψ + nθ) rn + O rn−1 ,

(236)

uniformly in θ. Likewise, for y ∈ R, u(0, y) = Re cn (iy)n + O |y|n−1 ≤ |cn ||y|n + O |y|n−1 .

(237)

Let θ∗ be any maximiser of θ 7→ cos(ψ + nθ), i.e. cos(ψ + nθ∗ ) = 1. Applying Eq. (235) with x = r cos θ∗ and y = r sin θ∗ and inserting Eqs. (236) and (237) gives |cn |rn + O rn−1 ≤ |cn | |sin θ∗ |n rn + O rn−1 .

57

(238)

Dividing by rn and letting r → ∞ yields |cn | ≤ |cn ||sin θ∗ |n , and since cn ̸= 0 and |sin θ∗ | ≤ 1 this forces |sin θ∗ | = 1, i.e. θ∗ ∈ {π/2, 3π/2}. But the equation cos(ψ + nθ) = 1 holds exactly when ψ + nθ ∈ 2π Z, so its solutions form an arithmetic progression of spacing 2π/n, of which exactly n lie in [0, 2π). For n ≥ 3 there are thus at least three maximisers, while only two values of θ in [0, 2π) satisfy |sin θ| = 1. This contradiction shows n ≤ 2. It remains to determine the coefficients. We have g(t) = c1 t + c2 t2 with g(0) = 0. From φY (−t) = φY (t) (Proposition 3.2 (iii)) we get exp(g(−t)) = exp(g(t)) for all real t, so the polynomial t 7→ g(−t) − g(t) takes values in 2πi Z; being continuous and vanishing at t = 0, it vanishes identically. Hence −c1 t + c2 t2 = c1 t + c2 t2 for all real t, whence c2 = c2 is real and c1 = −c1 is purely imaginary, say c1 = iµ with µ ∈ R. Then |φY (t)| = exp(c2 t2 ) ≤ 1 for all t ∈ R forces c2 ≤ 0, so c2 = − 21 σ 2 with σ 2 := −2c2 ∈ R≥0 . Therefore φY (t) = exp(iµt− 12 σ 2 t2 ), which by Notation 2.4 and Theorem 3.4 (i) means Y ∼ N (µ, σ 2 ).

A.5

Cramér’s decomposition theorem

The growth estimate needed for Cramér’s theorem is supplied by the following classical coefficient bound, which replaces an appeal to Hadamard’s factorisation theorem. Lemma A.5 (Borel–Carathéodory coefficient bound). Let h be an entire function with P h(0) = 0 and Taylor expansion h(z) = m≥1 am z m . Then for every r > 0 and every m ≥ 1, |am | rm ≤ 4 M (r),

M (r) := max Re h(z). |z|=r

(239)

Proof. Note first that M (r) ≥ Re h(0) = 0 by the mean value property applied to the harmonic function Re h. Write am = αm + iβm with αm , βm ∈ R. On the circle z = reiϑ , Re h reiϑ =

X

rm αm cos(mϑ) − βm sin(mϑ) ,

(240)

m≥1

the series converging uniformly in ϑ. Termwise integration gives the Fourier coefficients 1 π

Z 2π 0

while the mean value property gives sign, 1 ±r αm = π m

1 π

Re h reiϑ cos(mϑ) dϑ = rm αm ,

Z 2π 0

Re h re

iϑ

R 2π 0

Z 2π 0

Re h reiϑ sin(mϑ) dϑ = −rm βm ,

(241)

Re h(reiϑ ) dϑ = 2π Re h(0) = 0. Therefore, for either

M (r) 1 ± cos(mϑ) dϑ ≤ π

Z 2π 0

1 ± cos(mϑ) dϑ = 2M (r), (242)

where we used 1 ± cos(mϑ) ≥ 0 and Re h(reiϑ ) ≤ M (r). Hence |αm |rm ≤ 2M (r), and the same argument with 1 ± sin(mϑ) gives |βm |rm ≤ 2M (r). Thus |am |rm ≤ (|αm | + |βm |)rm ≤ 4M (r). Proof of Theorem 3.20. Let Y1 ⊥ ⊥ Y2 with S := Y1 + Y2 ∼ N (µ, σ 2 ). 2 The degenerate case. If σ = 0 then φS = φY1 φY2 has modulus 1 everywhere, and since |φYj | ≤ 1 this forces |φYj | ≡ 1. Then the symmetrised variable Yj −Yj′ , with Yj′ an independent copy of Yj , has characteristic function |φYj |2 ≡ 1, hence is almost surely 0 by Theorem 3.4 (i); so Yj is almost surely constant, i.e. degenerate Gaussian. Assume from now on σ > 0, and, replacing Y1 by Y1 − µ, that µ = 0.

58

Step 1: Gaussian tails. Let (Y1′ , Y2′ ) be an independent copy of (Y1 , Y2 ) and put Yej := Yj −Yj′ and Se := Ye1 + Ye2 ∼ N (0, 2σ 2 ). The variables Ye1 , Ye2 are independent and symmetric, so P[Ye2 ≥ 0] ≥ 12 and hence, for x > 0, P Se > x ≥ P Ye1 > x, Ye2 ≥ 0 = P Ye1 > x P Ye2 ≥ 0 ≥ 12 P Ye1 > x . 









 







(243)

A Chernoff bound for Se ∼ N (0, 2σ 2 ) gives P[Se > x] ≤ exp(−x2 /(4σ 2 )), so by symmetry 

P |Ye1 | > x ≤ 4 exp − 



x2  , 4σ 2

x > 0.

(244)

Let θ be a median of Y1 , so that P[Y1′ ≤ θ] ≥ 21 and P[Y1′ ≥ θ] ≥ 12 . Then P Ye1 > x ≥ P Y1 > θ + x, Y1′ ≤ θ ≥ 21 P Y1 > θ + x , 











(245)

and symmetrically for the left tail, so Eq. (244) yields P[|Y1 − θ| > x] ≤ 16 exp(−x2 /(4σ 2 )). Consequently E[exp(a|Y1 |)] < ∞ for every a > 0, and the same holds for Y2 . Step 2: entire extensions. By Theorem 3.13 both φY1 and φY2 extend to entire functions, given by φYj (z) = E[eizYj ]. The entire functions z 7→ φY1 (z)φY2 (z) and z 7→ exp(− 12 σ 2 z 2 ) agree on R, hence on C by the identity theorem: 



φY1 (z) φY2 (z) = exp − 12 σ 2 z 2 ,

z ∈ C.

(246)

In particular neither factor has a zero in C. Step 3: growth. Fix j ∈ {1, 2} and let l be the other index. By Proposition 3.15, for z = x + iy, |φYj (z)| ≤ φYj (iy) = E[e−yYj ]. By Eq. (246) at z = iy, φYj (iy) φYl (iy) = exp( 12 σ 2 y 2 ), and Jensen’s inequality gives φYl (iy) = E[e−yYl ] ≥ exp(−y E[Yl ]). Hence, with C := |E[Yl ]|,     exp 12 σ 2 y 2 φYj (z) ≤ ≤ exp 12 σ 2 y 2 + C|y| ≤ exp 21 σ 2 |z|2 + C|z| . φYl (iy)

(247)

Step 4: conclusion. Since φYj is entire and nowhere zero on the simply connected domain C, there is an entire hj with φYj = ehj , and we may normalise hj (0) = 0 because φYj (0) = 1. By Eq. (247), Re hj (z) = log φYj (z) ≤ 21 σ 2 |z|2 + C|z|,

(248)

so in the notation of Lemma A.5, M (r) ≤ 12 σ 2 r2 + Cr. For the Taylor coefficients am of hj and every m ≥ 3, Eq. (239) gives |am | ≤

4M (r) 2σ 2 r2 + 4Cr ≤ −−−→ 0, r→∞ rm rm

(249)

so am = 0 for all m ≥ 3 and hj is a polynomial of degree at most 2. Thus φYj = exp(hj ) with a polynomial hj , and Marcinkiewicz’ Theorem 3.17 shows that Yj is Gaussian, possibly degenerate.

B

A Proof of the Kagan–Linnik–Rao Theorem

This appendix proves Theorem 4.2, the one result that Sections 4 to 7 take as given. Its first part – the equality of images and ranks – was already proved in Proposition 4.4; what remains is the column dichotomy Theorem 4.2 (1)–(2). The strategy is the classical one that goes back to Marcinkiewicz (1939), Darmois (1953) and Skitovich (1954), and that is developed systematically in Kagan et al. (1973, Chapters 3 and 10) and Linnik and Ostrovskii (1977). It has three steps. 59

(1) Take distinguished logarithms in the identity between the two characteristic functions. Because each representation is an independent sum of one-dimensional sources evaluated along the columns of the mixing matrix, the identity becomes a linear relation among ridge functions t 7→ ψ(aT t) on Rp . (2) Kill all but one of these ridge functions with finite-difference operators: a ridge function in direction a is annihilated by the difference operator in any direction h ⊥ a, and the hypothesis that the columns are pairwise non-proportional is exactly what guarantees that each unwanted direction admits such an h that does not annihilate the term we want to keep. (3) What survives is a one-dimensional function all of whose iterated differences vanish. By Fréchet’s functional equation such a function is a polynomial, and Marcinkiewicz’ Theorem 3.17 converts “polynomial cumulant generating function” into “Gaussian”. Everything below uses only the elementary properties of characteristic functions collected in Proposition 3.2 – in particular Eqs. (8) and (9) – the distinguished logarithm of Lemma 3.10, Marcinkiewicz’ Theorem 3.17, Proposition 4.4 for the first part, and elementary calculus.

B.1

Finite differences

Definition B.1 (Difference operator). For a function F defined on a subset of Rp and a vector h ∈ Rp we write ∆h F (t) := F (t + h) − F (t),

(250)

whenever both t and t + h lie in the domain of F . Difference operators in different directions commute, and iterating gives the inclusion–exclusion formula ∆h1 · · · ∆hm F (t) =

X



(−1)m−|S| F t +



hr ,

(251)

r∈S

S⊆{1,...,m}

which requires F to be defined at the 2m points t +

X

r∈S hr .

P

A function of the form t 7→ ψ(aT t), constant along the hyperplane a⊥ , is called a ridge function in the direction a. The next two lemmas are the whole geometric input. Lemma B.2 (Ridge functions are annihilated by orthogonal differences). Let a ∈ Rp , let ψ be a function on a subset of R and put F (t) := ψ(aT t). If h ∈ Rp satisfies aT h = 0, then ∆h F = 0 on the domain of F . Proof. F (t + h) = ψ aT t + aT h = ψ(aT t) = F (t).

Lemma B.3 (Separating two non-proportional directions). Let b, d ∈ Rp \ {0} be nonproportional. Then there exists h ∈ Rp with dT h = 0,

bT h ̸= 0.

(252)

Moreover the set of such h is open in the hyperplane d⊥ , so h may be chosen with ∥h∥ arbitrarily small, and sh is admissible for every s ̸= 0. Proof. Suppose no such h exists, i.e. d⊥ ⊆ b⊥ . Taking orthogonal complements reverses the
⊥
⊥ inclusion, so span{b} = b⊥ ⊆ d⊥ = span{d}. As b = ̸ 0 this makes b a non-zero multiple of d, contradicting non-proportionality. The set {h ∈ d⊥ : bT h ̸= 0} is the complement in d⊥ of the closed subspace b⊥ ∩ d⊥ , hence open in d⊥ ; and it is invariant under h 7→ sh for s ̸= 0, so it contains vectors of arbitrarily small norm. 60

The analytic input is the local form of Fréchet’s functional equation. We prove it from scratch; the argument is a smoothing bootstrap followed by one differentiation. Lemma B.4 (Fréchet’s functional equation, local form). Let m ∈ N and let ρ > σ > 0 and ε > 0 satisfy σ + mε ≤ ρ. Let F : (−ρ, ρ) → C be continuous and assume ∆c1 · · · ∆cm F (u) = 0

whenever |c1 |, . . . , |cm | < ε and |u| < σ.

(253)

Then F coincides on − σ4 , σ4 with a polynomial of degree at most m − 1.

Proof. Note first that the condition σ + mε ≤ ρ makes Eq. (253) meaningful: by Eq. (251) P the left hand side only involves the values of F at the points u + r∈S cr , all of which lie in (−ρ, ρ). Put 

κ := 12 min ε,

σ  > 0, 2m2

(254)

and define the symmetric averaging operator 1 Mκ G (u) := 2κ

1 G(u + s) ds = 2κ −κ

Z κ

Z u+κ

G(v) dv, u−κ

(255)

defined for |u| < A − κ whenever G is continuous on (−A, A). Step 1: Mκ raises smoothness. If G is continuous on (−A, A) then, by the fundamental theorem of calculus, Mκ G is continuously differentiable on (−A + κ, A − κ) with
′

Mκ G (u) =

1 G(u + κ) − G(u − κ) . 2κ

(256)

Consequently, if G ∈ C r on (−A, A) then Mκ G ∈ C r+1 on (−A + κ, A − κ). Step 2: the hypothesis says (I − Mκ )m F = 0. We have F (u) − Mκ F (u) = −

1 2κ

Z κ −κ

∆s F (u) ds,

(257)

and translations – hence all the operators ∆s and Mκ – commute with one another, so Fubini’s theorem gives (I − Mκ )m F (u) =





(−1)m (2κ)m

Z κ

Z κ

···

−κ

−κ

∆s1 · · · ∆sm F (u) ds1 · · · dsm = 0

(258)

for |u| < σ: every sr occurring in the integral satisfies |sr | ≤ κ < ε, so the integrand vanishes identically by Eq. (253). Step 3: bootstrap. Expanding the binomial in Eq. (258) and solving for the n = 0 term, F =

m X

!

(−1)

n=1

n+1

m Mκn F n

on (−σ, σ).

(259)

Every term on the right applies Mκ at least once. Hence, by Step 1, if F ∈ C r on (−A, A) for some A ≤ ρ, then the right
hand side of Eq. (259) is C r+1 on (−A + mκ, A − mκ), and therefore F ∈ C r+1 on −A′ , A′ with A′ := min(σ, A − mκ). Starting from F ∈ C 0 on (−ρ, ρ) and iterating m times, and using ρ − mκ ≥ ρ − mε ≥ σ at the first step, we obtain F ∈ C m on (−Am , Am ) with Am ≥ σ − (m − 1)mκ ≥ σ − m2 κ ≥ 43 σ, where the last inequality uses κ ≤ σ/(4m2 ). 61

Step 4: differentiate. Fix u with |u| < σ/4 and let c1 , . . . , cm range over (−δ, δ), where
δ := 12 min ε, σ/(2m) , so that |u| + mδ < 34 σ ≤ Am and δ < ε. By Eq. (251) and Eq. (253), 

(−1)m−|S| F u +

X

X



cr = 0,

(260)

r∈S

S⊆{1,...,m}

identically in (c1 , . . . , cm ) ∈ (−δ, δ)m . All terms are C m functions of (c1 , . . . , cm ) by Step 3, so we may apply ∂ m /∂c1 · · · ∂cm . A summand indexed by S = ̸ {1, . . . , m} does not depend on cr0 for any r0 ∈ / S and is therefore annihilated; the summand indexed by S = {1, . . . , m} contributes F (m) (u + c1 + · · · + cm ). Hence F (m) (u + c1 + · · · + cm ) = 0, and setting
c1 = · · · = cm = 0 gives F (m) (u) = 0. As u was arbitrary in the interval − σ4 , σ4 , which is connected, F is there a polynomial of degree at most m − 1. Remark B.5 (What the hypotheses of Lemma B.4 do and do not say). Two features of Eq. (253) are essential. First, the steps c1 , . . . , cm are allowed to vary over a neighbourhood of the origin, which is what makes the averaging in Eq. (258) legitimate; for a single fixed step vector the conclusion is false, as the 1-periodic function F (u) = sin(2πu) with m = 1, c1 = 1 shows. Second, only continuity of F is assumed; no differentiability is available a priori, and Step 3 is precisely what manufactures it. The global version of the statement – a continuous F on R with ∆c1 · · · ∆cm F ≡ 0 for all c1 , . . . , cm ∈ R is a polynomial of degree < m – is due to Fréchet (1909), and the equation carries his name. Restricted-domain versions, in which the increments and the base point are confined to a subset, have been studied since; Ger (1994) shows that a solution on a sufficiently regular domain always extends to a polynomial function on the whole space. We have given a direct proof of the one-dimensional continuous case, which is all that is needed below and is short enough not to warrant a detour.

B.2

Proof of the column dichotomy

Proof of Theorem 4.2 (1)–(2). Write ki := k (i) , let a1 , . . . , ak1 be the columns of A(1) and b1 , . . . , bk2 the columns of A(2) , and put ϕj := φZ (1) ,

γl := φZ (2) ,

j

c := µ(2) − µ(1) .

l

(261)

By Eqs. (8) and (9), the two representations Eq. (54) give, for every t ∈ Rp , exp i tT µ(1)

k1
Y

ϕj aTj t = φX (t) = exp i tT µ(2)

j=1

k2
Y

(262)

γl bTl t .

l=1

Step 1: the logarithmic form. By Proposition 3.2 (i)–(ii) each ϕj and each γl is continuous and equals 1 at the origin, so there is a δ0 > 0 such that all of them are non-vanishing on (−δ0 , δ0 ), and Lemma 3.10 supplies the cumulant generating functions ψj := ψZ (1) and j

χl := ψZ (2) , continuous on (−δ0 , δ0 ) with ψj (0) = χl (0) = 0 and ϕj = exp(ψj ), γl = exp(χl ) l

there. Choose R > 0 so small that |aTj t| < δ0 and |bTl t| < δ0 for all j, l whenever ∥t∥ < R, and define on the ball BR (0) ⊆ Rp Φ(t) :=

k1 X

ψj aTj t −

j=1

k2 X

χl bTl t − i tT c.

(263)

l=1

Then Φ is continuous, Φ(0) = 0, and by Eq. (262) exp(Φ) ≡ 1 on BR (0). Hence Φ/(2πi) is a continuous, integer-valued function on the connected set BR (0) vanishing at the origin, so it 62

vanishes identically: k1 X

ψj aTj t −

j=1

k2 X

χl bTl t = i tT c,

t ∈ BR (0).

(264)

l=1

This is a linear relation among ridge functions, and it is the only consequence of the hypotheses that we shall use. Step 2: choosing the directions. We treat both parts at once. Fix the index l0 ∈ {1, . . . , k2 } of the column under consideration and let (

e :=

bl0 , a j0 ,

in case (1), in case (2), where bl0 = λaj0 ,

(265)

and let D be the multiset of the remaining column directions, (

D :=

{a1 , . . . , ak1 } ⊎ {bl : l ̸= l0 }, {aj : j ̸= j0 } ⊎ {bl : l ̸= l0 },

in case (1), in case (2),

(266)

where ⊎ is the disjoint union, so that |D| = k1 + k2 − 1 in case (1) and k1 + k2 − 2 in case (2) even if the same vector happens to be a column of both matrices. We claim that no d ∈ D is proportional to e. In case (1) this holds for d = aj by the hypothesis of part (1) and for d = bl , l ̸= l0 , because the columns of A(2) are pairwise non-proportional. In case (2) it holds for d = aj , j ̸= j0 , because the columns of A(1) are pairwise non-proportional, and for d = bl , l ̸= l0 , because bl ∥ aj0 would give bl ∥ λaj0 = bl0 , again contradicting pairwise non-proportionality. All columns are non-zero by hypothesis (i), so e ̸= 0 and d ̸= 0. Let m0 := |D| and set m := max(m0 , 2). For r = 1, . . . , m0 apply Lemma B.3 to the pair (e, dr ): it yields hr ∈ Rp with dTr hr = 0 and αr := eT hr ̸= 0, and hr may be taken with ∥hr ∥ < R/(2m). For r = m0 + 1, . . . , m (there are at most two such r, and only when m0 < 2) simply put hr := τ e/∥e∥2 with τ > 0 small, so that αr = τ ̸= 0 and ∥hr ∥ < R/(2m) as well. Note that each hr may be replaced by sr hr for any sr ∈ (−1, 1) \ {0} without disturbing any of these properties, and that this replaces αr by sr αr . Step 3: annihilating all but one term. Apply the commuting operators ∆h1 , . . . , ∆hm to Eq. (264). This is legitimate on BR/2 (0), because ∥t∥ < R/2 and ∥hr ∥ < R/(2m) force P ∥t + r∈S hr ∥ < R for every S ⊆ {1, . . . , m}. Each ridge function indexed by a direction dr ∈ D is annihilated by ∆hr (Lemma B.2), hence by the whole composition; and the affine term i tT c is annihilated by any two of the operators, of which there are m ≥ 2. What survives is (

∆h1 · · · ∆hm Ξ e t 

T


= 0,

∥t∥ < R/2,

where

Ξ :=

−χl0 , ψj0 − χl0 (λ · ),

in case (1), in case (2). (267)

In case (2) we used that bTl0 t = λ aTj0 t = λ eT t, so that the two surviving ridge functions share the direction e and combine into a single one. Step 4: Fréchet. Set ρ := ∥e∥R. Every u with |u| < ρ is of the form u = eT t for some t ∈ BR (0), and for such t the choice of R puts every aTj t and every bTl t inside (−δ0 , δ0 ). In case (1), Ξ = −χl0 is evaluated at u = bTl0 t; in case (2) it is built from ψj0 evaluated at u = aTj0 t and from χl0 evaluated at λu = bTl0 t. In either case all arguments lie in (−δ0 , δ0 ), so Ξ is defined and continuous on (−ρ, ρ). Because e = ̸ 0, the linear functional t 7→ eT t maps BR/2 (0)

63

onto (−σ, σ) with σ := ∥e∥R/2 = ρ/2, and by Eq. (251) the left hand side of Eq. (267) equals
∆α1 · · · ∆αm Ξ (eT t). Hence ∆α1 · · · ∆αm Ξ (u) = 0,

(268)

|u| < σ.

Replacing each hr by sr hr with sr ∈ (−1, 1), as permitted in Step 2, replaces αr by sr αr without affecting anything else, and as sr ranges over (−1, 1) the number sr αr ranges over the whole interval (−|αr |, |αr |). So Eq. (268) holds with αr replaced by any cr with 0 < |cr | < ε := minr |αr |; and for cr = 0 it holds trivially, because ∆0 = 0. This is exactly the hypothesis Eq. (253). Finally |αr | = |eT hr | ≤ ∥e∥∥hr ∥ < ∥e∥R/(2m) = ρ/(2m), so σ + mε < ρ/2 + ρ/2 = ρ and Lemma B.4 applies. It makes Ξ a polynomial of degree at most
σ σ m − 1 on − 4 , 4 , a neighbourhood of the origin in R. Step 5: conclusion. In case (1) we have obtained that χl0 = −Ξ is a polynomial near the origin, so φZ (2) (t) = exp χl0 (t)

(269)

l0

(2)

near the origin with a polynomial exponent. By Marcinkiewicz’ Theorem 3.17, Zl0 is a (possibly degenerate) Gaussian random variable, which is assertion (1). In case (2), Ξ(t) = ψj0 (t) − χl0 (λt) is a polynomial near the origin; writing g := −Ξ and exponentiating, φZ (2) (λt) = exp χl0 (λt) = exp ψj0 (t) exp g(t) = φZ (1) (t) · exp g(t)

j0

l0

(270)

near the origin, which is Eq. (56). For the final assertion of (2), note that this is exactly the situation of Remark 3.18 (i) with the factor λ present, which is why that remark was stated (2) (1) with the factor: it gives directly that Zl0 is Gaussian if and only if Zj0 is. (The factor cannot (2)

be dropped: by Eq. (8) the left hand side is φλZ (2) (t), and it is only because λ = ̸ 0 that λZl0 l0

(2)

being Gaussian is equivalent to Zl0 being Gaussian; cf. Footnote 1.) Remark B.6 (Reading the proof). Three points are worth isolating. Where each hypothesis enters. That the columns are non-zero is used twice: in Step 2, to make Lemma B.3 applicable to the pairs (e, dr ), and in Step 4, where it makes eT t surject onto an interval, so that a statement about t ∈ Rp becomes a statement about u ∈ R. Pairwise non-proportionality is used only in Step 2, to make Lemma B.3 applicable, and only for those columns that have to be separated from e: in case (2) for the columns of both matrices, in case (1) for the columns of A(2) only – there the separation of bl0 from the columns of A(1) is the hypothesis of part (1) itself, and the columns of A(1) may perfectly well be proportional to one another. Mutual independence of the components, hypothesis (iii)(a), is used once, to factorise the characteristic functions in Eq. (262). Non-constancy of the sources, hypothesis (iii)(b), is not needed for (1)–(2) at all; it enters only through Proposition 4.4, i.e. in the first part of Theorem 4.2. Why the degree bound is irrelevant. Step 4 produces a polynomial of degree at most m − 1, where m = max(k1 + k2 − 1, 2) in case (1) and m = max(k1 + k2 − 2, 2) in case (2); this bound grows with the number of sources and is of no use by itself. It is Marcinkiewicz’ theorem that collapses it to 2 – which is precisely why Theorem 3.17 was called the single most important analytic input of these notes in Section 3.4. Why the argument is local. Everything happens on a ball around the origin, because Lemma 3.10 produces the cumulant generating functions only there. This is not a defect: by Remark 3.8 a local identity between characteristic functions carries no information about 64

the laws in general, and the reason it does here is that the local conclusion is “exp of a polynomial”, a form rigid enough for Theorem 3.17 to upgrade it to a global statement about the distribution. Remark B.7 (Relation to the Darmois–Skitovich theorem). The classical theorem of Darmois (1953) and Skitovich (1954) says this: if Y1 , . . . , Yk are independent and the two linear forms P P L1 = j αj Yj and L2 = j βj Yj are independent of each other, then every Yj with αj βj ̸= 0 is Gaussian. It is the case p = 2 of Theorem 4.2. Indeed, put X := [L1 , L2 ]T ∈ R2 and compare the two representations "

I2

L1 L2

#

"

= X =

α1 · · · β1 · · · (2)

#

αk Y, βk

(271) (2)

so that k (1) = 2 with A(1) = I2 , and k (2) = k with aj = [αj , βj ]T . A non-zero column aj is proportional to a column of I2 precisely when αj βj = 0 (zero columns having been removed by Remark 4.1), so for every j with αj βj ̸= 0 part (1) declares Yj Gaussian. Two caveats: (2) the columns aj need not be pairwise non-proportional, so one first passes to a normalised representation as in Remark 4.1, which merges the sources belonging to proportional columns; and one then recovers the individual Yj from the merged sum by finitely many applications of Cramér’s Theorem 3.20, the case of an a.s. constant merged source being covered by the convention that constants are degenerate Gaussians (Notation 2.4). The two-representation formulation of Theorem 4.2 is what makes the result directly usable for identifiability, since it compares two models rather than two forms, and it is proved by the same finite-difference argument. See Kagan et al. (1973, Chapter 3) for the classical statement and Ghurye and Olkin (1962) for a multivariate version.

References S. Amari, A. Cichocki, and H. H. Yang. A new learning algorithm for blind signal separation. In D. S. Touretzky, M. C. Mozer, and M. E. Hasselmo, editors, Advances in Neural Information Processing Systems 8 (NIPS 1995), pages 757–763. MIT Press, Cambridge, MA, 1996. https://papers.nips.cc/paper/1115-a-new-learning-algorithm-for-blindsignal-separation. S. Amari, T.-P. Chen, and A. Cichocki. Stability analysis of learning algorithms for blind source separation. Neural Networks, 10(8):1345–1351, 1997. doi: 10.1016/S0893-6080(97)000397. S. Amari and J.-F. Cardoso. Blind source separation—semiparametric statistical approach. IEEE Transactions on Signal Processing, 45(11):2692–2700, 1997. doi: 10.1109/78.650095. S. Amari. Natural gradient works efficiently in learning. Neural Computation, 10(2):251–276, 1998. doi: 10.1162/089976698300017746. K. P. Balanda and H. L. MacGillivray. Kurtosis: a critical review. The American Statistician, 42(2):111–119, 1988. doi: 10.1080/00031305.1988.10475539. A. J. Bell and T. J. Sejnowski. An information-maximization approach to blind separation and blind deconvolution. Neural Computation, 7(6):1129–1159, 1995. doi: 10.1162/neco. 1995.7.6.1129. J.-F. Cardoso and B. H. Laheld. Equivariant adaptive source separation. IEEE Transactions on Signal Processing, 44(12):3017–3030, 1996. doi: 10.1109/78.553476. J.-F. Cardoso. Infomax and maximum likelihood for blind source separation. IEEE Signal Processing Letters, 4(4):112–114, 1997. doi: 10.1109/97.566704. 65

J.-F. Cardoso. Multidimensional independent component analysis. In Proceedings of the 1998 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP ’98), volume 4, pages 1941–1944, Seattle, WA, 1998a. IEEE. doi: 10.1109/ICASSP.1998.681443. J.-F. Cardoso. Blind signal separation: statistical principles. Proceedings of the IEEE, 86(10):2009–2025, 1998b. doi: 10.1109/5.720250. P. Comon. Independent component analysis, a new concept? Signal Processing, 36(3):287–314, 1994. doi: 10.1016/0165-1684(94)90029-9. P. Comon and C. Jutten, editors. Handbook of Blind Source Separation: Independent Component Analysis and Applications. Academic Press (Elsevier), Oxford, 1st edition, 2010. ISBN 978-0-12-374726-6 (print), 978-0-08-088494-3 (e-book). H. Cramér. Über eine Eigenschaft der normalen Verteilungsfunktion. Mathematische Zeitschrift, 41(1):405–414, 1936. doi: 10.1007/BF01180430. R. B. Darlington. Is kurtosis really “peakedness?” The American Statistician, 24(2):19–22, 1970. doi: 10.1080/00031305.1970.10478885. G. Darmois. Analyse générale des liaisons stochastiques: étude particulière de l’analyse factorielle linéaire. Revue de l’Institut International de Statistique / Review of the International Statistical Institute, 21(1/2):2–8, 1953. doi: 10.2307/1401511. J. Eriksson and V. Koivunen. Identifiability, separability, and uniqueness of linear ICA models. IEEE Signal Processing Letters, 11(7):601–604, 2004. doi: 10.1109/LSP.2004.830118. J. Eriksson and V. Koivunen. Complex random vectors and ICA models: identifiability, uniqueness, and separability. IEEE Transactions on Information Theory, 52(3):1017–1029, 2006. doi: 10.1109/TIT.2005.864440. Preprint: https://arxiv.org/abs/cs/0512063. W. Feller. An Introduction to Probability Theory and Its Applications, volume II. Wiley Series in Probability and Mathematical Statistics. John Wiley & Sons, New York, 2nd edition, 1971. 669 pp. ISBN 978-0-471-25709-7. M. Fréchet. Une définition fonctionnelle des polynomes. Nouvelles annales de mathématiques, 4e série, 9:145–162, 1909. https://www.numdam.org/item/NAM_1909_4_9__145_0/. R. Ger. On extensions of polynomial functions. Results in Mathematics, 26(3–4):281–289, 1994. doi: 10.1007/BF03323050. S. G. Ghurye and I. Olkin. A characterization of the multivariate normal distribution. The Annals of Mathematical Statistics, 33(2):533–541, 1962. doi: 10.1214/aoms/1177704579. L. Gresele, J. von Kügelgen, V. Stimper, B. Schölkopf, and M. Besserve. Independent mechanism analysis, a new concept? In Advances in Neural Information Processing Systems 34 (NeurIPS 2021), pages 28233–28248. Curran Associates, 2021. https://proceedings. neurips.cc/paper/2021/hash/edc27f139c3b4e4bb29d1cdbc45663f9- Abstract.html. Preprint: https://arxiv.org/abs/2106.05200. P. O. Hoyer, D. Janzing, J. M. Mooij, J. Peters, and B. Schölkopf. Nonlinear causal discovery with additive noise models. In Advances in Neural Information Processing Systems 21 (NIPS 2008), pages 689–696, 2008. https://proceedings.neurips.cc/paper/2008/ hash/f7664060cc52bc6f3d620bcedc94a4b6-Abstract.html. A. Hyvärinen and H. Morioka. Unsupervised feature extraction by time-contrastive learning and nonlinear ICA. In Advances in Neural Information Processing Systems 29 (NIPS 2016), pages 3765–3773. Curran Associates, 2016. https://proceedings.neurips.cc/paper/ 2016/hash/d305281faf947ca7acade9ad5c8c818c-Abstract.html. Preprint: https:// arxiv.org/abs/1605.06336. A. Hyvärinen and H. Morioka. Nonlinear ICA of temporally dependent stationary sources. 66

In A. Singh and J. Zhu, editors, Proceedings of the 20th International Conference on Artificial Intelligence and Statistics (AISTATS), volume 54 of Proceedings of Machine Learning Research, pages 460–469. PMLR, 2017. https://proceedings.mlr.press/v54/ hyvarinen17a.html. A. Hyvärinen and E. Oja. Independent component analysis: algorithms and applications. Neural Networks, 13(4–5):411–430, 2000. doi: 10.1016/S0893-6080(00)00026-5. A. Hyvärinen and P. Pajunen. Nonlinear independent component analysis: existence and uniqueness results. Neural Networks, 12(3):429–439, 1999. doi: 10.1016/S08936080(98)00140-3. A. Hyvärinen, J. Karhunen, and E. Oja. Independent Component Analysis. Wiley Series on Adaptive and Learning Systems for Signal Processing, Communications, and Control. John Wiley & Sons, New York, 2001. 504 pp. ISBN 978-0-471-40540-5. doi: 10.1002/0471221317. A. M. Kagan, Yu. V. Linnik, and C. R. Rao. Characterization Problems in Mathematical Statistics. Wiley Series in Probability and Mathematical Statistics. John Wiley & Sons, New York, 1973. xii+499 pp. ISBN 978-0-471-45421-2. Translated from the Russian by B. Ramachandran; Russian original: Kharakterizatsionnye zadachi matematicheskoi statistiki, Nauka, Moscow, 1972. O. Kallenberg. Foundations of Modern Probability, volume 99 of Probability Theory and Stochastic Modelling. Springer, Cham, 3rd edition, 2021. xii+946 pp. ISBN 978-3-03061870-4. doi: 10.1007/978-3-030-61871-1. I. Kaplansky. A common error concerning kurtosis. Journal of the American Statistical Association, 40(230):259, 1945. doi: 10.1080/01621459.1945.10501856. I. Khemakhem, D. P. Kingma, R. P. Monti, and A. Hyvärinen. Variational autoencoders and nonlinear ICA: a unifying framework. In S. Chiappa and R. Calandra, editors, Proceedings of the 23rd International Conference on Artificial Intelligence and Statistics (AISTATS), volume 108 of Proceedings of Machine Learning Research, pages 2207–2217. PMLR, 2020. https://proceedings.mlr.press/v108/khemakhem20a.html. Preprint: https://arxiv. org/abs/1907.04809. T. Kim, T. Eltoft, and T.-W. Lee. Independent vector analysis: an extension of ICA to multivariate components. In Independent Component Analysis and Blind Signal Separation (ICA 2006), volume 3889 of Lecture Notes in Computer Science, pages 165–172. Springer, Berlin, Heidelberg, 2006. ISBN 978-3-540-32630-4. doi: 10.1007/11679363_21. A. Klenke. Probability Theory: A Comprehensive Course. Universitext. Springer, Cham, 3rd edition, 2020. ISBN 978-3-030-56401-8. doi: 10.1007/978-3-030-56402-5. Lévy’s continuity theorem is Section 15.3. T.-W. Lee, M. Girolami, and T. J. Sejnowski. Independent component analysis using an extended infomax algorithm for mixed subgaussian and supergaussian sources. Neural Computation, 11(2):417–441, 1999. doi: 10.1162/089976699300016719. Yu. V. Linnik and I. V. Ostrovskii. Decomposition of Random Variables and Vectors, volume 48 of Translations of Mathematical Monographs. American Mathematical Society, Providence, RI, 1977. ix+380 pp. ISBN 978-0-8218-1598-4. Translated from the Russian; translation edited by J. Rosenblatt. doi: 10.1090/mmono/048. E. Lukacs. Characteristic Functions. Charles Griffin & Company Limited, London, 2nd, revised and enlarged edition, 1970. x+350 pp. ISBN 0-85264-170-2. LCCN 70-513840. D. J. C. MacKay. Maximum likelihood and covariant algorithms for independent component analysis. Unpublished report, Cavendish Laboratory, University of Cambridge, 1996. https://www.inference.org.uk/mackay/ica.pdf. Version 3.8, 8 January 1999, with 67

minor corrections of 8 October 2002. D. J. C. MacKay. Information Theory, Inference, and Learning Algorithms. Cambridge University Press, Cambridge, 2003. xii+628 pp. ISBN 978-0-521-64298-9. Chapter 34, “Independent Component Analysis and Latent Variable Modelling”, pp. 437–444. https: //www.inference.org.uk/mackay/itila/. J. Marcinkiewicz. Sur une propriété de la loi de Gauss. Mathematische Zeitschrift, 44(1):612– 618, 1939. doi: 10.1007/BF01210677. J. J. A. Moors. The meaning of kurtosis: Darlington reexamined. The American Statistician, 40(4):283–284, 1986. doi: 10.1080/00031305.1986.10475415. T. Pandeva and P. Forré. Multi-view independent component analysis with shared and individual sources. In R. J. Evans and I. Shpitser, editors, Proceedings of the 39th Conference on Uncertainty in Artificial Intelligence (UAI), volume 216 of Proceedings of Machine Learning Research, pages 1639–1650. PMLR, 2023a. https://proceedings.mlr.press/ v216/pandeva23a.html. Preprint: https://arxiv.org/abs/2210.02083. T. Pandeva and P. Forré. Multi-view independent component analysis for omics data integration. ICLR 2023 Workshop on Machine Learning and Global Health, 2023b. https://openreview.net/forum?id=r5KL-AfXt75. T. Pandeva, M. J. Jonker, L. Hamoen, J. Mooij, and P. Forré. Robust multi-view coexpression network inference. In B. Huang and M. Drton, editors, Proceedings of the 4th Conference on Causal Learning and Reasoning (CLeaR), volume 275 of Proceedings of Machine Learning Research, pages 490–513. PMLR, 2025. https://proceedings.mlr. press/v275/pandeva25a.html. Preprint: https://arxiv.org/abs/2409.19991. G. Pólya. Remarks on characteristic functions. In J. Neyman, editor, Proceedings of the Berkeley Symposium on Mathematical Statistics and Probability, pages 115–123. University of California Press, Berkeley and Los Angeles, 1949. https://digitalassets.lib.berkeley. edu/math/ucb/text/math_s1_article-08.pdf. H. Robbins and S. Monro. A stochastic approximation method. The Annals of Mathematical Statistics, 22(3):400–407, 1951. doi: 10.1214/aoms/1177729586. S. Shimizu, P. O. Hoyer, A. Hyvärinen, and A. Kerminen. A linear non-Gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7:2003–2030, 2006. https://jmlr.org/papers/v7/shimizu06a.html. S. Shimizu, T. Inazumi, Y. Sogawa, A. Hyvärinen, Y. Kawahara, T. Washio, P. O. Hoyer, and K. Bollen. DirectLiNGAM: a direct method for learning a linear non-Gaussian structural equation model. Journal of Machine Learning Research, 12:1225–1248, 2011. https: //jmlr.org/papers/v12/shimizu11a.html. V. P. Skitovich. Linear forms of independent random variables and the normal distribution law. Izvestiya Akademii Nauk SSSR, Seriya Matematicheskaya, 18(2):185–200, 1954. In Russian. https://www.mathnet.ru/eng/im3497. Announced in Doklady Akademii Nauk SSSR (N.S.), 89:217–219, 1953. P. H. Westfall. Kurtosis as peakedness, 1905–2014. R.I.P. The American Statistician, 68(3):191– 195, 2014. doi: 10.1080/00031305.2014.917055.

68

Related documents

• Watermarking of audio signals using independent component analysis#230081
• Evaluating multi-level membership inference risk in federated EEG learning.#340419
• Singer-Informed Vocal Source Separation for Multi-Singer Music Mixtures#609674
• On theoretical guarantees and a blessing of dimensionality for nonconvex sampling#635556
• On theoretical guarantees and a blessing of dimensionality for nonconvex sampling#633207