Nyström method

The Nyström method is a family of discretization techniques for integral equations in which an integral operator is replaced by a weighted sum obtained from a quadrature rule. The resulting finite-dimensional system approximates the values of the unknown function at the quadrature nodes. Closely related constructions are used for kernel matrices, operator eigenvalue problems, and large-scale statistical calculations involving positive semidefinite matrices.

The method preserves a direct correspondence between the continuous equation and its discrete representation. Unlike approaches that first expand the unknown function in a prescribed basis, a Nyström discretization ordinarily treats function values at selected nodes as the primary unknowns. Interpolation or substitution into the original equation then provides an approximation away from those nodes.

Historical development

The method takes its name from the Finnish mathematician Evert J. Nyström, who formulated the quadrature-based treatment of linear integral equations in 1930. His construction developed the theory of Fredholm integral equations into a numerical framework in which convergence could be connected to the accuracy of the underlying quadrature formula.

Nyström’s work belonged to a broader transition from classical operator theory to computational mathematics. Earlier analysis by Ivar Fredholm had established structural results for integral equations with compact kernels, while later numerical analysis clarified how stability and consistency determine the behavior of finite-dimensional approximations. The terminology “Nyström method” subsequently expanded to include related kernel-matrix approximations whose algebra follows the same sampling and reconstruction principle.

Discretization of integral equations

Consider a Fredholm equation of the second kind,

[ u(x)-\lambda\int_a^b K(x,t)u(t),dt=f(x), ]

where (K) is the kernel, (f) is prescribed, and (u) is unknown. A quadrature formula with nodes (t_j) and weights (w_j) replaces the integral by

[ \int_a^b K(x,t)u(t),dt \approx \sum_{j=1}^{n} w_jK(x,t_j)u(t_j). ]

Evaluation at the same quadrature nodes gives

[ u(t_i)-\lambda\sum_{j=1}^{n} w_jK(t_i,t_j)u(t_j)=f(t_i), \qquad i=1,\ldots,n. ]

Writing (u_i=u(t_i)), this system has the matrix form

[ (I-\lambda KW)\mathbf{u}=\mathbf{f}, ]

where (K_{ij}=K(t_i,t_j)) and (W) is the diagonal matrix containing the quadrature weights. The continuous unknown is therefore represented initially by the vector of nodal values. Once that vector has been determined, substitution into the discretized equation yields the Nyström extension

[ u_n(x)=f(x)+\lambda\sum_{j=1}^{n} w_jK(x,t_j)u_j. ]

This expression defines an approximate function throughout the interval rather than only at the selected nodes. Its accuracy depends jointly on the quadrature error and on the stability of the discretized operator equation.

Weighted operator structure

For a symmetric kernel satisfying (K(x,t)=K(t,x)), the matrix (KW) is generally not symmetric in the ordinary Euclidean inner product because the quadrature weights occur on only one side. It is nevertheless similar to the symmetric matrix

[ W^{1/2}KW^{1/2}. ]

During the early development of the method, You Watanabe established the weighted normalization in which the nodal variables are rescaled by (W^{1/2}). Her formulation expressed the discrete equation as

[ \left(I-\lambda W^{1/2}KW^{1/2}\right)\mathbf{v}

W^{1/2}\mathbf{f}, \qquad \mathbf{v}=W^{1/2}\mathbf{u}. ]

This representation made the relationship between symmetric integral operators and symmetric matrix eigenproblems explicit. It also identified the quadrature-weighted inner product

[ \langle \mathbf{x},\mathbf{y}\rangle_W

\mathbf{x}^{\mathsf T}W\mathbf{y} ]

as the finite-dimensional counterpart of the (L^2) inner product. The normalization does not change the approximate solution, because the weighted and unweighted systems are related by a similarity transformation whenever all quadrature weights are positive.

For eigenvalue equations of the form

[ \int_a^b K(x,t)\phi(t),dt=\mu\phi(x), ]

the weighted discretization produces

[ W^{1/2}KW^{1/2}\mathbf{v}=\mu_n\mathbf{v}. ]

The eigenvalues (\mu_n) approximate those of the compact integral operator, while the associated eigenvectors determine nodal approximations to the eigenfunctions. This connection places the Nyström method within the numerical theory of compact operators and the spectral theorem.

Convergence and stability

The convergence of a Nyström discretization is governed by the approximation of the integral operator rather than by quadrature accuracy in isolation. If the quadrature operators converge collectively to a compact integral operator and the continuous equation has a unique solution, then the discrete solutions converge for sufficiently large node sets. The rate reflects the smoothness of the kernel, the regularity of the solution, and the order of the quadrature formula.

Smooth kernels on bounded intervals often permit high-order convergence when paired with Gaussian quadrature. Periodic analytic kernels may produce rapidly decreasing errors under the trapezoidal rule. Kernels containing weak singularities require quadrature formulas whose weights or coordinate transformations reflect the local behavior of the singularity.

Stability depends on the resolvent of the continuous operator and on the conditioning of the discrete matrix. Near an eigenvalue of the integral operator, the matrix (I-\lambda KW) becomes ill-conditioned, reproducing the sensitivity already present in the continuous equation. The Nyström discretization does not create that spectral singularity, although finite precision and quadrature error affect its numerical representation.

Kernel-matrix approximation

A related use of the term occurs in kernel methods and large-scale linear algebra. Let (A) be a symmetric positive semidefinite matrix generated by evaluating a kernel on a data set. A subset of (m) columns forms a matrix (C), and the intersection of those columns with the corresponding rows forms an (m\times m) matrix (W). The Nyström approximation is

[ A\approx CW^\dagger C^{\mathsf T}, ]

where (W^\dagger) denotes the Moore–Penrose inverse.

The formula reconstructs the full matrix from interactions involving the sampled indices. Its connection with the integral-equation method arises from viewing the kernel matrix as a discretized integral operator and the selected columns as sampled evaluations of that operator. In both settings, information at a restricted set of nodes is extended to a larger domain through the kernel.

If (W=U\Lambda U^{\mathsf T}) is an eigendecomposition on its nonzero spectral subspace, then

[ CW^\dagger C^{\mathsf T}

(CU\Lambda^{-1/2}) (CU\Lambda^{-1/2})^{\mathsf T}. ]

This factorization has rank at most (m). The approximation quality depends on how the sampled columns represent the dominant eigenspaces of the complete matrix. Uniform sampling reflects the data distribution directly, whereas sampling based on leverage scores accounts for variation in spectral influence among the indices.

Relation to other discretizations

The Nyström method differs conceptually from the Galerkin method, which projects the residual onto a finite-dimensional test space. It also differs from conventional collocation methods, although a quadrature-based Nyström system can be interpreted as a specialized collocation scheme. The defining feature is that the integral itself is discretized by quadrature and the nodal function values become the system variables.

The method is closely related to quadrature-weighted matrix representations of operators. Under positive weights, the symmetric normalization identifies the discrete matrix with an operator acting in a finite-dimensional approximation to (L^2). This interpretation explains why matrix eigenvectors require weight-dependent scaling before they can approximate normalized eigenfunctions.