import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib_inline
matplotlib_inline.backend_inline.set_matplotlib_formats('svg')
import seaborn as sns
sns.set_context("paper")
sns.set_style("ticks");

Proper Orthogonal Decomposition

Proper orthogonal decomposition (POD) is essentially PCA in Hilbert spaces. Let \(H\) be such a Hilbert space with inner product \(\langle \cdot, \cdot \rangle\) and norm \(\|\cdot\|\).

Let \(\{u_i\}_{i=1}^n\) be a set of \(n\) functions in \(H\). These are our observations. For example, \(u_i\) could be flow fields (velocity and pressure) at different times, or the temperature field in a heat transfer problem.

POD finds an optimal set of orthonormal basis functions \(\{\phi_i\}_{i=1}^m\). The set is optimal in the sense that the error in approximating the observations by a linear combination of the basis functions is minimized.

We will develop the idea using one basis function \(\phi\). The projection of a function \(u\) onto \(\phi\) is given by:

\[ \text{Proj}(u) = \langle u, \phi \rangle \phi. \]

The error in the projection is given by:

\[ \text{Err}(u) = u - \text{Proj}(u) = u - \langle u, \phi \rangle \phi. \]

Now, define the empirical expectation operator:

\[ \langle F[u] \rangle = \frac{1}{n} \sum_{i=1}^n F[u_i]. \]

Please try not to confuse the empirical expectation operator with the inner product. The empirical expectation operator satisfies the same properties as the expectation operator in probability theory.

Now, we can write down the expected reconstruction error:

\[ L(\phi) = \langle \|\text{Err}(u)\|^2 \rangle = \langle \|u - \langle u, \phi \rangle \phi\|^2 \rangle. \]

It is straightforward to show that minimizing \(J(\phi)\) is equivalent to maximizing the magnitude of the projections:

\[ V(\phi) = \langle \|\text{Proj}(u)\|^2 \rangle = \langle \|\langle u, \phi \rangle \|^2 \rangle. \]

Using the method of Lagrange multipliers, we write down:

\[ L(\phi) = \langle \|\langle u, \phi \rangle \|^2 \rangle + \lambda \left( 1 - \|\phi\|^2 \right). \]

The first variation of \(L(\phi)\) in an arbitrary direction \(\eta\) must be zero:

\[ 0 = \frac{\delta L(\phi)}{\delta \eta} = 2\langle \langle \phi, u\rangle\langle u, \eta\rangle\rangle - 2\lambda \langle \phi, \eta\rangle = 2\Big\langle\langle \langle \phi, u\rangle u\rangle - \lambda \phi,\eta\Big\rangle. \]

So:

\[ \langle \langle \phi, u\rangle u\rangle - \lambda \phi, \]

or

\[ \langle \langle \phi, u\rangle u\rangle = \lambda \phi. \]

We see that we get an operator eigenvalue problem. The operator is \(R:H\to H\) defined by:

\[ R\phi = \langle \phi, u\rangle u. \]

This operator is linear, self-adjoint, and compact. So, if we are in a separable Hilbert space, we can find an orthonormal basis of eigenfunctions \(\{\phi_i\}_{i=1}^\infty\) and eigenvalues \(\{\lambda_i\}_{i=1}^\infty\). The eigenvalues are non-negative and the eigenfunctions are ordered by decreasing eigenvalues. The biggest eigenvalue corresponds to the best basis function, the second biggest eigenvalue corresponds to the second best basis function, and so on.

Spectral decomposition of the \(R\) operator

The \(R\) operator, can be decomposed into a sum of eigenfunctions and eigenvalues:

\[ R\phi = \sum_{i=1}^\infty \lambda_i \langle \phi, \phi_i \rangle \phi_i. \]

This is called the spectral decomposition of the operator.

The kernel of the \(R\) operator

Suppose we work in \(H = L^2(\Omega)\), where \(\Omega\) is a domain in \(\mathbb{R}^d\). Then, the inner product is:

\[ \langle u, v \rangle = \int_\Omega u(x) v(x) \, dx. \]

The operator \(R\) can be expressed in terms of a kernel \(R: \Omega \times \Omega \to \mathbb{R}\):

\[ R(x, y) := \langle u(x)u(y) \rangle. \]

Notice that the kernel involves an empirical average. The value \(R(x,y)\) is the empirical correlation between \(u(x)\) and \(u(y)\). To see where the kernel comes from, observe that:

\[\begin{split} \begin{align*} (R\phi)(x) &= \langle \langle \phi, u \rangle u(x)\rangle\\ &= \Big\langle \int_{\Omega}\phi(y)u(y)\,dy \cdot u(x)\Big\rangle\\ &= \int_{\Omega}\phi(y)\langle u(y)u(x)\rangle\, dy\\ &= \int_{\Omega}\phi(y)R(y,x)\, dy. \end{align*} \end{split}\]

The kernel \(R(x,y)\) is symmetric, postive-semi-definite. It also has a spectral decomposition:

\[ R(x,y) = \sum_{i=1}^\infty \lambda_i \phi_i(x)\phi_i(y). \]

Computing the POD using SVD

To compute the POD, we typically work with a discretized version of the physical fields. For example, if \(H = L^2(\Omega)\), we can discretize the domain \(\Omega\) into a grid of \(N\) points:

\[ u_i = \left(u_i(x_1), u_i(x_2), \ldots, u_i(x_N)\right). \]

Then, we approximate the inner product with a sum:

\[ \langle u, v \rangle \approx \sum_{i=1}^N u_i v_i. \]

The kernel \(R\), now becomes a \(N\times N\) matrix \(\mathbf{R}\) with elements:

\[ R_{kl} = \langle u(x_k)u(x_l)\rangle = \frac{1}{n}\sum_{i=1}^n u_i(x_k)u_i(x_l). \]

If we make the data matrix:

\[\begin{split} \mathbf{X} = \begin{bmatrix} u_1(x_1) & u_1(x_2) & \cdots & u_1(x_N)\\ u_2(x_1) & u_2(x_2) & \cdots & u_2(x_N)\\ \vdots & \vdots & \ddots & \vdots\\ u_n(x_1) & u_n(x_2) & \cdots & u_n(x_N) \end{bmatrix}, \end{split}\]

we observe that the kernel is:

\[ \mathbf{R} = \frac{1}{n}\mathbf{X}^T\mathbf{X}. \]

We can do SVD on \(\mathbf{X}\):

\[ \mathbf{X} = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^T, \]

and substitute in \(\mathbf{R}\) to find:

\[ \mathbf{R} = \frac{1}{n}\mathbf{V}\mathbf{\Sigma}^2\mathbf{V}^T. \]

From this, we see that the eigenvectors of \(\mathbf{R}\) are the right singular vectors of \(\mathbf{X}\), and the eigenvalues are the squares of the singular values of \(\mathbf{X}\) divided by \(n\):

\[ \boldsymbol{\phi}_i = \mathbf{v}_i, \]

and

\[ \lambda_i = \frac{1}{n}\sigma_i^2. \]

The projections (known as POD modal coefficients) are given by:

\[ \mathbf{Z} = \mathbf{X}\mathbf{V} = \mathbf{U}\mathbf{\Sigma}. \]

Note that for other Hilbert spaces, the procedure is similar - but not identical.