Spectral pollution

It is possible for the Rayleigh–Ritz method to produce values which do not converge to actual values in the spectrum of the operator as the truncation gets large. These values are known as spectral pollution.^[4][6][9] In some cases (such as for the Schrödinger equation), there is no approximation which both includes all eigenvalues of the equation, and contains no pollution.^[10]

The spectrum of the compression (and thus pollution) is bounded by the numerical range of the operator; in many cases it is bounded by a subset of the numerical range known as the essential numerical range.^[11][12]

Example

The matrix $${\displaystyle A={\begin{bmatrix}2&0&0\\0&2&1\\0&1&2\end{bmatrix}}}$$ has eigenvalues ${\displaystyle 1,2,3}$ and the corresponding eigenvectors $${\displaystyle \mathbf {x} _{\lambda =1}={\begin{bmatrix}0\\1\\-1\end{bmatrix}},\quad \mathbf {x} _{\lambda =2}={\begin{bmatrix}1\\0\\0\end{bmatrix}},\quad \mathbf {x} _{\lambda =3}={\begin{bmatrix}0\\1\\1\end{bmatrix}}.}$$ Let us take $${\displaystyle V={\begin{bmatrix}0&0\\1&0\\0&1\end{bmatrix}},}$$ then $${\displaystyle V^{*}AV={\begin{bmatrix}2&1\\1&2\end{bmatrix}}}$$ with eigenvalues ${\displaystyle 1,3}$ and the corresponding eigenvectors $${\displaystyle \mathbf {y} _{\mu =1}={\begin{bmatrix}1\\-1\end{bmatrix}},\quad \mathbf {y} _{\mu =3}={\begin{bmatrix}1\\1\end{bmatrix}},}$$ so that the Ritz values are ${\displaystyle 1,3}$ and the Ritz vectors are $${\displaystyle \mathbf {\tilde {x}} _{{\tilde {\lambda }}=1}={\begin{bmatrix}0\\1\\-1\end{bmatrix}},\quad \mathbf {\tilde {x}} _{{\tilde {\lambda }}=3}={\begin{bmatrix}0\\1\\1\end{bmatrix}}.}$$ We observe that each one of the Ritz vectors is exactly one of the eigenvectors of ${\displaystyle A}$ for the given ${\displaystyle V}$ as well as the Ritz values give exactly two of the three eigenvalues of ${\displaystyle A}$. A mathematical explanation for the exact approximation is based on the fact that the column space of the matrix ${\displaystyle V}$ happens to be exactly the same as the subspace spanned by the two eigenvectors ${\displaystyle \mathbf {x} _{\lambda =1}}$ and ${\displaystyle \mathbf {x} _{\lambda =3}}$ in this example.

Using the normal matrix

The definition of the singular value ${\displaystyle \sigma }$ and the corresponding left and right singular vectors is ${\displaystyle Mv=\sigma u}$ and ${\displaystyle M^{*}u=\sigma v}$. Having found one set (left of right) of approximate singular vectors and singular values by applying naively the Rayleigh–Ritz method to the Hermitian normal matrix ${\displaystyle M^{*}M\in \mathbb {C} ^{N\times N}}$ or ${\displaystyle MM^{*}\in \mathbb {C} ^{M\times M}}$, whichever one is smaller size, one could determine the other set of left of right singular vectors simply by dividing by the singular values, i.e., ${\displaystyle u=Mv/\sigma }$ and ${\displaystyle v=M^{*}u/\sigma }$. However, the division is unstable or fails for small or zero singular values.

An alternative approach, e.g., defining the normal matrix as ${\displaystyle A=M^{*}M\in \mathbb {C} ^{N\times N}}$ of size ${\displaystyle N\times N}$, takes advantage of the fact that for a given ${\displaystyle N\times m}$ matrix ${\displaystyle W\in \mathbb {C} ^{N\times m}}$ with orthonormal columns the eigenvalue problem of the Rayleigh–Ritz method for the ${\displaystyle m\times m}$ matrix $${\displaystyle W^{*}AW=W^{*}M^{*}MW=(MW)^{*}MW}$$ can be interpreted as a singular value problem for the ${\displaystyle N\times m}$ matrix ${\displaystyle MW}$. This interpretation allows simple simultaneous calculation of both left and right approximate singular vectors as follows.

1. Compute the ${\displaystyle N\times m}$ matrix ${\displaystyle MW}$.
2. Compute the thin, or economy-sized, SVD ${\displaystyle MW=\mathbf {U} \Sigma \mathbf {V} _{h},}$ with ${\displaystyle N\times m}$ matrix ${\displaystyle \mathbf {U} }$, ${\displaystyle m\times m}$ diagonal matrix ${\displaystyle \Sigma }$, and ${\displaystyle m\times m}$ matrix ${\displaystyle \mathbf {V} _{h}}$.
3. Compute the matrices of the Ritz left ${\displaystyle U=\mathbf {U} }$ and right ${\displaystyle V_{h}=\mathbf {V} _{h}W^{*}}$ singular vectors.
4. Output approximations ${\displaystyle U,\Sigma ,V_{h}}$, called the Ritz singular triplets, to selected singular values and the corresponding left and right singular vectors of the original matrix ${\displaystyle M}$ representing an approximate Truncated singular value decomposition (SVD) with left singular vectors restricted to the column-space of the matrix ${\displaystyle W}$.

The algorithm can be used as a post-processing step where the matrix ${\displaystyle W}$ is an output of an eigenvalue solver, e.g., such as LOBPCG, approximating numerically selected eigenvectors of the normal matrix ${\displaystyle A=M^{*}M}$.

In quantum physics

In quantum physics, where the spectrum of the Hamiltonian is the set of discrete energy levels allowed by a quantum mechanical system, the Rayleigh–Ritz method is used to approximate the energy states and wavefunctions of a complicated atomic or nuclear system.^[8] In fact, for any system more complicated than a single hydrogen atom, there is no known exact solution for the spectrum of the Hamiltonian.^[7]

In this case, a trial wave function, ${\displaystyle \Psi }$, is tested on the system. This trial function is selected to meet boundary conditions (and any other physical constraints). The exact function is not known; the trial function contains one or more adjustable parameters, which are varied to find a lowest energy configuration.

It can be shown that the ground state energy, ${\displaystyle E_{0}}$, satisfies an inequality: $${\displaystyle E_{0}\leq {\frac {\langle \Psi |{\hat {H}}|\Psi \rangle }{\langle \Psi |\Psi \rangle }}.}$$

That is, the ground-state energy is less than this value. The trial wave-function will always give an expectation value larger than or equal to the ground-energy.

If the trial wave function is known to be orthogonal to the ground state, then it will provide a boundary for the energy of some excited state.

The Ritz ansatz function is a linear combination of N known basis functions ${\displaystyle \left\lbrace \Psi _{i}\right\rbrace }$, parametrized by unknown coefficients: $${\displaystyle \Psi =\sum _{i=1}^{N}c_{i}\Psi _{i}.}$$

With a known Hamiltonian, we can write its expected value as $${\displaystyle \varepsilon ={\frac {\left\langle \displaystyle \sum _{i=1}^{N}c_{i}\Psi _{i}\right|{\hat {H}}\left|\displaystyle \sum _{i=1}^{N}c_{i}\Psi _{i}\right\rangle }{\left\langle \left.\displaystyle \sum _{i=1}^{N}c_{i}\Psi _{i}\right|\displaystyle \sum _{i=1}^{N}c_{i}\Psi _{i}\right\rangle }}={\frac {\displaystyle \sum _{i=1}^{N}\displaystyle \sum _{j=1}^{N}c_{i}^{*}c_{j}H_{ij}}{\displaystyle \sum _{i=1}^{N}\displaystyle \sum _{j=1}^{N}c_{i}^{*}c_{j}S_{ij}}}\equiv {\frac {A}{B}}.}$$

The basis functions are usually not orthogonal, so that the overlap matrix S has nonzero nondiagonal elements. Either ${\displaystyle \left\lbrace c_{i}\right\rbrace }$ or ${\displaystyle \left\lbrace c_{i}^{*}\right\rbrace }$ (the conjugation of the first) can be used to minimize the expectation value. For instance, by making the partial derivatives of ${\displaystyle \varepsilon }$ over ${\displaystyle \left\lbrace c_{i}^{*}\right\rbrace }$ zero, the following equality is obtained for every k = 1, 2, ..., N: $${\displaystyle {\frac {\partial \varepsilon }{\partial c_{k}^{*}}}={\frac {\displaystyle \sum _{j=1}^{N}c_{j}(H_{kj}-\varepsilon S_{kj})}{B}}=0,}$$ which leads to a set of N secular equations: $${\displaystyle \sum _{j=1}^{N}c_{j}\left(H_{kj}-\varepsilon S_{kj}\right)=0\quad {\text{for}}\quad k=1,2,\dots ,N.}$$

In the above equations, energy ${\displaystyle \varepsilon }$ and the coefficients ${\displaystyle \left\lbrace c_{j}\right\rbrace }$ are unknown. With respect to c, this is a homogeneous set of linear equations, which has a solution when the determinant of the coefficients to these unknowns is zero: $${\displaystyle \det \left(H-\varepsilon S\right)=0,}$$ which in turn is true only for N values of ${\displaystyle \varepsilon }$. Furthermore, since the Hamiltonian is a hermitian operator, the H matrix is also hermitian and the values of ${\displaystyle \varepsilon _{i}}$ will be real. The lowest value among ${\displaystyle \varepsilon _{i}}$ (i=1,2,..,N), ${\displaystyle \varepsilon _{0}}$, will be the best approximation to the ground state for the basis functions used. The remaining N-1 energies are estimates of excited state energies. An approximation for the wave function of state i can be obtained by finding the coefficients ${\displaystyle \left\lbrace c_{j}\right\rbrace }$ from the corresponding secular equation.

In mechanical engineering

The Rayleigh–Ritz method is often used in mechanical engineering for finding the approximate real resonant frequencies of multi degree of freedom systems, such as spring mass systems or flywheels on a shaft with varying cross section. It is an extension of Rayleigh's method. It can also be used for finding buckling loads and post-buckling behaviour for columns.

Consider the case whereby we want to find the resonant frequency of oscillation of a system. First, write the oscillation in the form, $${\displaystyle y(x,t)=Y(x)\cos \omega t}$$ with an unknown mode shape ${\displaystyle Y(x)}$. Next, find the total energy of the system, consisting of a kinetic energy term and a potential energy term. The kinetic energy term involves the square of the time derivative of ${\displaystyle y(x,t)}$ and thus gains a factor of ${\displaystyle \omega ^{2}}$. Thus, we can calculate the total energy of the system and express it in the following form: $${\displaystyle E=T+V\equiv A[Y(x)]\omega ^{2}\sin ^{2}\omega t+B[Y(x)]\cos ^{2}\omega t}$$

By conservation of energy, the average kinetic energy must be equal to the average potential energy. Thus, $${\displaystyle \omega ^{2}={\frac {B[Y(x)]}{A[Y(x)]}}=R[Y(x)]}$$ which is also known as the Rayleigh quotient. Thus, if we knew the mode shape ${\displaystyle Y(x)}$, we would be able to calculate ${\displaystyle A[Y(x)]}$ and ${\displaystyle B[Y(x)]}$, and in turn get the eigenfrequency. However, we do not yet know the mode shape. In order to find this, we can approximate ${\displaystyle Y(x)}$ as a combination of a few approximating functions ${\displaystyle Y_{i}(x)}$ $${\displaystyle Y(x)=\sum _{i=1}^{N}c_{i}Y_{i}(x)}$$ where ${\displaystyle c_{1},c_{2},\cdots ,c_{N}}$ are constants to be determined. In general, if we choose a random set of ${\displaystyle c_{1},c_{2},\cdots ,c_{N}}$, it will describe a superposition of the actual eigenmodes of the system. However, if we seek ${\displaystyle c_{1},c_{2},\cdots ,c_{N}}$ such that the eigenfrequency ${\displaystyle \omega ^{2}}$ is minimised, then the mode described by this set of ${\displaystyle c_{1},c_{2},\cdots ,c_{N}}$ will be close to the lowest possible actual eigenmode of the system. Thus, this finds the lowest eigenfrequency. If we find eigenmodes orthogonal to this approximated lowest eigenmode, we can approximately find the next few eigenfrequencies as well.

In general, we can express ${\displaystyle A[Y(x)]}$ and ${\displaystyle B[Y(x)]}$ as a collection of terms quadratic in the coefficients ${\displaystyle c_{i}}$: $${\displaystyle B[Y(x)]=\sum _{i}\sum _{j}c_{i}c_{j}K_{ij}=\mathbf {c} ^{\mathsf {T}}K\mathbf {c} }$$ $${\displaystyle A[Y(x)]=\sum _{i}\sum _{j}c_{i}c_{j}M_{ij}=\mathbf {c} ^{\mathsf {T}}M\mathbf {c} }$$ where ${\displaystyle K}$ and ${\displaystyle M}$ are the stiffness matrix and mass matrix of a discrete system respectively.

The minimization of ${\displaystyle \omega ^{2}}$ becomes: $${\displaystyle {\frac {\partial \omega ^{2}}{\partial c_{i}}}={\frac {\partial }{\partial c_{i}}}{\frac {\mathbf {c} ^{\mathsf {T}}K\mathbf {c} }{\mathbf {c} ^{\mathsf {T}}M\mathbf {c} }}=0}$$

Solving this, $${\displaystyle \mathbf {c} ^{\mathsf {T}}M\mathbf {c} {\frac {\partial \mathbf {c} ^{\mathsf {T}}K\mathbf {c} }{\partial \mathbf {c} }}-\mathbf {c} ^{\mathsf {T}}K\mathbf {c} {\frac {\partial \mathbf {c} ^{\mathsf {T}}M\mathbf {c} }{\partial \mathbf {c} }}=0}$$ $${\displaystyle K\mathbf {c} -{\frac {\mathbf {c} ^{\mathsf {T}}K\mathbf {c} }{\mathbf {c} ^{\mathsf {T}}M\mathbf {c} }}M\mathbf {c} =\mathbf {0} }$$ $${\displaystyle K\mathbf {c} -\omega ^{2}M\mathbf {c} =\mathbf {0} }$$

For a non-trivial solution of c, we require determinant of the matrix coefficient of c to be zero. $${\displaystyle \det(K-\omega ^{2}M)=0}$$

This gives a solution for the first N eigenfrequencies and eigenmodes of the system, with N being the number of approximating functions.

Simple case of double spring-mass system

The following discussion uses the simplest case, where the system has two lumped springs and two lumped masses, and only two mode shapes are assumed. Hence M = [m₁, m₂] and K = [k₁, k₂].

A mode shape is assumed for the system, with two terms, one of which is weighted by a factor B, e.g. Y = [1, 1] + B[1, −1]. Simple harmonic motion theory says that the velocity at the time when deflection is zero, is the angular frequency ${\displaystyle \omega }$ times the deflection (y) at time of maximum deflection. In this example the kinetic energy (KE) for each mass is ${\textstyle {\frac {1}{2}}\omega ^{2}Y_{1}^{2}m_{1}}$ etc., and the potential energy (PE) for each spring is ${\textstyle {\frac {1}{2}}k_{1}Y_{1}^{2}}$ etc.

We also know that without damping, the maximal KE equals the maximal PE. Thus, $${\displaystyle \sum _{i=1}^{2}\left({\frac {1}{2}}\omega ^{2}Y_{i}^{2}M_{i}\right)=\sum _{i=1}^{2}\left({\frac {1}{2}}K_{i}Y_{i}^{2}\right)}$$

The overall amplitude of the mode shape cancels out from each side, always. That is, the actual size of the assumed deflection does not matter, just the mode shape.

Mathematical manipulations then obtain an expression for ${\displaystyle \omega }$, in terms of B, which can be differentiated with respect to B, to find the minimum, i.e. when ${\displaystyle d\omega /dB=0}$. This gives the value of B for which ${\displaystyle \omega }$ is lowest. This is an upper bound solution for ${\displaystyle \omega }$ if ${\displaystyle \omega }$ is hoped to be the predicted fundamental frequency of the system because the mode shape is assumed, but we have found the lowest value of that upper bound, given our assumptions, because B is used to find the optimal 'mix' of the two assumed mode shape functions.

There are many tricks with this method, the most important is to try and choose realistic assumed mode shapes. For example, in the case of beam deflection problems it is wise to use a deformed shape that is analytically similar to the expected solution. A quartic may fit most of the easy problems of simply linked beams even if the order of the deformed solution may be lower. The springs and masses do not have to be discrete, they can be continuous (or a mixture), and this method can be easily used in a spreadsheet to find the natural frequencies of quite complex distributed systems, if you can describe the distributed KE and PE terms easily, or else break the continuous elements up into discrete parts.

This method could be used iteratively, adding additional mode shapes to the previous best solution, or you can build up a long expression with many Bs and many mode shapes, and then differentiate them partially.

In dynamical systems

The Koopman operator allows a finite-dimensional nonlinear system to be encoded as an infinite-dimensional linear system. In general, both of these problems are difficult to solve, but for the latter we can use the Ritz-Galerkin method to approximate a solution.^[14]