K-means, Gaussian mixtures, and soft assignments

Put $H:=A_1-A_0$ and $A_s:=A_0+sH$ for $0\le s\le 1$. For every nonzero $v$, $v^{\mathsf{T}}{A}_{s}v=(1-s)v^{\mathsf{T}}{A}_{0}v+s{v}^{\mathsf{T}}{A}_{1}v>0$, so each $A_s$ is symmetric positive definite. By [F1], choose an orthogonal eigenbasis $A_s=U\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_p)U^{\mathsf{T}}$ with every ${\lambda }_j>0$. The determinant laws [F9], including the triangular formula for the diagonal matrix, give $\det A_s=(\det U{)}^2{\prod }_j{\lambda }_j$ and $(\det U{)}^2=\det (U^{\mathsf{T}}U)=\det I=1$, so $\det A_s={\prod }_j{\lambda }_j>0$. Since the determinant is a polynomial by [F2] and $\log$ is differentiable, hence continuous, on $(0,\mathrm{\infty })$ by [F4], $g(s):=\log \det A_s$ is real and continuous on $[0,1]$.

On $0<s<1$, Jacobi's formula [F3], the scalar derivative [F4], and the chain rule [F6] give $g^{\mathrm{\prime }}(s)=\mathrm{tr}(A_s^{-1}H)$.

Differentiate the expression in step 2.1 along the same affine path. The inverse derivative [F5] and linearity of the finite trace [F7] give $g^{\mathrm{\prime }\mathrm{\prime }}(s)=-\mathrm{tr}(A_s^{-1}H{A}_s^{-1}H)$.

By [F1], form the symmetric positive-definite inverse square root $A_s^{-1\mathrm{/}2}$, and set $B:=A_s^{-1\mathrm{/}2}H{A}_s^{-1\mathrm{/}2}$. Since $H$ is symmetric, $B$ is symmetric; if $A_0\mathrm{\ne }{A}_1$, then $H\mathrm{\ne }0$ and invertibility of $A_s^{-1\mathrm{/}2}$ gives $B\mathrm{\ne }0$. Expanding the finite trace sum and interchanging its indices gives $\mathrm{tr}(XY)=\mathrm{tr}(YX)$; applying this identity yields $\mathrm{tr}(A_s^{-1}H{A}_s^{-1}H)=\mathrm{tr}(B^2)={\sum }_{i,j}{B}_{ij}^2>0$. Thus $g^{\mathrm{\prime }\mathrm{\prime }}(s)<0$ for every $0<s<1$. For any $0<a<b<1$, $g^{\mathrm{\prime }}$ is continuous on $[a,b]$ and differentiable on $(a,b)$ because $g$ is twice differentiable on $(0,1)$. By [F8] applied to $g^{\mathrm{\prime }}$, there is $c\in (a,b)$ with $g^{\mathrm{\prime }}(b)-g^{\mathrm{\prime }}(a)=g^{\mathrm{\prime }\mathrm{\prime }}(c)(b-a)<0$. Hence $g^{\mathrm{\prime }}$ is strictly decreasing on $(0,1)$.

For $0<t<1$, step 1.1 shows that $g$ is continuous on $[0,1]$ and differentiable on $(0,1)$, so [F8] applies on both $[0,t]$ and $[t,1]$. Their secant slopes are $g^{\mathrm{\prime }}(a)$ and $g^{\mathrm{\prime }}(b)$ for some $a\in (0,t)$ and $b\in (t,1)$. Step 4.1 makes $g^{\mathrm{\prime }}$ strictly decreasing, so the first slope is strictly larger than the second; rearranging yields $g(t)>(1-t)g(0)+tg(1)$. If $A_0=A_1$, the path is constant and equality holds.

Fix $i$ and write $R_i(\sigma ):={\sum }_{k\in M_i}{\pi }_k+{\sum }_{k\mathrm{\notin }{M}_i}{\pi }_k\exp (-(d_{ik}-d_i^{\ast })\mathrm{/}(2{\sigma }^2))$. Every difference in the second sum is strictly positive, so each exponential there tends to zero by [F1], since its exponent tends to $-\mathrm{\infty }$ as $\sigma \downarrow 0$; the sum is finite. Hence $R_i(\sigma )\to {\sum }_{k\in M_i}{\pi }_k>0$.

Factoring $\exp (-d_i^{\ast }\mathrm{/}(2{\sigma }^2))$ from the denominator of ${\rho }_{ik}^{(\sigma )}$ gives ${\rho }_{ik}^{(\sigma )}={\pi }_k\exp (-(d_{ik}-d_i^{\ast })\mathrm{/}(2{\sigma }^2))\mathrm{/}{R}_i(\sigma )$. If $k\in M_i$, its exponential is one and the limit is ${\pi }_k\mathrm{/}{\sum }_{j\in M_i}{\pi }_j$; if $k\mathrm{\notin }{M}_i$, its positive distance gap makes the numerator tend to zero while the denominator tends to a positive number. This proves the stated responsibility limits.

The same factorization gives $-2{\sigma }^2\log {\sum }_k{\pi }_k{e}^{-d_{ik}\mathrm{/}(2{\sigma }^2)}=d_i^{\ast }-2{\sigma }^2\log R_i(\sigma )$. The nearest terms give $R_i(\sigma )\ge {\sum }_{k\in M_i}{\pi }_k>0$, while $d_{ik}-d_i^{\ast }\ge 0$ and ${\sum }_k{\pi }_k=1$ give $R_i(\sigma )\le 1$. Thus $\log R_i(\sigma )$ stays between two finite constants, so the second term tends to zero. Summing over the finite set of observations gives ${\mathrm{\Psi }}_{\sigma }\to {\sum }_i{d}_i^{\ast }$.

For each labeling $z\in \{1,\dots ,K{\}}^n$, $d_{i,z_i}\ge d_i^{\ast }$ for every $i$, so ${\sum }_i{d}_{i,z_i}\ge {\sum }_i{d}_i^{\ast }$. Define $z_i:=\min M_i$, which exists because $M_i$ is a nonempty finite subset of $\{1,\dots ,K\}$. This labeling attains equality, proving the displayed minimum without any choice principle.

Set $A:={\mathrm{\Sigma }}^{1\mathrm{/}2}$, the symmetric positive-definite square root in the Statement, so $A^2=\mathrm{\Sigma }$. By [F8], write $A=U\mathrm{diag}(a_1,\dots ,a_p)U^{\mathsf{T}}$ with $U^{\mathsf{T}}U=I$ and $a_j>0$. The determinant laws [F11], including the triangular formula for the diagonal matrix, give $(\det U{)}^2=\det (U^{\mathsf{T}}U)=\det I=1$ and $\det A=(\det U{)}^2{\prod }_j{a}_j={\prod }_j{a}_j>0$. Also $\det \mathrm{\Sigma }=\det (A^2)=(\det A{)}^2$, so $\det A=\sqrt{\det \mathrm{\Sigma }}$. In the same basis, ${\mathrm{\Sigma }}^{-1}=A^{-2}=U\mathrm{diag}(a_1^{-2},\dots ,a_p^{-2})U^{\mathsf{T}}$ and $A^{\mathsf{T}}{\mathrm{\Sigma }}^{-1}A=I$.

By [F2] and [F4], the joint law of $Z$ is the product of its standard-normal marginals. Define $g(y):={\prod }_{j=1}^p{\varphi }_1(y_j)=(2\pi {)}^{-p\mathrm{/}2}{e}^{-\mathrm{\parallel }y{\mathrm{\parallel }}_2^2\mathrm{/}2}$. Tonelli [F5] and the product-measure identity [F6] give ${\int }_{{\mathbb{R}}^p}gd{\lambda }_p={\prod }_{j=1}^p{\int }_{\mathbb{R}}{\varphi }_{1}d{\lambda }_1=1$ and, for every Borel rectangle $R={\prod }_{j=1}^p{B}_j$, ${\int }_{R}gd{\lambda }_p={\prod }_{j=1}^p{\int }_{B_j}{\varphi }_{1}d{\lambda }_1$. Thus $gd{\lambda }_p$ and the product law of $Z$ are probability measures agreeing on rectangles. The rational half-open rectangles form a generating pi-system; $(-N,N{]}^p$ is an increasing exhaustion from it, and the common rectangle values make their measures equal and finite on every such box. Hence [F7] identifies the measures. AC supplies countable choice for the product-Lebesgue identity [F6].

Let $T(y)=\mu +Ay$. By step 1.1, $A$ is invertible, so $T$ is a $C^1$ diffeomorphism of ${\mathbb{R}}^p$ with $\mathrm{\mid }\det DT\mathrm{\mid }=\det A$. The linear measure-scaling result [F10] confirms this is the Lebesgue volume factor. The displayed function ${\varphi }_p$ is continuous and positive, so $f(x)={\mathbf{1}}_E(x){\varphi }_p(x;\mu ,\mathrm{\Sigma })$ is nonnegative Borel for every Borel $E\subseteq {\mathbb{R}}^p$. AC supplies the countable-choice assumption of [F9] by [F12]. Apply [F9] to $f$ and the change of variables by $T$ to obtain ${\int }_E{\varphi }_p(x;\mu ,\mathrm{\Sigma })dx={\int }_{T^{-1}(E)}{\varphi }_p(T(y);\mu ,\mathrm{\Sigma })\det Ady\mathrm{.}$ For $x=T(y)$, $x-\mu =Ay$ and $A^{\mathsf{T}}{\mathrm{\Sigma }}^{-1}A=A{A}^{-2}A=I$. Thus the exponent in [F1] becomes $-\mathrm{\parallel }y{\mathrm{\parallel }}_2^2\mathrm{/}2$, while $(\det \mathrm{\Sigma }{)}^{-1\mathrm{/}2}\det A=1$ by step 1.1. Hence ${\varphi }_p(T(y);\mu ,\mathrm{\Sigma })\det A=g(y)$. The last integral equals ${\int }_{T^{-1}(E)}g(y)dy=\mathbb{P}(X\in E)$ by step 1.2 and the definition of $X$. Thus ${\varphi }_p$ is the Lebesgue density of $X$.

By [F3], $X$ has law $N_p(\mu ,\mathrm{\Sigma })$. Step 2.1 identifies that law on every Borel set with integration against ${\varphi }_p$; taking $E={\mathbb{R}}^p$ gives ${\int }_{{\mathbb{R}}^p}{\varphi }_p(x;\mu ,\mathrm{\Sigma })dx=\mathbb{P}(X\in {\mathbb{R}}^p)=1$. Nonnegativity follows directly from [F1] and $\det \mathrm{\Sigma }>0$ in step 1.1.

Fix the current centers. For each observation, the old label is one of the available centers, so choosing a nearest center cannot increase its squared distance. If the label changes, its new distance is strictly smaller: an equal-distance tie would retain the old label by [F2]. Summing the pointwise inequalities gives nonincrease of $J$, strict when any label changes.

Fix a nonempty cluster $C$ and set $T(u)=(u{)}_{i\in C}$ from ${\mathbb{R}}^p$ to $({\mathbb{R}}^p{)}^C$, with product inner product, and $b=(x_i{)}_{i\in C}$. For $u,v\in {\mathbb{R}}^p$, $⟨Tu,Tv⟩=\mathrm{\mid }C\mathrm{\mid }⟨u,v⟩$, so $T^{\ast }T=\mathrm{\mid }C\mathrm{\mid }I$; also $⟨Tu,b⟩=⟨u,{\sum }_{i\in C}{x}_i⟩$, so $T^{\ast }b={\sum }_{i\in C}{x}_i$. Minimizing $\mathrm{\parallel }Tu-b\mathrm{\parallel }$ is equivalent to minimizing its square, since squaring is increasing on nonnegative reals. Thus [F3] gives the unique minimizer ${\bar{x}}_C=\mathrm{\mid }C{\mathrm{\mid }}^{-1}{\sum }_{i\in C}{x}_i$; uniqueness follows because $T$ is injective for nonempty $C$. Replacing a center by this mean cannot increase the cluster contribution and strictly decreases it if the old center differed. Empty clusters contribute zero and retain their centers.

A nonterminal iteration changes either the labels or at least one nonempty center. By steps 1.1 and 1.2, its full assignment-and-center update therefore strictly decreases $J$. For every $t\ge 1$, the state $(z^{(t)},{\mu }^{(t)})$ follows a center update: its nonempty centers are the means for $z^{(t)}$, so its objective value is determined by that labeling alone; empty centers do not affect it. Consequently, among these post-update states, a labeling cannot recur across nonterminal iterations, since its objective value would repeat despite strict decrease. There are only $K^n$ labelings, so after the first update there are at most $K^n-1$ further nonterminal iterations before a terminal iteration occurs. This proves finite termination without requiring the arbitrary initial centers to be cluster means or excluding $z^{(1)}=z^{(0)}$.

At termination the labels and centers are unchanged by the next iteration. The assignment rule therefore certifies that each label is nearest under the retention convention, while the center update certifies that every nonempty-cluster center is its sample mean. This is the stated coordinatewise fixed point.

By [F1] and [F2], every numerator and denominator in ${\gamma }_{ik}$ is finite and strictly positive. Hence ${\gamma }_{ik}>0$. Summing the responsibility formula over $k$ gives ${\sum }_k{\gamma }_{ik}=1$, so $N_k>0$ and ${\sum }_k{N}_k=n$. Thus all displayed weighted means and scatters are defined.

For each nonzero $v\in {\mathbb{R}}^p$, $v^{\mathsf{T}}{S}_{k}v=N_k^{-1}{\sum }_i{\gamma }_{ik}(v^{\mathsf{T}}(x_i-{\bar{x}}_k){)}^2$. If it were zero, positivity of each ${\gamma }_{ik}$ would force $v^{\mathsf{T}}{x}_i=v^{\mathsf{T}}{\bar{x}}_k$ for every $i$. That would put all observations in the proper affine hyperplane $\{x:v^{\mathsf{T}}x=v^{\mathsf{T}}{\bar{x}}_k\}$, contradicting their full affine span. Since $S_k$ is visibly symmetric, it is positive definite.

Expanding the component density in [F1,F2] in the definition of $Q$ gives $Q(\vartheta \mid \theta )={\sum }_k{N}_k\log {\rho }_k-\frac{np}{2}\log (2\pi )-\frac{1}{2}{\sum }_k(N_k\log \det V_k+{\sum }_i{\gamma }_{ik}(x_i-m_k{)}^{\mathsf{T}}{V}_k^{-1}(x_i-m_k))\mathrm{.}$ Here [F13] supplies $\log (\exp u)=u$ and the product and reciprocal laws. For each positive $a$, $\log (a^{-1\mathrm{/}2})=-\frac{1}{2}\log a$ follows by applying those laws to $(a^{-1\mathrm{/}2}{)}^2=1\mathrm{/}a$; applying the product law $p$ times gives $\log ((2\pi {)}^{-p\mathrm{/}2})=-\frac{p}{2}\log (2\pi )$. All logarithms are finite: weights are positive, and [F2] gives finite positive component densities for each positive-definite covariance, so their determinants are positive.

Define ${\alpha }_k=N_k\mathrm{/}n$, and for each subset $B\subseteq \{1,\dots ,K\}$ set $P(B):={\sum }_{k\in B}{\alpha }_k$ and $R(B):={\sum }_{k\in B}{\rho }_k$. Finite sums make these set functions additive on disjoint subsets; a disjoint sequence has at most $K$ nonempty members, so they are countably additive. Step 1.1 and the candidate weight constraint give total mass one, hence [F3] makes $P,R$ probability measures. Each ${\rho }_k>0$, so $P\ll R$. For every subset $B$, ${\sum }_{k\in B}({\alpha }_k\mathrm{/}{\rho }_k){\rho }_k=P(B)$, hence $dP\mathrm{/}dR(k)={\alpha }_k\mathrm{/}{\rho }_k$. The KL definition [F4] therefore gives $D(P\mathrm{\parallel }R)={\sum }_k{\alpha }_k\log ({\alpha }_k\mathrm{/}{\rho }_k)\ge 0$ by [A1,F5]. The weight-dependent part of $Q$ at $\alpha$ minus that at $\rho$ is ${\sum }_k{N}_k\log ({\alpha }_k\mathrm{/}{\rho }_k)=nD(P\mathrm{\parallel }R)\ge 0$, with equality only when $P=R$, hence ${\rho }_k={\alpha }_k$ for every $k$ by [F5].