Gauss quadrature approximation by Lanczos algorithm

一、 Gauss Quadrature

1.1 Gauss quadrature with weight function

For
$$I(f)={\int }_a^b\rho (x)f(x)dx\text{\hspace{1em}(1)}$$
If the following quadrature formula[2,p220]
$${\int }_a^b\rho (x)f(x)dx\approx \sum _{k=0}^n{A}_{k}f(x_k)\text{\hspace{1em}(2)}$$
with the $(2n+1)$ algebraic degree of precision, that is, when the degree of polynomial $f(x)$ less than or equal to $2n+1$, the left and right sides of equation (2) are strictly equal. Thus, it’s called the Gauss-type quadrature formula with weights $\rho (x)$.
The nodes $x_k$ are called Gauss points (nodes), $A_k$s are called quadrature coefficients about $\rho (x)$ (weights or coefficients).

1.2 The caculation of the weights and nodes

(Note: This section is mainly about how to caculate the Gauss points and coefficents. The theoretical foundation can refer to [2, p223]. )

• 1)Undeterminated coefficient method: As the Gaussian quadrature formula in equation (2) is accurately for any polynomial with degree at most $2n+1$. Thus, for the polynomial sequence { 1 , x , x 2 , ⋯   , x 2 n + 1 } \{1,x,x^{2},\cdots,x^{2n+1}\} {1,x,x2,⋯,x2n+1}, equation (2) are all exactly equal.
  Based on the above analysis, we can construct the following $(2n+2)$ equations about quadrature nodes { x k } k = 0 n \{x_{k}\}_{k=0}^{n} {xk​}k=0n​ and quadrature coefficient { A k } k = 0 n \{A_{k}\}_{k=0}^{n} {Ak​}k=0n​ :
  $$\{\begin{array}{l}{\sum }_{k=0}^n{A}_k={\int }_a^{b}1\mathrm{d}x\\ {\sum }_{k=0}^n{x}_k{A}_k={\int }_a^{b}x\mathrm{d}x\\ {\sum }_{k=0}^n{x}_k^2{A}_k={\int }_a^b{x}^2\mathrm{d}x\\ ⋮\\ {\sum }_{k=0}^n{x}_k^{2n+1}{A}_k={\int }_a^b{x}^{2n+1}\mathrm{d}x\end{array}$$

Note: when $n\ge 2$, it is difficult to solve the above system of equations directly.

• 2)Orthogonal polynomials method: The basic idea mainly based on the following two theorem.

[2, Theorem 5.5.1, p221] Let $I(f)={\int }_a^{b}f(x)\mathrm{d}x$ and its interpolation quadrature formula is $$I_n(f)=\sum _{k=0}^n{A}_{k}f(x_k),$$ and denote $$W_{n+1}(x)=(x-x_0)(x-x_1)\cdots (x-x_n).$$
Then { x k } k = 0 n \{x_{k}\}_{k=0}^{n} {xk​}k=0n​ are Gaussian nodes if and only if $W_{n+1}(x)$ is orthogonal to any polynomial $p(x)$ with degree at most $n$, that is $${\int }_a^{b}p(x)W_{n+1}(x)\mathrm{d}x=0.$$

[2, Theorem 5.5.3, p223] Let { q n ( x ) } n = 0 ∞ \{q_{n}(x)\}_{n=0}^{\infty} {qn​(x)}n=0∞​ is an orthogonal polynomilas sequence on $[a,b]$, thus $q_n(x)$ has $n$ different zeros on $[a,b]$.

` As the orthogonal polynomials sequence { q k ( x ) } k = 0 n + 1 \{q_{k}(x)\}_{k=0}^{n+1} {qk​(x)}k=0n+1​ form a basis for a polynomial space of degree n+1, and $q_n(x)$ is orthogoanl to any polynomial of degree at most n on the interval [a,b].
Thus, for a given weight function $\rho (x)$ (can be $\rho (x)=1$), if we can find a orthogonal polynomials sequence { q 0 ( x ) , q 1 ( x ) , ⋯   , q n ( x ) , q n + 1 ( x ) } \{q_{0}(x),q_{1}(x),\cdots,q_{n}(x),q_{n+1}(x)\} {q0​(x),q1​(x),⋯,qn​(x),qn+1​(x)} related to the weight funciton $\rho (x)$.
Where $q_{n+1}(x)$ is exact degree $n+1$, suppose its zeros are $t_0,\cdots ,t_n$, then they are also the Gauss points { x k } k = 0 n \{x_{k}\}_{k=0}^{n} {xk​}k=0n​, their coresponding coefficients { A k } k = 0 n \{A_{k}\}_{k=0}^{n} {Ak​}k=0n​ are derived by
$$A_k={\int }_a^b\rho (x)l_k(x)dx,\text{\hspace{1em}(3)}$$
where $l_k(t)$ are base function of interpolation
$$l_k(t)=\prod _{j=0,j\ne k}^n\frac{t-t_j}{t_k-t_j}.$$

Note: Using the zeros of orthogoanl polynomials to construct the Gauss quadrature formula, which is only effective for some special weight functions, such as $1$(Legendre), $\frac{1}{\sqrt{1-x^2}}$(Chebyshev 1), $\sqrt{1-x^2}$(Chebyshev 2) and $e^{-x^2}$(Hermite), et.al.
However, for the general weight function, the above undertermined coefficient method is usually adopted.

---

二、Lanczos algorithm

2.1 Krylov subspace

• Krylov subspace: $A$ is real matrix (unsymmetric) of order $n$. Let $v$ be a given vector and
  $$K_k=\left(v,Av,\cdots ,A^{k-1}v\right)$$
  be the Krylov matrix of dimension $n\times k$. The subspace spanned by the columns of matrix $K_k$ is called a Krylov subspace and denoted by $K_k(A,v)$ or $K(A,v)$.

Note: The natural basis of the Krylov subspace $K(A,v)$ given by the columns of the Krylov matirx $K_k$ is badly conditioned when $k$ is large.

2.2 Arnoldi algorithm

Arnoldi algorithm constructs an orthonormal basis of the Krylov subspace $K(A,v)$ by applying a variant of the Gram-Schmidt orthogonalization process.

• Arnoldi algorithm: Set vectors $v^{(j+1)}=A{v}^{(j)}$ with $v^{(1)}=v$, then $K(A,v)$ is spanned by the vectors { v ( j ) } j = 1 k \{v^{(j)}\}_{j=1}^{k} {v(j)}j=1k​.For constucting orthogonal basis vectors $v^j$, instead of orthogonalizing $A^{j}v$ against the previous vectors, orthogonalize $A{v}^j$.
  Starting from $v^1=v$(normalized,that is $v^1=v/\Vert v\Vert$),the $(j+1)$st vector of the basis is computed by using the previous vectors,
• Projection: $h_{i,j}=(A{v}^j,v^i),i=1,\cdots ,j,$
• Orthogonalize: ${\tilde{v}}^{j+1}=A{v}^j-{\sum }_{i=1}^j{h}_{i,j}{v}^i,$
• Length: $h_{j+1,j}=\Vert {\tilde{v}}^{j+1}\Vert$,(if $h_{j+1,j}=0$ then stop)
• Normalize: $v^{j+1}=\frac{{\tilde{v}}^{j+1}}{h_{j+1,j}}$

As $h_{j+1,j}{v}^{j+1}=A{v}^j-{\sum }_{i=1}^j{h}_{i,j}{v}^i$, we have
$$A{v}^j=\sum _{i=1}^j{h}_{i,j}{v}^i+h_{j+1,j}{v}^{j+1}.$$
If we collect the vectors $v^j,j=1,\cdots ,k$ in a matrix $V_k$, that is $V_k=[v_1,v_2,\cdots ,v_k]$ (thus $V_k$ is column orthogonal matrix), the relations defining the vectors $v^{k+1}$ can be written in the matrix form as
$$A{V}_k=V_k{H}_k+h_{k+1,k}{v}^{k+1}{\left(e^k\right)}^T,\text{\hspace{1em}(4)}$$
where $H_k$ is an upper Hessenberg matrix with elements $h_{i,j}$, $h_{i,j}=0,j=1,\cdots ,i-2,i>2$.
$$H_k=\left[\begin{array}{lllll}{h}_{1,1}& h_{1,2}& h_{1,3}& \cdots & h_{1,k}\\ h_{2,1}& h_{2,2}& h_{2,3}& \cdots & h_{2,k}\\ 0& h_{3,2}& h_{3,3}& \cdots & h_{3,k}\\ 0& 0& h_{4,3}& \cdots & h_{4,k}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& h_{k,k-1}& h_{k,k}\end{array}\right]$$

2.3 Lanczos algorithm

Multiplying equation (4) by $V_k^T$ and using orthogonality, we have
$$H_k=V_k^{T}A{V}_k.$$
If matrix $A$ is symmetric, then the Hessenberg matrix $H_k$ is tridiagonal and denoted by $J_k$, that is $J_{i,j}=0,j=i+2,\cdots ,k.$ (or $h_{i,j}=0,i=1,\cdots ,j-2,j+2,\cdots ,k$). This implies that the new vector $v^{k+1}$ can be computed by using only the two previous $v^k$ and $v^{k-1}$.
$$h_{j+1,j}{v}^{j+1}=A{v}^j-\sum _{i=1}^j{h}_{i,j}{v}^i=A{v}^j-h_{j-1,j}{v}^{j-1}-h_{j,j}{v}^j.$$
In matrix form, set ${\eta }_i$ denote the nonzero off-diagonal entries of $J_k$,
$$A{V}_k=V_k{J}_k+{\eta }_k{v}^{k+1}{\left(e^k\right)}^T.\text{\hspace{1em}(5)}$$

Equation (5) describes in matrix form the elegant Lanczos algorithm.
For simple notation system, starting from a nonzero vector $v^1=v/\Vert v\Vert ,{\alpha }_1=(A{v}^1,v^1),{\tilde{v}}^2=A{v}^1-{\alpha }_1{v}^1$, and then for $k=2,3,\cdots ,$
$$\begin{array}{c}{\eta }_{k-1}=\Vert {\tilde{v}}^k\Vert ,\\ v^k=\frac{{\tilde{v}}^k}{{\eta }_{k-1}},\\ {\alpha }_k=\left(v^k,A{v}^k\right)={\left(v^k\right)}^{T}A{v}^k,\\ {\tilde{v}}^{k+1}=A{v}^k-{\alpha }_k{v}^k-{\eta }_{k-1}{v}^{k-1}.\end{array}$$
$$J_k=\left[\begin{array}{lllll}{\alpha }_1& {\eta }_1& 0& \cdots & 0\\ {\eta }_1& {\alpha }_2& {\eta }_2& \cdots & 0\\ 0& {\eta }_2& {\alpha }_3& \cdots & 0\\ 0& 0& {\eta }_3& \cdots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& {\eta }_{k-1}& {\alpha }_k\end{array}\right]$$

三、Two important theory

3.1 Theorem 1

Theorem1 [1,Th.4.1]: Let ${\chi }_k(\lambda )$ be the determinant of $J_k-\lambda I$(which is a monic polynomial); then
$$v^k=p_k(A)v^1,p_k(\lambda )=(-1{)}^{k-1}\frac{{\chi }_{k-1}(\lambda )}{{\eta }_1\cdots {\eta }_{k-1}},k>1,p_1\equiv 1$$
The polynomials $p_k$ of degree $k-1$ are called the normalized Lanczos polynomials.

This theorem describes the most important property of Lanczos alg orithm:The Lanczos vectors $v^k$ are given as a polynomial in matrix $A$ applied to the initial vector $v^1$.

It is easy to verify that （expand on the last row of $J_{k+1}$）
$$\det (J_{k+1})={\alpha }_{k+1}\det (J_k)-{\eta }_k^2\det (J_{k-1}),\text{\hspace{1em}(6)}$$
and
$$\det (J_{k+1}-\lambda I)=({\alpha }_{k+1}-\lambda )\det (J_k-\lambda I)-{\eta }_k^2\det (J_{k-1}-\lambda I).\text{\hspace{1em}(7)}$$
From the expression of polynomial $p_k(\lambda )$ in Theorem 1, we know that the Lanczos polynomials satisfy a scalar three-term recurrence,
$${\eta }_k{p}_{k+1}(\lambda )=\left(\lambda -{\alpha }_k\right)p_k(\lambda )-{\eta }_{k-1}{p}_{k-1}(\lambda ),k=1,2,\dots ,\text{\hspace{1em}(8)}$$
with initial conditions $p_0\equiv 0,p_1\equiv 1$.

3.2 Theorem 2

Theorem2 [1, Th.4.2]:
Consider the Lanczos vectors $v^k$. There exists a measure $\alpha (\lambda )$ such that
$$\left(v^k,v^l\right)=⟨p_k,p_l⟩={\int }_a^b{p}_k(\lambda )p_l(\lambda )d\alpha (\lambda ),\text{\hspace{1em}(9)}$$
where $\alpha \le {\lambda }_1={\lambda }_{min}$ and $b\ge {\lambda }_n={\lambda }_{max}$, ${\lambda }_{min}$ and ${\lambda }_{max}$ being the smallest and largest eigenvalues of $A$, and $p_i$ are the Lanczos polynomials associated with $A$ and $v^1$.
$$\alpha (\lambda )=\{\begin{array}{ll}0,& \text{ if }\lambda <{\lambda }_1\\ {\sum }_{j=1}^i{\left[{\hat{v}}_j\right]}^2,& \text{ if }{\lambda }_i\le \lambda <{\lambda }_{i+1}\\ {\sum }_{j=1}^n{\left[{\hat{v}}_j\right]}^2,& \text{ if }{\lambda }_n\le \lambda \end{array}\text{\hspace{1em}(10)}$$
where $\hat{v}=Q^T{v}^1$,$v^1=v/\Vert v\Vert$, and the spectral decomposition of $A$ is $A=Q\Lambda Q^T$.

That is, normalized Lanczos polynomials sequence { p k ( λ ) } \{p_{k}(\lambda)\} {pk​(λ)} is orthogonal in the sense of measure $\alpha (\lambda )$.

Note: For the sake of simplicity, here suppose that the eigenvalues of $A$ are distinct.

四、Gauss quadrature by Lanczos algorithm

For symetric matrix $A$ with eigenvalue decomposition $A=Q\Lambda Q^T$, $f$ is smooth funciton, $v$ is random unit vector, then
$$v^{T}f(A)v=v^{T}Qf(\Lambda )Q^{T}v$$
set $\hat{v}=Q^{T}v$, then
$$v^{T}f(A)v={\hat{v}}^{T}f(\Lambda )v=\sum _{j=1}^{m}f({\lambda }_i)[{\hat{v}}_j{]}^2$$
Consider the above sum as Riemann-Stieltjes integral
$$\sum _{j=1}^{m}f({\lambda }_i)[{\hat{v}}_j{]}^2={\int }_a^{b}f(t)d\alpha (t).\text{\hspace{1em}(11)}$$
where the measure $\alpha (t)$ is defined as:
$$\alpha (t)=\{\begin{array}{ll}0,& \text{ if }t<{\lambda }_1=a\\ {\sum }_{j=1}^i{\left[{\hat{v}}_j\right]}^2,& \text{ if }{\lambda }_i\le t<{\lambda }_{i+1}\\ {\sum }_{j=1}^m{\left[{\hat{v}}_j\right]}^2,& \text{ if }b={\lambda }_m\le t.\end{array}$$

The integral in Equation (11) can be estimated using the Gauss quadrature,
$${\int }_a^{b}f(t)d\alpha (t)\approx \sum _{k=1}^m{A}_{k}f(x_k).\text{\hspace{1em}(12)}$$
As stated in the Theorem2, { p k } k = 1 m + 1 \{p_{k}\}_{k=1}^{m+1} {pk​}k=1m+1​ are orthonormal polynomial and the measure funciton is also $\alpha (t)$.
Thus, the Gaussian quadrature nodes in (12) can be derived by caculating the zeros of polynomial $p_{m+1}(x)$ with degree of $m$.
We can rewirte the three-term recurrence in Equation (8) in matrix form. Precisely, let $\mathbf{p}(\lambda )=[p_1(\lambda ),p_2(\lambda ),\cdots ,p_k(\lambda ){]}^T$ and $p_0(\lambda )=0$, then
$$\lambda \mathbf{p}(\lambda )=\mathbf{p}(\lambda )J_k+{\eta }_k{p}_{k+1}(\lambda )e_k^T.\text{\hspace{1em}(13)}$$
Then the zeros of $p_{k+1}(\lambda )$ are equal to the eigenvalues of $J_k$.
Thus, the Gaussian quadrature nodes in (12) can be derived by caculating the eigenvalues of tridiagonal Jacobi matrix $J_m$, where $J_m$ is derived by m-steps Lanczos algorithm associate with $A$ and $v$.

Note:$p_{k+1}(\lambda )$ is polynomial of degree k.

In conclusion, the Gauss point prolbem in Equation(12) can be get by the eigenvalues of $J_k$ which derived by Lanczos algorithm for matrix $A$ and $v$.

---

References:
[1] Golub, G. H. and Meurant, G. Matrices, Moments and Quadrature with Applications, Princeton University Press, 2010.
[2] 孙志忠，袁慰平，闻震初. 数值分析(第3版)，东南大学出版社.