The ICA BookICA 学习书 Yuda Bi

Appendix A附录 A

Mathematical Background数学基础

Facts used throughout: orthogonal and permutation matrices, matrix square roots, the polar decomposition, Kronecker products and tensors, matrix calculus, optimization on the sphere and the orthogonal group, and the asymptotic statistics of Part IV (delta method, sandwich covariance, Fisher information).

全书用到的数学工具:正交矩阵与置换矩阵、矩阵平方根、极分解、Kronecker 积与张量、矩阵微积分、球面与正交群上的优化,以及第四部分用到的渐近统计(delta 方法、三明治协方差、Fisher 信息)。

A.1 Conventions约定

This appendix collects the matrix and statistical facts used throughout the book. Short and instructive proofs are given in full; for the rest we point to standard references. All matrices are real. For matrices of the same size, ⟨A,B⟩=tr⁡(A⊤B)=∑ijaijbij\langle\bA, \bB\rangle = \tr(\bA^\top\bB) = \sum_{ij}a_{ij}b_{ij} is the Frobenius inner product and ∥A∥F=⟨A,A⟩1/2\|\bA\|_F = \langle\bA, \bA\rangle^{1/2} the Frobenius norm; ∥x∥\|\bx\| is the Euclidean norm of a vector. The identity tr⁡(AB)=tr⁡(BA)\tr(\bA\bB) = \tr(\bB\bA) is used constantly. For symmetric matrices, A⪰B\bA \succeq \bB (respectively A≻B\bA \succ \bB) means that A−B\bA - \bB is positive semidefinite (positive definite).

本附录汇集全书用到的矩阵与统计知识。简短而有启发的证明完整给出,其余给出标准参考文献。所有矩阵都是实矩阵。对同样大小的矩阵,⟨A,B⟩=tr⁡(A⊤B)=∑ijaijbij\langle\bA, \bB\rangle = \tr(\bA^\top\bB) = \sum_{ij}a_{ij}b_{ij} 是 Frobenius 内积,∥A∥F=⟨A,A⟩1/2\|\bA\|_F = \langle\bA, \bA\rangle^{1/2} 是 Frobenius 范数;∥x∥\|\bx\| 是向量的欧氏范数。恒等式 tr⁡(AB)=tr⁡(BA)\tr(\bA\bB) = \tr(\bB\bA) 会反复用到。对对称矩阵,A⪰B\bA \succeq \bB(或 A≻B\bA \succ \bB)表示 A−B\bA - \bB 半正定(或正定)。

Every square matrix splits as M=sym⁡(M)+skew⁡(M)\bM = \operatorname{sym}(\bM) + \operatorname{skew}(\bM), with symmetric part sym⁡(M)=(M+M⊤)/2\operatorname{sym}(\bM) = (\bM + \bM^\top)/2 and skew-symmetric part skew⁡(M)=(M−M⊤)/2\operatorname{skew}(\bM) = (\bM - \bM^\top)/2. The two parts are orthogonal: if S\bS is symmetric and K\bK is skew-symmetric, then tr⁡(SK)=tr⁡((SK)⊤)=−tr⁡(KS)=−tr⁡(SK)\tr(\bS\bK) = \tr\bigl((\bS\bK)^\top\bigr) = -\tr(\bK\bS) = -\tr(\bS\bK), so ⟨S,K⟩=0\langle\bS, \bK\rangle = 0. This simple fact underlies the gradients on the orthogonal group in Section A.7第 A.7 节.

每个方阵都可以分解为 M=sym⁡(M)+skew⁡(M)\bM = \operatorname{sym}(\bM) + \operatorname{skew}(\bM),其中对称部分为 sym⁡(M)=(M+M⊤)/2\operatorname{sym}(\bM) = (\bM + \bM^\top)/2,反对称部分为 skew⁡(M)=(M−M⊤)/2\operatorname{skew}(\bM) = (\bM - \bM^\top)/2。两部分相互正交:若 S\bS 对称、K\bK 反对称,则 tr⁡(SK)=tr⁡((SK)⊤)=−tr⁡(KS)=−tr⁡(SK)\tr(\bS\bK) = \tr\bigl((\bS\bK)^\top\bigr) = -\tr(\bK\bS) = -\tr(\bS\bK),所以 ⟨S,K⟩=0\langle\bS, \bK\rangle = 0。Section A.7第 A.7 节中正交群上的梯度就建立在这个简单事实之上。

A.2 Orthogonal and permutation matrices正交矩阵与置换矩阵

A square matrix Q\bQ is orthogonal if Q⊤Q=I\bQ^\top\bQ = \bI, equivalently QQ⊤=I\bQ\bQ^\top = \bI or Q−1=Q⊤\bQ^{-1} = \bQ^\top. Orthogonal matrices preserve inner products, (Qu)⊤(Qv)=u⊤v(\bQ\bu)^\top(\bQ\bv) = \bu^\top\bv, and hence lengths, angles and Frobenius norms: ∥QM∥F=∥MQ∥F=∥M∥F\|\bQ\bM\|_F = \|\bM\bQ\|_F = \|\bM\|_F. They form a group, the orthogonal group O(n)O(n). Since (det⁡Q)2=det⁡(Q⊤Q)=1(\det\bQ)^2 = \det(\bQ^\top\bQ) = 1, every orthogonal matrix has determinant ±1\pm 1. Those with determinant +1+1 form the group SO(n)SO(n) of rotations; the others include reflections such as diag⁡(1,…,1,−1)\diag(1, \dots, 1, -1).

若方阵 Q\bQ 满足 Q⊤Q=I\bQ^\top\bQ = \bI(等价地,QQ⊤=I\bQ\bQ^\top = \bI,或 Q−1=Q⊤\bQ^{-1} = \bQ^\top),则称它是正交矩阵。正交矩阵保持内积 (Qu)⊤(Qv)=u⊤v(\bQ\bu)^\top(\bQ\bv) = \bu^\top\bv,因而保持长度、夹角与 Frobenius 范数:∥QM∥F=∥MQ∥F=∥M∥F\|\bQ\bM\|_F = \|\bM\bQ\|_F = \|\bM\|_F。它们构成一个群,即正交群 O(n)O(n)。由于 (det⁡Q)2=det⁡(Q⊤Q)=1(\det\bQ)^2 = \det(\bQ^\top\bQ) = 1,正交矩阵的行列式为 ±1\pm 1。行列式为 +1+1 的正交矩阵构成旋转群 SO(n)SO(n);其余的包括 diag⁡(1,…,1,−1)\diag(1, \dots, 1, -1) 这样的反射。

Proposition A.1 (dimension of the orthogonal group)命题 A.1(正交群的维数)

O(n)O(n) is a smooth manifold of dimension n(n−1)/2n(n-1)/2. Its tangent space at Q\bQ is {ΩQ:Ω⊤=−Ω}\{\boldsymbol{\Omega}\bQ : \boldsymbol{\Omega}^\top = -\boldsymbol{\Omega}\}, and for every skew-symmetric Ω\boldsymbol{\Omega} the curve t↦exp⁡(tΩ)Qt \mapsto \exp(t\boldsymbol{\Omega})\bQ lies in O(n)O(n) and passes through Q\bQ with velocity ΩQ\boldsymbol{\Omega}\bQ.

O(n)O(n) 是维数为 n(n−1)/2n(n-1)/2 的光滑流形。它在 Q\bQ 处的切空间为 {ΩQ:Ω⊤=−Ω}\{\boldsymbol{\Omega}\bQ : \boldsymbol{\Omega}^\top = -\boldsymbol{\Omega}\};对每个反对称矩阵 Ω\boldsymbol{\Omega},曲线 t↦exp⁡(tΩ)Qt \mapsto \exp(t\boldsymbol{\Omega})\bQ 落在 O(n)O(n) 中,并以速度 ΩQ\boldsymbol{\Omega}\bQ 经过 Q\bQ。

Proof.证明

O(n)O(n) is the level set {Q:QQ⊤=I}\{\bQ : \bQ\bQ^\top = \bI\} of the map Q↦QQ⊤\bQ \mapsto \bQ\bQ^\top, whose derivative Q˙↦Q˙Q⊤+QQ˙⊤\dot\bQ \mapsto \dot\bQ\bQ^\top + \bQ\dot\bQ^\top maps onto the symmetric matrices (given a symmetric S\bS, take Q˙=12SQ\dot\bQ = \tfrac12\bS\bQ). By the regular value theorem, O(n)O(n) is a smooth manifold whose tangent space at Q\bQ is the kernel of this derivative. Writing Q˙=ΩQ\dot\bQ = \boldsymbol{\Omega}\bQ, the kernel condition becomes Ω+Ω⊤=0\boldsymbol{\Omega} + \boldsymbol{\Omega}^\top = \bzero. A skew-symmetric matrix is determined by its n(n−1)/2n(n-1)/2 entries above the diagonal, which gives the dimension. Finally, exp⁡(tΩ)⊤exp⁡(tΩ)=exp⁡(−tΩ)exp⁡(tΩ)=I\exp(t\boldsymbol{\Omega})^\top\exp(t\boldsymbol{\Omega}) = \exp(-t\boldsymbol{\Omega})\exp(t\boldsymbol{\Omega}) = \bI, and the derivative at t=0t = 0 is ΩQ\boldsymbol{\Omega}\bQ.

O(n)O(n) 是映射 Q↦QQ⊤\bQ \mapsto \bQ\bQ^\top 的水平集 {Q:QQ⊤=I}\{\bQ : \bQ\bQ^\top = \bI\},该映射的导数 Q˙↦Q˙Q⊤+QQ˙⊤\dot\bQ \mapsto \dot\bQ\bQ^\top + \bQ\dot\bQ^\top 映满全体对称矩阵(对给定的对称矩阵 S\bS,取 Q˙=12SQ\dot\bQ = \tfrac12\bS\bQ 即可)。由正则值定理,O(n)O(n) 是光滑流形,它在 Q\bQ 处的切空间就是这个导数的核。记 Q˙=ΩQ\dot\bQ = \boldsymbol{\Omega}\bQ,核的条件就变成 Ω+Ω⊤=0\boldsymbol{\Omega} + \boldsymbol{\Omega}^\top = \bzero。反对称矩阵由对角线以上的 n(n−1)/2n(n-1)/2 个元素确定,由此得到维数。最后,exp⁡(tΩ)⊤exp⁡(tΩ)=exp⁡(−tΩ)exp⁡(tΩ)=I\exp(t\boldsymbol{\Omega})^\top\exp(t\boldsymbol{\Omega}) = \exp(-t\boldsymbol{\Omega})\exp(t\boldsymbol{\Omega}) = \bI,且它在 t=0t = 0 处的导数为 ΩQ\boldsymbol{\Omega}\bQ。

For n=2n = 2, exp⁡(θ(0−110))=(cos⁡θ−sin⁡θsin⁡θcos⁡θ)\exp\bigl(\theta\bigl(\begin{smallmatrix} 0 & -1 \\ 1 & 0 \end{smallmatrix}\bigr)\bigr) = \bigl(\begin{smallmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{smallmatrix}\bigr), the rotation by θ\theta: one parameter, as n(n−1)/2=1n(n-1)/2 = 1 predicts. In general every rotation is the exponential of a skew-symmetric matrix. An orthogonal change of basis B\bB brings any orthogonal Q\bQ to a block-diagonal form with 2×22 \times 2 blocks exp⁡(θkJ)\exp(\theta_k\bJ), J=(0−110)\bJ = \bigl(\begin{smallmatrix} 0 & -1 \\ 1 & 0 \end{smallmatrix}\bigr), and diagonal entries ±1\pm 1 (Horn & Johnson, 2013). When det⁡Q=1\det\bQ = 1 the entries −1-1 come in pairs, each pair forming a block exp⁡(πJ)=−I2\exp(\pi\bJ) = -\bI_2, so Q=exp⁡(Ω)\bQ = \exp(\boldsymbol{\Omega}) with the skew-symmetric matrix Ω=B diag⁡(θ1J,…,θkJ,0,…,0) B⊤\boldsymbol{\Omega} = \bB\,\diag(\theta_1\bJ, \dots, \theta_k\bJ, 0, \dots, 0)\,\bB^\top. Conversely, det⁡exp⁡(Ω)=etr⁡Ω=1\det\exp(\boldsymbol{\Omega}) = e^{\tr\boldsymbol{\Omega}} = 1.

当 n=2n = 2 时,exp⁡(θ(0−110))=(cos⁡θ−sin⁡θsin⁡θcos⁡θ)\exp\bigl(\theta\bigl(\begin{smallmatrix} 0 & -1 \\ 1 & 0 \end{smallmatrix}\bigr)\bigr) = \bigl(\begin{smallmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{smallmatrix}\bigr),即旋转 θ\theta 角:只有一个参数,与 n(n−1)/2=1n(n-1)/2 = 1 一致。一般地,每个旋转都是某个反对称矩阵的指数。通过正交基变换 B\bB,任何正交矩阵 Q\bQ 都可以化为分块对角形式,其中的 2×22 \times 2 块为 exp⁡(θkJ)\exp(\theta_k\bJ),J=(0−110)\bJ = \bigl(\begin{smallmatrix} 0 & -1 \\ 1 & 0 \end{smallmatrix}\bigr),其余对角元为 ±1\pm 1(Horn & Johnson, 2013)。当 det⁡Q=1\det\bQ = 1 时,对角元 −1-1 成对出现,每一对构成一个块 exp⁡(πJ)=−I2\exp(\pi\bJ) = -\bI_2,于是 Q=exp⁡(Ω)\bQ = \exp(\boldsymbol{\Omega}),其中 Ω=B diag⁡(θ1J,…,θkJ,0,…,0) B⊤\boldsymbol{\Omega} = \bB\,\diag(\theta_1\bJ, \dots, \theta_k\bJ, 0, \dots, 0)\,\bB^\top 是反对称矩阵。反之,det⁡exp⁡(Ω)=etr⁡Ω=1\det\exp(\boldsymbol{\Omega}) = e^{\tr\boldsymbol{\Omega}} = 1。

Givens rotations. The Givens rotation G(i,j,θ)\bG(i, j, \theta) equals the identity except for the entries gii=gjj=cos⁡θg_{ii} = g_{jj} = \cos\theta, gij=−sin⁡θg_{ij} = -\sin\theta and gji=sin⁡θg_{ji} = \sin\theta; it rotates the (i,j)(i, j) coordinate plane by θ\theta and leaves the other coordinates unchanged. There are n(n−1)/2n(n-1)/2 such planes, matching the dimension of O(n)O(n). Every orthogonal matrix is a product of at most n(n−1)/2n(n-1)/2 Givens rotations and a diagonal matrix of signs: applying Givens rotations to zero the entries below the diagonal of Q\bQ column by column, as in the Givens QR factorization (Golub & Van Loan, 2013), leaves an upper triangular orthogonal matrix, which must be diagonal with entries ±1\pm 1. Jacobi-type algorithms such as JADE (Chapter 9第 9 章) search over O(n)O(n) by sweeping through the planes and optimizing one angle at a time.

Givens 旋转。Givens 旋转 G(i,j,θ)\bG(i, j, \theta) 除 gii=gjj=cos⁡θg_{ii} = g_{jj} = \cos\theta、gij=−sin⁡θg_{ij} = -\sin\theta、gji=sin⁡θg_{ji} = \sin\theta 之外与单位阵相同;它把 (i,j)(i, j) 坐标平面旋转 θ\theta 角,而保持其余坐标不变。这样的平面共有 n(n−1)/2n(n-1)/2 个,恰好等于 O(n)O(n) 的维数。每个正交矩阵都是至多 n(n−1)/2n(n-1)/2 个 Givens 旋转与一个符号对角矩阵的乘积:像用 Givens 旋转做 QR 分解那样(Golub & Van Loan, 2013),逐列消去 Q\bQ 对角线以下的元素,剩下一个上三角的正交矩阵,它必然是对角元为 ±1\pm 1 的对角矩阵。JADE(Chapter 9第 9 章)等 Jacobi 型算法在 O(n)O(n) 中搜索的方式,就是依次扫过各个平面,每次只优化一个角度。

Permutation matrices. For a permutation π\pi of {1,…,n}\{1, \dots, n\}, the permutation matrix Π\bPi has entries πij=1\pi_{ij} = 1 if i=π(j)i = \pi(j) and 00 otherwise, so (Πs)π(j)=sj(\bPi\bs)_{\pi(j)} = s_j: it moves component jj to position π(j)\pi(j). Permutation matrices are orthogonal, and so are the signed permutation matrices ΠD\bPi\bD with D=diag⁡(±1)\bD = \diag(\pm 1). These are the only orthogonal scaled permutations.

置换矩阵。对 {1,…,n}\{1, \dots, n\} 的置换 π\pi,置换矩阵 Π\bPi 的元素为:当 i=π(j)i = \pi(j) 时 πij=1\pi_{ij} = 1,否则为 00。于是 (Πs)π(j)=sj(\bPi\bs)_{\pi(j)} = s_j:它把第 jj 个分量移到第 π(j)\pi(j) 个位置。置换矩阵是正交的,带符号的置换矩阵 ΠD\bPi\bD(D=diag⁡(±1)\bD = \diag(\pm 1))也是正交的,而且它们是仅有的正交的带尺度置换矩阵。

Lemma A.2 (orthogonal scaled permutations)引理 A.2(正交的带尺度置换矩阵)

Let Π\bPi be a permutation matrix and D\bD a diagonal matrix. Then ΠD\bPi\bD is orthogonal if and only if D=diag⁡(±1)\bD = \diag(\pm 1). The signed permutation matrices form a subgroup of O(n)O(n) with 2nn!2^n n! elements.

设 Π\bPi 为置换矩阵,D\bD 为对角矩阵。则 ΠD\bPi\bD 正交当且仅当 D=diag⁡(±1)\bD = \diag(\pm 1)。带符号的置换矩阵构成 O(n)O(n) 的一个子群,共有 2nn!2^n n! 个元素。

Proof.证明

(ΠD)⊤(ΠD)=DΠ⊤ΠD=D2(\bPi\bD)^\top(\bPi\bD) = \bD\bPi^\top\bPi\bD = \bD^2, which equals I\bI exactly when every di2=1d_i^2 = 1. For two signed permutations, DΠ′=Π′(Π′⊤DΠ′)\bD\bPi' = \bPi'(\bPi'^\top\bD\bPi') and Π′⊤DΠ′\bPi'^\top\bD\bPi' is again a diagonal sign matrix, so the product (ΠD)(Π′D′)(\bPi\bD)(\bPi'\bD') is a signed permutation; so is the inverse (ΠD)⊤=Π⊤(ΠDΠ⊤)(\bPi\bD)^\top = \bPi^\top(\bPi\bD\bPi^\top). Distinct pairs (Π,D)(\bPi, \bD) give distinct matrices, because the position of the nonzero entry in each column determines Π\bPi and its sign determines D\bD; hence there are n!⋅2nn!\cdot 2^n of them.

(ΠD)⊤(ΠD)=DΠ⊤ΠD=D2(\bPi\bD)^\top(\bPi\bD) = \bD\bPi^\top\bPi\bD = \bD^2,它等于 I\bI 当且仅当每个 di2=1d_i^2 = 1。对两个带符号置换矩阵,DΠ′=Π′(Π′⊤DΠ′)\bD\bPi' = \bPi'(\bPi'^\top\bD\bPi'),而 Π′⊤DΠ′\bPi'^\top\bD\bPi' 仍是对角符号矩阵,所以乘积 (ΠD)(Π′D′)(\bPi\bD)(\bPi'\bD') 仍是带符号置换矩阵;逆矩阵 (ΠD)⊤=Π⊤(ΠDΠ⊤)(\bPi\bD)^\top = \bPi^\top(\bPi\bD\bPi^\top) 也是。不同的 (Π,D)(\bPi, \bD) 给出不同的矩阵,因为每列非零元的位置确定 Π\bPi,其符号确定 D\bD;因此共有 n!⋅2nn!\cdot 2^n 个。

In ICA, separation means that the global matrix is a scaled permutation ΠD\bPi\bD (Proposition 1.3命题 1.3). After whitening the global matrix is orthogonal (Proposition 3.4命题 3.4), and the lemma shows that the scaling is then reduced to signs: the solutions form the finite group of signed permutations.

在 ICA 中,分离意味着全局矩阵是带尺度的置换矩阵 ΠD\bPi\bD(Proposition 1.3命题 1.3)。白化之后全局矩阵是正交的(Proposition 3.4命题 3.4),由引理可知,此时尺度只剩下符号:所有解构成带符号置换矩阵这个有限群。

A.3 Symmetric matrices and matrix functions对称矩阵与矩阵函数

Spectral theorem. Every symmetric C∈Rn×n\bC \in \R^{n \times n} can be written C=EΛE⊤=∑iλieiei⊤\bC = \mathbf{E}\bLambda\mathbf{E}^\top = \sum_i \lambda_i\be_i\be_i^\top, with E\mathbf{E} orthogonal with columns ei\be_i and Λ=diag⁡(λ1,…,λn)\bLambda = \diag(\lambda_1, \dots, \lambda_n) real (Horn & Johnson, 2013). The columns are orthonormal eigenvectors, Cei=λiei\bC\be_i = \lambda_i\be_i. If the eigenvalues are distinct, each eigenvector is unique up to sign; for a repeated eigenvalue, any orthonormal basis of its eigenspace can be used.

谱定理。每个对称矩阵 C∈Rn×n\bC \in \R^{n \times n} 都可以写成 C=EΛE⊤=∑iλieiei⊤\bC = \mathbf{E}\bLambda\mathbf{E}^\top = \sum_i \lambda_i\be_i\be_i^\top,其中 E\mathbf{E} 是以 ei\be_i 为列的正交矩阵,Λ=diag⁡(λ1,…,λn)\bLambda = \diag(\lambda_1, \dots, \lambda_n) 为实对角矩阵(Horn & Johnson, 2013)。各列是标准正交的特征向量,Cei=λiei\bC\be_i = \lambda_i\be_i。若特征值互不相同,则每个特征向量在相差一个符号的意义下唯一;对重特征值,可以取其特征子空间的任意一组标准正交基。

Definiteness. C\bC is positive semidefinite, C⪰0\bC \succeq \bzero, if v⊤Cv≥0\bv^\top\bC\bv \ge 0 for all v\bv, equivalently if all λi≥0\lambda_i \ge 0; it is positive definite, C≻0\bC \succ \bzero, if v⊤Cv>0\bv^\top\bC\bv > 0 for all v≠0\bv \ne \bzero, equivalently if all λi>0\lambda_i > 0. Covariance matrices are positive semidefinite, since v⊤Cxv=Var⁡(v⊤x)\bv^\top\bC_{\bx}\bv = \Var(\bv^\top\bx). A positive definite matrix has a unique Cholesky factorization C=LL⊤\bC = \mathbf{L}\mathbf{L}^\top with L\mathbf{L} lower triangular with a positive diagonal (Golub & Van Loan, 2013).

正定性。若对所有 v\bv 都有 v⊤Cv≥0\bv^\top\bC\bv \ge 0,等价地,所有 λi≥0\lambda_i \ge 0,则称 C\bC 半正定,记作 C⪰0\bC \succeq \bzero;若对所有 v≠0\bv \ne \bzero 都有 v⊤Cv>0\bv^\top\bC\bv > 0,等价地,所有 λi>0\lambda_i > 0,则称 C\bC 正定,记作 C≻0\bC \succ \bzero。协方差矩阵是半正定的,因为 v⊤Cxv=Var⁡(v⊤x)\bv^\top\bC_{\bx}\bv = \Var(\bv^\top\bx)。正定矩阵有唯一的 Cholesky 分解 C=LL⊤\bC = \mathbf{L}\mathbf{L}^\top,其中 L\mathbf{L} 是对角元为正的下三角矩阵(Golub & Van Loan, 2013)。

Rayleigh quotient. For a unit vector w\bw, w⊤Cw=∑iλi(ei⊤w)2\bw^\top\bC\bw = \sum_i \lambda_i(\be_i^\top\bw)^2 is an average of the eigenvalues with weights (ei⊤w)2(\be_i^\top\bw)^2 that sum to one. Hence max⁡∥w∥=1w⊤Cw=λmax⁡\max_{\|\bw\| = 1}\bw^\top\bC\bw = \lambda_{\max} and min⁡∥w∥=1w⊤Cw=λmin⁡\min_{\|\bw\| = 1}\bw^\top\bC\bw = \lambda_{\min}, attained at the corresponding eigenvectors. This variational characterization is the basis of PCA, and Section A.7第 A.7 节 revisits it as an optimization problem on the sphere.

Rayleigh 商。对单位向量 w\bw,w⊤Cw=∑iλi(ei⊤w)2\bw^\top\bC\bw = \sum_i \lambda_i(\be_i^\top\bw)^2 是各特征值以 (ei⊤w)2(\be_i^\top\bw)^2 为权重的加权平均,权重之和为一。因此 max⁡∥w∥=1w⊤Cw=λmax⁡\max_{\|\bw\| = 1}\bw^\top\bC\bw = \lambda_{\max},min⁡∥w∥=1w⊤Cw=λmin⁡\min_{\|\bw\| = 1}\bw^\top\bC\bw = \lambda_{\min},分别在相应的特征向量处取到。这一变分刻画是 PCA 的基础,Section A.7第 A.7 节还会把它作为球面上的优化问题重新讨论。

Matrix functions. For a function ff defined on the eigenvalues, f(C)=Ef(Λ)E⊤=∑if(λi)eiei⊤f(\bC) = \mathbf{E}f(\bLambda)\mathbf{E}^\top = \sum_i f(\lambda_i)\be_i\be_i^\top. The result does not depend on the choice of eigenvectors, because f(C)=p(C)f(\bC) = p(\bC) for any polynomial pp with p(λi)=f(λi)p(\lambda_i) = f(\lambda_i) for all ii. Examples are powers, exp⁡(C)\exp(\bC), log⁡C\log\bC for C≻0\bC \succ \bzero, and square roots. Matrix functions commute with orthogonal similarity: f(QCQ⊤)=Qf(C)Q⊤f(\bQ\bC\bQ^\top) = \bQ f(\bC)\bQ^\top.

矩阵函数。对定义在各特征值上的函数 ff,令 f(C)=Ef(Λ)E⊤=∑if(λi)eiei⊤f(\bC) = \mathbf{E}f(\bLambda)\mathbf{E}^\top = \sum_i f(\lambda_i)\be_i\be_i^\top。这个定义与特征向量的选取无关,因为对任何满足 p(λi)=f(λi)p(\lambda_i) = f(\lambda_i)(对所有 ii)的多项式 pp,都有 f(C)=p(C)f(\bC) = p(\bC)。例子有幂、exp⁡(C)\exp(\bC)、C≻0\bC \succ \bzero 时的 log⁡C\log\bC,以及平方根。矩阵函数与正交相似变换可交换:f(QCQ⊤)=Qf(C)Q⊤f(\bQ\bC\bQ^\top) = \bQ f(\bC)\bQ^\top。

Proposition A.3 (matrix square roots)命题 A.3(矩阵平方根)

Every positive semidefinite C\bC has a unique positive semidefinite square root, that is, a matrix B⪰0\bB \succeq \bzero with B2=C\bB^2 = \bC, namely C1/2=EΛ1/2E⊤\bC^{1/2} = \mathbf{E}\bLambda^{1/2}\mathbf{E}^\top. If C≻0\bC \succ \bzero, then C−1/2=(C1/2)−1=EΛ−1/2E⊤\bC^{-1/2} = (\bC^{1/2})^{-1} = \mathbf{E}\bLambda^{-1/2}\mathbf{E}^\top and C−1/2C C−1/2=I\bC^{-1/2}\bC\,\bC^{-1/2} = \bI.

每个半正定矩阵 C\bC 都有唯一的半正定平方根,即满足 B2=C\bB^2 = \bC 的矩阵 B⪰0\bB \succeq \bzero,它就是 C1/2=EΛ1/2E⊤\bC^{1/2} = \mathbf{E}\bLambda^{1/2}\mathbf{E}^\top。若 C≻0\bC \succ \bzero,则 C−1/2=(C1/2)−1=EΛ−1/2E⊤\bC^{-1/2} = (\bC^{1/2})^{-1} = \mathbf{E}\bLambda^{-1/2}\mathbf{E}^\top,且 C−1/2C C−1/2=I\bC^{-1/2}\bC\,\bC^{-1/2} = \bI。

Proof.证明

EΛ1/2E⊤\mathbf{E}\bLambda^{1/2}\mathbf{E}^\top is positive semidefinite and its square is EΛE⊤=C\mathbf{E}\bLambda\mathbf{E}^\top = \bC. For uniqueness, let B⪰0\bB \succeq \bzero with B2=C\bB^2 = \bC. Then B\bB commutes with C=B2\bC = \bB^2, so it maps each eigenspace of C\bC into itself: if Cv=λv\bC\bv = \lambda\bv, then C(Bv)=BCv=λBv\bC(\bB\bv) = \bB\bC\bv = \lambda\bB\bv. Restricted to the eigenspace of λ\lambda, B\bB is symmetric and positive semidefinite with square λI\lambda\bI, so its eigenvalues are nonnegative square roots of λ\lambda, and it equals λ I\sqrt\lambda\,\bI there. Hence B=∑iλi eiei⊤\bB = \sum_i \sqrt{\lambda_i}\,\be_i\be_i^\top. The statements about C−1/2\bC^{-1/2} follow from the eigendecomposition.

EΛ1/2E⊤\mathbf{E}\bLambda^{1/2}\mathbf{E}^\top 是半正定的,且其平方为 EΛE⊤=C\mathbf{E}\bLambda\mathbf{E}^\top = \bC。再证唯一性:设 B⪰0\bB \succeq \bzero 且 B2=C\bB^2 = \bC。则 B\bB 与 C=B2\bC = \bB^2 可交换,所以 B\bB 把 C\bC 的每个特征子空间映到自身:若 Cv=λv\bC\bv = \lambda\bv,则 C(Bv)=BCv=λBv\bC(\bB\bv) = \bB\bC\bv = \lambda\bB\bv。把 B\bB 限制在特征值 λ\lambda 的特征子空间上,它是对称半正定的,且平方为 λI\lambda\bI,所以它的特征值都是 λ\lambda 的非负平方根,从而它在该子空间上等于 λ I\sqrt\lambda\,\bI。因此 B=∑iλi eiei⊤\bB = \sum_i \sqrt{\lambda_i}\,\be_i\be_i^\top。关于 C−1/2\bC^{-1/2} 的结论由特征分解直接得到。

Non-symmetric square roots, that is, matrices B\bB with BB⊤=C\bB\bB^\top = \bC such as the Cholesky factor, are not unique. If C≻0\bC \succ \bzero and BB⊤=B′B′⊤=C\bB\bB^\top = \bB'\bB'^\top = \bC, then Q=B−1B′\bQ = \bB^{-1}\bB' satisfies QQ⊤=B−1CB−⊤=I\bQ\bQ^\top = \bB^{-1}\bC\bB^{-\top} = \bI, so B′=BQ\bB' = \bB\bQ with Q\bQ orthogonal, and every such BQ\bB\bQ works. Inverting gives the whitening matrices of Proposition 3.2命题 3.2.

非对称的平方根,即满足 BB⊤=C\bB\bB^\top = \bC 的矩阵 B\bB(例如 Cholesky 因子),并不唯一。若 C≻0\bC \succ \bzero 且 BB⊤=B′B′⊤=C\bB\bB^\top = \bB'\bB'^\top = \bC,则 Q=B−1B′\bQ = \bB^{-1}\bB' 满足 QQ⊤=B−1CB−⊤=I\bQ\bQ^\top = \bB^{-1}\bC\bB^{-\top} = \bI,所以 B′=BQ\bB' = \bB\bQ,Q\bQ 正交;反之每个这样的 BQ\bB\bQ 都满足要求。取逆即得Proposition 3.2命题 3.2中的白化矩阵。

A.4 The SVD and the polar decomposition奇异值分解与极分解

Every W∈Rn×m\bW \in \R^{n \times m} has a singular value decomposition (SVD) W=LΣR⊤\bW = \mathbf{L}\bSigma\bR^\top, with L∈O(n)\mathbf{L} \in O(n), R∈O(m)\bR \in O(m), and Σ∈Rn×m\bSigma \in \R^{n \times m} zero except for the diagonal entries σ1≥σ2≥⋯≥0\sigma_1 \ge \sigma_2 \ge \dots \ge 0, the singular values (Golub & Van Loan, 2013). We call the factors L\mathbf{L} and R\bR, for left and right, to keep U\bU and V\bV for the rotation and the whitening matrix of the main text. The σi2\sigma_i^2 are the eigenvalues of WW⊤\bW\bW^\top and of W⊤W\bW^\top\bW, whose eigenvectors are the columns of L\mathbf{L} and of R\bR. A square W\bW is invertible exactly when all its singular values are positive.

每个 W∈Rn×m\bW \in \R^{n \times m} 都有奇异值分解(SVD)W=LΣR⊤\bW = \mathbf{L}\bSigma\bR^\top,其中 L∈O(n)\mathbf{L} \in O(n),R∈O(m)\bR \in O(m),Σ∈Rn×m\bSigma \in \R^{n \times m} 除对角元 σ1≥σ2≥⋯≥0\sigma_1 \ge \sigma_2 \ge \dots \ge 0(即奇异值)外全为零(Golub & Van Loan, 2013)。我们把两个因子记作 L\mathbf{L} 与 R\bR(取 left 与 right 之意),以便把 U\bU 与 V\bV 留给正文中的旋转与白化矩阵。σi2\sigma_i^2 是 WW⊤\bW\bW^\top 与 W⊤W\bW^\top\bW 的特征值,相应的特征向量分别是 L\mathbf{L} 与 R\bR 的各列。方阵 W\bW 可逆,当且仅当它的奇异值全为正。

Polar decomposition. An invertible square W\bW factors as

极分解。可逆方阵 W\bW 可以分解为

W=(WW⊤)1/2 QW,QW=(WW⊤)−1/2 W=LR⊤,(A.1)\bW = (\bW\bW^\top)^{1/2}\,\bQ_{\bW}, \qquad \bQ_{\bW} = (\bW\bW^\top)^{-1/2}\,\bW = \mathbf{L}\bR^\top,\tag{A.1}

where (WW⊤)1/2=LΣL⊤(\bW\bW^\top)^{1/2} = \mathbf{L}\bSigma\mathbf{L}^\top is positive definite and QW\bQ_{\bW} is orthogonal; indeed (WW⊤)−1/2W=LΣ−1L⊤LΣR⊤=LR⊤(\bW\bW^\top)^{-1/2}\bW = \mathbf{L}\bSigma^{-1}\mathbf{L}^\top\mathbf{L}\bSigma\bR^\top = \mathbf{L}\bR^\top. This is the matrix analogue of writing a nonzero complex number as a modulus times a phase. The same orthogonal factor appears on the other side, W=QW(W⊤W)1/2\bW = \bQ_{\bW}(\bW^\top\bW)^{1/2}. A singular W\bW also has polar decompositions, but its orthogonal factor is not unique.

其中 (WW⊤)1/2=LΣL⊤(\bW\bW^\top)^{1/2} = \mathbf{L}\bSigma\mathbf{L}^\top 正定,QW\bQ_{\bW} 正交;事实上 (WW⊤)−1/2W=LΣ−1L⊤LΣR⊤=LR⊤(\bW\bW^\top)^{-1/2}\bW = \mathbf{L}\bSigma^{-1}\mathbf{L}^\top\mathbf{L}\bSigma\bR^\top = \mathbf{L}\bR^\top。这相当于把非零复数写成模乘以相位的矩阵版本。同一个正交因子也出现在另一侧:W=QW(W⊤W)1/2\bW = \bQ_{\bW}(\bW^\top\bW)^{1/2}。奇异的 W\bW 也有极分解,但其正交因子不唯一。

Proposition A.4 (the closest orthogonal matrix)命题 A.4(最近的正交矩阵)

Let W∈Rn×n\bW \in \R^{n \times n} be invertible. Then (WW⊤)−1/2W=LR⊤(\bW\bW^\top)^{-1/2}\bW = \mathbf{L}\bR^\top is the unique orthogonal matrix closest to W\bW in Frobenius norm: ∥W−Q∥F>∥W−(WW⊤)−1/2W∥F\|\bW - \bQ\|_F > \|\bW - (\bW\bW^\top)^{-1/2}\bW\|_F for every other orthogonal Q\bQ.

设 W∈Rn×n\bW \in \R^{n \times n} 可逆。则 (WW⊤)−1/2W=LR⊤(\bW\bW^\top)^{-1/2}\bW = \mathbf{L}\bR^\top 是在 Frobenius 范数下离 W\bW 最近的唯一正交矩阵:对其他任何正交矩阵 Q\bQ,∥W−Q∥F>∥W−(WW⊤)−1/2W∥F\|\bW - \bQ\|_F > \|\bW - (\bW\bW^\top)^{-1/2}\bW\|_F。

Proof.证明

For orthogonal Q\bQ, ∥W−Q∥F2=∥W∥F2+n−2tr⁡(Q⊤W)\|\bW - \bQ\|_F^2 = \|\bW\|_F^2 + n - 2\tr(\bQ^\top\bW), so we must maximize tr⁡(Q⊤W)=tr⁡(Q⊤LΣR⊤)=tr⁡(MΣ)\tr(\bQ^\top\bW) = \tr(\bQ^\top\mathbf{L}\bSigma\bR^\top) = \tr(\bM\bSigma) with M=R⊤Q⊤L\bM = \bR^\top\bQ^\top\mathbf{L}, which is orthogonal. Then tr⁡(MΣ)=∑imiiσi≤∑iσi\tr(\bM\bSigma) = \sum_i m_{ii}\sigma_i \le \sum_i \sigma_i, because each column of M\bM is a unit vector and so ∣mii∣≤1|m_{ii}| \le 1. Since every σi>0\sigma_i > 0, equality holds if and only if mii=1m_{ii} = 1 for all ii; a unit column with ii-th entry 1 is the ii-th column of the identity, so this means M=I\bM = \bI, that is, Q=LR⊤\bQ = \mathbf{L}\bR^\top. By (A.1), LR⊤=(WW⊤)−1/2W\mathbf{L}\bR^\top = (\bW\bW^\top)^{-1/2}\bW.

对正交矩阵 Q\bQ,∥W−Q∥F2=∥W∥F2+n−2tr⁡(Q⊤W)\|\bW - \bQ\|_F^2 = \|\bW\|_F^2 + n - 2\tr(\bQ^\top\bW),所以问题化为最大化 tr⁡(Q⊤W)=tr⁡(Q⊤LΣR⊤)=tr⁡(MΣ)\tr(\bQ^\top\bW) = \tr(\bQ^\top\mathbf{L}\bSigma\bR^\top) = \tr(\bM\bSigma),其中 M=R⊤Q⊤L\bM = \bR^\top\bQ^\top\mathbf{L} 是正交矩阵。于是 tr⁡(MΣ)=∑imiiσi≤∑iσi\tr(\bM\bSigma) = \sum_i m_{ii}\sigma_i \le \sum_i \sigma_i,因为 M\bM 的每一列都是单位向量,所以 ∣mii∣≤1|m_{ii}| \le 1。由于所有 σi>0\sigma_i > 0,等号成立当且仅当对所有 ii 都有 mii=1m_{ii} = 1;第 ii 个分量为 1 的单位列向量就是单位阵的第 ii 列,所以这意味着 M=I\bM = \bI,即 Q=LR⊤\bQ = \mathbf{L}\bR^\top。由 (A.1),LR⊤=(WW⊤)−1/2W\mathbf{L}\bR^\top = (\bW\bW^\top)^{-1/2}\bW。

Uses in ICA. Symmetric FastICA orthogonalizes its current estimate by W←(WW⊤)−1/2W\bW \leftarrow (\bW\bW^\top)^{-1/2}\bW ((8.5)), which by the proposition is the orthogonal matrix closest to W\bW. It treats all rows alike: since matrix functions commute with orthogonal similarity, replacing W\bW by ΠW\bPi\bW replaces the result by Π(WW⊤)−1/2W\bPi(\bW\bW^\top)^{-1/2}\bW, whereas Gram–Schmidt orthogonalization depends on the order of the rows. ZCA whitening produces the whitened mixing matrix (AA⊤)−1/2A(\bA\bA^\top)^{-1/2}\bA, the polar factor of A\bA (Section 3.4第 3.4 节). Numerically, QW\bQ_{\bW} is computed from the SVD, or by iterations such as Newton's iteration X←(X+X−⊤)/2\bX \leftarrow (\bX + \bX^{-\top})/2 started at X=W\bX = \bW (Higham, 1986). Applied to BA⊤\bB\bA^\top in place of W\bW, the same argument solves the orthogonal Procrustes problem min⁡Q∈O(n)∥QA−B∥F\min_{\bQ \in O(n)}\|\bQ\bA - \bB\|_F.

在 ICA 中的应用。对称 FastICA 用 W←(WW⊤)−1/2W\bW \leftarrow (\bW\bW^\top)^{-1/2}\bW 对当前估计作正交化((8.5)),由命题可知,这正是离 W\bW 最近的正交矩阵。它对各行一视同仁:由于矩阵函数与正交相似变换可交换,把 W\bW 换成 ΠW\bPi\bW,结果就变成 Π(WW⊤)−1/2W\bPi(\bW\bW^\top)^{-1/2}\bW;而 Gram–Schmidt 正交化的结果依赖于各行的次序。ZCA 白化给出的白化混合矩阵 (AA⊤)−1/2A(\bA\bA^\top)^{-1/2}\bA 正是 A\bA 的极因子(Section 3.4第 3.4 节)。数值上,QW\bQ_{\bW} 可由 SVD 求得,也可以用迭代法计算,例如从 X=W\bX = \bW 出发的 Newton 迭代 X←(X+X−⊤)/2\bX \leftarrow (\bX + \bX^{-\top})/2(Higham, 1986)。把 W\bW 换成 BA⊤\bB\bA^\top,同样的论证还可以解决正交 Procrustes 问题 min⁡Q∈O(n)∥QA−B∥F\min_{\bQ \in O(n)}\|\bQ\bA - \bB\|_F。

A.5 Kronecker products, vec and tensorsKronecker 积、vec 与张量

Kronecker products. For A∈Rp×q\bA \in \R^{p \times q} and B∈Rr×s\bB \in \R^{r \times s}, the Kronecker product A⊗B∈Rpr×qs\bA \otimes \bB \in \R^{pr \times qs} is the block matrix whose (i,j)(i, j) block is aijBa_{ij}\bB. It satisfies (A⊗B)(C⊗D)=AC⊗BD(\bA \otimes \bB)(\bC \otimes \bD) = \bA\bC \otimes \bB\bD whenever the products exist, (A⊗B)⊤=A⊤⊗B⊤(\bA \otimes \bB)^\top = \bA^\top \otimes \bB^\top, and (A⊗B)−1=A−1⊗B−1(\bA \otimes \bB)^{-1} = \bA^{-1} \otimes \bB^{-1}. For square A\bA (p×pp \times p) and B\bB (r×rr \times r), the eigenvalues of A⊗B\bA \otimes \bB are the products of an eigenvalue of A\bA and an eigenvalue of B\bB, so tr⁡(A⊗B)=tr⁡A tr⁡B\tr(\bA \otimes \bB) = \tr\bA\,\tr\bB and det⁡(A⊗B)=(det⁡A)r(det⁡B)p\det(\bA \otimes \bB) = (\det\bA)^r(\det\bB)^p. By the first two rules, the Kronecker product of orthogonal matrices is orthogonal.

Kronecker 积。对 A∈Rp×q\bA \in \R^{p \times q} 与 B∈Rr×s\bB \in \R^{r \times s},Kronecker 积 A⊗B∈Rpr×qs\bA \otimes \bB \in \R^{pr \times qs} 是第 (i,j)(i, j) 块为 aijBa_{ij}\bB 的分块矩阵。它满足:只要乘积有意义,就有 (A⊗B)(C⊗D)=AC⊗BD(\bA \otimes \bB)(\bC \otimes \bD) = \bA\bC \otimes \bB\bD;(A⊗B)⊤=A⊤⊗B⊤(\bA \otimes \bB)^\top = \bA^\top \otimes \bB^\top;(A⊗B)−1=A−1⊗B−1(\bA \otimes \bB)^{-1} = \bA^{-1} \otimes \bB^{-1}。对方阵 A\bA(p×pp \times p)与 B\bB(r×rr \times r),A⊗B\bA \otimes \bB 的特征值是 A\bA 的一个特征值与 B\bB 的一个特征值之积,所以 tr⁡(A⊗B)=tr⁡A tr⁡B\tr(\bA \otimes \bB) = \tr\bA\,\tr\bB,det⁡(A⊗B)=(det⁡A)r(det⁡B)p\det(\bA \otimes \bB) = (\det\bA)^r(\det\bB)^p。由前两条规则,正交矩阵的 Kronecker 积仍是正交矩阵。

The vec operator. vec⁡(X)\vecop(\bX) stacks the columns of X\bX into one long vector. The key identity is

vec 算子。vec⁡(X)\vecop(\bX) 把 X\bX 的各列依次叠成一个长向量。关键的恒等式是

vec⁡(AXB)=(B⊤⊗A) vec⁡(X).(A.2)\vecop(\bA\bX\bB) = (\bB^\top \otimes \bA)\,\vecop(\bX).\tag{A.2}

To see it, note that column kk of AXB\bA\bX\bB is ∑jbjkAxj\sum_j b_{jk}\bA\bx_j, where xj\bx_j is column jj of X\bX; this is block kk of (B⊤⊗A)vec⁡(X)(\bB^\top \otimes \bA)\vecop(\bX), whose (k,j)(k, j) block is bjkAb_{jk}\bA. Special cases are vec⁡(ab⊤)=b⊗a\vecop(\ba\bb^\top) = \bb \otimes \ba and tr⁡(A⊤B)=vec⁡(A)⊤vec⁡(B)\tr(\bA^\top\bB) = \vecop(\bA)^\top\vecop(\bB). In NumPy, which stores arrays by rows, vec⁡(X)\vecop(\bX) is X.reshape(-1, order="F"); flattening by rows instead gives the identity with A⊗B⊤\bA \otimes \bB^\top in place of B⊤⊗A\bB^\top \otimes \bA.

验证如下:AXB\bA\bX\bB 的第 kk 列为 ∑jbjkAxj\sum_j b_{jk}\bA\bx_j,其中 xj\bx_j 是 X\bX 的第 jj 列;这正是 (B⊤⊗A)vec⁡(X)(\bB^\top \otimes \bA)\vecop(\bX) 的第 kk 块,因为 B⊤⊗A\bB^\top \otimes \bA 的第 (k,j)(k, j) 块为 bjkAb_{jk}\bA。特例有 vec⁡(ab⊤)=b⊗a\vecop(\ba\bb^\top) = \bb \otimes \ba 与 tr⁡(A⊤B)=vec⁡(A)⊤vec⁡(B)\tr(\bA^\top\bB) = \vecop(\bA)^\top\vecop(\bB)。NumPy 按行存储数组,vec⁡(X)\vecop(\bX) 对应 X.reshape(-1, order="F");若按行展平,恒等式中的 B⊤⊗A\bB^\top \otimes \bA 要换成 A⊗B⊤\bA \otimes \bB^\top。

Tensors. An order-kk tensor is a kk-way array T=(ti1…ik)\mathcal{T} = (t_{i_1 \dots i_k}). The outer product of vectors a,b,c,d\ba, \bb, \mathbf{c}, \mathbf{d} is the order-4 tensor a∘b∘c∘d\ba \circ \bb \circ \mathbf{c} \circ \mathbf{d} with entries aibjckdla_ib_jc_kd_l; tensors of this form have rank one. A tensor is symmetric if its entries are invariant under permutations of the indices. A CP (canonical polyadic) decomposition writes T=∑r=1Rλr ar∘br∘cr∘dr\mathcal{T} = \sum_{r=1}^R \lambda_r\,\ba_r \circ \bb_r \circ \mathbf{c}_r \circ \mathbf{d}_r, and the smallest such RR is the CP rank. Tensor rank behaves differently from matrix rank. It can exceed every dimension of the array, computing it is NP-hard in general (Hillar & Lim, 2013), and the CP decomposition is often essentially unique, whereas a low-rank matrix factorization is unique only up to an invertible linear transformation (Kolda & Bader, 2009). Uniqueness conditions, such as Kruskal's, are discussed in Section 12.4第 12.4 节.

张量。kk 阶张量是一个 kk 维数组 T=(ti1…ik)\mathcal{T} = (t_{i_1 \dots i_k})。向量 a,b,c,d\ba, \bb, \mathbf{c}, \mathbf{d} 的外积是四阶张量 a∘b∘c∘d\ba \circ \bb \circ \mathbf{c} \circ \mathbf{d},其元素为 aibjckdla_ib_jc_kd_l;这种形式的张量秩为一。若张量的元素在下标的任意置换下不变,则称它是对称的。CP(canonical polyadic)分解把张量写成 T=∑r=1Rλr ar∘br∘cr∘dr\mathcal{T} = \sum_{r=1}^R \lambda_r\,\ba_r \circ \bb_r \circ \mathbf{c}_r \circ \mathbf{d}_r,这样的 RR 的最小值称为 CP 秩。张量秩与矩阵秩的性质很不一样:它可以超过数组的每一维,一般情形下计算它是 NP 难的(Hillar & Lim, 2013),而且 CP 分解常常在本质上唯一,而低秩矩阵分解只在相差一个可逆线性变换的意义下唯一(Kolda & Bader, 2009)。Kruskal 条件等唯一性条件见Section 12.4第 12.4 节。

Multilinear transformation. For a matrix M\bM with nn columns and an order-4 tensor T\mathcal{T} with all dimensions equal to nn, define

多重线性变换。对有 nn 列的矩阵 M\bM 与各维都等于 nn 的四阶张量 T\mathcal{T},定义

(T×1M×2M×3M×4M)ijkl=∑p,q,r,tmip mjq mkr mlt tpqrt.(\mathcal{T} \times_1 \bM \times_2 \bM \times_3 \bM \times_4 \bM)_{ijkl} = \sum_{p,q,r,t} m_{ip}\,m_{jq}\,m_{kr}\,m_{lt}\,t_{pqrt}.

It maps the rank-one tensor a∘a∘a∘a\ba \circ \ba \circ \ba \circ \ba to (Ma)∘(Ma)∘(Ma)∘(Ma)(\bM\ba) \circ (\bM\ba) \circ (\bM\ba) \circ (\bM\ba), and for orthogonal M\bM it preserves the Frobenius norm ∥T∥F2=∑ijkltijkl2\|\mathcal{T}\|_F^2 = \sum_{ijkl}t_{ijkl}^2 (Exercise A.3习题 A.3). The fourth-order cumulant tensor Kx\mathcal{K}_{\bx} of a random vector, with entries cum⁡(xi,xj,xk,xl)\cum(x_i, x_j, x_k, x_l), is symmetric, and by multilinearity (Proposition 2.8命题 2.8) KMx=Kx×1M×2M×3M×4M\mathcal{K}_{\bM\bx} = \mathcal{K}_{\bx} \times_1 \bM \times_2 \bM \times_3 \bM \times_4 \bM. For independent sources, Ks\mathcal{K}_{\bs} is diagonal, its only nonzero entries being (Ks)pppp=κp(\mathcal{K}_{\bs})_{pppp} = \kappa_p. Hence, for x=As\bx = \bA\bs,

它把秩一张量 a∘a∘a∘a\ba \circ \ba \circ \ba \circ \ba 映为 (Ma)∘(Ma)∘(Ma)∘(Ma)(\bM\ba) \circ (\bM\ba) \circ (\bM\ba) \circ (\bM\ba);当 M\bM 正交时,它保持 Frobenius 范数 ∥T∥F2=∑ijkltijkl2\|\mathcal{T}\|_F^2 = \sum_{ijkl}t_{ijkl}^2 不变(Exercise A.3习题 A.3)。随机向量的四阶累积量张量 Kx\mathcal{K}_{\bx} 以 cum⁡(xi,xj,xk,xl)\cum(x_i, x_j, x_k, x_l) 为元素,它是对称张量;由多重线性(Proposition 2.8命题 2.8),KMx=Kx×1M×2M×3M×4M\mathcal{K}_{\bM\bx} = \mathcal{K}_{\bx} \times_1 \bM \times_2 \bM \times_3 \bM \times_4 \bM。对独立的源,Ks\mathcal{K}_{\bs} 是对角张量,唯一的非零元为 (Ks)pppp=κp(\mathcal{K}_{\bs})_{pppp} = \kappa_p。因此对 x=As\bx = \bA\bs,

Kx=∑p=1nκp ap∘ap∘ap∘ap,(A.3)\mathcal{K}_{\bx} = \sum_{p=1}^n \kappa_p\,\ba_p \circ \ba_p \circ \ba_p \circ \ba_p,\tag{A.3}

a symmetric CP decomposition of rank at most nn whose factors are the columns of A\bA; this is (2.3) in tensor form. Arranging the entries as an m2×m2m^2 \times m^2 matrix, with rows indexed by (i,j)(i, j) and columns by (k,l)(k, l), gives ∑pκp vec⁡(apap⊤)vec⁡(apap⊤)⊤\sum_p \kappa_p\,\vecop(\ba_p\ba_p^\top)\vecop(\ba_p\ba_p^\top)^\top, a symmetric matrix of rank at most nn. For whitened data the factors are orthogonal, which JADE exploits (Chapter 9第 9 章); with more sources than sensors the rank exceeds the dimension, which only tensors allow (Chapter 12第 12 章).

这是一个秩至多为 nn 的对称 CP 分解,其因子就是 A\bA 的各列;它正是 (2.3) 的张量形式。把各元素排成 m2×m2m^2 \times m^2 的矩阵,行以 (i,j)(i, j) 为指标、列以 (k,l)(k, l) 为指标,得到 ∑pκp vec⁡(apap⊤)vec⁡(apap⊤)⊤\sum_p \kappa_p\,\vecop(\ba_p\ba_p^\top)\vecop(\ba_p\ba_p^\top)^\top,这是秩至多为 nn 的对称矩阵。对白化后的数据,这些因子相互正交,JADE 正是利用了这一点(Chapter 9第 9 章);当源多于传感器时,秩超过了维数,而这只有张量才允许(Chapter 12第 12 章)。

A.6 Matrix calculus矩阵微积分

For a scalar function ff of a matrix W∈Rn×m\bW \in \R^{n \times m}, the gradient ∂f/∂W\partial f/\partial\bW is the n×mn \times m matrix of partial derivatives ∂f/∂wij\partial f/\partial w_{ij}. The most reliable way to compute it is through the differential: if df=tr⁡(G⊤dW)df = \tr(\bG^\top d\bW) for every direction dWd\bW, then ∂f/∂W=G\partial f/\partial\bW = \bG (Magnus & Neudecker, 2019). The basic rules are

对矩阵 W∈Rn×m\bW \in \R^{n \times m} 的标量函数 ff,梯度 ∂f/∂W\partial f/\partial\bW 是由偏导数 ∂f/∂wij\partial f/\partial w_{ij} 构成的 n×mn \times m 矩阵。最可靠的计算方法是通过微分:若对每个方向 dWd\bW 都有 df=tr⁡(G⊤dW)df = \tr(\bG^\top d\bW),则 ∂f/∂W=G\partial f/\partial\bW = \bG(Magnus & Neudecker, 2019)。基本规则如下:

d tr⁡(AW)=tr⁡(A dW),∂ tr⁡(AW)∂W=A⊤,d(a⊤Wb)=tr⁡(ba⊤dW),∂ (a⊤Wb)∂W=ab⊤,d tr⁡(W⊤MW)=tr⁡(W⊤(M+M⊤) dW),∂ tr⁡(W⊤MW)∂W=(M+M⊤)W,d(W−1)=−W−1(dW) W−1.\begin{aligned} d\,\tr(\bA\bW) &= \tr(\bA\,d\bW), & \frac{\partial\,\tr(\bA\bW)}{\partial\bW} &= \bA^\top, \\[6pt] d(\ba^\top\bW\bb) &= \tr(\bb\ba^\top d\bW), & \frac{\partial\,(\ba^\top\bW\bb)}{\partial\bW} &= \ba\bb^\top, \\[6pt] d\,\tr(\bW^\top\bM\bW) &= \tr\bigl(\bW^\top(\bM + \bM^\top)\,d\bW\bigr), & \frac{\partial\,\tr(\bW^\top\bM\bW)}{\partial\bW} &= (\bM + \bM^\top)\bW, \\[6pt] d(\bW^{-1}) &= -\bW^{-1}(d\bW)\,\bW^{-1}. \end{aligned}

The last rule follows by differentiating WW−1=I\bW\bW^{-1} = \bI. The determinant, which brings the Jacobian of y=Wx\by = \bW\bx into every likelihood and entropy of this book, has a particularly simple gradient.

最后一条由对 WW−1=I\bW\bW^{-1} = \bI 求微分得到。行列式把变换 y=Wx\by = \bW\bx 的 Jacobi 因子带进了本书中的每个似然与熵,而它的梯度格外简单。

Proposition A.5 (gradient of the log-determinant)命题 A.5(对数行列式的梯度)

For invertible W\bW, dlog⁡∣det⁡W∣=tr⁡(W−1dW)d\log|\det\bW| = \tr(\bW^{-1}d\bW), hence ∂log⁡∣det⁡W∣/∂W=W−⊤\partial\log|\det\bW|/\partial\bW = \bW^{-\top}.

对可逆的 W\bW,dlog⁡∣det⁡W∣=tr⁡(W−1dW)d\log|\det\bW| = \tr(\bW^{-1}d\bW),因此 ∂log⁡∣det⁡W∣/∂W=W−⊤\partial\log|\det\bW|/\partial\bW = \bW^{-\top}。

Proof.证明

det⁡(W+dW)=det⁡W det⁡(I+W−1dW)\det(\bW + d\bW) = \det\bW\,\det(\bI + \bW^{-1}d\bW), and det⁡(I+Δ)=1+tr⁡Δ+O(∥Δ∥2)\det(\bI + \boldsymbol{\Delta}) = 1 + \tr\boldsymbol{\Delta} + O(\|\boldsymbol{\Delta}\|^2): in the expansion of the determinant, the product of the diagonal entries contributes 1+tr⁡Δ+O(∥Δ∥2)1 + \tr\boldsymbol{\Delta} + O(\|\boldsymbol{\Delta}\|^2), and every other term contains at least two off-diagonal entries of Δ\boldsymbol{\Delta}. Hence log⁡∣det⁡(W+dW)∣=log⁡∣det⁡W∣+tr⁡(W−1dW)+O(∥dW∥2)\log|\det(\bW + d\bW)| = \log|\det\bW| + \tr(\bW^{-1}d\bW) + O(\|d\bW\|^2). Writing tr⁡(W−1dW)=tr⁡((W−⊤)⊤dW)\tr(\bW^{-1}d\bW) = \tr\bigl((\bW^{-\top})^\top d\bW\bigr) identifies the gradient.

det⁡(W+dW)=det⁡W det⁡(I+W−1dW)\det(\bW + d\bW) = \det\bW\,\det(\bI + \bW^{-1}d\bW),而 det⁡(I+Δ)=1+tr⁡Δ+O(∥Δ∥2)\det(\bI + \boldsymbol{\Delta}) = 1 + \tr\boldsymbol{\Delta} + O(\|\boldsymbol{\Delta}\|^2):在行列式的展开式中,对角元之积贡献 1+tr⁡Δ+O(∥Δ∥2)1 + \tr\boldsymbol{\Delta} + O(\|\boldsymbol{\Delta}\|^2),其余每一项都至少含有 Δ\boldsymbol{\Delta} 的两个非对角元。因此 log⁡∣det⁡(W+dW)∣=log⁡∣det⁡W∣+tr⁡(W−1dW)+O(∥dW∥2)\log|\det(\bW + d\bW)| = \log|\det\bW| + \tr(\bW^{-1}d\bW) + O(\|d\bW\|^2)。把 tr⁡(W−1dW)\tr(\bW^{-1}d\bW) 写成 tr⁡((W−⊤)⊤dW)\tr\bigl((\bW^{-\top})^\top d\bW\bigr),即得梯度。

Example A.6 (gradient of the ICA log-likelihood)例 A.6(ICA 对数似然的梯度)

Let f(W)=log⁡∣det⁡W∣+∑iE[log⁡pi(wi⊤x)]f(\bW) = \log|\det\bW| + \sum_i \E[\log p_i(\bw_i^\top\bx)], where wi⊤\bw_i^\top is row ii of W\bW, y=Wx\by = \bW\bx, and ψi=−pi′/pi\psi_i = -p_i'/p_i is the score function of source ii. Then df=tr⁡(W−1dW)−E[∑iψi(yi) dwi⊤x]=tr⁡(W−1dW)−tr⁡(E[xψ(y)⊤] dW)df = \tr(\bW^{-1}d\bW) - \E\bigl[\sum_i \psi_i(y_i)\,d\bw_i^\top\bx\bigr] = \tr(\bW^{-1}d\bW) - \tr\bigl(\E[\bx\bpsi(\by)^\top]\,d\bW\bigr), so

设 f(W)=log⁡∣det⁡W∣+∑iE[log⁡pi(wi⊤x)]f(\bW) = \log|\det\bW| + \sum_i \E[\log p_i(\bw_i^\top\bx)],其中 wi⊤\bw_i^\top 是 W\bW 的第 ii 行,y=Wx\by = \bW\bx,ψi=−pi′/pi\psi_i = -p_i'/p_i 是源 ii 的得分函数。则 df=tr⁡(W−1dW)−E[∑iψi(yi) dwi⊤x]=tr⁡(W−1dW)−tr⁡(E[xψ(y)⊤] dW)df = \tr(\bW^{-1}d\bW) - \E\bigl[\sum_i \psi_i(y_i)\,d\bw_i^\top\bx\bigr] = \tr(\bW^{-1}d\bW) - \tr\bigl(\E[\bx\bpsi(\by)^\top]\,d\bW\bigr),所以

∂f∂W=W−⊤−E[ψ(y) x⊤],∂f∂W W⊤W=(I−E[ψ(y) y⊤])W.\frac{\partial f}{\partial\bW} = \bW^{-\top} - \E\bigl[\bpsi(\by)\,\bx^\top\bigr], \qquad \frac{\partial f}{\partial\bW}\,\bW^\top\bW = \bigl(\bI - \E[\bpsi(\by)\,\by^\top]\bigr)\bW.

The second expression, the gradient multiplied on the right by the positive definite matrix W⊤W\bW^\top\bW, is the relative or natural gradient of Chapter 7第 7 章; it is still an ascent direction, because ⟨G,GW⊤W⟩=∥GW⊤∥F2≥0\langle\bG, \bG\bW^\top\bW\rangle = \|\bG\bW^\top\|_F^2 \ge 0 for G=∂f/∂W\bG = \partial f/\partial\bW.

第二个表达式,即右乘正定矩阵 W⊤W\bW^\top\bW 之后的梯度,就是Chapter 7第 7 章中的相对梯度(自然梯度);它仍是上升方向,因为对 G=∂f/∂W\bG = \partial f/\partial\bW 有 ⟨G,GW⊤W⟩=∥GW⊤∥F2≥0\langle\bG, \bG\bW^\top\bW\rangle = \|\bG\bW^\top\|_F^2 \ge 0。

Formulas like these are easy to get wrong by a transpose, and easy to check. The script below verifies (A.2), Proposition A.5命题 A.5 by central differences, and Proposition A.4命题 A.4 against 20,000 random orthogonal matrices drawn uniformly from O(4)O(4).

这类公式很容易在转置上出错,也很容易检验。下面的脚本验证 (A.2),用中心差分验证Proposition A.5命题 A.5,并用从 O(4)O(4) 中均匀抽取的 20000 个随机正交矩阵检验Proposition A.4命题 A.4。

import numpy as np

rng = np.random.default_rng(0)

# vec(AXB) = (B^T kron A) vec(X), with vec stacking columns
A, X, B = (rng.standard_normal(s) for s in [(3, 4), (4, 5), (5, 2)])
vec = lambda M: M.reshape(-1, order="F")
print(np.allclose(vec(A @ X @ B), np.kron(B.T, A) @ vec(X)))

# gradient of log|det W| by central differences
W, h = rng.standard_normal((4, 4)), 1e-6
logdet = lambda M: np.linalg.slogdet(M)[1]
G = np.zeros((4, 4))
for i in range(4):
    for j in range(4):
        D = np.zeros((4, 4)); D[i, j] = h
        G[i, j] = (logdet(W + D) - logdet(W - D)) / (2 * h)
print(f"{np.abs(G - np.linalg.inv(W).T).max():.1e}")

# closest orthogonal matrix: L R^T = (W W^T)^(-1/2) W
L, sv, Rt = np.linalg.svd(W)
d, E = np.linalg.eigh(W @ W.T)
print(np.allclose(L @ Rt, E @ np.diag(d ** -0.5) @ E.T @ W))
def random_orthogonal(n):                 # uniformly distributed on O(n)
    Q, R = np.linalg.qr(rng.standard_normal((n, n)))
    return Q * np.sign(np.diag(R))
dist = min(np.linalg.norm(W - random_orthogonal(4)) for _ in range(20000))
print(f"{np.linalg.norm(W - L @ Rt):.3f} {dist:.3f}")

It prints True, a finite-difference error of 5.7e-08, True, and the distances 2.886 from W\bW to its polar factor and 3.047 to the best of the random orthogonal matrices.

程序依次输出 True、有限差分误差 5.7e-08、True,以及 W\bW 到其极因子的距离 2.886 和到随机正交矩阵中最近者的距离 3.047。

A.7 Optimization on the sphere and the orthogonal group球面与正交群上的优化

Many ICA contrasts are optimized under constraints. One-unit methods search over unit vectors w\bw (Chapter 5第 5 章, Chapter 8第 8 章), and methods for whitened data search over orthogonal matrices. Both constraint sets are smooth manifolds, and the same three ingredients recur: first-order (Lagrange) conditions, the Riemannian gradient, and retractions, which map a step taken in the tangent space back onto the manifold (Edelman, Arias & Smith, 1998; Absil, Mahony & Sepulchre, 2008).

许多 ICA 对比函数都是在约束下优化的。单单元方法在单位向量 w\bw 中搜索(Chapter 5第 5 章,Chapter 8第 8 章),针对白化数据的方法则在正交矩阵中搜索。这两个约束集都是光滑流形,反复出现的是同样三样东西:一阶(拉格朗日)条件、黎曼梯度,以及把在切空间中迈出的一步拉回流形的收缩映射(retraction)(Edelman, Arias & Smith, 1998; Absil, Mahony & Sepulchre, 2008)。

Proposition A.7 (optimality conditions on the sphere)命题 A.7(球面上的最优性条件)

Let ff be twice continuously differentiable on Rn\R^n, let w\bw be a unit vector, and write λ=w⊤∇f(w)\lambda = \bw^\top\nabla f(\bw). (a) If w\bw is a local maximum or minimum of ff on the unit sphere, then ∇f(w)=λw\nabla f(\bw) = \lambda\bw; equivalently, the Riemannian gradient (I−ww⊤)∇f(w)(\bI - \bw\bw^\top)\nabla f(\bw) vanishes. (b) At such a critical point, if v⊤(∇2f(w)−λI)v<0\bv^\top(\nabla^2 f(\bw) - \lambda\bI)\bv < 0 for every nonzero v\bv with w⊤v=0\bw^\top\bv = 0, then w\bw is a strict local maximum on the sphere; if v⊤(∇2f(w)−λI)v>0\bv^\top(\nabla^2 f(\bw) - \lambda\bI)\bv > 0 for some such v\bv, it is not a local maximum.

设 ff 在 Rn\R^n 上二阶连续可微,w\bw 是单位向量,记 λ=w⊤∇f(w)\lambda = \bw^\top\nabla f(\bw)。(a)若 w\bw 是 ff 在单位球面上的局部极大点或极小点,则 ∇f(w)=λw\nabla f(\bw) = \lambda\bw;等价地,黎曼梯度 (I−ww⊤)∇f(w)(\bI - \bw\bw^\top)\nabla f(\bw) 为零。(b)在这样的临界点处,若对每个满足 w⊤v=0\bw^\top\bv = 0 的非零 v\bv 都有 v⊤(∇2f(w)−λI)v<0\bv^\top(\nabla^2 f(\bw) - \lambda\bI)\bv < 0,则 w\bw 是球面上的严格局部极大点;若对某个这样的 v\bv 有 v⊤(∇2f(w)−λI)v>0\bv^\top(\nabla^2 f(\bw) - \lambda\bI)\bv > 0,则它不是局部极大点。

Proof.证明

For a unit vector v⊥w\bv \perp \bw, the curve w(t)=cos⁡t w+sin⁡t v\bw(t) = \cos t\,\bw + \sin t\,\bv stays on the sphere, with w(0)=w\bw(0) = \bw, w′(0)=v\bw'(0) = \bv and w′′(0)=−w\bw''(0) = -\bw. Hence the first derivative of f(w(t))f(\bw(t)) at t=0t = 0 is ∇f(w)⊤v\nabla f(\bw)^\top\bv, and the second is v⊤∇2f(w)v−∇f(w)⊤w\bv^\top\nabla^2 f(\bw)\bv - \nabla f(\bw)^\top\bw. At a local extremum the first derivative vanishes for every tangent v\bv, so ∇f(w)\nabla f(\bw) is parallel to w\bw, which is (a). At a critical point the second derivative equals v⊤(∇2f(w)−λI)v\bv^\top(\nabla^2 f(\bw) - \lambda\bI)\bv. If it is positive for some v\bv, ff increases along that curve, so w\bw is not a local maximum. Every point of the sphere near w\bw has the form w(t)\bw(t) for a small tt and some unit v⊥w\bv \perp \bw, and by Taylor's theorem f(w(t))=f(w)+t22v⊤(∇2f(w)−λI)v+o(t2)f(\bw(t)) = f(\bw) + \tfrac{t^2}{2}\bv^\top(\nabla^2 f(\bw) - \lambda\bI)\bv + o(t^2) uniformly in v\bv, because the unit tangent vectors form a compact set; a negative definite form therefore gives f(w(t))<f(w)f(\bw(t)) < f(\bw) for small t≠0t \ne 0.

对与 w\bw 正交的单位向量 v\bv,曲线 w(t)=cos⁡t w+sin⁡t v\bw(t) = \cos t\,\bw + \sin t\,\bv 始终在球面上,且 w(0)=w\bw(0) = \bw,w′(0)=v\bw'(0) = \bv,w′′(0)=−w\bw''(0) = -\bw。因此 f(w(t))f(\bw(t)) 在 t=0t = 0 处的一阶导数为 ∇f(w)⊤v\nabla f(\bw)^\top\bv,二阶导数为 v⊤∇2f(w)v−∇f(w)⊤w\bv^\top\nabla^2 f(\bw)\bv - \nabla f(\bw)^\top\bw。在局部极值点,一阶导数对每个切向量 v\bv 都为零,所以 ∇f(w)\nabla f(\bw) 与 w\bw 平行,这就是(a)。在临界点处,二阶导数等于 v⊤(∇2f(w)−λI)v\bv^\top(\nabla^2 f(\bw) - \lambda\bI)\bv。若它对某个 v\bv 为正,则 ff 沿这条曲线增大,所以 w\bw 不是局部极大点。球面上靠近 w\bw 的每一点都可以写成 w(t)\bw(t) 的形式(tt 很小,v⊥w\bv \perp \bw 为单位向量);由于单位切向量构成紧集,由 Taylor 定理,f(w(t))=f(w)+t22v⊤(∇2f(w)−λI)v+o(t2)f(\bw(t)) = f(\bw) + \tfrac{t^2}{2}\bv^\top(\nabla^2 f(\bw) - \lambda\bI)\bv + o(t^2) 关于 v\bv 一致成立,因此当这个二次型负定时,对很小的 t≠0t \ne 0 有 f(w(t))<f(w)f(\bw(t)) < f(\bw)。

For example, f(w)=w⊤Cwf(\bw) = \bw^\top\bC\bw with C\bC symmetric has ∇f=2Cw\nabla f = 2\bC\bw and ∇2f=2C\nabla^2 f = 2\bC, so its critical points on the sphere are the unit eigenvectors, with λ=2λk\lambda = 2\lambda_k at ek\be_k. Suppose λ1>λ2≥⋯≥λn\lambda_1 > \lambda_2 \ge \dots \ge \lambda_n. At e1\be_1, the form 2v⊤(C−λ1I)v2\bv^\top(\bC - \lambda_1\bI)\bv is negative for every nonzero v⊥e1\bv \perp \be_1, so e1\be_1 is a strict local maximum. At ek\be_k with k≥2k \ge 2, the tangent direction v=e1\bv = \be_1 gives 2(λ1−λk)>02(\lambda_1 - \lambda_k) > 0, so ek\be_k is not a local maximum. The only local maxima are ±e1\pm\be_1: PCA has no spurious local optima. Chapter 5第 5 章 carries out the same analysis for kurtosis (Theorem 5.3定理 5.3).

举例来说,对对称矩阵 C\bC,f(w)=w⊤Cwf(\bw) = \bw^\top\bC\bw 有 ∇f=2Cw\nabla f = 2\bC\bw 与 ∇2f=2C\nabla^2 f = 2\bC,所以它在球面上的临界点就是单位特征向量,在 ek\be_k 处 λ=2λk\lambda = 2\lambda_k。设 λ1>λ2≥⋯≥λn\lambda_1 > \lambda_2 \ge \dots \ge \lambda_n。在 e1\be_1 处,二次型 2v⊤(C−λ1I)v2\bv^\top(\bC - \lambda_1\bI)\bv 对每个与 e1\be_1 正交的非零 v\bv 都为负,所以 e1\be_1 是严格局部极大点。在 ek\be_k(k≥2k \ge 2)处,取切方向 v=e1\bv = \be_1 得 2(λ1−λk)>02(\lambda_1 - \lambda_k) > 0,所以 ek\be_k 不是局部极大点。唯一的局部极大点是 ±e1\pm\be_1:PCA 没有虚假的局部最优。Chapter 5第 5 章对峭度作了同样的分析(Theorem 5.3定理 5.3)。

The vector (I−ww⊤)∇f(w)(\bI - \bw\bw^\top)\nabla f(\bw), the projection of the gradient onto the tangent space {v:w⊤v=0}\{\bv : \bw^\top\bv = 0\}, is the Riemannian gradient of ff on the sphere. Gradient ascent takes a step μv\mu\bv along it and maps back by normalization, w←(w+μv)/∥w+μv∥\bw \leftarrow (\bw + \mu\bv)/\|\bw + \mu\bv\|, the simplest retraction; the normalization step of FastICA has the same form.

向量 (I−ww⊤)∇f(w)(\bI - \bw\bw^\top)\nabla f(\bw) 是梯度在切空间 {v:w⊤v=0}\{\bv : \bw^\top\bv = 0\} 上的投影,称为 ff 在球面上的黎曼梯度。梯度上升沿它走一步 μv\mu\bv,再通过归一化拉回球面:w←(w+μv)/∥w+μv∥\bw \leftarrow (\bw + \mu\bv)/\|\bw + \mu\bv\|,这是最简单的收缩映射;FastICA 的归一化步骤也是这种形式。

The orthogonal group. Let ff be a smooth function of U∈Rn×n\bU \in \R^{n \times n} with Euclidean gradient G=∂f/∂U\bG = \partial f/\partial\bU. By Proposition A.1命题 A.1, the directions tangent to O(n)O(n) at U\bU are ΩU\boldsymbol{\Omega}\bU with Ω\boldsymbol{\Omega} skew-symmetric, and the derivative of ff along t↦exp⁡(tΩ)Ut \mapsto \exp(t\boldsymbol{\Omega})\bU at t=0t = 0 is ⟨G,ΩU⟩=⟨GU⊤,Ω⟩=⟨skew⁡(GU⊤),Ω⟩\langle\bG, \boldsymbol{\Omega}\bU\rangle = \langle\bG\bU^\top, \boldsymbol{\Omega}\rangle = \langle\operatorname{skew}(\bG\bU^\top), \boldsymbol{\Omega}\rangle, because the symmetric part of GU⊤\bG\bU^\top is orthogonal to Ω\boldsymbol{\Omega} (Section A.1第 A.1 节). Therefore

正交群。设 ff 是 U∈Rn×n\bU \in \R^{n \times n} 的光滑函数,欧氏梯度为 G=∂f/∂U\bG = \partial f/\partial\bU。由Proposition A.1命题 A.1,O(n)O(n) 在 U\bU 处的切方向为 ΩU\boldsymbol{\Omega}\bU(Ω\boldsymbol{\Omega} 反对称),ff 沿 t↦exp⁡(tΩ)Ut \mapsto \exp(t\boldsymbol{\Omega})\bU 在 t=0t = 0 处的导数为 ⟨G,ΩU⟩=⟨GU⊤,Ω⟩=⟨skew⁡(GU⊤),Ω⟩\langle\bG, \boldsymbol{\Omega}\bU\rangle = \langle\bG\bU^\top, \boldsymbol{\Omega}\rangle = \langle\operatorname{skew}(\bG\bU^\top), \boldsymbol{\Omega}\rangle,因为 GU⊤\bG\bU^\top 的对称部分与 Ω\boldsymbol{\Omega} 正交(Section A.1第 A.1 节)。因此

grad⁡f(U)=skew⁡(GU⊤) U(A.4)\operatorname{grad} f(\bU) = \operatorname{skew}(\bG\bU^\top)\,\bU\tag{A.4}

is the Riemannian gradient, and U\bU is a critical point of ff on O(n)O(n) if and only if GU⊤\bG\bU^\top is symmetric. This is the Lagrange condition for the constraint UU⊤=I\bU\bU^\top = \bI, whose multiplier is a symmetric matrix Λ\bLambda with G=ΛU\bG = \bLambda\bU. Retractions include the matrix exponential U←exp⁡(μΩ)U\bU \leftarrow \exp(\mu\boldsymbol{\Omega})\bU, the Cayley transform U←(I−μ2Ω)−1(I+μ2Ω)U\bU \leftarrow (\bI - \tfrac{\mu}{2}\boldsymbol{\Omega})^{-1}(\bI + \tfrac{\mu}{2}\boldsymbol{\Omega})\bU (Exercise A.1习题 A.1), and the polar factor of U+μΩU\bU + \mu\boldsymbol{\Omega}\bU, that is, symmetric orthogonalization (Section A.4第 A.4 节). In two dimensions each of them simply changes the rotation angle.

就是黎曼梯度;U\bU 是 ff 在 O(n)O(n) 上的临界点,当且仅当 GU⊤\bG\bU^\top 对称。这正是约束 UU⊤=I\bU\bU^\top = \bI 下的拉格朗日条件,其乘子是满足 G=ΛU\bG = \bLambda\bU 的对称矩阵 Λ\bLambda。常用的收缩映射有矩阵指数 U←exp⁡(μΩ)U\bU \leftarrow \exp(\mu\boldsymbol{\Omega})\bU、Cayley 变换 U←(I−μ2Ω)−1(I+μ2Ω)U\bU \leftarrow (\bI - \tfrac{\mu}{2}\boldsymbol{\Omega})^{-1}(\bI + \tfrac{\mu}{2}\boldsymbol{\Omega})\bU(Exercise A.1习题 A.1),以及 U+μΩU\bU + \mu\boldsymbol{\Omega}\bU 的极因子,即对称正交化(Section A.4第 A.4 节)。在二维情形,它们都只是改变旋转角。

A.8 Asymptotic statistics渐近统计

Part IV and several other chapters make statements about large samples. We write XT=OP(rT)X_T = O_P(r_T) if XT/rTX_T/r_T is bounded in probability, that is, for every ε>0\varepsilon > 0 there is an MM with P(∣XT/rT∣>M)<εP(|X_T/r_T| > M) < \varepsilon for all TT, and XT=oP(rT)X_T = o_P(r_T) if XT/rT→0X_T/r_T \to 0 in probability. By the central limit theorem, the average of TT i.i.d. random vectors with mean μ\bmu and covariance Σ\bSigma satisfies T(xˉ−μ)→N(0,Σ)\sqrt T(\bar\bx - \bmu) \to \Ncal(\bzero, \bSigma) in distribution, so xˉ−μ=OP(T−1/2)\bar\bx - \bmu = O_P(T^{-1/2}). Slutsky's lemma allows us to add oP(1)o_P(1) terms to a sequence that converges in distribution, or to multiply it by a sequence that converges in probability to a constant, without changing the limit beyond the obvious (van der Vaart, 1998).

第四部分以及其他若干章都涉及大样本下的结论。若 XT/rTX_T/r_T 依概率有界,即对任意 ε>0\varepsilon > 0,存在 MM 使对所有 TT 都有 P(∣XT/rT∣>M)<εP(|X_T/r_T| > M) < \varepsilon,就记 XT=OP(rT)X_T = O_P(r_T);若 XT/rTX_T/r_T 依概率收敛到 0,就记 XT=oP(rT)X_T = o_P(r_T)。由中心极限定理,TT 个均值为 μ\bmu、协方差为 Σ\bSigma 的独立同分布随机向量的平均满足 T(xˉ−μ)→N(0,Σ)\sqrt T(\bar\bx - \bmu) \to \Ncal(\bzero, \bSigma)(依分布),所以 xˉ−μ=OP(T−1/2)\bar\bx - \bmu = O_P(T^{-1/2})。Slutsky 引理允许我们给依分布收敛的序列加上 oP(1)o_P(1) 项,或乘以依概率收敛到常数的序列,而极限只发生显然的变化(van der Vaart, 1998)。

Proposition A.8 (delta method)命题 A.8(delta 方法)

Let T(θ^T−θ)→N(0,Σ)\sqrt T(\hat\btheta_T - \btheta) \to \Ncal(\bzero, \bSigma) in distribution, and let gg be a real function differentiable at θ\btheta. Then T(g(θ^T)−g(θ))→N(0, ∇g(θ)⊤Σ ∇g(θ))\sqrt T\bigl(g(\hat\btheta_T) - g(\btheta)\bigr) \to \Ncal\bigl(0,\ \nabla g(\btheta)^\top\bSigma\,\nabla g(\btheta)\bigr). For a vector-valued gg, replace ∇g(θ)⊤\nabla g(\btheta)^\top by the Jacobian matrix of gg at θ\btheta.

设 T(θ^T−θ)\sqrt T(\hat\btheta_T - \btheta) 依分布收敛到 N(0,Σ)\Ncal(\bzero, \bSigma),gg 是在 θ\btheta 处可微的实函数。则 T(g(θ^T)−g(θ))→N(0, ∇g(θ)⊤Σ ∇g(θ))\sqrt T\bigl(g(\hat\btheta_T) - g(\btheta)\bigr) \to \Ncal\bigl(0,\ \nabla g(\btheta)^\top\bSigma\,\nabla g(\btheta)\bigr)。若 gg 取向量值,把 ∇g(θ)⊤\nabla g(\btheta)^\top 换成 gg 在 θ\btheta 处的 Jacobi 矩阵即可。

Proof.证明

Differentiability gives g(θ^T)−g(θ)=∇g(θ)⊤(θ^T−θ)+oP(∥θ^T−θ∥)g(\hat\btheta_T) - g(\btheta) = \nabla g(\btheta)^\top(\hat\btheta_T - \btheta) + o_P(\|\hat\btheta_T - \btheta\|), and ∥θ^T−θ∥=OP(T−1/2)\|\hat\btheta_T - \btheta\| = O_P(T^{-1/2}). Multiply by T\sqrt T: the remainder becomes oP(1)o_P(1), and the linear term converges to N(0,∇g⊤Σ∇g)\Ncal(0, \nabla g^\top\bSigma\nabla g). Slutsky's lemma completes the proof.

由可微性,g(θ^T)−g(θ)=∇g(θ)⊤(θ^T−θ)+oP(∥θ^T−θ∥)g(\hat\btheta_T) - g(\btheta) = \nabla g(\btheta)^\top(\hat\btheta_T - \btheta) + o_P(\|\hat\btheta_T - \btheta\|),且 ∥θ^T−θ∥=OP(T−1/2)\|\hat\btheta_T - \btheta\| = O_P(T^{-1/2})。乘以 T\sqrt T 后,余项变为 oP(1)o_P(1),线性项收敛到 N(0,∇g⊤Σ∇g)\Ncal(0, \nabla g^\top\bSigma\nabla g)。由 Slutsky 引理即得结论。

Example A.9 (variance of the sample kurtosis)例 A.9(样本峭度的方差)

Let yy have zero mean and moments μj=E[yj]\mu_j = \E[y^j] up to order 8, and let k^=m^4−3m^22\hat k = \hat m_4 - 3\hat m_2^2 with m^j=1T∑ty(t)j\hat m_j = \frac1T\sum_t y(t)^j. By the central limit theorem, T(m^2−μ2,m^4−μ4)\sqrt T(\hat m_2 - \mu_2, \hat m_4 - \mu_4) is asymptotically normal with variances μ4−μ22\mu_4 - \mu_2^2 and μ8−μ42\mu_8 - \mu_4^2 and covariance μ6−μ2μ4\mu_6 - \mu_2\mu_4, and g(a,b)=b−3a2g(a, b) = b - 3a^2 has gradient (−6μ2,1)(-6\mu_2, 1). The delta method gives T(k^−kurt⁡(y))→N(0,v)\sqrt T(\hat k - \kurt(y)) \to \Ncal(0, v) with

设 yy 均值为零,直到 8 阶的矩 μj=E[yj]\mu_j = \E[y^j] 都有限,令 k^=m^4−3m^22\hat k = \hat m_4 - 3\hat m_2^2,其中 m^j=1T∑ty(t)j\hat m_j = \frac1T\sum_t y(t)^j。由中心极限定理,T(m^2−μ2,m^4−μ4)\sqrt T(\hat m_2 - \mu_2, \hat m_4 - \mu_4) 渐近正态,两个分量的方差分别为 μ4−μ22\mu_4 - \mu_2^2 与 μ8−μ42\mu_8 - \mu_4^2,协方差为 μ6−μ2μ4\mu_6 - \mu_2\mu_4;而 g(a,b)=b−3a2g(a, b) = b - 3a^2 的梯度为 (−6μ2,1)(-6\mu_2, 1)。由 delta 方法,T(k^−kurt⁡(y))→N(0,v)\sqrt T(\hat k - \kurt(y)) \to \Ncal(0, v),其中

v=μ8−μ42−12 μ2(μ6−μ2μ4)+36 μ22(μ4−μ22).(A.5)v = \mu_8 - \mu_4^2 - 12\,\mu_2(\mu_6 - \mu_2\mu_4) + 36\,\mu_2^2(\mu_4 - \mu_2^2).\tag{A.5}

For a standard Gaussian, (μ2,μ4,μ6,μ8)=(1,3,15,105)(\mu_2, \mu_4, \mu_6, \mu_8) = (1, 3, 15, 105) and v=24v = 24, the 24/T\sqrt{24/T} used in Exercise 1.5习题 1.5. For the unit-variance Laplace law, (μ2,μ4,μ6,μ8)=(1,6,90,2520)(\mu_2, \mu_4, \mu_6, \mu_8) = (1, 6, 90, 2520) and v=1656v = 1656: kurtosis estimates of super-Gaussian sources are far noisier (Section 5.6第 5.6 节). For symmetric laws, estimating the mean as well does not change vv: to first order, centering adds −4μ3yˉ-4\mu_3\bar y to k^\hat k, which changes vv by 64μ2μ32−8μ3μ564\mu_2\mu_3^2 - 8\mu_3\mu_5, and μ3=0\mu_3 = 0 for symmetric laws. A simulation with T=2000T = 2000 and 20 00020\,000 repetitions, subtracting the sample mean, gives TVar⁡(k^)≈23.5T\Var(\hat k) \approx 23.5, 9.859.85 and 16551655 for Gaussian, uniform and Laplace data, with Monte Carlo standard errors of about 0.250.25, 0.100.10 and 4040, against the predicted 2424, 9.879.87 (Exercise A.4习题 A.4) and 16561656.

对标准高斯分布,(μ2,μ4,μ6,μ8)=(1,3,15,105)(\mu_2, \mu_4, \mu_6, \mu_8) = (1, 3, 15, 105),v=24v = 24,这正是Exercise 1.5习题 1.5中用到的 24/T\sqrt{24/T}。对单位方差的拉普拉斯分布,(μ2,μ4,μ6,μ8)=(1,6,90,2520)(\mu_2, \mu_4, \mu_6, \mu_8) = (1, 6, 90, 2520),v=1656v = 1656:超高斯源的峭度估计要嘈杂得多(Section 5.6第 5.6 节)。对对称分布,即使同时估计均值,vv 也不变:在一阶近似下,中心化给 k^\hat k 加上 −4μ3yˉ-4\mu_3\bar y,使 vv 改变 64μ2μ32−8μ3μ564\mu_2\mu_3^2 - 8\mu_3\mu_5,而对称分布的 μ3=0\mu_3 = 0。用 T=2000T = 2000、重复 20 00020\,000 次并减去样本均值的模拟,对高斯、均匀与拉普拉斯数据得到 TVar⁡(k^)≈23.5T\Var(\hat k) \approx 23.5、9.859.85 与 16551655(蒙特卡罗标准误分别约为 0.250.25、0.100.10 与 4040),预言值分别为 2424、9.879.87(Exercise A.4习题 A.4)与 16561656。

M-estimators. Most ICA estimators are defined as solutions of estimating equations, such as E[φi(yi) yj]=0\E[\varphi_i(y_i)\,y_j] = 0 for i≠ji \ne j in Chapter 6第 6 章, with the expectation replaced by a sample average. Their asymptotic covariance has a standard form.

M 估计量。大多数 ICA 估计量都定义为估计方程的解,例如Chapter 6第 6 章中的 E[φi(yi) yj]=0\E[\varphi_i(y_i)\,y_j] = 0(i≠ji \ne j),其中期望换成样本平均。它们的渐近协方差有标准的形式。

Proposition A.10 (asymptotic normality of M-estimators)命题 A.10(M 估计量的渐近正态性)

Let x(1),…,x(T)\bx(1), \dots, \bx(T) be i.i.d., and let θ^T\hat\btheta_T solve 1T∑t=1Th(x(t),θ)=0\frac1T\sum_{t=1}^T \mathbf{h}(\bx(t), \btheta) = \bzero, where E[h(x,θ0)]=0\E[\mathbf{h}(\bx, \btheta_0)] = \bzero. Assume that θ^T→θ0\hat\btheta_T \to \btheta_0 in probability, that h\mathbf{h} is twice continuously differentiable in θ\btheta with first and second derivatives dominated near θ0\btheta_0 by an integrable function of x\bx, that H=E[∂h(x,θ0)/∂θ⊤]\bH = \E[\partial\mathbf{h}(\bx, \btheta_0)/\partial\btheta^\top] is invertible, and that B=E[h(x,θ0) h(x,θ0)⊤]\bB = \E[\mathbf{h}(\bx, \btheta_0)\,\mathbf{h}(\bx, \btheta_0)^\top] is finite. Then T(θ^T−θ0)→N(0, H−1BH−⊤)\sqrt T(\hat\btheta_T - \btheta_0) \to \Ncal(\bzero,\ \bH^{-1}\bB\bH^{-\top}).

设 x(1),…,x(T)\bx(1), \dots, \bx(T) 独立同分布,θ^T\hat\btheta_T 是方程 1T∑t=1Th(x(t),θ)=0\frac1T\sum_{t=1}^T \mathbf{h}(\bx(t), \btheta) = \bzero 的解,且 E[h(x,θ0)]=0\E[\mathbf{h}(\bx, \btheta_0)] = \bzero。假设 θ^T\hat\btheta_T 依概率收敛到 θ0\btheta_0;h\mathbf{h} 关于 θ\btheta 二阶连续可微,且其一阶与二阶导数在 θ0\btheta_0 附近被某个关于 x\bx 可积的函数控制;H=E[∂h(x,θ0)/∂θ⊤]\bH = \E[\partial\mathbf{h}(\bx, \btheta_0)/\partial\btheta^\top] 可逆;B=E[h(x,θ0) h(x,θ0)⊤]\bB = \E[\mathbf{h}(\bx, \btheta_0)\,\mathbf{h}(\bx, \btheta_0)^\top] 有限。则 T(θ^T−θ0)→N(0, H−1BH−⊤)\sqrt T(\hat\btheta_T - \btheta_0) \to \Ncal(\bzero,\ \bH^{-1}\bB\bH^{-\top})。

Proof.证明

(Sketch.) A Taylor expansion of the estimating equations around θ0\btheta_0 gives 0=1T∑th(x(t),θ0)+H^T(θ^T−θ0)\bzero = \frac1T\sum_t \mathbf{h}(\bx(t), \btheta_0) + \hat\bH_T(\hat\btheta_T - \btheta_0), where H^T\hat\bH_T is an average of derivatives evaluated at points between θ^T\hat\btheta_T and θ0\btheta_0; by the law of large numbers and the domination condition, H^T→H\hat\bH_T \to \bH in probability. Hence T(θ^T−θ0)=−H^T−1 T−1/2∑th(x(t),θ0)\sqrt T(\hat\btheta_T - \btheta_0) = -\hat\bH_T^{-1}\,T^{-1/2}\sum_t \mathbf{h}(\bx(t), \btheta_0) with probability tending to one, and the central limit theorem and Slutsky's lemma give the limit. Van der Vaart (1998, Chapter 5) gives the details and weaker conditions.

(证明概要。)把估计方程在 θ0\btheta_0 附近作 Taylor 展开,得 0=1T∑th(x(t),θ0)+H^T(θ^T−θ0)\bzero = \frac1T\sum_t \mathbf{h}(\bx(t), \btheta_0) + \hat\bH_T(\hat\btheta_T - \btheta_0),其中 H^T\hat\bH_T 是在 θ^T\hat\btheta_T 与 θ0\btheta_0 之间的点处取值的导数的平均;由大数定律与控制条件,H^T\hat\bH_T 依概率收敛到 H\bH。因此,在概率趋于一的事件上,T(θ^T−θ0)=−H^T−1 T−1/2∑th(x(t),θ0)\sqrt T(\hat\btheta_T - \btheta_0) = -\hat\bH_T^{-1}\,T^{-1/2}\sum_t \mathbf{h}(\bx(t), \btheta_0),再由中心极限定理与 Slutsky 引理即得极限。细节以及更弱的条件见 van der Vaart(1998)第 5 章。

The covariance H−1BH−⊤\bH^{-1}\bB\bH^{-\top} is called the sandwich. Chapter 10第 10 章 evaluates it for ICA estimators (Theorem 10.8定理 10.8), and Exercise A.5习题 A.5 works out the simplest case, the estimation of a location parameter.

协方差 H−1BH−⊤\bH^{-1}\bB\bH^{-\top} 称为三明治(sandwich)协方差。Chapter 10第 10 章对 ICA 估计量计算了它(Theorem 10.8定理 10.8),Exercise A.5习题 A.5则处理最简单的情形,即位置参数的估计。

Fisher information and the Cramér–Rao bound. For a family of densities p(x;θ)p(\bx; \btheta), the Fisher information is I(θ)=E[∇θlog⁡p(x;θ) ∇θlog⁡p(x;θ)⊤]\mathcal{I}(\btheta) = \E\bigl[\nabla_{\btheta}\log p(\bx; \btheta)\,\nabla_{\btheta}\log p(\bx; \btheta)^\top\bigr], which equals −E[∇θ2log⁡p(x;θ)]-\E[\nabla^2_{\btheta}\log p(\bx; \btheta)] under regularity conditions. The Cramér–Rao bound states that, under regularity conditions, every unbiased estimator based on TT i.i.d. observations satisfies Cov⁡(θ^)⪰(T I(θ))−1\Cov(\hat\btheta) \succeq (T\,\mathcal{I}(\btheta))^{-1}; the maximum-likelihood estimator is asymptotically normal with covariance I(θ)−1/T\mathcal{I}(\btheta)^{-1}/T, so it attains the bound asymptotically (Lehmann & Casella, 1998; van der Vaart, 1998). In the sandwich formula, maximum likelihood takes h\mathbf{h} to be the score, so that H=−I\bH = -\mathcal{I}, B=I\bB = \mathcal{I} and H−1BH−⊤=I−1\bH^{-1}\bB\bH^{-\top} = \mathcal{I}^{-1}. For ICA, Chapter 10第 10 章 derives the bound for the gain matrix (Theorem 10.10定理 10.10).

Fisher 信息与 Cramér–Rao 界。对密度族 p(x;θ)p(\bx; \btheta),Fisher 信息定义为 I(θ)=E[∇θlog⁡p(x;θ) ∇θlog⁡p(x;θ)⊤]\mathcal{I}(\btheta) = \E\bigl[\nabla_{\btheta}\log p(\bx; \btheta)\,\nabla_{\btheta}\log p(\bx; \btheta)^\top\bigr],在正则条件下它等于 −E[∇θ2log⁡p(x;θ)]-\E[\nabla^2_{\btheta}\log p(\bx; \btheta)]。Cramér–Rao 界是说:在正则条件下,基于 TT 个独立同分布观测的任何无偏估计量都满足 Cov⁡(θ^)⪰(T I(θ))−1\Cov(\hat\btheta) \succeq (T\,\mathcal{I}(\btheta))^{-1};最大似然估计量渐近正态,协方差为 I(θ)−1/T\mathcal{I}(\btheta)^{-1}/T,所以它渐近地达到这个界(Lehmann & Casella, 1998; van der Vaart, 1998)。在三明治公式中,最大似然取 h\mathbf{h} 为得分,于是 H=−I\bH = -\mathcal{I},B=I\bB = \mathcal{I},H−1BH−⊤=I−1\bH^{-1}\bB\bH^{-\top} = \mathcal{I}^{-1}。对 ICA,Chapter 10第 10 章推导了增益矩阵的这个界(Theorem 10.10定理 10.10)。

For a location family p(x−θ)p(x - \theta), the score with respect to θ\theta is −p′(x−θ)/p(x−θ)=ψ(x−θ)-p'(x - \theta)/p(x - \theta) = \psi(x - \theta), where ψ=−p′/p\psi = -p'/p is the score function in the book's notation, and the Fisher information is J=E[ψ(s)2]\mathcal{J} = \E[\psi(s)^2]. For the standard Gaussian, ψ(u)=u\psi(u) = u and J=1\mathcal{J} = 1; for the unit-variance Laplace law, ψ(u)=2 sign⁡(u)\psi(u) = \sqrt2\,\sign(u) and J=2\mathcal{J} = 2. Among laws of a given variance, the Gaussian carries the least information.

对位置族 p(x−θ)p(x - \theta),关于 θ\theta 的得分为 −p′(x−θ)/p(x−θ)=ψ(x−θ)-p'(x - \theta)/p(x - \theta) = \psi(x - \theta),其中 ψ=−p′/p\psi = -p'/p 即本书记号中的得分函数,Fisher 信息为 J=E[ψ(s)2]\mathcal{J} = \E[\psi(s)^2]。对标准高斯分布,ψ(u)=u\psi(u) = u,J=1\mathcal{J} = 1;对单位方差的拉普拉斯分布,ψ(u)=2 sign⁡(u)\psi(u) = \sqrt2\,\sign(u),J=2\mathcal{J} = 2。在方差给定的分布中,高斯分布携带的信息最少。

Proposition A.11 (Fisher information at unit variance)命题 A.11(单位方差下的 Fisher 信息)

Let ss have zero mean, unit variance and a continuously differentiable density p>0p > 0 with score ψ=−p′/p\psi = -p'/p, J=E[ψ(s)2]<∞\mathcal{J} = \E[\psi(s)^2] < \infty, and u p(u)→0u\,p(u) \to 0 as ∣u∣→∞|u| \to \infty. Then J≥1\mathcal{J} \ge 1, with equality if and only if s∼N(0,1)s \sim \Ncal(0, 1).

设 ss 均值为零、方差为一,具有连续可微的密度 p>0p > 0,其得分为 ψ=−p′/p\psi = -p'/p,J=E[ψ(s)2]<∞\mathcal{J} = \E[\psi(s)^2] < \infty,并且当 ∣u∣→∞|u| \to \infty 时 u p(u)→0u\,p(u) \to 0。则 J≥1\mathcal{J} \ge 1,等号成立当且仅当 s∼N(0,1)s \sim \Ncal(0, 1)。

Proof.证明

Integrating by parts, E[s ψ(s)]=−∫u p′(u) du=[−u p(u)]−∞∞+∫p(u) du=1\E[s\,\psi(s)] = -\int u\,p'(u)\,du = \bigl[-u\,p(u)\bigr]_{-\infty}^{\infty} + \int p(u)\,du = 1; the integral converges absolutely because E∣s ψ(s)∣≤(E[s2] J)1/2\E|s\,\psi(s)| \le (\E[s^2]\,\mathcal{J})^{1/2}. By the Cauchy–Schwarz inequality, 1=(E[s ψ(s)])2≤E[s2] E[ψ(s)2]=J1 = (\E[s\,\psi(s)])^2 \le \E[s^2]\,\E[\psi(s)^2] = \mathcal{J}. Equality holds if and only if ψ(s)=c s\psi(s) = c\,s almost surely for a constant cc; since ψ\psi is continuous and p>0p > 0, this means ψ(u)=cu\psi(u) = cu for all uu, so (log⁡p)′=−cu(\log p)' = -cu and pp is a Gaussian density with variance 1/c1/c. Unit variance forces c=1c = 1.

分部积分得 E[s ψ(s)]=−∫u p′(u) du=[−u p(u)]−∞∞+∫p(u) du=1\E[s\,\psi(s)] = -\int u\,p'(u)\,du = \bigl[-u\,p(u)\bigr]_{-\infty}^{\infty} + \int p(u)\,du = 1;由于 E∣s ψ(s)∣≤(E[s2] J)1/2\E|s\,\psi(s)| \le (\E[s^2]\,\mathcal{J})^{1/2},这个积分绝对收敛。由 Cauchy–Schwarz 不等式,1=(E[s ψ(s)])2≤E[s2] E[ψ(s)2]=J1 = (\E[s\,\psi(s)])^2 \le \E[s^2]\,\E[\psi(s)^2] = \mathcal{J}。等号成立当且仅当存在常数 cc 使 ψ(s)=c s\psi(s) = c\,s 几乎必然成立;由于 ψ\psi 连续且 p>0p > 0,这意味着对所有 uu 有 ψ(u)=cu\psi(u) = cu,于是 (log⁡p)′=−cu(\log p)' = -cu,pp 是方差为 1/c1/c 的高斯密度。单位方差迫使 c=1c = 1。

The same integration by parts applies to absolutely continuous densities such as the Laplace density, so the inequality holds for them too. Consequently, for two unit-variance sources JiJj≥1\mathcal{J}_i\mathcal{J}_j \ge 1, with equality only when both are Gaussian, which is why the conditions JiJj>1\mathcal{J}_i\mathcal{J}_j > 1 in Chapter 7第 7 章 and Chapter 10第 10 章 fail only for pairs of Gaussian sources.

同样的分部积分也适用于拉普拉斯密度这样的绝对连续密度,所以不等式对它们同样成立。因此,对两个单位方差的源有 JiJj≥1\mathcal{J}_i\mathcal{J}_j \ge 1,仅当二者都是高斯的时取等号;这就是为什么Chapter 7第 7 章与Chapter 10第 10 章中的条件 JiJj>1\mathcal{J}_i\mathcal{J}_j > 1 只对一对高斯源才不成立。

Key points

本章要点

  • O(n)O(n) has dimension n(n−1)/2n(n-1)/2, and its tangent directions at Q\bQ are ΩQ\boldsymbol{\Omega}\bQ with Ω\boldsymbol{\Omega} skew-symmetric. Givens rotations generate SO(n)SO(n), and together with a diagonal sign matrix they give every element of O(n)O(n). A scaled permutation is orthogonal only when its scales are ±1\pm 1.
  • A positive semidefinite matrix has a unique positive semidefinite square root, C1/2=EΛ1/2E⊤\bC^{1/2} = \mathbf{E}\bLambda^{1/2}\mathbf{E}^\top; non-symmetric square roots differ by orthogonal factors.
  • The polar factor (WW⊤)−1/2W=LR⊤(\bW\bW^\top)^{-1/2}\bW = \mathbf{L}\bR^\top is the orthogonal matrix closest to W\bW; it is what symmetric orthogonalization computes.
  • vec⁡(AXB)=(B⊤⊗A)vec⁡(X)\vecop(\bA\bX\bB) = (\bB^\top \otimes \bA)\vecop(\bX). Cumulant tensors transform multilinearly, and for x=As\bx = \bA\bs the fourth-order cumulant tensor is a symmetric CP decomposition along the columns of A\bA.
  • Compute gradients through differentials; in particular ∂log⁡∣det⁡W∣/∂W=W−⊤\partial\log|\det\bW|/\partial\bW = \bW^{-\top}.
  • On the sphere, critical points satisfy ∇f(w)=λw\nabla f(\bw) = \lambda\bw and the Riemannian gradient is (I−ww⊤)∇f(\bI - \bw\bw^\top)\nabla f; on O(n)O(n) it is skew⁡(GU⊤)U\operatorname{skew}(\bG\bU^\top)\bU.
  • The delta method, the sandwich covariance of M-estimators and the Cramér–Rao bound are the tools of Part IV; a unit-variance law has Fisher information J≥1\mathcal{J} \ge 1, with equality only for the Gaussian.
  • O(n)O(n) 的维数为 n(n−1)/2n(n-1)/2,在 Q\bQ 处的切方向为 ΩQ\boldsymbol{\Omega}\bQ(Ω\boldsymbol{\Omega} 反对称)。Givens 旋转生成 SO(n)SO(n),再乘上一个符号对角矩阵即可得到 O(n)O(n) 的任意元素。带尺度的置换矩阵仅当尺度为 ±1\pm 1 时才是正交的。
  • 半正定矩阵有唯一的半正定平方根 C1/2=EΛ1/2E⊤\bC^{1/2} = \mathbf{E}\bLambda^{1/2}\mathbf{E}^\top;非对称的平方根之间相差正交因子。
  • 极因子 (WW⊤)−1/2W=LR⊤(\bW\bW^\top)^{-1/2}\bW = \mathbf{L}\bR^\top 是离 W\bW 最近的正交矩阵,对称正交化计算的正是它。
  • vec⁡(AXB)=(B⊤⊗A)vec⁡(X)\vecop(\bA\bX\bB) = (\bB^\top \otimes \bA)\vecop(\bX)。累积量张量按多重线性方式变换;对 x=As\bx = \bA\bs,四阶累积量张量是沿 A\bA 各列的对称 CP 分解。
  • 通过微分来计算梯度;特别地,∂log⁡∣det⁡W∣/∂W=W−⊤\partial\log|\det\bW|/\partial\bW = \bW^{-\top}。
  • 在球面上,临界点满足 ∇f(w)=λw\nabla f(\bw) = \lambda\bw,黎曼梯度为 (I−ww⊤)∇f(\bI - \bw\bw^\top)\nabla f;在 O(n)O(n) 上则为 skew⁡(GU⊤)U\operatorname{skew}(\bG\bU^\top)\bU。
  • delta 方法、M 估计量的三明治协方差与 Cramér–Rao 界是第四部分的工具;单位方差分布的 Fisher 信息 J≥1\mathcal{J} \ge 1,仅对高斯分布取等号。

A.9 Exercises习题

Exercise A.1习题 A.1 basic基础

Let Ω\boldsymbol{\Omega} be skew-symmetric. Show that I−Ω\bI - \boldsymbol{\Omega} is invertible and that the Cayley transform Q=(I−Ω)−1(I+Ω)\bQ = (\bI - \boldsymbol{\Omega})^{-1}(\bI + \boldsymbol{\Omega}) is orthogonal. For n=2n = 2 and Ω=θ2(0−110)\boldsymbol{\Omega} = \tfrac{\theta}{2}\bigl(\begin{smallmatrix} 0 & -1 \\ 1 & 0 \end{smallmatrix}\bigr), which rotation is Q\bQ?

设 Ω\boldsymbol{\Omega} 反对称。证明 I−Ω\bI - \boldsymbol{\Omega} 可逆,并且 Cayley 变换 Q=(I−Ω)−1(I+Ω)\bQ = (\bI - \boldsymbol{\Omega})^{-1}(\bI + \boldsymbol{\Omega}) 是正交矩阵。当 n=2n = 2、Ω=θ2(0−110)\boldsymbol{\Omega} = \tfrac{\theta}{2}\bigl(\begin{smallmatrix} 0 & -1 \\ 1 & 0 \end{smallmatrix}\bigr) 时,Q\bQ 是哪个旋转?

Show solution查看解答

For every v\bv, v⊤Ωv\bv^\top\boldsymbol{\Omega}\bv is a scalar equal to its transpose v⊤Ω⊤v=−v⊤Ωv\bv^\top\boldsymbol{\Omega}^\top\bv = -\bv^\top\boldsymbol{\Omega}\bv, so it is zero. Hence v⊤(I∓Ω)v=∥v∥2>0\bv^\top(\bI \mp \boldsymbol{\Omega})\bv = \|\bv\|^2 > 0 for v≠0\bv \ne \bzero, and both I−Ω\bI - \boldsymbol{\Omega} and I+Ω\bI + \boldsymbol{\Omega} are invertible. Since (I−Ω)⊤=I+Ω(\bI - \boldsymbol{\Omega})^\top = \bI + \boldsymbol{\Omega} and (I+Ω)⊤=I−Ω(\bI + \boldsymbol{\Omega})^\top = \bI - \boldsymbol{\Omega}, we get Q⊤=(I−Ω)(I+Ω)−1\bQ^\top = (\bI - \boldsymbol{\Omega})(\bI + \boldsymbol{\Omega})^{-1}. The matrices (I+Ω)−1(\bI + \boldsymbol{\Omega})^{-1} and (I−Ω)−1(\bI - \boldsymbol{\Omega})^{-1} commute, being functions of the same matrix, so Q⊤Q=(I−Ω)(I+Ω)−1(I−Ω)−1(I+Ω)=(I−Ω)(I−Ω)−1(I+Ω)−1(I+Ω)=I\bQ^\top\bQ = (\bI - \boldsymbol{\Omega})(\bI + \boldsymbol{\Omega})^{-1}(\bI - \boldsymbol{\Omega})^{-1}(\bI + \boldsymbol{\Omega}) = (\bI - \boldsymbol{\Omega})(\bI - \boldsymbol{\Omega})^{-1}(\bI + \boldsymbol{\Omega})^{-1}(\bI + \boldsymbol{\Omega}) = \bI. For n=2n = 2, write a=θ/2a = \theta/2; then Q=11+a2(1−a2−2a2a1−a2)\bQ = \frac{1}{1 + a^2}\bigl(\begin{smallmatrix} 1 - a^2 & -2a \\ 2a & 1 - a^2 \end{smallmatrix}\bigr), the rotation by the angle 2arctan⁡(θ/2)2\arctan(\theta/2), since cos⁡ϕ=(1−a2)/(1+a2)\cos\phi = (1 - a^2)/(1 + a^2) and sin⁡ϕ=2a/(1+a2)\sin\phi = 2a/(1 + a^2) for ϕ=2arctan⁡a\phi = 2\arctan a. For small θ\theta the angle is θ+O(θ3)\theta + O(\theta^3), so the Cayley transform agrees with the exponential to first order, as a retraction should.

对任意 v\bv,标量 v⊤Ωv\bv^\top\boldsymbol{\Omega}\bv 等于它的转置 v⊤Ω⊤v=−v⊤Ωv\bv^\top\boldsymbol{\Omega}^\top\bv = -\bv^\top\boldsymbol{\Omega}\bv,所以为零。于是对 v≠0\bv \ne \bzero 有 v⊤(I∓Ω)v=∥v∥2>0\bv^\top(\bI \mp \boldsymbol{\Omega})\bv = \|\bv\|^2 > 0,I−Ω\bI - \boldsymbol{\Omega} 与 I+Ω\bI + \boldsymbol{\Omega} 都可逆。由于 (I−Ω)⊤=I+Ω(\bI - \boldsymbol{\Omega})^\top = \bI + \boldsymbol{\Omega},(I+Ω)⊤=I−Ω(\bI + \boldsymbol{\Omega})^\top = \bI - \boldsymbol{\Omega},得 Q⊤=(I−Ω)(I+Ω)−1\bQ^\top = (\bI - \boldsymbol{\Omega})(\bI + \boldsymbol{\Omega})^{-1}。(I+Ω)−1(\bI + \boldsymbol{\Omega})^{-1} 与 (I−Ω)−1(\bI - \boldsymbol{\Omega})^{-1} 都是同一个矩阵的函数,因而可交换,所以 Q⊤Q=(I−Ω)(I+Ω)−1(I−Ω)−1(I+Ω)=(I−Ω)(I−Ω)−1(I+Ω)−1(I+Ω)=I\bQ^\top\bQ = (\bI - \boldsymbol{\Omega})(\bI + \boldsymbol{\Omega})^{-1}(\bI - \boldsymbol{\Omega})^{-1}(\bI + \boldsymbol{\Omega}) = (\bI - \boldsymbol{\Omega})(\bI - \boldsymbol{\Omega})^{-1}(\bI + \boldsymbol{\Omega})^{-1}(\bI + \boldsymbol{\Omega}) = \bI。当 n=2n = 2 时,记 a=θ/2a = \theta/2,则 Q=11+a2(1−a2−2a2a1−a2)\bQ = \frac{1}{1 + a^2}\bigl(\begin{smallmatrix} 1 - a^2 & -2a \\ 2a & 1 - a^2 \end{smallmatrix}\bigr);由于对 ϕ=2arctan⁡a\phi = 2\arctan a 有 cos⁡ϕ=(1−a2)/(1+a2)\cos\phi = (1 - a^2)/(1 + a^2)、sin⁡ϕ=2a/(1+a2)\sin\phi = 2a/(1 + a^2),它是旋转 2arctan⁡(θ/2)2\arctan(\theta/2) 角的矩阵。当 θ\theta 很小时这个角为 θ+O(θ3)\theta + O(\theta^3),所以 Cayley 变换与矩阵指数在一阶上一致,这正是收缩映射应有的性质。

Exercise A.2习题 A.2 basic基础

Let z\bz be white, let GG be a smooth function with derivative gg, and let f(w)=E[G(w⊤z)]f(\bw) = \E[G(\bw^\top\bz)]; assume that the expectations exist and that differentiation under the expectation is allowed. Show that the critical points of ff on the unit sphere satisfy E[z g(w⊤z)]=βw\E[\bz\,g(\bw^\top\bz)] = \beta\bw with β=E[(w⊤z) g(w⊤z)]\beta = \E[(\bw^\top\bz)\,g(\bw^\top\bz)]. What happens for G(u)=u2G(u) = u^2?

设 z\bz 是白的,GG 是光滑函数、导数为 gg,令 f(w)=E[G(w⊤z)]f(\bw) = \E[G(\bw^\top\bz)];假设各期望存在,且可以在期望号下求导。证明 ff 在单位球面上的临界点满足 E[z g(w⊤z)]=βw\E[\bz\,g(\bw^\top\bz)] = \beta\bw,其中 β=E[(w⊤z) g(w⊤z)]\beta = \E[(\bw^\top\bz)\,g(\bw^\top\bz)]。当 G(u)=u2G(u) = u^2 时会怎样?

Show solution查看解答

Differentiating under the expectation, ∇f(w)=E[z g(w⊤z)]\nabla f(\bw) = \E[\bz\,g(\bw^\top\bz)]. By Proposition A.7命题 A.7(a), a critical point satisfies ∇f(w)=λw\nabla f(\bw) = \lambda\bw with λ=w⊤∇f(w)=E[(w⊤z) g(w⊤z)]=β\lambda = \bw^\top\nabla f(\bw) = \E[(\bw^\top\bz)\,g(\bw^\top\bz)] = \beta. This is the equation from which FastICA starts (Chapter 8第 8 章). For G(u)=u2G(u) = u^2, g(u)=2ug(u) = 2u and ∇f(w)=2E[zz⊤]w=2w\nabla f(\bw) = 2\E[\bz\bz^\top]\bw = 2\bw because z\bz is white, so every unit vector is a critical point: f≡1f \equiv 1 on the sphere. Second-order statistics carry no information about directions in whitened data, as Section 3.5第 3.5 节 explains.

在期望号下求导得 ∇f(w)=E[z g(w⊤z)]\nabla f(\bw) = \E[\bz\,g(\bw^\top\bz)]。由Proposition A.7命题 A.7(a),临界点满足 ∇f(w)=λw\nabla f(\bw) = \lambda\bw,其中 λ=w⊤∇f(w)=E[(w⊤z) g(w⊤z)]=β\lambda = \bw^\top\nabla f(\bw) = \E[(\bw^\top\bz)\,g(\bw^\top\bz)] = \beta。这正是 FastICA 的出发点(Chapter 8第 8 章)。当 G(u)=u2G(u) = u^2 时,g(u)=2ug(u) = 2u,由于 z\bz 是白的,∇f(w)=2E[zz⊤]w=2w\nabla f(\bw) = 2\E[\bz\bz^\top]\bw = 2\bw,所以每个单位向量都是临界点:ff 在球面上恒等于 1。白化数据的二阶统计量不携带任何关于方向的信息,正如Section 3.5第 3.5 节所解释的那样。

Exercise A.3习题 A.3 intermediate进阶

Show that for orthogonal Q\bQ the multilinear transform T↦T×1Q×2Q×3Q×4Q\mathcal{T} \mapsto \mathcal{T} \times_1 \bQ \times_2 \bQ \times_3 \bQ \times_4 \bQ preserves the Frobenius norm. Deduce that if z\bz is white and y=Uz\by = \bU\bz with U\bU orthogonal, then ∑ijklcum⁡(yi,yj,yk,yl)2\sum_{ijkl}\cum(y_i, y_j, y_k, y_l)^2 does not depend on U\bU, so that maximizing ∑ikurt⁡(yi)2\sum_i\kurt(y_i)^2 over U\bU is the same as minimizing the sum of squares of all cross-cumulants (Comon, 1994).

证明:对正交矩阵 Q\bQ,多重线性变换 T↦T×1Q×2Q×3Q×4Q\mathcal{T} \mapsto \mathcal{T} \times_1 \bQ \times_2 \bQ \times_3 \bQ \times_4 \bQ 保持 Frobenius 范数不变。由此推出:若 z\bz 是白的,y=Uz\by = \bU\bz 且 U\bU 正交,则 ∑ijklcum⁡(yi,yj,yk,yl)2\sum_{ijkl}\cum(y_i, y_j, y_k, y_l)^2 与 U\bU 无关;因而在 U\bU 上最大化 ∑ikurt⁡(yi)2\sum_i\kurt(y_i)^2,等价于最小化所有交叉累积量的平方和(Comon, 1994)。

Show solution查看解答

Arrange the entries of T\mathcal{T} in a vector t\bt, with the first index running fastest. The transformed entries ∑pqrtqipqjqqkrqlt tpqrt\sum_{pqrt} q_{ip}q_{jq}q_{kr}q_{lt}\,t_{pqrt} then form the vector (Q⊗Q⊗Q⊗Q)t(\bQ \otimes \bQ \otimes \bQ \otimes \bQ)\bt. By the rules of Section A.5第 A.5 节, a Kronecker product of orthogonal matrices is orthogonal, so the Euclidean norm of t\bt, which is the Frobenius norm of T\mathcal{T}, is preserved. For y=Uz\by = \bU\bz, multilinearity gives Ky=Kz×1U×2U×3U×4U\mathcal{K}_{\by} = \mathcal{K}_{\bz} \times_1 \bU \times_2 \bU \times_3 \bU \times_4 \bU, so ∥Ky∥F=∥Kz∥F\|\mathcal{K}_{\by}\|_F = \|\mathcal{K}_{\bz}\|_F for every U\bU. Split ∥Ky∥F2\|\mathcal{K}_{\by}\|_F^2 into the diagonal entries, cum⁡(yi,yi,yi,yi)2=kurt⁡(yi)2\cum(y_i, y_i, y_i, y_i)^2 = \kurt(y_i)^2, and the squared cross-cumulants. Their total is fixed, so maximizing the first part is the same as minimizing the second. Independent outputs have zero cross-cumulants (Proposition 2.8命题 2.8), so at a separating rotation the whole norm sits on the diagonal and ∑ikurt⁡(yi)2\sum_i\kurt(y_i)^2 attains its largest possible value ∥Kz∥F2\|\mathcal{K}_{\bz}\|_F^2.

把 T\mathcal{T} 的元素按第一个下标变化最快的次序排成向量 t\bt。变换后的元素 ∑pqrtqipqjqqkrqlt tpqrt\sum_{pqrt} q_{ip}q_{jq}q_{kr}q_{lt}\,t_{pqrt} 构成向量 (Q⊗Q⊗Q⊗Q)t(\bQ \otimes \bQ \otimes \bQ \otimes \bQ)\bt。由Section A.5第 A.5 节中的规则,正交矩阵的 Kronecker 积仍是正交矩阵,所以 t\bt 的欧氏范数,也就是 T\mathcal{T} 的 Frobenius 范数,保持不变。对 y=Uz\by = \bU\bz,由多重线性得 Ky=Kz×1U×2U×3U×4U\mathcal{K}_{\by} = \mathcal{K}_{\bz} \times_1 \bU \times_2 \bU \times_3 \bU \times_4 \bU,所以对每个 U\bU 都有 ∥Ky∥F=∥Kz∥F\|\mathcal{K}_{\by}\|_F = \|\mathcal{K}_{\bz}\|_F。把 ∥Ky∥F2\|\mathcal{K}_{\by}\|_F^2 分成对角元部分 cum⁡(yi,yi,yi,yi)2=kurt⁡(yi)2\cum(y_i, y_i, y_i, y_i)^2 = \kurt(y_i)^2 与交叉累积量的平方两部分。二者之和固定,所以最大化前者等价于最小化后者。独立的输出没有交叉累积量(Proposition 2.8命题 2.8),所以在实现分离的旋转处,全部范数都集中在对角元上,∑ikurt⁡(yi)2\sum_i\kurt(y_i)^2 达到可能的最大值 ∥Kz∥F2\|\mathcal{K}_{\bz}\|_F^2。

Exercise A.4习题 A.4 intermediate进阶

Using (A.5), compute the asymptotic variance vv of the sample kurtosis for the unit-variance uniform law, whose moments are μ2=1\mu_2 = 1, μ4=9/5\mu_4 = 9/5, μ6=27/7\mu_6 = 27/7 and μ8=9\mu_8 = 9. Compare it with the Gaussian and Laplace values, and find the sample size at which the standard error of k^\hat k equals 0.10.1 for each of the three laws.

利用 (A.5),计算单位方差均匀分布的样本峭度的渐近方差 vv,其各阶矩为 μ2=1\mu_2 = 1,μ4=9/5\mu_4 = 9/5,μ6=27/7\mu_6 = 27/7,μ8=9\mu_8 = 9。把它与高斯分布和拉普拉斯分布的值比较,并对这三个分布分别求出使 k^\hat k 的标准误等于 0.10.1 的样本量。

Show solution查看解答

v=9−8125−12(277−95)+36(95−1)=14425−86435+1445=1008−4320+5040175=1728175≈9.87v = 9 - \tfrac{81}{25} - 12\bigl(\tfrac{27}{7} - \tfrac95\bigr) + 36\bigl(\tfrac95 - 1\bigr) = \tfrac{144}{25} - \tfrac{864}{35} + \tfrac{144}{5} = \tfrac{1008 - 4320 + 5040}{175} = \tfrac{1728}{175} \approx 9.87. So vv is about 9.879.87 for the uniform law, 2424 for the Gaussian and 16561656 for the Laplace law: bounded, sub-Gaussian laws give precise kurtosis estimates, and heavy tails give very imprecise ones. The standard error v/T\sqrt{v/T} equals 0.10.1 when T=100 vT = 100\,v, that is, at about 990 samples for the uniform law, 2400 for the Gaussian and 165 600 for the Laplace law.

v=9−8125−12(277−95)+36(95−1)=14425−86435+1445=1008−4320+5040175=1728175≈9.87v = 9 - \tfrac{81}{25} - 12\bigl(\tfrac{27}{7} - \tfrac95\bigr) + 36\bigl(\tfrac95 - 1\bigr) = \tfrac{144}{25} - \tfrac{864}{35} + \tfrac{144}{5} = \tfrac{1008 - 4320 + 5040}{175} = \tfrac{1728}{175} \approx 9.87。所以均匀分布的 vv 约为 9.879.87,高斯分布为 2424,拉普拉斯分布为 16561656:有界的亚高斯分布给出精确的峭度估计,重尾分布给出的估计则很不精确。标准误 v/T\sqrt{v/T} 等于 0.10.1 当且仅当 T=100 vT = 100\,v,即均匀分布约需 990 个样本,高斯分布需 2400 个,拉普拉斯分布需 165 600 个。

Exercise A.5习题 A.5 advanced硬核

Let x=θ+sx = \theta + s, where ss has a symmetric, continuously differentiable density p>0p > 0 with score ψ=−p′/p\psi = -p'/p and Fisher information J=E[ψ(s)2]<∞\mathcal{J} = \E[\psi(s)^2] < \infty. Let φ\varphi be odd and twice continuously differentiable with φ′\varphi' and φ′′\varphi'' bounded, with E[φ(s)2]<∞\E[\varphi(s)^2] < \infty, E[φ′(s)]≠0\E[\varphi'(s)] \ne 0 and φ(u)p(u)→0\varphi(u)p(u) \to 0 as ∣u∣→∞|u| \to \infty, and let θ^\hat\theta be a consistent solution of ∑tφ(x(t)−θ)=0\sum_t \varphi(x(t) - \theta) = 0. (a) Use Proposition A.10命题 A.10 to show that the asymptotic variance of T(θ^−θ)\sqrt T(\hat\theta - \theta) is E[φ(s)2]/(E[φ′(s)])2\E[\varphi(s)^2]/(\E[\varphi'(s)])^2. (b) Show that E[φ′(s)]=E[φ(s)ψ(s)]\E[\varphi'(s)] = \E[\varphi(s)\psi(s)], and deduce that this variance is at least 1/J1/\mathcal{J}, with equality if and only if φ\varphi is proportional to ψ\psi. (c) Evaluate the bound and the variance for the unit-variance Laplace law and φ=tanh⁡\varphi = \tanh.

设 x=θ+sx = \theta + s,其中 ss 具有对称、连续可微的密度 p>0p > 0,得分为 ψ=−p′/p\psi = -p'/p,Fisher 信息 J=E[ψ(s)2]<∞\mathcal{J} = \E[\psi(s)^2] < \infty。设 φ\varphi 是二阶连续可微的奇函数,φ′\varphi' 与 φ′′\varphi'' 有界,并满足 E[φ(s)2]<∞\E[\varphi(s)^2] < \infty、E[φ′(s)]≠0\E[\varphi'(s)] \ne 0,且当 ∣u∣→∞|u| \to \infty 时 φ(u)p(u)→0\varphi(u)p(u) \to 0;令 θ^\hat\theta 是 ∑tφ(x(t)−θ)=0\sum_t \varphi(x(t) - \theta) = 0 的相合解。(a)利用Proposition A.10命题 A.10证明 T(θ^−θ)\sqrt T(\hat\theta - \theta) 的渐近方差为 E[φ(s)2]/(E[φ′(s)])2\E[\varphi(s)^2]/(\E[\varphi'(s)])^2。(b)证明 E[φ′(s)]=E[φ(s)ψ(s)]\E[\varphi'(s)] = \E[\varphi(s)\psi(s)],并由此推出这个方差不小于 1/J1/\mathcal{J},等号成立当且仅当 φ\varphi 与 ψ\psi 成比例。(c)对单位方差的拉普拉斯分布与 φ=tanh⁡\varphi = \tanh,计算这个下界与方差。

Show solution查看解答

(a) Here h(x,θ)=φ(x−θ)\mathbf{h}(x, \theta) = \varphi(x - \theta), and E[φ(s)]=0\E[\varphi(s)] = 0 because φ\varphi is odd and pp symmetric, so the true θ\theta solves the population equation. Since ∂h/∂θ=−φ′(x−θ)\partial\mathbf{h}/\partial\theta = -\varphi'(x - \theta), we have H=−E[φ′(s)]≠0\bH = -\E[\varphi'(s)] \ne 0 and B=E[φ(s)2]<∞\bB = \E[\varphi(s)^2] < \infty. The first and second derivatives of h\mathbf{h}, −φ′(x−θ)-\varphi'(x - \theta) and φ′′(x−θ)\varphi''(x - \theta), are bounded by constants, so all the conditions of Proposition A.10命题 A.10 hold, and the sandwich is B/H2=E[φ(s)2]/(E[φ′(s)])2\bB/\bH^2 = \E[\varphi(s)^2]/(\E[\varphi'(s)])^2. (b) Integrating by parts and using φp→0\varphi p \to 0 at ±∞\pm\infty, E[φ′(s)]=∫φ′(u)p(u) du=−∫φ(u)p′(u) du=∫φ(u)ψ(u)p(u) du=E[φ(s)ψ(s)]\E[\varphi'(s)] = \int\varphi'(u)p(u)\,du = -\int\varphi(u)p'(u)\,du = \int\varphi(u)\psi(u)p(u)\,du = \E[\varphi(s)\psi(s)]. By the Cauchy–Schwarz inequality, (E[φψ])2≤E[φ2] E[ψ2]=E[φ2] J(\E[\varphi\psi])^2 \le \E[\varphi^2]\,\E[\psi^2] = \E[\varphi^2]\,\mathcal{J}, so E[φ2]/(E[φ′])2≥1/J\E[\varphi^2]/(\E[\varphi'])^2 \ge 1/\mathcal{J}, with equality if and only if φ=cψ\varphi = c\psi. The maximum-likelihood choice φ=ψ\varphi = \psi attains the Cramér–Rao bound 1/(TJ)1/(T\mathcal{J}). (c) The Laplace density is not differentiable at 0, but it is absolutely continuous and the integration by parts in (b) still holds. The function tanh⁡\tanh is odd and bounded, and its derivatives 1−tanh⁡21 - \tanh^2 and −2tanh⁡ (1−tanh⁡2)-2\tanh\,(1 - \tanh^2) are bounded, so it satisfies the hypotheses. With b=1/2b = 1/\sqrt2, J=1/b2=2\mathcal{J} = 1/b^2 = 2, so the bound is 1/21/2. Numerical integration gives E[tanh⁡2(s)]≈0.3235\E[\tanh^2(s)] \approx 0.3235 and E[tanh⁡′(s)]=E[1−tanh⁡2(s)]≈0.6765\E[\tanh'(s)] = \E[1 - \tanh^2(s)] \approx 0.6765, so the asymptotic variance of the tanh estimator is about 0.3235/0.67652≈0.7070.3235/0.6765^2 \approx 0.707, some 41% above the bound. Chapter 10第 10 章 carries out the same comparison for ICA (Proposition 10.9命题 10.9).

(a)这里 h(x,θ)=φ(x−θ)\mathbf{h}(x, \theta) = \varphi(x - \theta);由于 φ\varphi 是奇函数且 pp 对称,E[φ(s)]=0\E[\varphi(s)] = 0,所以真值 θ\theta 满足总体方程。由 ∂h/∂θ=−φ′(x−θ)\partial\mathbf{h}/\partial\theta = -\varphi'(x - \theta) 得 H=−E[φ′(s)]≠0\bH = -\E[\varphi'(s)] \ne 0,B=E[φ(s)2]<∞\bB = \E[\varphi(s)^2] < \infty。h\mathbf{h} 的一阶与二阶导数 −φ′(x−θ)-\varphi'(x - \theta) 与 φ′′(x−θ)\varphi''(x - \theta) 都被常数控制,所以Proposition A.10命题 A.10的条件全部成立,三明治协方差为 B/H2=E[φ(s)2]/(E[φ′(s)])2\bB/\bH^2 = \E[\varphi(s)^2]/(\E[\varphi'(s)])^2。(b)分部积分并利用 φp\varphi p 在 ±∞\pm\infty 处趋于零,得 E[φ′(s)]=∫φ′(u)p(u) du=−∫φ(u)p′(u) du=∫φ(u)ψ(u)p(u) du=E[φ(s)ψ(s)]\E[\varphi'(s)] = \int\varphi'(u)p(u)\,du = -\int\varphi(u)p'(u)\,du = \int\varphi(u)\psi(u)p(u)\,du = \E[\varphi(s)\psi(s)]。由 Cauchy–Schwarz 不等式,(E[φψ])2≤E[φ2] E[ψ2]=E[φ2] J(\E[\varphi\psi])^2 \le \E[\varphi^2]\,\E[\psi^2] = \E[\varphi^2]\,\mathcal{J},所以 E[φ2]/(E[φ′])2≥1/J\E[\varphi^2]/(\E[\varphi'])^2 \ge 1/\mathcal{J},等号成立当且仅当 φ=cψ\varphi = c\psi。最大似然的选择 φ=ψ\varphi = \psi 达到 Cramér–Rao 界 1/(TJ)1/(T\mathcal{J})。(c)拉普拉斯密度在 0 处不可微,但它绝对连续,(b)中的分部积分仍然成立。tanh⁡\tanh 是有界的奇函数,其导数 1−tanh⁡21 - \tanh^2 与 −2tanh⁡ (1−tanh⁡2)-2\tanh\,(1 - \tanh^2) 也有界,所以它满足题设条件。取 b=1/2b = 1/\sqrt2,J=1/b2=2\mathcal{J} = 1/b^2 = 2,下界为 1/21/2。数值积分给出 E[tanh⁡2(s)]≈0.3235\E[\tanh^2(s)] \approx 0.3235,E[tanh⁡′(s)]=E[1−tanh⁡2(s)]≈0.6765\E[\tanh'(s)] = \E[1 - \tanh^2(s)] \approx 0.6765,所以 tanh 估计量的渐近方差约为 0.3235/0.67652≈0.7070.3235/0.6765^2 \approx 0.707,比下界高约 41%。Chapter 10第 10 章对 ICA 作了同样的比较(Proposition 10.9命题 10.9)。

A.10 References参考文献

  1. Absil, P.-A., Mahony, R., & Sepulchre, R. (2008). Optimization Algorithms on Matrix Manifolds. Princeton University Press.
  2. Comon, P. (1994). Independent component analysis, a new concept? Signal Processing, 36(3), 287–314.
  3. Edelman, A., Arias, T. A., & Smith, S. T. (1998). The geometry of algorithms with orthogonality constraints. SIAM Journal on Matrix Analysis and Applications, 20(2), 303–353.
  4. Golub, G. H., & Van Loan, C. F. (2013). Matrix Computations (4th ed.). Johns Hopkins University Press.
  5. Higham, N. J. (1986). Computing the polar decomposition—with applications. SIAM Journal on Scientific and Statistical Computing, 7(4), 1160–1174.
  6. Hillar, C. J., & Lim, L.-H. (2013). Most tensor problems are NP-hard. Journal of the ACM, 60(6), Article 45.
  7. Horn, R. A., & Johnson, C. R. (2013). Matrix Analysis (2nd ed.). Cambridge University Press.
  8. Kolda, T. G., & Bader, B. W. (2009). Tensor decompositions and applications. SIAM Review, 51(3), 455–500.
  9. Lehmann, E. L., & Casella, G. (1998). Theory of Point Estimation (2nd ed.). Springer.
  10. Magnus, J. R., & Neudecker, H. (2019). Matrix Differential Calculus with Applications in Statistics and Econometrics (3rd ed.). Wiley.
  11. van der Vaart, A. W. (1998). Asymptotic Statistics. Cambridge University Press.