Latent factors, ICA, NMF, and representation ambiguity

By [F1], bilinearity of covariance gives $\mathrm{Cov}(X)=A\mathrm{Cov}(S)A^T+A\mathrm{Cov}(S,\epsilon )+\mathrm{Cov}(\epsilon ,S)A^T+\mathrm{Cov}(\epsilon )=A{A}^T+D$.

The transformed factors have mean zero, covariance $R{I}_q{R}^T=I_q$, and cross-covariance $\mathrm{Cov}(RS,\epsilon )=R\mathrm{Cov}(S,\epsilon )=0$. Also $A_R{S}_R=A{R}^{T}RS=AS$, and $A_R{A}_R^T+D=A{R}^{T}R{A}^T+D=A{A}^T+D$. Hence the transformed pair is another covariance factor model with the same covariance.

By [F1], $C=\mathrm{Cov}(X)=A{I}_p{A}^T=A{A}^T$. For every nonzero $v\in {\mathbb{R}}^p$, $v^{T}Cv=\parallel A^{T}v{\parallel }^2>0$ because $A$ is invertible; hence $C$ is symmetric positive definite.

By [F2], write $A=U\mathrm{\Sigma }{V}^T$ with $U,V$ orthogonal. Since $A$ is invertible, every diagonal entry of $\mathrm{\Sigma }$ is positive. Thus $C=U{\mathrm{\Sigma }}^2{U}^T$, and $B=U{\mathrm{\Sigma }}^{-1}{U}^T$ is symmetric positive definite with $B^2=C^{-1}$.

If $K$ is another symmetric positive-definite matrix with $K^2=C^{-1}$, then $K$ commutes with $C^{-1}$ and preserves each eigenspace $E_{\lambda }$ of $C^{-1}$. On each such subspace, [F3] gives an orthonormal eigenbasis for $K{\mathrm{\mid }}_{E_{\lambda }}$. Its eigenvalues $\mu$ are positive and satisfy ${\mu }^2=\lambda$, so $K{\mathrm{\mid }}_{E_{\lambda }}=\sqrt{\lambda }I$. This determines $K$ uniquely and agrees with $B$.

Therefore $C^{-1\mathrm{/}2}=B$ and $C^{-1\mathrm{/}2}A=U{\mathrm{\Sigma }}^{-1}{U}^{T}U\mathrm{\Sigma }{V}^T=U{V}^T$. Since $(U{V}^T{)}^T(U{V}^T)=I_p$, the whitened mixing is orthogonal and .

For the rescaled factors, each entry is $(W^{\mathrm{\prime }}{H}^{\mathrm{\prime }}{)}_{ij}={\sum }_{k=1}^r(c_k{W}_{ik})(c_k^{-1}{H}_{kj})={\sum }_{k=1}^r{W}_{ik}{H}_{kj}=(WH{)}_{ij}$. Thus every rate is unchanged.

For the permuted factors, $(W^{\mathrm{\prime }\mathrm{\prime }}{H}^{\mathrm{\prime }\mathrm{\prime }}{)}_{ij}={\sum }_{k=1}^r{W}_{i,\pi (k)}{H}_{\pi (k),j}=(WH{)}_{ij}$ because $\pi$ reorders the finite sum. Equal rate matrices give the same independent Poisson law and likelihood by [F1].

When all rates are positive, [F1] and independence give joint likelihood $L(W,H)={\prod }_{i,j}{e}^{-{\lambda }_{ij}}{\lambda }_{ij}^{x_{ij}}\mathrm{/}{x}_{ij}!$. Taking logarithms and using $\log \prod a_{ij}=\sum \log a_{ij}$ gives the displayed formula.

If ${\lambda }_{ij}=0$, [F1] gives cell likelihood $1$ when $x_{ij}=0$ and $0$ when $x_{ij}>0$; their log-likelihoods are respectively $0$ and $-\mathrm{\infty }$, exactly the stated conventions. The factorial term depends only on the fixed data, so removing it leaves the same maximizing objective in $W,H$.

By [F9], substitution $x=\sqrt{2}u$ in [F6] gives $\int g(x)dx=1$, so $g$ defines a probability law $\gamma$. These integrals can be read as Lebesgue integrals: on each $[-m,m]$, the uniform dyadic lower and upper sums of the continuous integrand converge to the same Riemann integral and bound its Lebesgue integral; increasing truncations then give the whole-line equality by [F3]. For each integer $k\ge 0$, $\mathrm{\mid }x{\mathrm{\mid }}^{k}g(x)$ is integrable: its tail is bounded by a constant times $e^{-x^2\mathrm{/}4}$, since $\mathrm{\mid }x{\mathrm{\mid }}^k{e}^{-x^2\mathrm{/}4}$ is bounded (differentiate it on $x>0$). If $G$ has law $\gamma$, integration by parts [F8] on $[-m,m]$, followed by [F3], using $g^{\mathrm{\prime }}(x)=-xg(x)$ and vanishing polynomial times Gaussian boundary terms, gives $\mathbb{E}{G}^{2n}=(2n-1)\mathbb{E}{G}^{2n-2}$ and odd moments zero. Hence $m_{2n}:=\mathbb{E}{G}^{2n}=(2n)!\mathrm{/}(2^{n}n!)$. Also ${\phi }_G^{\mathrm{\prime }}(t)=\int ix{e}^{itx}g(x)dx=-t{\phi }_G(t)$ by the same integration by parts, justified by $\int \mathrm{\mid }x\mathrm{\mid }g(x)dx<\mathrm{\infty }$. By [F11], differentiation of $e^{t^2\mathrm{/}2}{\phi }_G(t)$ gives zero and its value at zero is one: ${\phi }_G(t)=e^{-t^2\mathrm{/}2}$ for every real $t$. All integrations by parts first concern smooth functions on compact intervals, whose derivatives are continuous and therefore integrable by [F10]; the same dyadic comparison identifies their Riemann and Lebesgue integrals.

Fix an integer $n\ge 1$ and define $H_{n,t}(x)=(2\sin (tx\mathrm{/}2)\mathrm{/}t{)}^{2n}$ for $t\mathrm{\ne }0$. It is nonnegative and bounded for fixed $t$. The binomial theorem and $2i\sin (v\mathrm{/}2)=e^{iv\mathrm{/}2}-e^{-iv\mathrm{/}2}$ give $H_{n,t}(x)=(-1{)}^n{t}^{-2n}{\sum }_{k=0}^{2n}(-1{)}^k\left(\genfrac{}{}{0px}{}{2n}{k}\right)e^{i(n-k)tx}\mathrm{.}$ For $0<\mathrm{\mid }t\mathrm{\mid }<r\mathrm{/}n$, its expectation under $Y$ therefore equals its expectation under $G$, since every frequency $(n-k)t$ lies in $(-r,r)$. On the prescribed sequence $t_{\mathrm{\ell }}=r\mathrm{/}(2n\mathrm{\ell })$, $H_{n,t_{\mathrm{\ell }}}(x)\to x^{2n}$; moreover $H_{n,t}(x)\le \mathrm{\mid }x{\mathrm{\mid }}^{2n}$ because $\mathrm{\mid }\sin v\mathrm{\mid }\le \mathrm{\mid }v\mathrm{\mid }$. By [F3] applied to $G$, $\mathbb{E}{H}_{n,t_{\mathrm{\ell }}}(G)\to m_{2n}$. Fatou's lemma [F2] applied to $Y$ now gives $\mathbb{E}\mathrm{\mid }Y{\mathrm{\mid }}^{2n}\le m_{2n}<\mathrm{\infty }$. Applying [F3] again, now with the integrable majorant $\mathrm{\mid }Y{\mathrm{\mid }}^{2n}$, shows equality $\mathbb{E}{Y}^{2n}=m_{2n}$.

Every absolute moment of $Y$ is finite, since $\mathrm{\mid }y{\mathrm{\mid }}^k\le 1+\mathrm{\mid }y{\mathrm{\mid }}^{2k}$. By [F4], equality of $\phi$ and ${\phi }_G$ on a neighborhood of zero implies equality of their derivatives there, hence $\mathbb{E}{Y}^k=\mathbb{E}{G}^k$ for every $k\ge 0$. For every $a\ge 0$, the nonnegative partial sums of $\cosh (aY)$ increase to it; [F3] and step 2.1 give $\mathbb{E}\cosh (aY)={\sum }_{n=0}^{\mathrm{\infty }}\frac{a^{2n}{m}_{2n}}{(2n)!}=e^{a^2\mathrm{/}2}\mathrm{.}$ Since $e^{a\mathrm{\mid }y\mathrm{\mid }}\le 2\cosh (ay)$, $\mathbb{E}{e}^{a\mathrm{\mid }Y\mathrm{\mid }}<\mathrm{\infty }$. The identical calculation applies to $G$. For every real $t$, the partial sums of the exponential series for $e^{itY}$ have modulus at most $e^{\mathrm{\mid }t\mathrm{\mid }\mathrm{\mid }Y\mathrm{\mid }}$, so [F3] and moment equality yield $\phi (t)={\sum }_{k=0}^{\mathrm{\infty }}\frac{(it{)}^k\mathbb{E}{Y}^k}{k!}={\sum }_{k=0}^{\mathrm{\infty }}\frac{(it{)}^k\mathbb{E}{G}^k}{k!}=e^{-t^2\mathrm{/}2}\mathrm{.}$ This proves global equality without presupposing analyticity of $\phi$.

For completeness identify the law directly, without importing a general Fourier uniqueness theorem. For $\epsilon >0$ set $g_{\epsilon }(y)={\epsilon }^{-1}g(y\mathrm{/}\epsilon )$. The Gaussian transform computed in step 1.1, with the substitution $s=\epsilon t$, gives $g_{\epsilon }(y)=\frac{1}{2\pi }{\int }_{\mathbb{R}}{e}^{-ity}{e}^{-{\epsilon }^2{t}^2\mathrm{/}2}dt\mathrm{.}$ If $\mu$ is the law of $Y$, then $h_{\mu ,\epsilon }(z)=\int g_{\epsilon }(z-x)\mu (dx)$ is a probability density by Tonelli. Substituting the displayed identity and applying Fubini gives $h_{\mu ,\epsilon }(z)=\frac{1}{2\pi }\int e^{-itz}{e}^{-{\epsilon }^2{t}^2\mathrm{/}2}\phi (t)dt\mathrm{.}$ Absolute integrability for this interchange follows from $\mu (\mathbb{R})=1$ and $\int e^{-{\epsilon }^2{t}^2\mathrm{/}2}dt<\mathrm{\infty }$. By step 3.1 this density equals $h_{\gamma ,\epsilon }$ for every $z$.

For a bounded continuous real $f$, Fubini and substitution give $\int f(z)h_{\mu ,\epsilon }(z)dz=\int \text{\hspace{0.17em}⁣}\int f(x+\epsilon u)g(u)du\mu (dx)\mathrm{.}$ The product-integrand is bounded in absolute integral by $\parallel f{\parallel }_{\mathrm{\infty }}$. Along $\epsilon =1\mathrm{/}\mathrm{\ell }$, dominated convergence gives limit $\int fd\mu$; the same holds for $\gamma$. Step 4.1 therefore implies $\int fd\mu =\int fd\gamma$ for every such $f$. For $a\in \mathbb{R}$ take the continuous functions $f_{a,n}(x)$ equal to $1$ on $x\le a$, equal to $1-n(x-a)$ on $a<x<a+1\mathrm{/}n$, and zero on $x\ge a+1\mathrm{/}n$. They are bounded by one and converge pointwise to ${\mathbf{1}}_{(-\mathrm{\infty },a]}$, so [F3] gives equality of both measures on every such half-line. These half-lines are a pi-system generating the Borel sets: open intervals are countable unions and differences of half-lines with rational endpoints, and every open subset of $\mathbb{R}$ is the union of the rational-endpoint intervals contained in it. Both measures have total mass one; [F7] proves $\mu =\gamma$. This proves the assertion.

Let ${\phi }_i(u)=\mathbb{E}{e}^{iu{S}_i}$ and ${\theta }_j(v)=\mathbb{E}{e}^{iv{T}_j}$. Each $T_j$ has a finite second moment, since $({\sum }_i{O}_{ji}{S}_i{)}^2\le ({\sum }_i{O}_{ji}^2)({\sum }_i{S}_i^2)$ pointwise. Thus [F2] makes all these functions $C^2$, and each is one at zero. Choose $\delta >0$ such that none vanishes on $(-\delta ,\delta )$; this is possible by continuity and finiteness of the families. Define ${\psi }_i(u)={\int }_0^u\frac{{\phi }_i^{\mathrm{\prime }}(s)}{{\phi }_i(s)}ds,{\chi }_j(v)={\int }_0^v\frac{{\theta }_j^{\mathrm{\prime }}(s)}{{\theta }_j(s)}ds\mathrm{.}$ By [F5] applied on every compact subinterval, these are complex-valued $C^2$ functions on that interval, equal to zero at zero. Differentiation shows $(e^{-{\psi }_i}{\phi }_i{)}^{\mathrm{\prime }}=0$, so $e^{{\psi }_i}={\phi }_i$; likewise $e^{{\chi }_j}={\theta }_j$. Moreover ${\psi }_i^{\mathrm{\prime }}(0)=0$ and ${\psi }_i^{\mathrm{\prime }\mathrm{\prime }}(0)={\phi }_i^{\mathrm{\prime }\mathrm{\prime }}(0)-{\phi }_i^{\mathrm{\prime }}(0{)}^2=-1$. This constructs the needed local logarithms without any global nonvanishing hypothesis.

For $\parallel t\parallel <\delta$, also $\parallel O^{T}t\parallel <\delta$. Factoring the characteristic function of $t^{T}T=(O^{T}t{)}^{T}S$ in both independent coordinate systems by [F3] gives ${\prod }_i{\phi }_i((O^{T}t{)}_i)={\prod }_j{\theta }_j(t_j)\mathrm{.}$ Consequently $D(t):={\sum }_i{\psi }_i((O^{T}t{)}_i)-{\sum }_j{\chi }_j(t_j)$ has $e^{D(t)}=1$. It is continuous on the ball, is zero at zero, and takes values in the discrete set $2\pi i\mathbb{Z}$. Along each line segment $s\mapsto st$ the intermediate value property of its imaginary part makes it constant; hence $D(t)=0$ throughout the ball. Differentiating this $C^2$ identity with respect to distinct coordinates $t_j,t_k$ gives ${\sum }_i{O}_{ji}{O}_{ki}{\psi }_i^{\mathrm{\prime }\mathrm{\prime }}((O^{T}t{)}_i)=0(j\mathrm{\ne }k)\mathrm{.}$

Fix $i$ and two distinct rows $j,k$ with $O_{ji}{O}_{ki}\mathrm{\ne }0$. Orthogonality maps the ball onto itself, so in the last identity one may set $O^{T}t=u$ for every $\parallel u\parallel <\delta$. First take $u=v{e}_i$, $\mathrm{\mid }v\mathrm{\mid }<\delta$, and then $u=0$. Subtraction cancels all summands except the $i$th, giving $O_{ji}{O}_{ki}({\psi }_i^{\mathrm{\prime }\mathrm{\prime }}(v)-{\psi }_i^{\mathrm{\prime }\mathrm{\prime }}(0))=0$. Thus ${\psi }_i^{\mathrm{\prime }\mathrm{\prime }}(v)=-1$ throughout $(-\delta ,\delta )$. With ${\psi }_i^{\mathrm{\prime }}(0)={\psi }_i(0)=0$, [F6] applied successively to ${\psi }_i\mathrm{\prime }(v)+v$ and ${\psi }_i(v)+v^2\mathrm{/}2$ gives ${\psi }_i(v)=-v^2\mathrm{/}2$ and ${\phi }_i(v)=e^{-v^2\mathrm{/}2}$ there. By [F4], $S_i$ is standard Gaussian.

The contrapositive of step 3.1 says that a non-Gaussian source column has at most one nonzero entry. Its squared norm is one because $O^{T}O=I$, so it has exactly one entry, whose square is one. This proves both assertions.

Work on the probability space of $S$ and define $T=A^{\mathrm{\prime }-1}AS$. The given equality implies by applying the fixed invertible linear map $A^{\mathrm{\prime }-1}$ to both laws. In particular the coordinates of $T$ are mutually independent, centered and have unit variance. Write $O=A^{\mathrm{\prime }-1}A$, so $T=OS$ on this probability space. Independence and centering give $\mathrm{Cov}(S)=\mathrm{Cov}(T)=I_p$, and covariance under this linear map yields $O{O}^T=I_p$. Since $O$ is square and invertible, also $O^{T}O=I_p$; it is orthogonal. This covariance calculation is the relative whitening reduction and requires neither a density nor moments beyond the given second moments.

By [F2], every non-Gaussian column of $O$ is a signed coordinate vector. Distinct such columns occupy distinct rows, because they are orthogonal. There are at least $p-1$ such columns. If there are $p$ they occupy all rows and $O$ is already a signed-permutation matrix. Otherwise there are exactly $p-1$, leaving one row $r$ and the one possible Gaussian column $i$. Its inner product with each occupied signed coordinate column is zero, so $O_{ji}=0$ on every occupied row $j$. Its norm is one, so $O_{ri}=1$ or $-1$. Thus $O$ is a signed-permutation matrix also in this case. For $p=1$ the same argument has zero occupied rows and its single norm-one column is $1$ or $-1$.

Put $Q=O^T$, again a signed-permutation matrix. From $O=A^{\mathrm{\prime }-1}A$ and $O^{-1}=O^T$ we get $A^{\mathrm{\prime }}=A{O}^T=AQ$. Also . These are exactly the asserted matrix and source-law conclusions.

Since $W,H$ are strictly positive, every current rate ${\lambda }_{ij}={\sum }_k{W}_{ik}{H}_{kj}$ is positive. For each row $i$, choose $j_0$ with $x_{i{j}_0}>0$ using [F1]. Then for every $k$, the numerator ${\sum }_j{H}_{kj}{x}_{ij}\mathrm{/}{\lambda }_{ij}$ has the positive term $H_{k{j}_0}{x}_{i{j}_0}\mathrm{/}{\lambda }_{i{j}_0}$, while its denominator ${\sum }_j{H}_{kj}$ is positive. Thus every updated $W_{ik}$ is positive.

The intermediate rates ${\tilde{\lambda }}_{ij}={\sum }_k{W}_{ik}^{\mathrm{n}\mathrm{e}\mathrm{w}}{H}_{kj}$ are positive. For each column $j$, choose $i_0$ with $x_{i_{0}j}>0$ using [F1]. For every $k$, the numerator ${\sum }_i{W}_{ik}^{\mathrm{n}\mathrm{e}\mathrm{w}}{x}_{ij}\mathrm{/}{\tilde{\lambda }}_{ij}$ has the positive term at $i_0$, and its denominator ${\sum }_i{W}_{ik}^{\mathrm{n}\mathrm{e}\mathrm{w}}$ is positive. Thus every updated $H_{kj}$ is positive. Repeating this argument inductively proves the claim for each finite sweep.