The fixed-point algorithm as an approximate Newton method: one-unit, deflation and symmetric versions, convergence rates, and the choice of nonlinearity.
The gradient algorithms of Chapter 7第 7 章 need a step size, and near a solution they converge only linearly, at a rate set by the smallest eigenvalue of their linearization. FastICA (Hyvärinen & Oja, 1997; Hyvärinen, 1999) takes another route. It works on whitened data, estimates one direction at a time on the unit sphere, and replaces the gradient step by a fixed-point iteration obtained from Newton's method. The result has no step size, converges quadratically or even cubically, and usually needs fewer than ten iterations. This chapter derives the iteration, proves its convergence in source coordinates, builds the versions that estimate several components, and discusses the choice of nonlinearity.
Chapter 7第 7 章中的梯度算法需要选择步长，而且在解附近只是线性收敛，收敛速度由线性化的最小特征值决定。FastICA（Hyvärinen & Oja, 1997; Hyvärinen, 1999）走的是另一条路：它在白化数据上工作，在单位球面上一次估计一个方向，并用由牛顿法导出的不动点迭代代替梯度步。这样得到的算法没有步长，二次甚至三次收敛，通常不到十次迭代就够了。本章推导这一迭代，在源坐标下证明其收敛性，构造同时估计多个成分的版本，并讨论非线性函数的选择。
The setting is that of Chapter 3第 3 章: the data are whitened, z=Vx=A~s with A~ orthogonal (Proposition 3.4命题 3.4), and the sources have zero mean and unit variance. A one-unit method looks for a unit vector w such that y=w⊤z is one of the sources. In source coordinates, q=A~⊤w is again a unit vector and y=q⊤s; the direction w recovers source i exactly when q=±ei, where ei is the i-th standard basis vector.
where ν is a standard Gaussian variable and G is a smooth function. Because z is white, every y=w⊤z with ∥w∥=1 has unit variance, so JG compares y with a Gaussian of the same variance. The algorithm uses the derivative g=G′ and its derivative g′. Table 8.1表 8.1 lists the three standard choices. For G(u)=u4/4, EG(y)−EG(ν)=kurt(y)/4, and JG is the squared kurtosis of Chapter 5第 5 章.
Table 8.1. Standard contrast functions of FastICA (1≤a≤2 is a constant, often a=1) and the quantity λ=E[sg(s)−g′(s)] that governs their behaviour, for unit-variance uniform and Laplace sources (computed with a=1 by numerical quadrature). For a Gaussian source λ=0 for every g.表 8.1 FastICA 的标准对比函数（1≤a≤2 为常数，常取 a=1），以及决定其行为的量 λ=E[sg(s)−g′(s)] 在单位方差均匀源和拉普拉斯源下的取值（取 a=1，用数值积分计算）。对高斯源，任何 g 都有 λ=0。
G(u)
g(u)
g′(u)
λ, uniformλ，均匀
λ, Laplaceλ，拉普拉斯
a1logcosh(au)
tanh(au)
a(1−tanh2(au))
0.126
−0.148
−e−u2/2
ue−u2/2
(1−u2)e−u2/2
0.217
−0.211
u4/4
u3
3u2
−1.2
3
The source directions are critical points of EG(w⊤z) on the sphere, and the sign of λi=E[sig(si)−g′(si)] decides their type.
Proposition 8.1 (critical points of the one-unit contrast)命题 8.1（单单元对比函数的临界点）
Let G be twice continuously differentiable, and let the sources be independent with zero mean, unit variance and finite moments of the orders needed below. On the unit sphere, f(q)=EG(q⊤s) has a critical point at every ±ei. It is a strict local maximum if λi>0 and a strict local minimum if λi<0 (Hyvärinen, 1999).
设 G 二阶连续可微，各源相互独立、均值为零、方差为一，且下面用到的各阶矩有限。在单位球面上，f(q)=EG(q⊤s) 在每个 ±ei 处取临界点。若 λi>0，它是严格局部极大值点；若 λi<0，它是严格局部极小值点（Hyvärinen, 1999）。
Proof.证明
The Euclidean gradient of f is E[sg(q⊤s)]. At q=ei its j-th component is E[sj]E[g(si)]=0 for j=i, so the gradient is parallel to ei and ei is a critical point on the sphere (the Lagrange condition). For the second order, take i=1 and move along the sphere: q=(1−∥δ∥2,δ) with δ∈Rn−1 small, so y=s1+h with h=(1−∥δ∥2−1)s1+∑j≥2δjsj. Expanding G(s1+h)=G(s1)+g(s1)h+21g′(s1)h2+o(h2) and using independence, zero means and unit variances,
The term with E[s1g(s1)] comes from E[g(s1)h], using 1−∥δ∥2−1=−21∥δ∥2+O(∥δ∥4) and E[g(s1)sj]=0 for j≥2; the term with E[g′(s1)] comes from 21E[g′(s1)h2], using E[g′(s1)(∑jδjsj)2]=E[g′(s1)]∥δ∥2, all other terms being of higher order. The quadratic form −21λ1∥δ∥2 is negative definite if λ1>0 and positive definite if λ1<0.
For the kurtosis contrast λi=κi, consistent with Theorem 5.3定理 5.3: Ey4 is locally maximal at super-Gaussian and minimal at sub-Gaussian sources. For logcosh and the Gaussian-shaped G the signs are reversed for the laws of Table 8.1表 8.1. Squaring in (8.1) removes the need to know the sign in advance: a critical point ei of f is a local maximum of JG when (EG(si)−EG(ν))λi>0. This is Theorem 1 of Hyvärinen (1999), stated there for even G. For kurtosis the product equals κi2/4, positive whenever κi=0; for the other two G of the table it is positive for the uniform, Laplace and logistic laws (we checked it by quadrature). As the next sections show, the FastICA iteration itself does not care about the sign at all.
8.3Derivation as an approximate Newton method作为近似牛顿法的推导
The optima of EG(w⊤z) under ∥w∥2=1 satisfy the Lagrange condition
EG(w⊤z) 在约束 ∥w∥2=1 下的最优点满足 Lagrange 条件
F(w)=E[zg(w⊤z)]−βw=0,β=E[w⋆⊤zg(w⋆⊤z)],(8.2)
where w⋆ is the optimum and the value of the multiplier β follows by multiplying the condition by w⋆⊤. Newton's method for (8.2) needs the Jacobian E[zz⊤g′(w⊤z)]−βI. Its first term is the expensive part, and the whiteness of z suggests the approximation
which ignores the dependence between zz⊤ and g′(w⊤z). At a source direction it is nevertheless exact in every direction orthogonal to w and wrong only along w itself, which the normalization takes care of (Exercise 8.5习题 8.5). With a diagonal Jacobian the Newton step is
这一近似忽略了 zz⊤ 与 g′(w⊤z) 之间的相依性。尽管如此，在源方向上，它在与 w 正交的所有方向上都是精确的，只在 w 本身的方向上有误差，而这个方向由归一化处理（Exercise 8.5习题 8.5）。雅可比矩阵变成对角阵后，牛顿步为
w+=w−E[g′(w⊤z)]−βE[zg(w⊤z)]−βw.
Multiplying by the scalar β−E[g′(w⊤z)] changes only the length and possibly the sign of w+, both of which the subsequent normalization ignores up to sign. The β terms then cancel and we obtain the FastICA update
There is no step size and no β to estimate. Because a factor of either sign was dropped, w may flip sign from one iteration to the next; convergence is therefore judged by ∣w+⊤w∣→1 for the normalized iterates, not by w+−w→0.
In source coordinates the population iteration reads q+=E[sg(q⊤s)]−E[g′(q⊤s)]q, followed by normalization, because A~⊤z=s and w⊤z=q⊤s. For the kurtosis nonlinearity g(u)=u3 it can be solved in closed form. Expanding (q⊤s)3, the moment E[sisjsksl] of independent, zero-mean, unit-variance sources vanishes unless the indices pair up, so
using ∥q∥=1. Since E[g′(q⊤s)]=3E[(q⊤s)2]=3, the update is simply qi+=κiqi3: each coordinate is cubed and weighted by the kurtosis of its source. The identity holds exactly, even for skewed sources; we checked it to machine precision with discrete sources whose expectations can be computed exactly.
(a) Kurtosis (Hyvärinen & Oja, 1997). For g(u)=u3 and sources with finite fourth moments, qi+=κiqi3 for every i. If the index i⋆ that maximizes ∣κi∣∣qi(0)∣ is unique, the iterates converge to ±ei⋆ up to sign, and for every j=i⋆ the ratios ρj(k)=∣qj(k)/qi⋆(k)∣ satisfy ρj(k+1)=(∣κj∣/∣κi⋆∣)(ρj(k))3: the convergence is cubic.
（a）峭度（Hyvärinen & Oja, 1997）。对 g(u)=u3 且源的四阶矩有限，对每个 i 都有 qi+=κiqi3。若使 ∣κi∣∣qi(0)∣ 最大的下标 i⋆ 唯一，则迭代在相差符号的意义下收敛到 ±ei⋆，并且对每个 j=i⋆，比值 ρj(k)=∣qj(k)/qi⋆(k)∣ 满足 ρj(k+1)=(∣κj∣/∣κi⋆∣)(ρj(k))3：收敛是三次的。
(b) General g (convergence rates: Hyvärinen, 1999; the constant cij is computed in the proof below). Let g be three times differentiable with bounded third derivative, let the sources have finite fourth moments, and let λi=E[sig(si)−g′(si)]=0. Then ei is a fixed point up to sign, and near it the ratios rj=qj/qi, j=i, satisfy
The convergence is therefore at least quadratic, and cubic if E[g′′(si)]E[sj3]=0 for all j=i, for example when g is odd and si is symmetric, or when all other sources are symmetric.
因此收敛至少是二次的；若对所有 j=i 都有 E[g′′(si)]E[sj3]=0，则是三次的，例如 g 为奇函数且 si 对称，或者其余各源都对称。
Proof.证明
(a) The identity was derived above. Put ui=∣κi∣qi. The unnormalized update gives ∣κi∣qi+=sign(κi)ui3, and normalization multiplies all coordinates by the same positive number, so ∣uj(k)/ui⋆(k)∣=∣uj(0)/ui⋆(0)∣3k. By uniqueness of i⋆ the base is below 1 (coordinates with κj=0 vanish after one step), hence all ratios tend to zero and q(k) approaches ±ei⋆; in particular κi⋆=0 and qi⋆(0)=0, so q+ never vanishes. The recursion for ρj is ∣κjqj3∣/∣κi⋆qi⋆3∣.
(b) At q=ei the update is E[sg(si)]−E[g′(si)]ei=λiei, because E[sjg(si)]=0 for j=i; so ei is fixed up to the sign of λi. Near it write q=qiei+δ with δ⊥ei and qi=1−∥δ∥2, and y=si+h with h=(qi−1)si+∑k=iδksk. Taylor's theorem with the bounded third derivative gives, for j=i,
where independence and zero means kill every other term (the factor qi−1 is of order ∥δ∥2). Similarly E[g′(y)]=E[g′(si)]+O(∥δ∥2), because E[g′′(si)h]=(qi−1)E[sig′′(si)]. The j-th coordinate of the unnormalized update is therefore uj=E[sjg(y)]−E[g′(y)]δj=21E[g′′(si)]E[sj3]δj2+O(∥δ∥3): the first-order terms cancel. The i-th coordinate is ui=λi+O(∥δ∥). Normalization does not change ratios, so rj+=uj/ui, and since δj=rj(1+O(∥r∥2)), (8.4) follows.
Three consequences are worth stating. First, the iteration converges to source directions whether they are maxima or minima of EG: only λi=0 matters, so FastICA finds sub- and super-Gaussian sources with the same nonlinearity. When λi<0 the vector flips sign at every step, which is why convergence is tested through ∣w+⊤w∣. Second, part (a) is a global statement: for a random start, ties in ∣κi∣∣qi(0)∣ have probability zero, so the kurtosis iteration converges from almost every starting point. Part (b) is local. Third, a Gaussian source has λi=0 (Stein's lemma) and is never a fixed point; with kurtosis its coordinate is set to zero in one step (Exercise 8.2习题 8.2).
Example 8.4 (the kurtosis iteration in two dimensions)例 8.4（二维情形下的峭度迭代）
Take a uniform source (κ1=−1.2) and a Laplace source (κ2=3), and start from q(0)=(0.8,0.6), which is closer to the first source. Since 1.2×0.8≈0.876<3×0.6≈1.039, part (a) predicts convergence to the second source, with ρ(k+1)=0.4(ρ(k))3 for ρ=∣q1/q2∣. The iterates are q(1)=(−0.688,0.726), q(2)=(0.323,0.946), q(3)=(−0.0159,0.99987) and q(4)=(1.6×10−6,1.0), and ρ runs through 1.33,0.948,0.341,0.0159,1.6×10−6,1.6×10−18. The number of correct digits roughly triples at each step once the iterate is close, and the first coordinate alternates in sign because κ1<0.
Figure 8.1. (a) Population iterations with exact expectations: (i) the kurtosis example above; (ii) g=tanh with a uniform and a Laplace source, cubic since the sources are symmetric; (iii) g=tanh with two centered exponential sources, quadratic, with ∣r(k+1)∣/(r(k))2→1.806=∣c12∣, where c12=−1.806 is the constant of (8.4). (b) Symmetric FastICA on the four-source example of Section 8.5第 8.5 节 (T=10000) with three nonlinearities; the curves show 1−mini∣wi(k)⊤wi(k−1)∣.图 8.1 （a）用精确期望计算的总体迭代：（i）上面的峭度例子；（ii）g=tanh，一个均匀源和一个拉普拉斯源，由于源对称而三次收敛；（iii）g=tanh，两个中心化的指数分布源，二次收敛，且 ∣r(k+1)∣/(r(k))2→1.806=∣c12∣，其中 c12=−1.806 是 (8.4) 中的常数。（b）在Section 8.5第 8.5 节的四源例子（T=10000）上用三种非线性函数运行对称 FastICA；曲线为 1−mini∣wi(k)⊤wi(k−1)∣。
With a finite sample the iteration uses E^ and converges to a fixed point of the sample map, which lies within O(T−1/2) of the population fixed point; its accuracy is the subject of Theorem 10.8定理 10.8. At that point the derivative of the sample map is no longer exactly zero but of the order of the sampling error. In our experiments with g=tanh and T from 103 to 106 its norm was of the order of 1/T, typically a few times 1/T, with a factor that depended on the sources and on the sample. The final iterations therefore converge linearly, but with a ratio so small that this is rarely visible. In the run of Figure 8.1图 8.1(b), with T=104, the measure drops by a factor of about 600 per iteration at the end, so the angle between successive iterates shrinks by a factor of about 25, a contraction of about 4/T.
Running the one-unit algorithm n times from random starts would often find the same source several times. Since the columns of A~ are orthonormal, the rows wi⊤ of the orthogonal unmixing matrix U must be orthonormal too, and FastICA enforces this in one of two ways. In this section, as in the FastICA literature, W denotes the matrix with rows wi⊤, which acts on whitened data and is the orthogonal matrix U of Table 1.1表 1.1; the unmixing matrix for x is WV.
从随机起点把单单元算法独立运行 n 次，往往会多次找到同一个源。由于 A~ 的列是标准正交的，正交解混矩阵 U 的各行 wi⊤ 也必须标准正交，FastICA 用以下两种方式之一来保证这一点。本节沿用 FastICA 文献的习惯，用 W 表示以 wi⊤ 为行、作用于白化数据的矩阵，即Table 1.1表 1.1中的正交矩阵 U；对 x 的解混矩阵是 WV。
Deflation. Estimate the vectors one after another. While iterating for wp, follow each update by a Gram–Schmidt step against the vectors already found, wp←wp−∑j<p(wp⊤wj)wj, and then normalize. If w1,…,wp−1 were exact source directions, the remaining sources would span exactly the orthogonal complement, and the iteration there is again a one-unit problem with fewer sources, so Theorem 8.3定理 8.3 applies. With estimated vectors, however, each error tilts the subspace in which all later vectors must lie, so errors accumulate along the sequence; the last vector has no freedom left at all, being the unit vector orthogonal to the others. The order in which sources come out depends on the starting points and the contrast, not on any property chosen by the user.
Let W be invertible with singular value decomposition W=LΣR⊤. Then (WW⊤)−1/2W=LR⊤. This matrix is orthogonal, it is the orthogonal matrix closest to W in Frobenius norm (the polar factor of W; Section A.4第 A.4 节, Proposition A.4命题 A.4), and it treats all rows alike: for every permutation matrix Π, the result for ΠW is Π times the result for W.
设 W 可逆，其奇异值分解为 W=LΣR⊤。则 (WW⊤)−1/2W=LR⊤。这个矩阵是正交的，是在 Frobenius 范数下离 W 最近的正交矩阵（即 W 的极分解因子；见Section A.4第 A.4 节与Proposition A.4命题 A.4），并且对各行一视同仁：对任意置换矩阵 Π，ΠW 的结果等于 Π 乘以 W 的结果。
Proof.证明
WW⊤=LΣ2L⊤, so (WW⊤)−1/2=LΣ−1L⊤ and (WW⊤)−1/2W=LΣ−1L⊤LΣR⊤=LR⊤, a product of orthogonal matrices. For the last claim, (ΠWW⊤Π⊤)−1/2=Π(WW⊤)−1/2Π⊤, since ΠM−1/2Π⊤ is a symmetric positive definite square root of (ΠMΠ⊤)−1; multiplying by ΠW gives Π(WW⊤)−1/2W.
Gram–Schmidt, by contrast, depends on the order of the rows (Exercise 8.3习题 8.3). The symmetric version estimates all components at once, so errors are not passed along a sequence, and the cost per iteration is that of n one-unit updates plus one n×n eigendecomposition. If an inverse square root is undesirable, the iteration W←23W−21WW⊤W (Hyvärinen, 1999), started after dividing W by ∥W∥F, converges quadratically to the same matrix (Exercise 8.4习题 8.4).
与之相比，Gram–Schmidt 正交化的结果依赖于各行的顺序（Exercise 8.3习题 8.3）。对称版本同时估计所有成分，误差不会沿着一个序列向后传递，每次迭代的代价是 n 次单单元更新再加一次 n×n 的特征分解。如果不想计算逆平方根，可以先把 W 除以 ∥W∥F，再迭代 W←23W−21WW⊤W（Hyvärinen, 1999），它二次收敛到同一个矩阵（Exercise 8.4习题 8.4）。
The script below is a complete implementation, applied to four sources of different types: uniform, Laplace, binary ±1 and Student-t with 5 degrees of freedom, all with unit variance.
import numpy as np
def whiten(X):
"""Center the rows of X; return Z = V X (identity sample covariance) and V."""
X = X - X.mean(axis=1, keepdims=True)
d, E = np.linalg.eigh(np.cov(X))
V = E @ np.diag(d ** -0.5) @ E.T
return V @ X, V
def sym_orth(W):
"""Symmetric orthogonalization W <- (W W^T)^{-1/2} W."""
d, E = np.linalg.eigh(W @ W.T)
return E @ np.diag(d ** -0.5) @ E.T @ W
def fastica_symmetric(X, g, dg, max_iter=100, tol=1e-12, seed=0):
Z, V = whiten(X)
n, T = Z.shape
W = sym_orth(np.random.default_rng(seed).standard_normal((n, n)))
for it in range(1, max_iter + 1):
Y = W @ Z
# row i: E[z g(y_i)] - E[g'(y_i)] w_i; then orthogonalize all rows together
W_new = sym_orth(g(Y) @ Z.T / T - dg(Y).mean(axis=1)[:, None] * W)
change = 1 - np.min(np.abs(np.sum(W_new * W, axis=1)))
W = W_new
print(f"iteration {it}: 1 - min|<w_new, w_old>| = {change:.1e}")
if change < tol:
break
return W @ V # unmixing matrix for the raw data
rng = np.random.default_rng(1)
T = 10000
S = np.vstack([rng.uniform(-np.sqrt(3), np.sqrt(3), T), # sub-Gaussian
rng.laplace(0, 1 / np.sqrt(2), T), # super-Gaussian
np.sign(rng.standard_normal(T)), # binary +-1
rng.standard_t(5, T) / np.sqrt(5 / 3)]) # heavy-tailed
A = rng.standard_normal((4, 4))
W = fastica_symmetric(A @ S, np.tanh, lambda y: 1 - np.tanh(y) ** 2)
P = W @ A
print(np.round(P / np.abs(P).max(axis=1, keepdims=True), 2) + 0.0)
Seven iterations suffice, and the global matrix, with each row scaled by its largest entry, is a signed permutation up to off-diagonal entries of about 0.01–0.02, the estimation error for T=10000. The same nonlinearity tanh has found the sub-Gaussian uniform and binary sources as well as the two super-Gaussian ones. Figure 8.1图 8.1(b) shows that the other two nonlinearities converge just as fast on these data. The convergence of symmetric FastICA, in which the units are coupled through the orthogonalization, has been analysed by Oja and Yuan (2006). In a population check with exact expectations and the kurtosis nonlinearity, the distance to the nearest signed permutation fell from 0.10 to 2.6×10−4 to 4.3×10−12, again a cubic pattern.
Robustness. Each observation enters the update through zg(w⊤z). For the kurtosis contrast g(u)=u3 is unbounded: a single outlier at y=10 weighs as much as a thousand typical samples at y=1, and the estimate inherits the fragility of the sample kurtosis discussed in Section 5.6第 5.6 节. For logcosh, ∣g∣≤1; for the Gaussian-shaped G, ∣g(u)∣≤e−1/2 and g(u)→0 as ∣u∣→∞, so extreme values are effectively ignored.
Accuracy. The asymptotic variance of the one-unit estimator (Theorem 10.8定理 10.8) has λi2=(E[sig(si)−g′(si)])2 in its denominator. A source with λi close to zero for the chosen g is estimated poorly and, with finite data, found slowly or not at all; a non-Gaussian source can even have λi=0 for one particular g (a source with zero kurtosis is invisible to g(u)=u3), and then another g must be used. The variance is smallest when g is an affine function of the score of the source (Proposition 10.9命题 10.9). Since one fixed g cannot be optimal for sources of different types, the EFICA algorithm of Koldovský, Tichavský and Oja (2006) adapts the nonlinearity to each component after a symmetric FastICA run; Chapter 10第 10 章 returns to it.
精度。单单元估计量的渐近方差（Theorem 10.8定理 10.8）的分母中含有 λi2=(E[sig(si)−g′(si)])2。对所选的 g，若某个源的 λi 接近零，它的估计就很差，在有限数据下找到它也很慢，甚至找不到；一个非高斯源对某个特定的 g 甚至可能恰好有 λi=0（峭度为零的源对 g(u)=u3 不可见），这时必须换用别的 g。当 g 是源得分函数的仿射函数时，方差最小（Proposition 10.9命题 10.9）。由于一个固定的 g 不可能对不同类型的源都最优，Koldovský、Tichavský 与 Oja（2006）的 EFICA 算法在对称 FastICA 运行之后为每个成分自适应地选择非线性函数；Chapter 10第 10 章会再回到这一点。
Speed. Kurtosis gives cubic convergence for any sources, the other two nonlinearities give cubic convergence for symmetric sources and quadratic convergence otherwise. In practice the difference amounts to at most an iteration or two (Figure 8.1图 8.1). Each iteration costs one pass over the data, of order nT operations per unit; evaluating tanh or an exponential is slower than a cube, but rarely enough to matter.
速度。峭度对任何源都给出三次收敛；另两种非线性函数对对称源三次收敛，否则二次收敛。实践中两者的差别至多一两次迭代（Figure 8.1图 8.1）。每次迭代需要遍历一次数据，每个单元的计算量为 nT 量级；计算 tanh 或指数函数比计算立方慢，但很少慢到值得在意的程度。
These considerations lead to the usual practice, recommended by Hyvärinen (1999): logcosh with a=1 as the default; the Gaussian-shaped G when the sources are strongly super-Gaussian or the data contain outliers; kurtosis mainly for sub-Gaussian sources in clean data, where its fourth-moment estimates are stable. Whatever the choice, the sign of λi need not be known in advance.
On whitened data, FastICA iterates w←E[zg(w⊤z)]−E[g′(w⊤z)]w followed by normalization: an approximate Newton step for the Lagrange condition, in which whiteness makes the Jacobian diagonal. There is no step size.
In source coordinates the kurtosis iteration is exactly qi+=κiqi3: it converges globally and cubically to the source that maximizes ∣κi∣∣qi(0)∣.
For a general g, every source with λi=E[sig(si)−g′(si)]=0 is a fixed point with at least quadratic convergence (cubic for symmetric sources). Sub- and super-Gaussian sources are found with the same g; Gaussian sources are never fixed points.
Several components are obtained by deflation, which accumulates errors along the sequence, or by symmetric orthogonalization W←(WW⊤)−1/2W, the closest orthogonal matrix, which treats all units alike.
Bounded nonlinearities (logcosh, Gaussian) are robust to outliers; the accuracy for each source is governed by λi, and the best g is the score of the source.
Here g(u)=u3 and g′(u)=3u2, so E[g′(w⊤z)]=3E[(w⊤z)2]=3w⊤E[zz⊤]w=3∥w∥2=3. Without whitening, E[(w⊤x)2]=w⊤Cxw is no longer 1 on the unit sphere, so the coefficient changes from point to point; worse, the derivation used E[zz⊤]=I to make the Jacobian diagonal, and the identity qi+=κiqi3 relied on A~ being orthogonal. On unwhitened data the iteration is no longer an approximate Newton method and its fixed points need not be the source directions.
Suppose source 1 is Gaussian and the others are not. (a) Using Stein's lemma (Exercise 7.4习题 7.4), show that λ1=0 for every g. (b) With the kurtosis nonlinearity, what happens to the coordinate q1 after one iteration? (c) Can the one-unit iteration converge to ±e1?
(a) For s1∼N(0,1), Stein's lemma gives E[s1g(s1)]=E[g′(s1)], so λ1=0. (b) κ1=0, so q1+=κ1q13=0: the Gaussian coordinate is removed in a single step, and the iteration continues among the non-Gaussian sources. (c) No. At q=e1 the unnormalized update is λ1e1=0, so e1 is not even a well-defined fixed point; for kurtosis, part (a) of Theorem 8.3定理 8.3 shows that the limit is the non-Gaussian source maximizing ∣κi∣∣qi(0)∣. This is how FastICA copes with one Gaussian source: it finds the others, and in the symmetric version the Gaussian one is what remains in the orthogonal complement.
Let W=(1011). (a) Orthonormalize its rows by Gram–Schmidt, first in the order (row 1, row 2), then in the order (row 2, row 1). (b) Compute (WW⊤)−1/2W and verify that it is orthogonal. (c) Explain, using Proposition 8.5命题 8.5, why (b) is the natural choice for symmetric FastICA.
(a) Order (1, 2): w1=(1,1)/2; then (0,1)−21(1,1)=(−21,21), normalized w2=(−1,1)/2. Order (2, 1): w2=(0,1); then (1,1)−(0,1)=(1,0)=w1. The two orders give different matrices; whichever row comes first is kept in its direction. (b) WW⊤=(2111). For a 2×2 matrix with positive determinant the polar factor is the rotation (cosθsinθ−sinθcosθ) that maximizes tr(R(θ)⊤W)=(w11+w22)cosθ+(w21−w12)sinθ, so (cosθ,sinθ)∝(w11+w22,w21−w12)=(2,−1) and (WW⊤)−1/2W=51(2−112), whose rows are orthonormal (a direct computation of LR⊤ from the singular value decomposition gives the same matrix). The first row turns by about 18.4∘ clockwise and the second by about 26.6∘ counterclockwise: both rows move, instead of one of them staying fixed as in Gram–Schmidt. (c) The result is the orthogonal matrix closest to W, and permuting the rows of W permutes the result in the same way: no unit is privileged, so each one is corrected by the least possible amount instead of absorbing all the constraints of the units processed before it.
Let W=LΣR⊤ be invertible and scaled so that its largest singular value is at most 1 (dividing by ∥W∥F achieves this). Show that the iteration W←23W−21WW⊤W keeps L and R and maps each singular value by f(σ)=(3σ−σ3)/2. Prove that f(k)(σ)→1 for every σ∈(0,1], with 1−f(σ)=21(1−σ)2(2+σ), and conclude that the iteration converges quadratically to LR⊤.
设 W=LΣR⊤ 可逆，并已缩放使其最大奇异值不超过 1（除以 ∥W∥F 即可做到）。证明迭代 W←23W−21WW⊤W 保持 L 和 R 不变，并把每个奇异值映为 f(σ)=(3σ−σ3)/2。证明对每个 σ∈(0,1] 都有 f(k)(σ)→1，且 1−f(σ)=21(1−σ)2(2+σ)，由此推出该迭代二次收敛到 LR⊤。
Show solution查看解答
WW⊤W=LΣ3R⊤, so 23W−21WW⊤W=Lf(Σ)R⊤ with f applied to each singular value. Expanding, (1−σ)2(2+σ)=σ3−3σ+2, which is 2(1−f(σ)). On (0,1], f(σ)−σ=21σ(1−σ2)≥0 and f(σ)≤1, so f(k)(σ) increases and stays in (0,1]; its limit is a fixed point of f in (0,1], and the only one is 1. The identity gives 1−f(σ)≤23(1−σ)2, quadratic convergence of every singular value, hence of Lf(k)(Σ)R⊤ to LR⊤=(WW⊤)−1/2W (which does not change when W is rescaled).
WW⊤W=LΣ3R⊤，所以 23W−21WW⊤W=Lf(Σ)R⊤，其中 f 作用于每个奇异值。展开得 (1−σ)2(2+σ)=σ3−3σ+2，它等于 2(1−f(σ))。在 (0,1] 上，f(σ)−σ=21σ(1−σ2)≥0 且 f(σ)≤1，所以 f(k)(σ) 单调递增并保持在 (0,1] 内；其极限是 f 在 (0,1] 中的不动点，而唯一的不动点是 1。由恒等式得 1−f(σ)≤23(1−σ)2，即每个奇异值都二次收敛，从而 Lf(k)(Σ)R⊤ 二次收敛到 LR⊤=(WW⊤)−1/2W（缩放 W 不改变这个矩阵）。
so that the approximation E[zz⊤g′(w⊤z)]≈E[g′(w⊤z)]I of the derivation is exact on the tangent space ei⊥. (b) Use (a) to show directly that the tangent components of the unnormalized update at q=ei+δ, δ⊥ei, have no first-order term in δ, and explain why this makes the convergence superlinear.
(a) The (j,k) entry is E[sjskg′(si)]. If j=k, at least one index differs from i and appears once, so independence and zero mean give 0 (if j=i=k: E[sig′(si)]E[sk]=0). If j=k=i: E[sj2]E[g′(si)]=E[g′(si)]. If j=k=i: E[si2g′(si)]. This is the stated matrix: E[g′(si)] on the diagonal except in position (i,i). On vectors orthogonal to ei it acts as E[g′(si)]I, and E[g′(q⊤s)]=E[g′(si)] at q=ei, so the approximation is exact there. (b) The unnormalized update is u(q)=E[sg(q⊤s)]−E[g′(q⊤s)]q. Its derivative at ei in a direction δ⊥ei is E[ss⊤g′(si)]δ−E[g′(si)]δ−(E[s⊤g′′(si)]δ)ei. By (a) the first two terms cancel, and the last is parallel to ei (in fact zero, since E[skg′′(si)]=0 for k=i). So the tangent components of u(ei+δ) are o(∥δ∥), while the radial component tends to λi=0. After normalization the new tangential error is o(∥δ∥): the map has zero derivative at its fixed point, which implies superlinear convergence; the second-order term computed in Theorem 8.3定理 8.3 makes it at least quadratic. In Newton terms, the method uses the exact Jacobian on exactly the directions that matter.
Hyvärinen, A. (1998). New approximations of differential entropy for independent component analysis and projection pursuit. In Advances in Neural Information Processing Systems 10 (pp. 273–279). MIT Press.
Hyvärinen, A. (1999). Fast and robust fixed-point algorithms for independent component analysis. IEEE Transactions on Neural Networks, 10(3), 626–634.
Hyvärinen, A., & Oja, E. (1997). A fast fixed-point algorithm for independent component analysis. Neural Computation, 9(7), 1483–1492.
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.
Oja, E., & Yuan, Z. (2006). The FastICA algorithm revisited: Convergence analysis. IEEE Transactions on Neural Networks, 17(6), 1370–1381.
