Tutorials › Statistical Methods for Single-Cell RNA Sequencing › Diffusion Maps from Eigenvectors

Pseudotime and Trajectory Inference · Tutorial 372 of 700

Diffusion Maps from Eigenvectors

Learn how to transform normalized graph eigenvectors into diffusion coordinates, scale them by eigenvalues over diffusion time, and check the resulting coordinates.

Advanced 11 min read

What You'll Learn

  • Define diffusion-map coordinates from nonstationary eigenvectors and eigenvalues of a normalized transition operator.
  • Relate eigenvectors of the symmetrically normalized operator to right eigenvectors of the transition matrix.
  • Calculate diffusion coordinates for small connected cell graphs at specified diffusion times.
  • Check graph connectivity, degree positivity, eigenvector normalization, and diffusion-time assumptions.
  • Recognize how eigenvector sign, repeated eigenvalues, and mode truncation affect coordinate interpretation.

From a Cell Graph to Coordinates

A cell graph represents cells as vertices and similarities between cells as edge weights. A transition operator converts those weights into one-step movement probabilities. Diffusion maps use eigenvectors of that operator to assign each cell coordinates: cells with similar patterns of transitions can have similar coordinates, even when they are not directly connected by a single edge.

The central calculation is compact, but its details matter. Coordinates depend on which operator’s eigenvectors are used, how those eigenvectors are normalized, which modes are retained, and the diffusion time. We will use one explicit convention throughout so that the numbers are reproducible.

Definition: Let \(W\) be a symmetric matrix of nonnegative edge weights for a graph with \(n\) cells. Define the degree of cell \(i\) as \(d_i=\sum_j W_{ij}\), let \(D\) be the diagonal matrix with entries \(d_i\), and let \(Z=\sum_i d_i\). Assume every \(d_i>0\). The row-stochastic transition matrix and its symmetrically normalized counterpart are $$ P=D^{-1}W,\qquad S=D^{-1/2}WD^{-1/2}. $$ The stationary weights are \(\pi_i=d_i/Z\). For an eigenvector \(u_\ell\) of \(S\), normalized to have Euclidean norm one, define \(\phi_\ell(i)=\sqrt{Z}\,u_\ell(i)/\sqrt{d_i}\). The diffusion-map coordinate for cell \(i\), mode \(\ell\), and integer diffusion time \(t\geq 1\) is \(x_{t,\ell}(i)=\lambda_\ell^t\phi_\ell(i)\), where \(S u_\ell=\lambda_\ell u_\ell\). The coordinate vector retains a chosen set of nonstationary modes, ordinarily \(\ell=1,\ldots,q\), excluding the stationary mode \(\ell=0\).

These definitions give the same coordinates as using right eigenvectors of \(P\) with a particular normalization. Indeed, set \(\phi_\ell=\sqrt{Z}D^{-1/2}u_\ell\). Since \(P=D^{-1/2}SD^{1/2}\), we have \(P\phi_\ell=\lambda_\ell\phi_\ell\). Also, \(\sum_i\pi_i\phi_\ell(i)\phi_m(i)=1\) if \(\ell=m\), and \(0\) otherwise. Thus the right eigenvectors are orthonormal under the stationary weights, not necessarily under the ordinary unweighted Euclidean inner product.

For a connected graph, the eigenvalue \(1\) is simple, and its right eigenvector under this convention is \(\phi_0(i)=1\) for every cell. This constant mode describes no differences among cells, so it is omitted from the coordinates. If the graph has self-loops with positive weight at every vertex, then the chain is aperiodic and all other eigenvalues have absolute value less than one. Self-loops are not required to define the coordinates, but on a bipartite graph without them an eigenvalue of \(-1\) can occur.

Why Eigenvalues Scale the Coordinates

The exponent \(t\) represents repeated transitions. If a mode of the transition operator has eigenvalue \(\lambda_\ell\), applying the operator \(t\) times multiplies that mode by \(\lambda_\ell^t\). The coordinate therefore combines the mode’s cell-specific eigenvector values with its contribution after \(t\) transitions. This is a coordinate construction; the later spectral treatment of diffusion distance is not needed to calculate it.

Derivation: The symmetric normalization and transition matrix are similar: $$ P=D^{-1}W=D^{-1/2}SD^{1/2}. $$ If \(S u_\ell=\lambda_\ell u_\ell\), then $$ P(D^{-1/2}u_\ell)=D^{-1/2}S u_\ell =\lambda_\ell(D^{-1/2}u_\ell). $$ Thus \(D^{-1/2}u_\ell\) is a right eigenvector of \(P\). Multiplying it by \(\sqrt{Z}\) gives the stationary-weight normalization \(\phi_\ell\). Repeated application gives \(P^t\phi_\ell=\lambda_\ell^t\phi_\ell\), which justifies the diffusion-time factor in each coordinate.

This argument also explains why symmetric normalization is useful here: \(S\) is a real symmetric matrix, so it has real eigenvalues and an orthonormal basis of eigenvectors. The transition matrix \(P\) need not be symmetric, but it has the corresponding right eigenvectors obtained by the degree adjustment above. Using eigenvectors of \(S\) directly as though they were right eigenvectors of \(P\) generally gives different, incorrectly scaled cell values.

Worked Example: Coordinates on a Three-Cell Path

Consider a hypothetical connected graph with three cells in a path and weight matrix $$ W=\begin{pmatrix}1&1&0\\1&2&1\\0&1&1\end{pmatrix}. $$ The diagonal entries are self-loop weights. The degrees are \(d=(2,4,2)\), so \(Z=8\), and the stationary weights are \(\pi=(1/4,1/2,1/4)\). The transition matrix is $$ P=\begin{pmatrix}1/2&1/2&0\\1/4&1/2&1/4\\0&1/2&1/2\end{pmatrix}. $$

The symmetric operator is $$ S=\begin{pmatrix}1/2&1/(2\sqrt{2})&0\\1/(2\sqrt{2})&1/2&1/(2\sqrt{2})\\0&1/(2\sqrt{2})&1/2\end{pmatrix}. $$ Its eigenvalues are \(1,1/2,0\), with corresponding unit eigenvectors $$ u_0=\begin{pmatrix}1/2\\1/\sqrt{2}\\1/2\end{pmatrix},\quad u_1=\begin{pmatrix}1/\sqrt{2}\\0\\-1/\sqrt{2}\end{pmatrix},\quad u_2=\begin{pmatrix}1/2\\-1/\sqrt{2}\\1/2\end{pmatrix}. $$ Each vector has squared entries summing to one. Transforming to the stationary-weight normalization gives $$ \phi_0=(1,1,1),\qquad \phi_1=(\sqrt{2},0,-\sqrt{2}),\qquad \phi_2=(1,-1,1). $$

For diffusion time \(t=1\), keep both nonstationary modes. Cell 1 has coordinates $$ x_{1,1}(1)=(1/2)(\sqrt{2})=\sqrt{2}/2\approx0.7071,\qquad x_{1,2}(1)=0(1)=0. $$ Cell 2 has coordinates \((0,0)\), and cell 3 has coordinates \((-0.7071,0)\). The zero second coordinate occurs because the corresponding eigenvalue is zero, not because the second eigenvector has zero entries at every cell. The graph is connected, all degrees are positive, and the self-loops ensure the nonstationary modes have absolute eigenvalues below one.

Time, Retained Modes, and Computation

For a fixed mode, changing \(t\) changes its contribution by a power of its eigenvalue. In the three-cell example, at \(t=2\), the first mode’s values are multiplied by \((1/2)^2=1/4\), while the second mode remains zero because \(0^2=0\). The coordinates become \((\sqrt{2}/4,0)\), \((0,0)\), and \((-\sqrt{2}/4,0)\), respectively.

In practice, one usually retains only a subset of nonstationary modes. The number \(q\) is a choice: all available modes can be retained for a small graph, while larger graphs require a computationally feasible subset. Each retained mode contributes one coordinate, so a cell’s map is the vector \((\lambda_1^t\phi_1(i),\ldots,\lambda_q^t\phi_q(i))\). The formula defines these coordinates; it does not by itself prescribe a universal choice of \(q\) or \(t\).

Worked Example: Changing Diffusion Time

Use the three-cell graph above, but retain only mode 1 and set \(t=2\). The eigenvalue is \(\lambda_1=1/2\), and the eigenfunction values are \((\sqrt{2},0,-\sqrt{2})\). Apply the coordinate formula separately to each cell: $$ x_{2,1}(1)=(1/2)^2\sqrt{2}=\sqrt{2}/4\approx0.3536,\quad x_{2,1}(2)=(1/2)^2(0)=0,\quad x_{2,1}(3)=(1/2)^2(-\sqrt{2})=-\sqrt{2}/4\approx-0.3536. $$

The calculation is consistent with the graph assumptions and the selected mode. Compared with \(t=1\), the nonzero coordinate values are halved, since \((1/2)^2/(1/2)=1/2\). This is a direct consequence of the eigenvalue power; it is not an additional rescaling based on the cell’s degree.

Worked Example: A Two-Cell Weighted Graph

Consider a hypothetical graph with $$ W=\begin{pmatrix}2&1\\1&2\end{pmatrix}. $$ Both degrees equal 3, so \(Z=6\), \(\pi=(1/2,1/2)\), and \(P=S=W/3\). The constant eigenvector has eigenvalue 1. The nonstationary unit eigenvector of \(S\) is \(u_1=(1,-1)/\sqrt{2}\), with eigenvalue \(\lambda_1=1/3\). Transforming gives \(\phi_1=\sqrt{6}\,u_1/\sqrt{3}=(1,-1)\).

At \(t=2\), the coordinates are $$ x_{2,1}(1)=(1/3)^2(1)=1/9\approx0.1111,\qquad x_{2,1}(2)=(1/3)^2(-1)=-1/9\approx-0.1111. $$ The graph is connected, its degrees are positive, and its self-loops ensure aperiodicity. The stationary mode \((1,1)\) is excluded, leaving one coordinate per cell.

Assumptions, Checks, and Failure Modes

Before interpreting computed coordinates, check that the graph and eigenvectors match the definitions. The weights should be nonnegative and symmetric for the stated symmetric-operator construction, and each cell must have positive degree. Connectivity matters: if the graph is disconnected, the eigenvalue \(1\) is not simple, and there is more than one stationary component. A single constant mode then cannot represent all component-wise stationary behavior.

Numerical calculations should also be checked for eigenvector residuals and normalization. For each computed pair, verify that \(S u_\ell\) is close to \(\lambda_\ell u_\ell\), that \(u_\ell^\mathsf{T}u_\ell\) is close to 1, and that transformed eigenvectors satisfy the chosen stationary-weight normalization. Small deviations can occur from numerical rounding; large deviations suggest a mismatched operator, an indexing error, or an inaccurate eigensolution.

  • Sign is arbitrary. Replacing \(u_\ell\) by \(-u_\ell\) also replaces every coordinate in that dimension by its negative. This does not change the underlying mode, but it reverses the plotted axis.
  • Repeated eigenvalues allow rotations. Within an eigenspace with the same eigenvalue, a different orthonormal basis can give different coordinates. Interpretations that depend on one specific axis in such a subspace may be unstable.
  • Mode truncation loses information. Keeping only \(q\) modes omits the others. A low-dimensional plot is therefore not the full set of coordinates.
  • Disconnected graphs need separate care. Multiple stationary modes can make coordinates from a single connected-graph interpretation misleading.
  • Negative eigenvalues can alternate in sign across times. If a nonstationary eigenvalue is negative, \(\lambda_\ell^t\) changes sign when integer \(t\) changes from even to odd. Do not treat such a sign change as an eigenvector error without checking the graph and operator.

Choosing the Eigenvector Representation

There are two equivalent computational routes under the assumptions above. One can diagonalize the symmetric matrix \(S\), then transform each unit eigenvector by \(\phi_\ell=\sqrt{Z}D^{-1/2}u_\ell\). Alternatively, one can calculate right eigenvectors of \(P\) and normalize them to have unit norm under the stationary weights \(\pi_i\). The first route takes advantage of symmetry and gives a standard Euclidean orthonormal eigenbasis before transformation. The second route works directly with the transition matrix but must respect its weighted normalization.

Neither route means that ordinary principal component coordinates or unscaled eigenvectors of \(S\) are automatically diffusion coordinates. Principal component analysis uses a different operator and construction. Unscaled eigenvectors of \(S\) also differ cell by cell from the right eigenvectors of \(P\), due to the degree adjustment. The coordinate formula requires both the correct transformed eigenvector and its eigenvalue power.

Worked Example: Detecting a Missing Degree Adjustment

In the three-cell example, the first unit eigenvector of \(S\) is \(u_1=(1/\sqrt{2},0,-1/\sqrt{2})\). If these entries are incorrectly used as the right eigenvector values, the first cell’s \(t=1\) coordinate would be reported as \((1/2)(1/\sqrt{2})\approx0.3536\). But the stated convention requires \(\phi_1(1)=\sqrt{2}\), so the correct coordinate is \((1/2)\sqrt{2}\approx0.7071\).

The discrepancy comes from omitting \(\sqrt{Z}/\sqrt{d_i}\). Here that factor is \(\sqrt{8}/\sqrt{2}=2\) for cell 1, so the correctly transformed eigenvector value is twice the untransformed value. Checking the right-eigenvector equation confirms the distinction: \(u_1\) is an eigenvector of \(S\), while \(\phi_1\) is the corresponding right eigenvector of \(P\).

Key takeaway: With this normalization, compute each retained diffusion coordinate as \(x_{t,\ell}(i)=\lambda_\ell^t\sqrt{Z}\,u_\ell(i)/\sqrt{d_i}\). Exclude the constant stationary mode, state the chosen diffusion time and retained modes, and check the graph assumptions and eigenvector normalization before interpreting the result.

Check Your Understanding

Work through the questions, then compare your results with the solutions.

  1. For the three-cell path example, what is the stationary coordinate mode, and why is it excluded from the diffusion map?
  2. For a mode with eigenvalue \(1/2\) and transformed eigenvector value \(3\) at a cell, calculate its coordinate at \(t=3\).
  3. Why can \(D^{-1/2}u_\ell\) be used as a right eigenvector of \(P\), and what extra factor gives the normalization used here?
  4. What issue arises if a graph has two disconnected components?

Solutions. (1) The stationary mode has \(\lambda_0=1\) and \(\phi_0(i)=1\) at every cell. It is excluded because it gives the same value to every cell and therefore does not distinguish cells. (2) \(x_{3,\ell}(i)=(1/2)^3(3)=3/8=0.375\). (3) Similarity gives \(P=D^{-1/2}SD^{1/2}\), so \(P(D^{-1/2}u_\ell)=\lambda_\ell D^{-1/2}u_\ell\). Multiplication by \(\sqrt{Z}\) gives the stationary-weight normalization \(\phi_\ell=\sqrt{Z}D^{-1/2}u_\ell\). (4) The eigenvalue \(1\) is no longer simple: each component contributes a stationary mode. The connected-graph interpretation that excludes one constant mode is therefore insufficient without handling the components.

Clinical Biostats