Part 7 - The Expectation Maximization Algorithm & Linear Dimensionality Reduction

Expectation Maximization (EM)

EM solves a maximum likelihood problem of the form 1,

L(θ)=i=1Mlogp(x(i);θ)=i=1Mlogz(i)=1Kp(x(i),z(i);θ)L(\mathbf{\theta}) = \sum_{i = 1}^M \log p(\mathbf{x}^{(i)}; \mathbf{\theta}) = \sum_{i = 1}^M \log \sum_{\mathbf{z}^{(i)} = 1}^K p(\mathbf{x}^{(i)}, \mathbf{z}^{(i)}; \mathbf{\theta})
  • θ\mathbf{\theta}: Parameters of the probalistic model we are trying to find.
  • {x(i)}i=1M\{x^{(i)}\}_{i = 1}^M: Observed training examples.
  • {z(i)}i=1M\{z^{(i)}\}_{i = 1}^M: Unobserved latent variables. (e.g., in GMM, z(i)z^{(i)} indicates which one of the KK clusters x(i)x^{(i)} belongs to, which is unobserved.)

Jensen’s Inequality

Theorem 1 (Jensen’s Inequality)

Suppose f:RRf : \mathbb{R} \mapsto \mathbb{R} is concave, then for all probability distributions pp, we have,

f(E[x])E[f(x)].f(\mathbb{E} [\mathbf{x}]) \geq \mathbb{E} [f(\mathbf{x})].

Where the expectation is taken with respect to the random variable x\mathbf{x} drawn from the probability distribution pp.

The equality holds if and only if

  1. x\mathbf{x} is a constant, or,
  2. ff is an affine function (i.e., f(x)=aTx+bf(\mathbf{x}) = a^T \mathbf{x} + b).

EM Derivation

Let’s derive the EM algorithm for the maximum likelihood problem.

L(θ)=i=1Mlogp(x(i);θ)=i=1Mlogz(i)=1Kp(x(i),z(i);θ)=i=1Mlogz(i)=1Kq(z(i))p(x(i),z(i);θ)q(z(i))=i=1MlogEz(i)q[p(x(i),z(i);θ)q(z(i))]\begin{aligned} L(\mathbf{\theta}) &= \sum_{i = 1}^M \log p(\mathbf{x}^{(i)}; \mathbf{\theta}) \newline &= \sum_{i = 1}^M \log \sum_{\mathbf{z}^{(i)} = 1}^K p(\mathbf{x}^{(i)}, \mathbf{z}^{(i)}; \mathbf{\theta}) \newline &= \sum_{i = 1}^M \log \sum_{\mathbf{z}^{(i)} = 1}^K q(\mathbf{z}^{(i)}) \frac{p(\mathbf{x}^{(i)}, \mathbf{z}^{(i)}; \mathbf{\theta})}{q(\mathbf{z}^{(i)})} \newline &= \sum_{i = 1}^M \log \mathbb{E}_{\mathbf{z}^{(i)} \sim q} \left[ \frac{p(\mathbf{x}^{(i)}, \mathbf{z}^{(i)}; \mathbf{\theta})}{q(\mathbf{z}^{(i)})} \right] \newline \end{aligned}

Now we can apply Jensen’s inequality to the above equation,

L(θ)i=1MEz(i)q[logp(x(i),z(i);θ)q(z(i))]\begin{aligned} L(\mathbf{\theta}) &\geq \sum_{i = 1}^M \mathbb{E}_{\mathbf{z}^{(i)} \sim q} \left[ \log \frac{p(\mathbf{x}^{(i)}, \mathbf{z}^{(i)}; \mathbf{\theta})}{q(\mathbf{z}^{(i)})} \right] \newline \end{aligned}

Let’s go back to sum notation,

L(θ)=i=1Mz(i)=1Kq(z(i))logp(x(i),z(i);θ)q(z(i))=i=1Mz(i)=1Kq(z(i))logp(x(i),z(i);θ)i=1Mz(i)=1Kq(z(i))logq(z(i))\begin{aligned} L(\mathbf{\theta}) &= \sum_{i = 1}^M \sum_{\mathbf{z}^{(i)} = 1}^K q(\mathbf{z}^{(i)}) \log \frac{p(\mathbf{x}^{(i)}, \mathbf{z}^{(i)}; \mathbf{\theta})}{q(\mathbf{z}^{(i)})} \newline &= \sum_{i = 1}^M \sum_{\mathbf{z}^{(i)} = 1}^K q(\mathbf{z}^{(i)}) \log p(\mathbf{x}^{(i)}, \mathbf{z}^{(i)}; \mathbf{\theta}) \newline &- \sum_{i = 1}^M \sum_{\mathbf{z}^{(i)} = 1}^K q(\mathbf{z}^{(i)}) \log q(\mathbf{z}^{(i)}) \newline \end{aligned}

Let’s call this last expression for ell(θ)ell(\mathbf{\theta}), this is a lower bound of the original objective L(θ)L(\mathbf{\theta}). The equality holds when p(x(i),z(i);θ)q(z(i))\frac{p(\mathbf{x}^{(i)}, \mathbf{z}^{(i)}; \mathbf{\theta})}{q(\mathbf{z}^{(i)})} is a constant.

This can be achieved for q(z(i))=p(z(i)x(i);θ)q(z^{(i)}) = p(z^{(i)} | x^{(i)}; \theta).

The EM algorithm aims to optimize the lower bound ell(θ)ell(\mathbf{\theta}),

θ=argmaxθ(θ)=argmaxθi=1Mz(i)=1Kq(z(i))logp(x(i),z(i);θ)q(z(i))\mathbf{\theta}^{\star} = \underset{\mathbf{\theta}}{\arg \max} \ell(\mathbf{\theta}) = \underset{\mathbf{\theta}}{\arg \max} \sum_{i = 1}^M \sum_{\mathbf{z}^{(i)} = 1}^K q(\mathbf{z}^{(i)}) \log \frac{p(\mathbf{x}^{(i)}, \mathbf{z}^{(i)}; \mathbf{\theta})}{q(\mathbf{z}^{(i)})}

EM repeatedly performs the following two steps until convergence.

At tt-th iteration,

  1. E-step: For each index ii, we compute,
q(t)(z(i))=p(z(i)x(i);θ(t))q^{(t)}(z^{(i)}) = p(z^{(i)} | x^{(i)}; \mathbf{\theta}^{(t)})
  1. M-step: Compute,
θ(t+1)=argmaxθi=1Mz(i)=1q(t)(z(i))logp(x(i),z(i);θ)\mathbf{\theta}^{(t + 1)} = \underset{\mathbf{\theta}}{\arg \max} \sum_{i = 1}^M \sum_{\mathbf{z}^{(i)} = 1} q^{(t)}(\mathbf{z}^{(i)}) \log p(\mathbf{x}^{(i)}, \mathbf{z}^{(i)}; \mathbf{\theta})

In the E-step, we do not fill in the unobserved z(i)z^{(i)} with hard values, but find a posterior distribution q(z(i))q(z^{(i)}), given x(i)x^{(i)} and θ(t)\theta^{(t)}, i.e.,

q(t)(z(i))=p(z(i)x(i);θ(t)).q^{(t)}(z^{(i)}) = p(z^{(i)} | x^{(i)}; \mathbf{\theta}^{(t)}).

In the M-step, we maximize the lower bound ell(θ)ell(\mathbf{\theta}), while holding q(t)(z(i))q^{(t)}(z^{(i)}) fixed, which is computed from the E-step.

The M-step optimization can be done efficiently in most cases. For example, in GMM, we have a closed-form solution for all parameters.

EM Convergence

Assuming θ(t)\mathbf{\theta}^{(t)} and θ(t+1)\mathbf{\theta}^{(t + 1)} are the parameters from two successive iterations of EM, we have,

L(θ(t))=(1)i=1Mlogp(x(i);θ(t))=(2)iMlogz(i)=1Kq(z(i))p(x(i),z(i);θ(t))q(z(i))=(3)i=1Mz(i)=1Kq(t)(z(i))logp(x(i),z(i);θ(t))q(t)(z(i))(4)i=1Mz(i)=1Kq(t)(z(i))logp(x(i),z(i);θ(t))(5)i=1Mlogz(i)=1Kq(t)(z(i))p(x(i),z(i);θ(t))q(t)(z(i))=(6)L(θ(t+1))\begin{aligned} L(\mathbf{\theta}^{(t)}) &\stackrel{(1)}{=} \sum_{i = 1}^M \log p(\mathbf{x}^{(i)}; \mathbf{\theta}^{(t)}) \stackrel{(2)}{=} \sum_i^M \log \sum_{\mathbf{z}^{(i)} = 1}^K q(\mathbf{z}^{(i)}) \frac{p(\mathbf{x}^{(i)}, \mathbf{z}^{(i)}; \mathbf{\theta}^{(t)})}{q^(\mathbf{z}^{(i)})} \newline & \stackrel{(3)}{=} \sum_{i = 1}^M \sum_{\mathbf{z}^{(i)} = 1}^K q^{(t)}(\mathbf{z}^{(i)}) \log \frac{p(\mathbf{x}^{(i)}, \mathbf{z}^{(i)}; \mathbf{\theta}^{(t)})}{q^{(t)}(\mathbf{z}^{(i)})} \newline & \stackrel{(4)}{\leq} \sum_{i = 1}^M \sum_{\mathbf{z}^{(i)} = 1}^K q^{(t)}(\mathbf{z}^{(i)}) \log p(\mathbf{x}^{(i)}, \mathbf{z}^{(i)}; \mathbf{\theta}^{(t)}) \newline & \stackrel{(5)}{\leq} \sum_{i = 1}^M \log \sum_{\mathbf{z}^{(i)} = 1}^K q^{(t)}(\mathbf{z}^{(i)}) \frac{p(\mathbf{x}^{(i)}, \mathbf{z}^{(i)}; \mathbf{\theta}^{(t)})}{q^{(t)}(\mathbf{z}^{(i)})} \stackrel{(6)}{=} L(\mathbf{\theta}^{(t + 1)}) \end{aligned}
  1. By definition, this is the (log) likelihood of the data.
  2. By marginalization over z(i)z^{(i)} and multiplication an arbitrary distribution q(z(i))q(z^{(i)}) to both numerator and denominator inside log.
  3. By Jensen’s inequality where equality condition is satisfied by setting q(t)(z(i))=p(z(i)x(i);θ(t))q^{(t)}(z^{(i)}) = p(z^{(i)} | x^{(i)}; \mathbf{\theta}^{(t)}).
  4. By M-step of EM, where we maximize (3), holding q(t)(z(i))q^{(t)}(z^{(i)}) fixed.
  5. By Jensen’s inequality (in reverse order). Note that we have already updated θ\mathbf{\theta} from θ(t)\mathbf{\theta}^{(t)} to θ(t+1)\mathbf{\theta}^{(t + 1)}, q(t)(z(i))q^{(t)}(z^{(i)}) may now not satisfy the equality condition.
  6. By definition, this is the lower bound of the likelihood.

Hence, EM causes the likelihood to increase monotonically.

Remark

Note (EM as coordinate ascent)

One remark that we have to do is, if we define the EM as,

J(q,θ)=i=1Mz(i)=1Kq(z(i))logp(x(i),z(i);θ)q(z(i))J(q, \mathbf{\theta}) = \sum_{i = 1}^M \sum_{\mathbf{z}^{(i)} = 1}^K q(\mathbf{z}^{(i)}) \log \frac{p(\mathbf{x}^{(i)}, \mathbf{z}^{(i)}; \mathbf{\theta})}{q(\mathbf{z}^{(i)})}

We can view this as coordinate ascent on JJ, in which the E-step maximizes JJ with respect to qq, and the M-step maximizes JJ with respect to θ\mathbf{\theta}. One can easily prove this using Lagrangian multipliers!

Clustering Summary

  • Clustering Task
    • Given a set of input vectors D={x(i)}i=1M\mathcal{D} = \{x^{(i)}\}_{i = 1}^M, with x(i)RNx^{(i)} \in \mathbb{R}^N, group similar x(i)x^{(i)} into clusters.
      • Estimate a cluster center, representing the data points in that cluster.
      • Predict the cluster for a new data point.
  • Exhaustive clustering
    • Cluster shape: Arbitrary shape.
    • Principle: Minimize an assumed clustering criterion with brute-force search.
    • Pros: Optimal under the given clustering criterion.
    • Cons: Impractical to construct the clustering criterion, prohibitive to compute.
  • KK-means
    • Cluster shape: Circular.
    • Principle: Minimize distance to cluster center.
    • Pros: Simple and scalable (MiniBatchKMeans).
    • Cons: Sensitive to initialization, could get bad solutions due to local minima, need to choose KK.
  • Gaussian Mixture Model (GMM)
    • Cluster shape: Elliptical.
    • Principle: Maximum likelihood using expectation maximization.
    • Pros: Elliptical cluster shapes.
    • Cons: Sensitive to initialization, could get bad solutions due to local minima, need to choose KK.
  • Feature normalization
    • Feature normalization is typically required for clustering.
    • E.g., algorithms based on Eucledian distance (KK-means).

Dimensionality Reduction

Transform high-dimensional vectors into low-dimensional vectors. Dimensions in the low-dimensional data represent co-occuring features in the high-dimensional data. Dimensions in the low-dimensional data may have semantic meaning.

For example, in document analysis.

  • High-dimensional data: Bag-of-word vectors of documents.
  • Low-dimensional data: Each dimension represents similarity to a topic.

Reasons for Dimensionality Reduction

  • Preprocessing, makes the dataset easier to use.
  • Reduce computational cost of running machine learning algorithms.
  • Can be used to “de-noise” data by projecting to lower-dimensional space and then projecting back to the original high-dimensional space.
  • Makes the results easier to understand (e.g., visualization).

Dimensionality Reduction VS. Feature Selection

The goal of feature selection is to remove features that are not informative with respect to the class label. This obviously reduces the dimensionality of the feature space.

Dimensionality reduction can be used to find a meaningful lower-dimensional feature space even when there is information in each feature dimension so that none can be discarded.

Another important property of dimensionality reduction is that it is unsupervised.

While dimensionality reduction can be seen as a simplistic form of (data) compression, it is not equivalent to it, as the goal of compression is to reduce the expected code length (which is lower bounded by entropy) of the representation not only the dimensionality.

For example, in lossless compression, arithmetic coding encodes the entire data into a single number, an arbitrary-precision fraction qq where 0.0q<1.00.0 \leq q < 1.0.

Linear Dimensionality Reduction

In linear dimensionality reduction we project the original data onto a lower-dimensional hyperplane (e.g., line, plane).

I.e., move and rotate the coordinate axis of the data, then we represent the data with coordinates in the new component space.

Mathematically, this can be written as,

x(i)=k=1Kzk(i)bk\mathbf{x}^{(i)} = \sum_{k = 1}^K z_k^{(i)} \mathbf{b}_k

where bk\mathbf{b}_k is a basis vector and zk(i)Rz_k^{(i)} \in \mathbb{R} is the corresponding weight.

Connection to Linear Regression

If we focus on the jj-th entry of x(i)\mathbf{x}^{(i)}, we have,

xj(i)=k=1Kzk(i)bjkx_j^{(i)} = \sum_{k = 1}^K z_k^{(i)} b_{jk}

This expression can be seen as linear regression.

  • xj(i)x_j^{(i)} is the target.
  • zk(i)z_k^{(i)} for each kk are the weights.
  • bjkb_{jk} for each kk are the features.

Alternatively, we may view zk(i)z_k^{(i)} as feature and bjkb_{jk} as weights.

Unlike linear regression, we only know “targets”. We must learn both features and weights.

Matrix Formulation

Let XRM×N\mathbf{X} \in \mathbb{R}^{M \times N} be the data matrix, with one data case x(i)RN\mathbf{x}^{(i)} \in \mathbb{R}^N per row.

X=[(x(1))T(x(2))T(x(M))T]\mathbf{X} = \begin{bmatrix} -- & (\mathbf{x}^{(1)})^T & -- \newline -- & (\mathbf{x}^{(2)})^T & -- \newline & \vdots & \newline -- & (\mathbf{x}^{(M)})^T & -- \end{bmatrix}

Let ZRM×K\mathbf{Z} \in \mathbb{R}^{M \times K} be the loading matrix and BRK×N\mathbf{B} \in \mathbb{R}^{K \times N} be the factor matrix.

Z=[z1(1)z2(1)zK(1)z1(2)z2(2)zK(2)z1(M)z2(M)zK(M)],B=[(b1)T(b2)T(bK)T]\mathbf{Z} = \begin{bmatrix} z_1^{(1)} & z_2^{(1)} & \ldots & z_K^{(1)} \newline z_1^{(2)} & z_2^{(2)} & \ldots & z_K^{(2)} \newline \vdots & \ddots & \ddots & \vdots \newline z_1^{(M)} & z_2^{(M)} & \ldots & z_K^{(M)} \end{bmatrix}, \mathbf{B} = \begin{bmatrix} -- & (\mathbf{b}_1)^T & -- \newline -- & (\mathbf{b}_2)^T & -- \newline & \vdots & \newline -- & (\mathbf{b}_K)^T & -- \end{bmatrix}

With this we can express X\mathbf{X} as follows,

X=ZB.X = ZB.

However, most real world data will be subject to noise. If we assume that ϵRM×N\mathbf{\epsilon} \in \mathbb{R}^{M \times N} is a matrix of noise values from some probability distribution, we have,

X=ZB+ϵ.X = ZB + \mathbf{\epsilon}.

Learning Criterion

The learning problem for linear dimensionality reduction is to estimate values for both Z\mathbf{Z} and B\mathbf{B} given only the noisy observations of X\mathbf{X}.

One possible learning criterion is to minimize the sum of squared errors when reconstructing X\mathbf{X} from Z\mathbf{Z} and B\mathbf{B}.

This leads to,

argminZ,BXZBF2\underset{\mathbf{Z}, \mathbf{B}}{\arg \min} \Vert \mathbf{X} - \mathbf{ZB} \Vert_F^2

Where F\Vert \cdot \Vert_F is the Frobenius norm of a matrix, defined as,

AF=ijAij2\Vert \mathbf{A} \Vert_F = \sqrt{\sum_{ij} A_{ij}^2}

Alternating Least Squares

By leveraging the OLS solution for linear regression, we can estimate Z\mathbf{Z} and B\mathbf{B} using Alternating Least Squares (ALS).

Starting from some random initialization, ALS iterates between two steps until covergence.

  1. Assume Z\mathbf{Z} is fixed and optimize B\mathbf{B},
B(ZTZ)1ZTX\mathbf{B} \leftarrow (\mathbf{Z}^T \mathbf{Z})^{-1} \mathbf{Z}^T \mathbf{X}
  1. Assume B\mathbf{B} is fixed and optimize Z\mathbf{Z},
ZT(BBT)1BXT\mathbf{Z}^T \leftarrow (\mathbf{B} \mathbf{B}^T)^{-1} \mathbf{B} \mathbf{X}^T

Lack of Uniqueness for Optimal Parameters

Suppose we run the ALS algorithm to convergence and obtain optimal parameters Z\mathbf{Z}^{\star} and B\mathbf{B}^{\star} such that,

=XZBF2\ell^{\star} = \Vert \mathbf{X} - \mathbf{Z}^{\star} \mathbf{B}^{\star} \Vert_F^2

Let RRK×K\mathbf{R} \in \mathbb{R}^{K \times K} be an arbitrary invertible matrix.

A K×KK \times K matrix R\mathbf{R} is invertible, if there exists a K×KK \times K matrix S\mathbf{S} such that RS=SR=I\mathbf{R} \mathbf{S} = \mathbf{S} \mathbf{R} = \mathbf{I}. Which we also can denote as R1=S\mathbf{R}^{-1} = \mathbf{S}.

We obtain a different set of parameters Z~=ZR\mathbf{\tilde{Z}} = \mathbf{Z}^{\star} \mathbf{R} and B~=R1B\mathbf{\tilde{B}} = \mathbf{R}^{-1} \mathbf{B}^{\star}, with the same optimal value,

=XZ(I)BF2=XZ(RR1)BF2=XZ~B~F2\ell^{\star} = \Vert \mathbf{X} - \mathbf{Z}^{\star} (\mathbf{I}) \mathbf{B}^{\star} \Vert_F^2 = \Vert \mathbf{X} - \mathbf{Z}^{\star} (\mathbf{R} \mathbf{R}^{-1}) \mathbf{B}^{\star} \Vert_F^2 = \Vert \mathbf{X} - \mathbf{\tilde{Z}} \mathbf{\tilde{B}} \Vert_F^2

We can obtain the global optimal solution(s) and make them unique by specifying additional criteria.

Singular Value Decomposition (SVD)

Let X\mathbf{X} be a M×NM \times N matrix, with MNM \geq N. It can be facotrized as,

X=U(Σ0)VT\mathbf{X} = \mathbf{U} \begin{pmatrix} \mathbf{\Sigma} \newline \mathbf{0} \end{pmatrix} \mathbf{V}^T

Where URM×M\mathbf{U} \in \mathbb{R}^{M \times M} and VRN×N\mathbf{V} \in \mathbb{R}^{N \times N} are orthogonal, i.e.,

UTU=UUT=IM,VTV=VVT=IN\mathbf{U}^T \mathbf{U} = \mathbf{U} \mathbf{U}^T = \mathbf{I}_M, \quad \mathbf{V}^T \mathbf{V} = \mathbf{V} \mathbf{V}^T = \mathbf{I}_N

Columns of U\mathbf{U} and V\mathbf{V} are called left and right singular vectors of X\mathbf{X}, respectively.

ΣRM×N\mathbf{\Sigma} \in \mathbb{R}^{M \times N} is a diagonal matrix,

Σ=diag(σ1,σ2,,σN),σ1σ2σN0\mathbf{\Sigma} = \text{diag}(\sigma_1, \sigma_2, \ldots, \sigma_N), \sigma_1 \geq \sigma_2 \geq \ldots \geq \sigma_N \geq 0

σi\sigma_i’s are called singular values of X\mathbf{X}.

SVD visualization of the matrices.
SVD visualization of the matrices.
SVD visualization of the different operations.
SVD visualization of the different operations.

In Figure 2, we see four different operations.

  • Upper Left: The unit disc with the two canonical unit vectors.
  • Upper Right: Transformed with M\mathbf{M}.
  • Lower Left: The action of VT\mathbf{V}^T. This is just a rotation.
  • Lower Right: The action of ΣVT\mathbf{\Sigma}\mathbf{V}^T. Σ\mathbf{\Sigma} scales vertically and horizontally.

Reduced-Form SVD

If only K<min{M,N}K < \min\{M, N\} singular values are non-zero, the SVD XRM×N\mathbf{X} \in \mathbb{R}^{M \times N} can be represented in reduced form as follows,

X=UΣKVT,\mathbf{X} = \mathbf{U} \mathbf{\Sigma}_K \mathbf{V}^T,

which we can write as,

X=k=1KσkukvkT.\mathbf{X} = \sum_{k = 1}^K \sigma_k \mathbf{u}_k \mathbf{v}_k^T.

Where

  • U=[u1,u2,,uM]RM×K\mathbf{U} = [\mathbf{u}_1, \mathbf{u}_2, \ldots, \mathbf{u}_M] \in \mathbb{R}^{M \times K}.
    • UTU=IK\mathbf{U}^T \mathbf{U} = \mathbf{I}_K.
  • V=[v1,v2,,vN]RN×K\mathbf{V} = [\mathbf{v}_1, \mathbf{v}_2, \ldots, \mathbf{v}_N] \in \mathbb{R}^{N \times K}.
    • VTV=IK\mathbf{V}^T \mathbf{V} = \mathbf{I}_K.
  • ΣK=diag(σ1,σ2,,σK)RK×K\mathbf{\Sigma}_K = \text{diag}(\sigma_1, \sigma_2, \ldots, \sigma_K) \in \mathbb{R}^{K \times K}.
    • σ1σ2σK0\sigma_1 \geq \sigma_2 \geq \ldots \geq \sigma_K \geq 0.
  • ukvkTRM×N\mathbf{u}_k \mathbf{v}_k^T \in \mathbb{R}^{M \times N} is the product of a column vector uk\mathbf{u}_k and a row vector vkT\mathbf{v}_k^T.
    • It has rank 1.
    • X\mathbf{X} is a weighted summation of KK rank-1 matrices.

Eckart-Young-Mirsky Theorem

Theorem 2 (Eckart-Young-Mirsky Theorem)

Given a M×NM \times N matrix X\mathbf{X} of rank Rmin{M,N}R \leq \min \{M, N\} and its singular value decomposition X=UΣRVT\mathbf{X} = \mathbf{U} \mathbf{\Sigma}_R \mathbf{V}^T with singular values σ1σ2σR>0\sigma_1 \geq \sigma_2 \geq \ldots \geq \sigma_R > 0 and the rest being zero, then among all M×NM \times N matrices of lower rank KRK \leq R, the best approximation is Y=UΣKVT\mathbf{Y}^{\star} = \mathbf{U} \mathbf{\Sigma}_K \mathbf{V}^T, where ΣK\mathbf{\Sigma}_K is the diagonal matrix with singular values σ1,σ2,,σK\sigma_1, \sigma_2, \ldots, \sigma_K in the sense that 2,

XYF2=min{XYF2;YRM×N,rankYK}\Vert \mathbf{X} - \mathbf{Y}^{\star} \Vert_F^2 = \min \{\Vert \mathbf{X} - \mathbf{Y} \Vert_F^2; \mathbf{Y} \in \mathbb{R}^{M \times N}, \text{rank} \mathbf{Y} \leq K\}

SVD provdies a unique solution to minimum Frobenius norm linear dimensionality reduction.

References

  1. Arthur P. Dempster, Nan M. Laird, and Donald B. Rubin (1977) “Maximum Likelihood from Incomplete Data via the EM Algorithm”. Journal of the Royal Statistical Society: Series B 39, 1-38. Link
  2. Carl Eckart and Gale Young (1936) “The Approximation of One Matrix by Another of Lower Rank”. Psychometrika 1, 211-218. Link