Asymptotic Performance and the Cramér–Rao Bound渐近性能与 Cramér–Rao 界
The gain matrix and interference-to-signal ratio, asymptotic variance of contrast-based estimators, the optimal nonlinearity, and the Cramér–Rao bound for ICA.
An ICA algorithm computes its estimate from a finite sample, so even when the model holds exactly, each output keeps a small amount of the other sources. The script of Section 1.5第 1.5 节 printed WA≈(0.03−1.001.000.02): separation succeeded, but the entries 0.03 and 0.02 are not zero, and they would change with a new sample. This chapter quantifies such residual errors. How large are they? How do they depend on the sample size T, on the distributions of the sources and on the algorithm? Is there a limit that no algorithm can beat?
The answers are asymptotic: they describe the leading behaviour of the errors as T→∞, and in practice they are accurate already for moderate T. Three facts organize the chapter. For the usual estimators the distribution of the errors does not depend on the mixing matrix (Section 10.2第 10.2 节). The errors shrink like 1/T, with constants that can be computed in closed form; for one-unit FastICA the constant is the asymptotic variance Vg of Hyvärinen (1997; 1999), and the best nonlinearity is the score function of the source (Section 10.5第 10.5 节). Finally, the Cramér–Rao bound (Tichavský, Koldovský & Oja, 2006) gives the smallest variance that any unbiased estimator can achieve; it is governed by the Fisher information of the sources and becomes infinite when two sources are Gaussian (Section 10.7第 10.7 节).
Throughout the chapter the model is square (m=n, A invertible), the sources satisfy Assumption 1.1假设 1.1 with zero mean and unit variance, and their densities pi are smooth enough for every expectation below to exist and for differentiation under the integral sign. The score function is ψi=−pi′/pi and the Fisher information is Ji=E[ψi2(si)] (Definition 6.1定义 6.1). One identity is used repeatedly. If f is differentiable and f(s)p(s)→0 as ∣s∣→∞, integration by parts gives ∫f′pds=−∫fp′ds=∫fψpds, that is,
With f=1, f(s)=s and f=ψ it yields E[ψ(s)]=0, E[sψ(s)]=1 and E[ψ′(s)]=J. For a standard Gaussian ψ(s)=s, so J=1; every other unit-variance law has J>1, because E[(ψ(s)−s)2]=J−2+1=J−1 vanishes only when ψ(s)=s.
10.2The gain matrix and the interference-to-signal ratio增益矩阵与干扰信号比
Let W^ be an estimate of the unmixing matrix computed from T samples. Its quality is judged not by W^ itself, which depends on the arbitrary scaling of the outputs, but by what it does to the sources.
设 W^ 是由 T 个样本算出的解混矩阵估计。评价它的好坏,不看 W^ 本身(它依赖于输出的任意尺度),而看它对源做了什么。
Definition 10.1 (gain matrix and its error)定义 10.1(增益矩阵及其误差)
The gain matrix of W^ is P^=W^A, so that the outputs are y^=P^s. Order the outputs so that output i estimates source i, and let D^=diag(P^11,…,P^nn). The gain error is the matrix ε with zero diagonal defined by
W^ 的增益矩阵是 P^=W^A,于是输出为 y^=P^s。把输出排好顺序,使第 i 个输出估计第 i 个源,并记 D^=diag(P^11,…,P^nn)。增益误差是对角元为零的矩阵 ε,定义为
D^−1P^=I+ε,εij=P^iiP^ij(i=j).(10.2)
Output i is then P^ii(si+∑j=iεijsj): after the unrecoverable scale P^ii has been divided out, εij is the amount of source j that leaks into the estimate of source i. Perfect separation means ε=0.
于是第 i 个输出等于 P^ii(si+∑j=iεijsj):除掉无法恢复的尺度 P^ii 之后,εij 就是源 j 漏进源 i 的估计中的分量。完全分离意味着 ε=0。
The interference-to-signal ratio of source j in output i is ISRij=E[εij2], i=j, where the expectation is over the random sample. The total ISR of output i is ∑j=iISRij.
源 j 在输出 i 中的干扰信号比(ISR)定义为 ISRij=E[εij2](i=j),期望对随机样本求取。输出 i 的总 ISR 为 ∑j=iISRij。
Because the sources have unit variance, εij2 is the ratio of the power that source j contributes to output i to the power of the wanted source in that output. ISRs are often quoted in decibels, 10log10ISR; the reciprocal is the signal-to-interference ratio (SIR). In a simulation A is known, so ε can be computed in every run and the ISR estimated by averaging εij2 over runs. For the estimators of this chapter, Tεij converges in distribution to a normal law with mean zero and some variance vij, the asymptotic variance. As is customary, we then use ISRij≈vij/T, replacing the exact mean square by its asymptotic value; Section 10.6第 10.6 节 checks this approximation numerically. The constant vij measures the statistical efficiency of a method independently of T.
由于源的方差为一,εij2 就是源 j 在输出 i 中贡献的功率与目标源功率之比。ISR 常以分贝表示,即 10log10ISR;其倒数是信干比(SIR)。做模拟时 A 已知,因此每次运行都能算出 ε,再对多次运行的 εij2 求平均即可估计 ISR。对本章讨论的估计量,Tεij 依分布收敛到均值为零、方差为某个 vij 的正态分布,vij 称为渐近方差。按惯例,我们用 ISRij≈vij/T,即以渐近值代替精确的均方;Section 10.6第 10.6 节用数值方法检验这一近似。常数 vij 衡量的是方法本身的统计效率,与 T 无关。
Why the mixing matrix does not matter为什么混合矩阵无关紧要
Most ICA estimators are defined by equations that involve the data only through the outputs. Maximum likelihood and Infomax solve E^[φ(y)y⊤]=I; the fixed points of FastICA satisfy whitening and fixed-point conditions on y; JADE optimizes cumulants of y. Such estimators are equivariant (Cardoso & Laheld, 1996; see also Proposition 7.1命题 7.1), and equivariance has a strong consequence for performance.
大多数 ICA 估计量由一组只通过输出依赖数据的方程定义。最大似然与 Infomax 求解 E^[φ(y)y⊤]=I;FastICA 的不动点满足关于 y 的白化条件和不动点条件;JADE 优化的是 y 的累积量。这样的估计量具有等变性(Cardoso & Laheld, 1996;另见Proposition 7.1命题 7.1),而等变性对性能有一个很强的推论。
Proposition 10.3 (performance does not depend on the mixing matrix)命题 10.3(性能与混合矩阵无关)
Suppose W^ solves equations H(W^X)=0 that depend on the data matrix X only through the outputs W^X, and that when several solutions exist, the one returned is chosen by a rule that also depends only on the outputs. Then the gain matrix P^=W^A obtained from X=AS equals the estimate obtained from the sources S themselves. In particular, the distributions of P^, ε and every ISR are the same for all invertible A.
Let WS be the solution returned for the data S. The matrix WSA−1 produces from X=AS exactly the same outputs, (WSA−1)(AS)=WSS, so it solves the equations for X and is the one returned. Hence W^(AS)=W^(S)A−1 and P^=W^(S), which does not involve A.
记 WS 为对数据 S 返回的解。矩阵 WSA−1 作用于 X=AS 时给出完全相同的输出,(WSA−1)(AS)=WSS,因此它是关于 X 的方程的解,而且就是被返回的那一个。于是 W^(AS)=W^(S)A−1,从而 P^=W^(S),与 A 无关。
Simulation studies can therefore take A=I without loss of generality, and an ill-conditioned mixing matrix does not by itself make an equivariant estimator less accurate. What can depend on A is the path of an iterative algorithm that is not itself equivariant, such as plain gradient ascent with a fixed step size, and the numerical round-off when A is nearly singular.
因此,模拟研究可以不失一般性地取 A=I;混合矩阵病态本身也不会降低等变估计量的精度。可能依赖于 A 的,是本身不具等变性的迭代算法的收敛路径(例如固定步长的普通梯度上升),以及 A 接近奇异时的数值舍入误差。
The ISR requires the outputs to be matched to the sources and the scales to be removed. Amari, Cichocki and Yang (1996) proposed a single number that needs no matching of outputs to sources. It does depend on the scales of the outputs, which are therefore normalized first (see the invariance property below). For an n×n gain matrix P=(pij),
Its main properties follow directly from the definition.
它的主要性质可以直接从定义得出。
Zero exactly at separation. Each bracket is a sum of ratios in [0,1], one of which equals 1, so it is nonnegative and vanishes exactly when its row (or column) has a single nonzero entry. Hence E(P)≥0, with equality if and only if P=ΠD as in Proposition 1.3命题 1.3.
Range. Each bracket is at most n−1, so 0≤E(P)≤2n(n−1); the normalized index E/(2n(n−1)) lies in [0,1].
Invariance.E is unchanged by permutations of the rows or of the columns. A row bracket is unchanged when its row is rescaled and a column bracket when its column is rescaled, but rescaling a row changes the column brackets. In practice E is computed after normalizing every output to unit variance, that is, every row of P^ to unit Euclidean norm (the sources having unit variance).
Relation to the ISR. With the rows normalized and the outputs matched to the sources, ∣P^ii∣=1+O(∥ε∥2) and E(P^)=2∑i=j∣εij∣+O(∥ε∥2). If Tεij is approximately normal, E∣εij∣≈2ISRij/π, so E[E(P^)]≈22/π∑i=jISRij. The Amari index decreases like 1/T, the ISR like 1/T: halving the index takes four times as many samples.
Example 10.4 (reading a gain matrix)例 10.4(解读一个增益矩阵)
Take the gain matrix printed in Section 1.5第 1.5 节, P^=(0.03−1.001.000.02), whose rows already have norm 1.00 to two decimals. The row brackets are 0.03 and 0.02, the column brackets 0.03 and 0.02, so E(P^)=0.10 and the normalized index is 0.025. Output 1 estimates s2 and output 2 estimates s1. Reordering the outputs accordingly gives ε12=0.02/(−1.00) (source 2 in the estimate of source 1) and ε21=0.03/1.00 (source 1 in the estimate of source 2). The realized interference-to-signal ratios are ε122=4×10−4 (−34.0 dB) and ε212=9×10−4 (−30.5 dB), and indeed 2(∣ε12∣+∣ε21∣)=0.10=E(P^).
10.4Estimating equations and asymptotic normality估计方程与渐近正态性
Most ICA estimators are solutions of estimating equationsT1∑t=1TF(x(t);θ)=0 whose population version E[F(x;θ)]=0 holds at the true parameter θ0; here θ collects the entries of W and θ0 is a separating matrix. The theory of such M-estimators (Section A.8第 A.8 节; van der Vaart, 1998) states that if θ0 is an isolated root, F is smooth enough and M=E[∂F/∂θ⊤] is invertible at θ0, then the root θ^ near θ0 is consistent and
The argument is a first-order expansion, 0=E^[F(x;θ^)]≈E^[F(x;θ0)]+M(θ^−θ0), followed by the central limit theorem for E^[F(x;θ0)]. The errors are therefore of order 1/T and the ISR of order 1/T. For ICA the expansion has a special structure: it splits into separate 2×2 problems, one for each pair of sources.
Proposition 10.5 (asymptotic variance of ICA estimating equations)命题 10.5(ICA 估计方程的渐近方差)
Let φ1,…,φn be differentiable nonlinearities, and let W^ solve
设 φ1,…,φn 是可微的非线性函数,W^ 满足
E^[φi(yi)yj]=0,i=j,(10.5)
together with n scale equations (for instance E^[yi2]=1) whose population solution is P=I up to signs. Let μi=E[siφi(si)], ρi=E[φi′(si)], βi=E[φi2(si)], and assume ρiρj=μiμj for all i=j. Then for the solution near separation, Tεij is asymptotically normal with mean zero and variance
以及 n 个尺度方程(例如 E^[yi2]=1),且这些方程的总体解在符号意义下为 P=I。记 μi=E[siφi(si)],ρi=E[φi′(si)],βi=E[φi2(si)],并假设对所有 i=j 有 ρiρj=μiμj。则对分离点附近的解,Tεij 渐近服从均值为零的正态分布,方差为
If the data are centered with the sample mean before the equations are solved, the same holds with βi replaced by Varφi(si)=βi−(E[φi(si)])2 for every i.
若先用样本均值把数据中心化、再求解这些方程,结论不变,只需对每个 i 把 βi 换成 Varφi(si)=βi−(E[φi(si)])2。
Proof.证明
Write P^ij=δij+eij with all eij=Op(T−1/2) (signs fixed). Since εij=eij/(1+eii)=eij+Op(T−1), it suffices to study eij. The outputs are yi=si+∑keiksk, and a Taylor expansion of φi gives
把增益矩阵写成 P^ij=δij+eij,其中所有 eij=Op(T−1/2)(符号已固定)。由于 εij=eij/(1+eii)=eij+Op(T−1),只需研究 eij。输出为 yi=si+∑keiksk,对 φi 作 Taylor 展开得
In the two sums the sample averages can be replaced by their expectations at a cost op(T−1/2). By independence and zero means, E[φi(si)sk] equals μi for k=i and 0 otherwise, and for j=i, E[φi′(si)sksj] equals ρi for k=j and 0 otherwise. With uij=E^[φi(si)sj], the equations for (i,j) and (j,i) become
The diagonal entries eii do not appear, so the scale equations play no role. Solving, eij=(μiuji−ρjuij)/(ρiρj−μiμj). The pair T(uij,uji) is an average of i.i.d. zero-mean terms; by the central limit theorem it is asymptotically normal with variances E[φi2(si)]E[sj2]=βi, βj and covariance E[φi(si)si]E[φj(sj)sj]=μiμj. The variance of ρjuij−μiuji is then (ρj2βi−2ρjμi2μj+μi2βj)/T. With centering, s is replaced by s−sˉ, where sˉ is the sample mean of the sources. To first order this only subtracts sˉjE[φi(si)] from uij, so φi(si) is replaced by φi(si)−E[φi(si)] in uij: the variance βi becomes Varφi(si), and the covariance μiμj is unchanged.
The pairs decouple: vij depends only on the laws of si and sj and on the two nonlinearities applied to them, and it does not change when φi is multiplied by a constant. With the normalization μi=1 of Chapter 7第 7 章, the denominator is (ρiρj−1)2, and ρiρj>1 is condition (iii) of Theorem 7.2定理 7.2 (ρi is the ki of that theorem): estimating equations whose solution is close to instability have a large variance. The most important case is the true score.
Corollary 10.6 (maximum likelihood with the true scores)推论 10.6(使用真实得分函数的最大似然)
If φi=ψi for every i (maximum likelihood with the true source densities), then vij=Jj/(JiJj−1) for all i=j.
若对每个 i 都取 φi=ψi(用真实源密度做最大似然),则对所有 i=j 有 vij=Jj/(JiJj−1)。
Proof.证明
By (10.1), μi=1 and ρi=βi=Ji. The numerator of (10.6) becomes Jj2Ji−2Jj+Jj=Jj(JiJj−1) and the denominator (JiJj−1)2. Since E[ψi(si)]=0, centering does not change this value.
Section 10.7第 10.7 节 shows that this value is the smallest possible: maximum likelihood with the true scores is asymptotically efficient. Pham and Garat (1997) studied the asymptotic performance of this family of quasi-maximum-likelihood estimators, in which the true scores are replaced by approximations.
Many algorithms whiten with the sample covariance and then rotate (symmetric FastICA, JADE). Their outputs are exactly decorrelated in the sample, E^[y^y^⊤]=I, and this constraint alone fixes half of the error.
Suppose the outputs satisfy E^[y^y^⊤]=I and, after their order and signs have been fixed, P^−I=Op(T−1/2). Then for i=j, εij+εji=−E^[sisj]+op(T−1/2), so T(εij+εji) converges in distribution to N(0,1). Consequently, if (Tεij,Tεji) converges in distribution to a limit with variances vij and vji, then vij+vji≥21; to first order, ISRij+ISRji≥1/(2T), whatever the source distributions and the rotation criterion.
We have E^[y^y^⊤]=P^C^sP^⊤ with C^s=E^[ss⊤]=I+Δ. For i=j, Δij is the average of the i.i.d. zero-mean terms si(t)sj(t), whose variance is E[si2sj2]=1, so TΔij converges in distribution to N(0,1); by the law of large numbers Δ→0 in probability. Write P^=(I+δ)(I+ε) with δ=D^−I diagonal; by hypothesis δ and ε are Op(T−1/2). Expanding (I+δ)(I+ε)(I+Δ)(I+ε)⊤(I+δ)=I gives 2δ+ε+ε⊤+Δ=op(T−1/2), since every other term is a product of a factor of order T−1/2 and a factor tending to zero. The off-diagonal part is the first claim. Finally, if (Tεij,Tεji) converges to (X,Y), then X+Y∼N(0,1), and 1=Var(X+Y)≤2Var(X)+2Var(Y) because Var(X+Y)+Var(X−Y)=2Var(X)+2Var(Y).
Proposition 3.6命题 3.6 states the same result for whitening with the sample covariance followed by a rotation. There the error ε=P^−I keeps its diagonal; dividing each row of P^ by its diagonal entry, as in (10.2), changes the off-diagonal entries only by Op(1/T), and centering changes C^s only by Op(1/T). The rotation step can only reduce the antisymmetric part εij−εji; the symmetric part is set by the sampling error of the covariance. Section 10.7第 10.7 节 shows that the best possible estimator does better than this whenever both sources are non-Gaussian.
10.5One-unit estimators and the optimal nonlinearity单单元估计量与最优非线性函数
One-unit FastICA (Algorithm 8.2算法 8.2) centers the data, whitens them with the sample covariance, z=V^x, and returns a unit vector w with E^[zg(w⊤z)]=λw for some scalar λ: the Lagrange condition for an extremum of E^[G(w⊤z)] on the sphere, which the update (8.3) solves. Because the sources have zero mean, we first analyse the estimator without the centering step and then state what centering changes. Suppose the solution estimates source i, and write y=w⊤z=q⊤s with q=(V^A)⊤w, the row of the gain matrix that produces y. The gain errors of this output are εk=qk/qi, k=i.
Theorem 10.8 (asymptotic variance of one-unit estimators)定理 10.8(单单元估计量的渐近方差)
Let g be differentiable, let the expectations below exist, and let E[sig(si)−g′(si)]=0. For the solution w of E^[zg(w⊤z)]=λw, ∥w∥=1, near the direction of source i, and for every k=i, Tεk is asymptotically normal with mean zero and variance
设 g 可微,下列期望存在,且 E[sig(si)−g′(si)]=0。对方程 E^[zg(w⊤z)]=λw、∥w∥=1 在源 i 方向附近的解 w,以及每个 k=i,Tεk 渐近服从均值为零的正态分布,方差为
The variance depends on the distribution of si alone: not on the other sources and, by Proposition 10.3命题 10.3, not on A (Hyvärinen, 1997; 1999). If the data are centered with the sample mean before whitening, as in Algorithm 8.2算法 8.2, the same holds with E[g2(si)] in the numerator replaced by Varg(si)=E[g2(si)]−(E[g(si)])2 (Wei, 2015).
Source coordinates. With M=V^A we have z=Ms and MC^sM⊤=I, where C^s=E^[ss⊤], because V^ whitens the sample. Since w=M−⊤q and M−1M−⊤=C^s, multiplying the fixed-point condition by M−1 gives E^[sg(y)]=λC^sq=λE^[sy], and the constraint becomes ∥w∥2=q⊤C^sq=E^[y2]=1.
Expansion. Take i=1 and write q=(1+δ)e1+∑k≥2ηkek with δ,ηk=Op(T−1/2), so that εk=ηk+Op(T−1). For k≥2, expand both sides of the k-th equation around y=s1 and replace the sample averages that multiply δ or ηl by their expectations; since E[skslg′(s1)]=δklE[g′(s1)] for k,l≥2 and E[sks1g′(s1)]=0,
Limit. The right side is an average of i.i.d. zero-mean terms with variance E[sk2]E[(g(s1)−μs1)2]=E[g2(s1)]−μ2, using E[s1g(s1)]=μ and E[s12]=1. The central limit theorem and division by μ−E[g′(s1)] give the result.
Centering. If the data are centered, the argument applies with s replaced by s−sˉ, where sˉ is the sample mean of the sources. To first order the only change is that E^[skg(s1)] gains the term −sˉkE[g(s1)]. Hence g(s1) is replaced by g(s1)−E[g(s1)] on the right side, and the variance becomes E[g2(s1)]−μ2−(E[g(s1)])2.
Three features of (10.7) matter in practice. The denominator is the quantity E[sg(s)−g′(s)] that decides whether the source direction is an extremum of the FastICA contrast (Chapter 8第 8 章); for a Gaussian source (10.1) with f=g gives E[sg(s)]=E[g′(s)], so the variance is infinite. Second, Vg does not change when g is multiplied by a nonzero constant or when a multiple of s is added to g: only the nonlinear part of g matters. Third, the numerator equals E[(g(s)−μs)2], the mean square of the part of g(s) that a linear regression on s does not explain. With centering the regression also contains a constant, so centering never increases Vg. The two versions agree when E[g(s)]=0, as for a symmetric source and an odd g; this is the case in all other examples of this chapter, including the simulations below, which center the data. For skewed sources they can differ appreciably: for the unit-variance exponential law shifted to zero mean (density e−(s+1) for s≥−1) and g=tanh, numerical quadrature gives Vg=3.75 without centering and 3.14 with it. Table 10.1表 10.1 lists Vg for the three standard nonlinearities of FastICA.
(10.7) 有三个在实践中要紧的特点。其一,分母中的 E[sg(s)−g′(s)] 正是决定源方向是否为 FastICA 对比函数极值点的那个量(Chapter 8第 8 章);对高斯源,在 (10.1) 中取 f=g 得 E[sg(s)]=E[g′(s)],所以方差为无穷大。其二,把 g 乘以非零常数,或给 g 加上 s 的倍数,Vg 都不变:起作用的只是 g 的非线性部分。其三,分子等于 E[(g(s)−μs)2],即 g(s) 中不能由 s 的线性回归解释的部分的均方。做中心化时,回归中还多了一个常数项,所以中心化绝不会增大 Vg。当 E[g(s)]=0 时(例如源分布对称而 g 为奇函数),两种情形的结果相同;本章其余的例子都属于这种情形,包括下文先做中心化的模拟。对偏斜的源,两者可能相差不小:对平移到零均值的单位方差指数分布(s≥−1 时密度为 e−(s+1))并取 g=tanh,数值积分给出不做中心化时 Vg=3.75,做中心化时 Vg=3.14。Table 10.1表 10.1列出了 FastICA 三个标准非线性函数的 Vg。
Table 10.1. Asymptotic variance Vg of one-unit FastICA for unit-variance sources, computed from (10.7) by numerical quadrature, and the lower bound 1/(J−1) of Proposition 10.9命题 10.9. The hyperbolic secant density is 21sech(πs/2); densities of the other laws are in Table 10.2表 10.2. The uniform density has bounded support, so J is not defined for it (see Table 10.2表 10.2).表 10.1 单位方差源下单单元 FastICA 的渐近方差 Vg(由 (10.7) 经数值积分算得),以及Proposition 10.9命题 10.9给出的下界 1/(J−1)。双曲正割分布的密度为 21sech(πs/2);其余分布的密度见Table 10.2表 10.2。均匀分布的支撑集有界,J 对它无定义(见Table 10.2表 10.2)。
Source law源的分布
J
g(u)=tanhu
g(u)=ue−u2/2
g(u)=u3
1/(J−1)
Laplace拉普拉斯
2
2.01
1.82
6.00
1.00
hyperbolic secant双曲正割
π2/8≈1.234
4.43
4.37
9.00
4.28
logistic逻辑斯谛
π2/9≈1.097
10.37
10.96
15.43
10.35
uniform均匀
—
0.69
0.71
0.43
—
No nonlinearity in the table is best for every law, and the cubic nonlinearity of kurtosis is the best for the uniform source but the worst for the three super-Gaussian ones. What is the best possible choice for a given source?
Proposition 10.9 (the score is the optimal nonlinearity)命题 10.9(得分函数是最优非线性函数)
Let s have unit variance and a density with score ψ and Fisher information 1<J<∞, and let g satisfy the conditions of Theorem 10.8定理 10.8 together with g(s)p(s)→0 as ∣s∣→∞. Then
设 s 方差为一,其密度的得分函数为 ψ,Fisher 信息满足 1<J<∞;又设 g 满足Theorem 10.8定理 10.8的条件,并且 ∣s∣→∞ 时 g(s)p(s)→0。则
Vg≥J−11,(10.8)
with equality if and only if g(s)=aψ(s)+bs for constants a=0 and b. In particular the score function itself is optimal, with Vψ=1/(J−1).
Let μ=E[sg(s)] and h(s)=g(s)−μs, so that the numerator of (10.7) is E[h2]. By (10.1), E[g′(s)]=E[ψ(s)g(s)], hence E[sg(s)−g′(s)]=E[(s−ψ(s))g(s)]=E[(s−ψ(s))h(s)], where the last step uses E[(s−ψ(s))s]=1−1=0. By the Cauchy–Schwarz inequality,
which is (10.8). Equality holds if and only if h is proportional to s−ψ, that is, g=μs+c(s−ψ) with c=0 (c=0 makes g linear and the denominator zero). For g=ψ: μ=1, the numerator is J−1 and the denominator (1−J)2. With centering, the numerator is E[h2] for h(s)=g(s)−μs−E[g(s)]; since E[s−ψ(s)]=0, the same argument gives the same bound, with equality if and only if g(s)=aψ(s)+bs+d for constants a=0, b and d, and Vψ is unchanged because E[ψ(s)]=0.
In terms of the contrast, g=G′=ψ means G=−logp up to scale and an added quadratic: the best one-unit contrast is minus the log-density of the source. The unit-variance hyperbolic secant law has score 2πtanh(πs/2), so tanh(au) with a=π/2 is optimal for it; with a=1 the variance is 4.43 instead of 4.28. The Laplace score is 2sign(s); tanh(au) approaches it as a grows, and quadrature gives Vg=2.01, 1.11 and 1.04 for a=1, 10 and 30, approaching the bound 1. For the kurtosis nonlinearity g(u)=u3, (10.7) becomes (E[s6]−(E[s4])2)/(E[s4]−3)2, which involves the sixth moment; for the Laplace law it is six times the optimum. This is the asymptotic form of the robustness argument of Section 5.6第 5.6 节. Since the density is unknown in practice, the optimum is approached either by robust general-purpose choices or by choosing for each component a nonlinearity adapted to its estimated distribution, as EFICA does (Section 10.7第 10.7 节).
The script below estimates T⋅ISR for one-unit FastICA on a Laplace source mixed with two uniform sources, for g=tanh and g(u)=u3, and compares it with Vg computed by Monte Carlo integration. Every run draws new sources and whitens with the sample covariance; the iteration starts at the true direction, so that the run measures the accuracy of the fixed point near source 1 rather than the chance of converging to another source.
where ± is two standard errors of the Monte Carlo average (the exact values are Vtanh=2.015 and Vu3=6). For tanh prediction and observation agree. For u3 the observed value is about 5% low, just outside the Monte Carlo interval: the asymptotic formula neglects higher-order terms, and these are larger for u3, whose fluctuations are driven by high moments of the source (the denominator involves the sample kurtosis, whose variance involves E[s8], which equals 2520 for the unit-variance Laplace law). Figure 10.1图 10.1 repeats the experiment over a range of T and for a uniform target. All six curves settle on their predicted constants; the cubic one for the Laplace source does so only for T near 104.
Figure 10.1.T⋅ISR of one-unit FastICA against the sample size, for a Laplace source mixed with a uniform one (a) and a uniform source mixed with a Laplace one (b), with three nonlinearities. Markers: averages over 1500 runs with ±2 standard errors; dashed lines: Vg from (10.7); dotted line in (a): the optimum 1/(J−1)=1 for the Laplace source.图 10.1 单单元 FastICA 的 T⋅ISR 随样本量的变化:(a)拉普拉斯源与一个均匀源混合,(b)均匀源与一个拉普拉斯源混合,各用三种非线性函数。标记:1500 次运行的平均值及 ±2 倍标准误;虚线:由 (10.7) 算出的 Vg;(a)中的点线:拉普拉斯源的最优值 1/(J−1)=1。
The variances computed so far belong to particular estimators. The Cramér–Rao bound limits all of them. Suppose the source densities are known, so that the unmixing matrix is the only unknown, and parametrize it around the truth as W=(I+E)A−1. The gain matrix is then P=I+E: the off-diagonal entries Eij are the separation errors and the diagonal entries are scale parameters. By the change of variables behind (6.1), one observation has log-density
Theorem 10.10 (Cramér–Rao bound for ICA)定理 10.10(ICA 的 Cramér–Rao 界)
Assume the standing assumptions, with source densities pi that are absolutely continuous and positive on R, have finite Fisher information Ji, and satisfy the usual regularity conditions of the Cramér–Rao inequality. At E=0 the Fisher information matrix of one observation is block diagonal: the diagonal (scale) parameters form one block, and each pair (Eij,Eji), i<j, forms the block
Since ∂log∣det(I+E)∣/∂Eij=[(I+E)−⊤]ij (Section A.6第 A.6 节) and ∂logpi(yi)/∂Eij=(pi′/pi)(yi)sj, the score for Eij at E=0 is Γij=δij−ψi(si)sj. For i=j and k=l, E[ΓijΓkl]=E[ψi(si)sjψk(sk)sl]. By independence this factorizes over the sources, and any factor of the form E[sl] or E[ψk(sk)] vanishes (the latter by (10.1) with f=1). The only nonzero cases are (k,l)=(i,j), giving E[ψi2(si)]E[sj2]=Ji, and (k,l)=(j,i), giving E[ψi(si)si]E[ψj(sj)sj]=1. In the same way E[ΓijΓkk]=−E[ψi(si)sj]+E[ψi(si)sjψk(sk)sk]=0 for i=j, because every term contains E[sj] or E[ψi(si)]. This proves the block structure. For T i.i.d. observations the information is T times larger, and the Cramér–Rao inequality bounds Var(E^ij) by the corresponding diagonal entry of the inverse, [(TFij)−1]11, where
Finally εij=E^ij/(1+E^ii) agrees with E^ij to first order.
最后,εij=E^ij/(1+E^ii) 在一阶近似下与 E^ij 相同。
Practical estimators are biased at finite T, so the bound is used in its asymptotic sense: under standard regularity conditions, no regular estimator that is consistent and asymptotically normal has an asymptotic variance vij below Jj/(JiJj−1) (van der Vaart, 1998). Comparing (10.10) with Corollary 10.6推论 10.6 shows that maximum likelihood with the true scores attains the bound.
实际的估计量在有限 T 下是有偏的,因此这个界按渐近意义使用:在标准的正则条件下,任何相合且渐近正态的正则估计量,其渐近方差 vij 都不会低于 Jj/(JiJj−1)(van der Vaart, 1998)。把 (10.10) 与Corollary 10.6推论 10.6比较可知,用真实得分函数的最大似然达到了这个界。
What the bound says这个界说明了什么
Two Gaussian sources. If Ji=Jj=1, the block (10.9) is singular and the bound is infinite: the data carry no information about the rotation within the pair, as Section 4.4第 4.4 节 shows from the identifiability side.
One Gaussian source. If Jj=1<Ji, then CRB(εij)=1/(T(Ji−1)) but CRB(εji)=Ji/(T(Ji−1))>1/T. The non-Gaussian source is located by its non-Gaussianity; the Gaussian output can only be found as the direction uncorrelated with it, and it inherits the sampling error of the covariance (Exercise 10.4习题 10.4).
Nearly Gaussian sources. For two sources with the same law the bound is J/(T(J2−1)). With T=1000 it is 6.7×10−4 (−31.8 dB) for Laplace sources (J=2) but 5.4×10−3 (−22.7 dB) for logistic sources (J≈1.097): identifiable is not the same as easy to estimate (Chapter 4第 4 章).
Strongly non-Gaussian sources. As Jj→∞, CRB(εij)→1/(TJi) and CRB(εji)→0.
Bounded sources. A uniform density jumps at the ends of its support, the regularity conditions fail, and the bound is not defined (Tichavský et al., 2006). Smooth densities that approach the uniform law, such as the generalized Gaussian with large exponent in Table 10.2表 10.2, have J→∞, which correctly suggests that such sources are easy to separate.
有界的源。均匀分布的密度在支撑集端点处跳跃,正则条件不成立,这个界没有定义(Tichavský et al., 2006)。逼近均匀分布的光滑密度(例如Table 10.2表 10.2中指数很大的广义高斯分布)的 J→∞,这正确地提示了这类源很容易分离。
Table 10.2. Fisher information J=E[ψ2(s)] of unit-variance laws (values checked by numerical quadrature). For the generalized Gaussian law, β=Γ(1/α)/Γ(3/α); J is infinite for α≤1/2 and grows like α/3 as α→∞. In this chapter ν denotes the degrees of freedom of the Student law, not the Gaussian variable of Table 1.1表 1.1.表 10.2 单位方差分布的 Fisher 信息 J=E[ψ2(s)](数值均经数值积分核对)。广义高斯分布中 β=Γ(1/α)/Γ(3/α);α≤1/2 时 J 为无穷大,α→∞ 时 J 像 α/3 一样增长。本章中 ν 表示 Student 分布的自由度,而不是Table 1.1表 1.1中的高斯变量。
Law分布
Density密度
J
Gaussian高斯
2π1e−s2/2
1
Laplace拉普拉斯
21e−2∣s∣
2
logistic逻辑斯谛
43πsech223πs
π2/9≈1.097
hyperbolic secant双曲正割
21sech2πs
π2/8≈1.234
Student tν, ν>2Student tν,ν>2
∝(1+ν−2s2)−(ν+1)/2
(ν−2)(ν+3)ν(ν+1); 1.25 for ν=5(ν−2)(ν+3)ν(ν+1);ν=5 时为 1.25
generalized Gaussian, α>1/2广义高斯,α>1/2
∝e−∣s/β∣α
Γ(1/α)2α2Γ(2−1/α)Γ(3/α); 1.37 for α=4Γ(1/α)2α2Γ(2−1/α)Γ(3/α);α=4 时为 1.37
uniform均匀
231 on [−3,3]在 [−3,3] 上为 231
not defined无定义
Which estimators attain the bound?哪些估计量能达到这个界?
Maximum likelihood with the true scores does (Corollary 10.6推论 10.6). The densities are unknown in practice, but each off-diagonal score ψi(si)sj is uncorrelated with every zero-mean function of a single source, so perturbing the marginal densities carries no first-order information about Eij; this is why the same bound remains the benchmark when the densities must be learned from the data (see Amari & Cardoso, 1997, for the semiparametric theory of ICA estimating functions). Other estimators fall short of the bound, for reasons that can be identified.
A one-unit contrast looks only at the distribution of its own output and ignores the non-Gaussianity of the interfering sources, which maximum likelihood exploits. Estimators with exactly decorrelated outputs pay a different price. By Proposition 10.7命题 10.7, T(εij+εji) is asymptotically N(0,1), so to first order the sum has variance 1/T, while the bound for the variance of the sum, (1,1)Fij−1(1,1)⊤/T, equals
since the difference of the two sides is proportional to (Ji−1)(Jj−1)≥0, with equality only when one of the sources is Gaussian. Hence no estimator with exactly decorrelated outputs is asymptotically efficient for a pair of non-Gaussian sources; this limitation of prewhitened (orthogonal) methods was analysed by Cardoso (1994). Symmetric FastICA (Algorithm 8.6算法 8.6) has both limitations, since it applies one nonlinearity g to all outputs and decorrelates them exactly. Tichavský, Koldovský and Oja (2006) derived its asymptotic accuracy (the result is restated as eq. (10) of Koldovský, Tichavský & Oja, 2006): for k=l,
with μk=E[skg(sk)], ρk=E[g′(sk)] and βk=E[g2(sk)]. For two Laplace sources and g=tanh this is 1.26, compared with 2.01 for the one-unit estimator and the bound 0.67. Exercise 10.5习题 10.5 derives (10.11) by linearizing the fixed-point conditions in the way of Proposition 10.5命题 10.5. Like (10.7), it assumes that the data are not centered with the sample mean, or that E[g(sk)]=E[g(sl)]=0; with centering, its numerator loses (E[g(sk)])2+(E[g(sl)])2 (Wei, 2015). EFICA (Koldovský, Tichavský & Oja, 2006) is designed to remove both limitations: it starts from symmetric FastICA, chooses for each component a nonlinearity adapted to its estimated distribution, and refines the estimate. Its authors proved that it is asymptotically efficient, attaining the Cramér–Rao bound, for generalized Gaussian sources with α>2, for which the nonlinearity it selects is the score function of the source, with the exponent estimated from the data. For super-Gaussian generalized Gaussian sources with 1≤α<2 its accuracy is close to the bound without attaining it, and for smaller α it is no longer close.
The gain error εij=P^ij/P^ii measures how much of source j leaks into the estimate of source i; ISRij=E[εij2]. Once the outputs are normalized, Amari's index summarizes P^ without matching outputs to sources; it decreases like 1/T.
For equivariant estimators the distribution of the errors does not depend on the mixing matrix, only on the sources and the method.
Estimators defined by estimating equations have errors of order 1/T; for ICA the asymptotic variance decouples into pairs of sources, and exact decorrelation makes T(εij+εji) asymptotically N(0,1).
One-unit FastICA has asymptotic variance Vg, which depends only on the estimated source (and, for skewed sources, on whether the data are centered); it is minimized by the score function, with value 1/(J−1).
The Cramér–Rao bound is Jj/(T(JiJj−1)): infinite for two Gaussian sources, large for nearly Gaussian ones, attained by maximum likelihood with the true scores, and in general not attained by one-unit or exactly decorrelating estimators.
Let P^=(0.1−1.02.00.05). (a) Compute the Amari index E(P^). (b) Multiply the first row by 10 and compute it again. (c) Normalize the rows to unit norm, match the outputs to the sources, and compute the gain errors, their squares in decibels, and 2∑i=j∣εij∣.
(a) Row brackets: (0.1+2.0)/2.0−1=0.05 and (1.0+0.05)/1.0−1=0.05. Column brackets: (0.1+1.0)/1.0−1=0.10 and (2.0+0.05)/2.0−1=0.025. Hence E=0.225. (b) The first row becomes (1,20); its bracket is still 0.05, but the columns (1,−1) and (20,0.05) give brackets 1 and 0.0025, so E=1.1025, although the separation is unchanged: the index is not invariant under rescaling of rows. (c) After normalization the rows are approximately (0.050,0.999) and (−0.999,0.050). Output 1 estimates s2 with ε21=0.1/2.0=0.05, and output 2 estimates s1 with ε12=0.05/(−1.0)=−0.05. Both realized ISRs are 2.5×10−3, that is, −26.0 dB. All four brackets now equal 0.05, so E=0.20=2(∣ε12∣+∣ε21∣).
Show that the unit-variance Laplace density 21e−2∣s∣ has Fisher information J=2. For two Laplace sources and T=1000, compute the Cramér–Rao bound for ε12 in decibels, and compare it with one-unit FastICA using the score and using g=tanh (Table 10.1表 10.1).
For s=0, logp(s)=−21log2−2∣s∣, so ψ(s)=2sign(s) and ψ2=2 almost everywhere; hence J=2. The bound is T1⋅4−12=3T2=6.7×10−4, or −31.8 dB. The one-unit estimator with the score has 1/(J−1)=1, giving 10−3 (−30.0 dB); with tanh, 2.015/T=2.0×10−3 (−27.0 dB). The gap between the first two numbers is the information in the non-Gaussianity of the interfering source, which a one-unit contrast ignores (Corollary 10.11推论 10.11).
Show that for g(u)=u3 the asymptotic variance (10.7) equals (E[s6]−(E[s4])2)/(E[s4]−3)2. Evaluate it for the unit-variance Laplace law (E[s4]=6, E[s6]=90) and uniform law (E[s4]=9/5, E[s6]=27/7). Then explain why a heavy tail can make the kurtosis estimator useless even when the kurtosis itself is large and finite.
With g(u)=u3: E[g2(s)]=E[s6], E[sg(s)]=E[s4] and E[g′(s)]=3E[s2]=3, which gives the formula. Laplace: (90−36)/(6−3)2=54/9=6. Uniform: the numerator is 27/7−81/25=108/175 and the denominator (9/5−3)2=36/25, so V=3/7≈0.43. The numerator is governed by the sixth moment, the denominator by the squared kurtosis. For a Student tν source (ν degrees of freedom) with 4<ν≤6 the kurtosis is finite but E[s6]=∞, so Vu3=∞: the kurtosis-based estimate does not even converge at the 1/T rate, whereas any bounded g, such as tanh, keeps Vg finite.
Let s1 be standard Gaussian and s2 unit-variance Laplace. Compute CRB(ε12) and CRB(ε21). Then show that maximum likelihood with the true scores attains both, by writing ε12 as the sum of two uncorrelated terms.
Here J1=1, J2=2 and J1J2−1=1, so CRB(ε12)=2/T and CRB(ε21)=1/T. In the proof of Proposition 10.5命题 10.5 take φ1(u)=u (the Gaussian score) and φ2=ψ2, so μ1=ρ1=β1=1, μ2=1 and ρ2=β2=2. Then e21=(μ2u12−ρ1u21)/(ρ1ρ2−μ1μ2)=u12−u21, with TVar(e21)=1−2+2=1, which attains the bound. The estimating equation E^[φ1(y1)y2]=E^[y1y2]=0 decorrelates the two outputs, and the computation of Proposition 10.7命题 10.7 gives e12=−u12−e21, with u12=E^[s1s2]. The two terms are uncorrelated: TCov(u12,e21)=TVar(u12)−TCov(u12,u21)=1−μ1μ2=0. Hence TVar(e12)=1+1=2. The Laplace output is located by its non-Gaussianity; the Gaussian output is found only as the direction uncorrelated with it, and its error adds the sampling error of the covariance to the error of the Laplace direction.
Symmetric FastICA (Algorithm 8.6算法 8.6) updates all rows of an orthogonal W by (8.3) and then orthogonalizes symmetrically, W←(W+W+⊤)−1/2W+. (a) Show that at a fixed point (up to sign changes of rows) the outputs satisfy E^[yy⊤]=I and E^[gk(yk)yl]=E^[gl(yl)yk] for k=l, where gk=dkg and dk=±1 is the sign of μk−ρk. (b) Derive (10.11).
(a) Whitened data and an orthogonal W give E^[yy⊤]=WW⊤=I. Let D=diag(dk). The iteration returns DW exactly when (W+W+⊤)−1/2W+=DW, that is, when W+W⊤D=(W+W+⊤)1/2 is symmetric positive definite. Since y=Wz and WW⊤=I, W+W⊤=E^[g(y)y⊤]−diag(E^[g′(y)])=:M. Symmetry of MD means Mkldl=Mlkdk; multiplying by dkdl gives dkE^[g(yk)yl]=dlE^[g(yl)yk]. Positivity of the diagonal, dk(E^[g(yk)yk]−E^[g′(yk)])>0, forces dk=sign(μk−ρk) for large T.
(b) Let μk′=dkμk and ρk′=dkρk be the constants of gk, so that τk=μk′−ρk′=∣μk−ρk∣, and let ukl=E^[gk(sk)sl] and Δkl=E^[sksl]. As in the proof of Proposition 10.5命题 10.5, E^[gk(yk)yl]=ukl+μk′elk+ρk′ekl+op(T−1/2), so the symmetry condition reads ukl+μk′elk+ρk′ekl=ulk+μl′ekl+ρl′elk. By Proposition 10.7命题 10.7, elk=−Δkl−ekl. Substituting and collecting terms,
The two averages are uncorrelated, since T times their covariance is E[sk(gk(sk)−μk′sk)]E[sl(gl(sl)−ρl′sl)]=0⋅(μl′−ρl′)=0. Their variances, times T, are E[(gk(sk)−μk′sk)2]=βk−μk2=γk and E[(gl(sl)−ρl′sl)2]=βl−2ρl′μl′+ρl′2=γl+τl2. Therefore TVar(ekl)→(γk+γl+τl2)/(τk+τl)2, and εkl=ekl+Op(T−1). If the data are first centered with the sample mean, s is replaced by s−sˉ. To first order this subtracts sˉlE[gk(sk)] from ukl and sˉkE[gl(sl)] from ulk, so in the two averages gk(sk) and gl(sl) are replaced by their deviations from their means. The averages remain uncorrelated, γk becomes γk−(E[g(sk)])2, γl+τl2 becomes γl+τl2−(E[g(sl)])2, and the numerator loses (E[g(sk)])2+(E[g(sl)])2.
这两个平均不相关,因为它们的协方差乘以 T 等于 E[sk(gk(sk)−μk′sk)]E[sl(gl(sl)−ρl′sl)]=0⋅(μl′−ρl′)=0。它们的方差乘以 T 分别为 E[(gk(sk)−μk′sk)2]=βk−μk2=γk 和 E[(gl(sl)−ρl′sl)2]=βl−2ρl′μl′+ρl′2=γl+τl2。因此 TVar(ekl)→(γk+γl+τl2)/(τk+τl)2,而 εkl=ekl+Op(T−1)。若先用样本均值把数据中心化,就要把 s 换成 s−sˉ。在一阶近似下,这使 ukl 减去 sˉlE[gk(sk)],ulk 减去 sˉkE[gl(sl)],相当于在两个平均中把 gk(sk) 和 gl(sl) 换成它们与各自均值之差。两个平均仍不相关,γk 变为 γk−(E[g(sk)])2,γl+τl2 变为 γl+τl2−(E[g(sl)])2,所以分子减少 (E[g(sk)])2+(E[g(sl)])2。
Amari, S., & Cardoso, J.-F. (1997). Blind source separation—semiparametric statistical approach. IEEE Transactions on Signal Processing, 45(11), 2692–2700.
Amari, S., Cichocki, A., & Yang, H. H. (1996). A new learning algorithm for blind signal separation. In Advances in Neural Information Processing Systems 8 (pp. 757–763). MIT Press.
Cardoso, J.-F. (1994). On the performance of orthogonal source separation algorithms. In Proceedings of EUSIPCO-94 (pp. 776–779). Edinburgh.
Cardoso, J.-F., & Laheld, B. H. (1996). Equivariant adaptive source separation. IEEE Transactions on Signal Processing, 44(12), 3017–3030.
Hyvärinen, A. (1997). One-unit contrast functions for independent component analysis: A statistical analysis. In Neural Networks for Signal Processing VII: Proceedings of the 1997 IEEE Workshop (pp. 388–397). IEEE.
Hyvärinen, A. (1999). Fast and robust fixed-point algorithms for independent component analysis. IEEE Transactions on Neural Networks, 10(3), 626–634.
Koldovský, Z., Tichavský, P., & Oja, E. (2006). Efficient variant of algorithm FastICA for independent component analysis attaining the Cramér–Rao lower bound. IEEE Transactions on Neural Networks, 17(5), 1265–1277.
Pham, D.-T., & Garat, P. (1997). Blind separation of mixture of independent sources through a quasi-maximum likelihood approach. IEEE Transactions on Signal Processing, 45(7), 1712–1725.
Tichavský, P., Koldovský, Z., & Oja, E. (2006). Performance analysis of the FastICA algorithm and Cramér–Rao bounds for linear independent component analysis. IEEE Transactions on Signal Processing, 54(4), 1189–1203.
van der Vaart, A. W. (1998). Asymptotic Statistics. Cambridge University Press.
Wei, T. (2015). An overview of the asymptotic performance of the family of the FastICA algorithms. In Latent Variable Analysis and Signal Separation (LVA/ICA 2015) (Lecture Notes in Computer Science, Vol. 9237, pp. 336–343). Springer.