Tridiagonalize \(A\) and approximate \(f(A)b \approx |b|_W V f(T) e_1\) via Lanczos algorithm. A is self-adjoint in the weights \( W\).
More...
|
| | UniversalLanczos () |
| | Allocate nothing, Call construct method before usage.
|
| |
| | UniversalLanczos (const ContainerType ©able, unsigned max_iterations) |
| | Allocate memory for the method.
|
| |
| template<class ... Params> |
| void | construct (Params &&...ps) |
| | Perfect forward parameters to one of the constructors.
|
| |
| void | set_max (unsigned new_max) |
| | Set the maximum number of iterations.
|
| |
| unsigned | get_max () const |
| | Get the current maximum number of iterations.
|
| |
| void | set_verbose (bool verbose) |
| | Set or unset debugging output during iterations.
|
| |
| double | get_bnorm () const |
| | Norm of b from last call to operator()
|
| |
| unsigned | get_iter () const |
| | Get the number of iterations in the last call to tridiag or solve (same as T.num_rows)
|
| |
| template<class MatrixType , class ContainerType0 , class ContainerType1 , class ContainerType2 , class FuncTe1 > |
| unsigned | solve (ContainerType0 &x, FuncTe1 f, MatrixType &&A, const ContainerType1 &b, const ContainerType2 &weights, double eps, double nrmb_correction=1., std::string error_norm="universal", double res_fac=1., unsigned q=1) |
| | \( x = f(A)b \approx ||b||_W V f(T) e_1 \) via Lanczos and matrix function computation. A is self-adjoint in the weights \( W\).
|
| |
| template<class MatrixType , class ContainerType0 , class ContainerType1 > |
| const dg::TriDiagonal< thrust::host_vector< double > > & | tridiag (MatrixType &&A, const ContainerType0 &b, const ContainerType1 &weights, double eps=1e-4, double nrmb_correction=1., std::string error_norm="compute_extreme_EV", double res_fac=1., unsigned q=1) |
| | Tridiagonalization of A using Lanczos method with \( f(x) = x^{-1} \).
|
| |
| template<class MatrixType , class ContainerType0 , class ContainerType1 , class ContainerType2 > |
| void | normMbVy (MatrixType &&A, const dg::TriDiagonal< thrust::host_vector< double > > &T, const ContainerType0 &y, ContainerType1 &x, const ContainerType2 &b, double bnorm) |
| | compute \( x = |b|_W V y \) from a given tridiagonal matrix T and in-place re-computation of V
|
| |
| template<class UnaryOp , class MatrixType , class ContainerType1 , class ContainerType2 > |
| const dg::TriDiagonal< thrust::host_vector< double > > & | tridiag (UnaryOp f, MatrixType &&A, const ContainerType1 &b, const ContainerType2 &weights, double eps, double nrmb_correction, std::string error_norm="residual", double res_fac=1., unsigned q=1) |
| | Tridiagonalization of A using Lanczos method.
|
| |
template<class ContainerType>
class dg::mat::UniversalLanczos< ContainerType >
Tridiagonalize \(A\) and approximate \(f(A)b \approx |b|_W V f(T) e_1\) via Lanczos algorithm. A is self-adjoint in the weights \( W\).
The M-Lanczos method is based on the paper Novel Numerical Methods for Solving the Time-Space Fractional Diffusion Equation in Two Dimensions by Q. Yang et al, but adopts a more efficient implementation similar to that in the PCG method. Further also the conventional Lanczos method can be found there and also in text books such as Iteratvie Methods for Sparse Linear Systems 2nd edition by Yousef Saad.
The Lanczos algorithm's main purpose is to find (approximate) Eigenvalues and Eigenvectors of a symmetric definite matrix. It does so by computing Eigenvalues of the (symmetric) tridiagonal matrix \( T\) (the "Ritz" or "Rayleigh-Ritz" values) that the Lanczos algorithm produces. The extreme values i.e. the maximum and minimum Eigenvalue of \( T\) converge fairly quickly to the extreme Eigenvalues of \(A\).
- Note
- We have several stopping criteria. The residual criterion stops when \( \tau ||r_i||_W = \tau ||\vec b||_W \beta_i (T^{-1})_{1m} \leq \epsilon_{rel} ||\vec b||_W + \epsilon_{aps} \) where \( \tau \) is a residual factor accounting for the condition of various matrix functions
-
The universal stopping criterion is based on the paper Estimating the error in matrix function approximations by Q. Eshghi, N. and Reichel L., The iteration stops when
\[
||\vec{e}_{f,m}||_W = ||\vec{b}||_W ||\left(\check{T} - f(\bar{T})\right)\vec{e}_1||_2 \leq \epsilon_{rel} ||\vec{b}||_W ||f(T)\vec e_1||_2 + \epsilon_{abs}
\]
with
\[
\bar{T} =
\begin{pmatrix}
T & \beta_m \vec{e}_m & & \\
\beta_m \vec{e}_m^T & \alpha_{m-1} & \beta_{m-2} & \\
& \beta_{m-2} & \alpha_{m-2} & \beta_{m-1-q} \\
& & \beta_{n-1-q} & \alpha_{n-q}
\end{pmatrix}
\]
\[
\check{T} =
\begin{pmatrix}
f(T) & 0 \\
0 & 0
\end{pmatrix}
\]
-
The "compute_extreme_EV" and "compute_max_EV" stopping criteria are based on tracking the convergence progress of \(\lambda^T_\min\) and \( \lambda^T_\max\) of the tridiagonal matrix \( T\).
\[
|\lambda^T_{\min, 10(i-1)} - \lambda^T_{\min , 10i}| < 10(\epsilon\lambda_\min + 10^{-12}\lambda_\max)\\
|\lambda^T_{\max, 10(i-1)} - \lambda^T_{\max , 10i}| < 10\epsilon\lambda_\max
\]
We test convergence only every 10th iteration for various reasons
- The Eigenvalue decomposition of \( T\) may become performance relevant for very large \( i\). Testing only every 10th iterations amortises the cost
- Lanczos is at risk of misconvergence which means the algorithm may "converge" to a \(\lambda^T\) that is not an extreme Eigenvalue because for several iterations no larger/smaller Eigenvalue appears (this happens not infrequently). Testing only every 10th iteration decreases this risk such that we do not see it in practise (but it does not vanish entirely)
Also here note that theoretically the minimum and maximum Eigenvalue should converge at the same rate, however the error constant for the minimum Eigenvalue is much worse (by the condition number \(\kappa\) of the matrix!). This is because all Eigenvalues that are found by \( T\) lie within the spectrum of \(A\), i.e.
\[
\lambda^A_\min < \lambda^T_\min < \lambda^T_\max < \lambda^T_\max
\\
\frac{\lambda^A_\max - \lambda^T_\max }{\lambda^A_\max} < 1
\\
\frac{\lambda^T_\min - \lambda^A_\min}{\lambda^A_\min} = \kappa\frac{\lambda^T_\min - \lambda^A_\min}{\lambda^A_\max} < \kappa
\]
-
THIS METHOD DOES NOT WORK (but is still instructive: On estimating the largest Eigenvalue with the Lanczos algorithm by B. N. Parlett, H. Simon and L. M. Stringer The proposed stopping criterion consists of checking if the found (approximate) Eigenvector is indeed an Eigenvector of \(A\). Suppose we find
\[
T_m s = \lambda s
\]
Then, we accept the corresponding vector \( y = V_m s\) as an Eigenvector of \( A\) if
\[
|| Ay - \lambda y ||_W \leq \epsilon || y||_W
\]
Since the \( V_m\) are normalised we esimate \( || y||_W > 0.9\) and we have \( A y= A V_m s = V_m T_m s + \beta_m \vec v_{m+1} \vec e_m^T s = \lambda V_m s + \beta_m \vec v_{m+1} s_m\) where \( s_m\) is the last component of \(s\). Thus, we have \(
A y - \lambda y = \beta_m \vec v_{m+1} s_m
\) and therefore
\[
\beta_m s_m \leq 0.9 \epsilon
\]
Unfortunately in practise \( \beta_m s_m \) is rather large and the convergence is extremely slow. Even if the Eigenvalues converge rather quickly, the Eigenvectors do not!
The common Lanczos method (and M-Lanczos) method are prone to loss of orthogonality for finite precision.
- Template Parameters
-
| ContainerType | Any class for which a specialization of TensorTraits exists and which fulfills the requirements of the there defined data and execution policies derived from AnyVectorTag and AnyPolicyTag. Among others
dg::HVec (serial), dg::DVec (cuda / omp), dg::MHVec (mpi + serial) or dg::MDVec (mpi + cuda / omp) std::vector<dg::DVec> (vector of shared device vectors), std::array<double, 4> (array of 4 doubles) or std::map < std::string, dg::DVec> ( a map of named vectors) double (scalar) and other primitive types ...
If there are several ContainerTypes in the argument list, then TensorTraits must exist for all of them |
- See also
- See The dg dispatch system for a detailed explanation of our type dispatch system
template<class ContainerType >
template<class MatrixType , class ContainerType0 , class ContainerType1 , class ContainerType2 , class FuncTe1 >
| unsigned dg::mat::UniversalLanczos< ContainerType >::solve |
( |
ContainerType0 & | x, |
|
|
FuncTe1 | f, |
|
|
MatrixType && | A, |
|
|
const ContainerType1 & | b, |
|
|
const ContainerType2 & | weights, |
|
|
double | eps, |
|
|
double | nrmb_correction = 1., |
|
|
std::string | error_norm = "universal", |
|
|
double | res_fac = 1., |
|
|
unsigned | q = 1 ) |
|
inline |
\( x = f(A)b \approx ||b||_W V f(T) e_1 \) via Lanczos and matrix function computation. A is self-adjoint in the weights \( W\).
Tridiagonalize \(A\) using Lanczos algorithm with a residual or universal stopping criterion
- Parameters
-
| x | output vector |
| f | the matrix function that is called like yH = f( T) where T is the tridiagonal matrix and returns the result of \( f(T)\vec e_1\) (for example dg::mat::make_FuncEigen_Te1( dg::SQRT<double>()) ) |
- See also
- Matrix-functions
- Parameters
-
| A | A self-adjoint, positive definit matrix |
| b | The initial vector that starts orthogonalization. |
- Attention
- If tridiagonalisation is used to determine Eigenvalues,
b should be initialised with random values to make misconvergence and breakdown due to beta = 0 less likely.
- Parameters
-
| weights | Weights that define the scalar product in which A is self-adjoint and in which the error norm is computed. |
| eps | relative accuracy of residual |
| nrmb_correction | the absolute error C in units of eps to be respected (ignored for "compute_*_EV" error norms) |
| error_norm | Either "residual" or "universal" or "compute_extreme_EV" or "compute_max_EV". In the latter two cases the iteration terminates if the extreme ( = maximum and minimum) /maximum Eigenvalue(s) of the tridiagonal matrix converge i.e. \( | \lambda_{\max,
i} - \lambda_{\max, i-1}|< \epsilon \lambda_{\max,i}\) (analogous for minimum). In the Eigenvalue cases nrmb_correction, res_fac and q are ignored. In the Eigenvalue case the algorithm trows if breakdown occurs. The error conditions is checked every iteration for the first 10 iterations and then only every 10th iteration. |
| res_fac | factor \( \tau\) in the "universal" and "resisdual" error_norm that is multiplied to the norm of the residual. Used to account for specific matrix function and operator in the convergence criterium |
| q | The q-number in the "universal" error_norm (ignored otherwise) |
- See also
- https://www.jstor.org/stable/2007471?seq=1
- Returns
- number of iterations of M-Lanczos routine
template<class ContainerType >
template<class UnaryOp , class MatrixType , class ContainerType1 , class ContainerType2 >
| const dg::TriDiagonal< thrust::host_vector< double > > & dg::mat::UniversalLanczos< ContainerType >::tridiag |
( |
UnaryOp | f, |
|
|
MatrixType && | A, |
|
|
const ContainerType1 & | b, |
|
|
const ContainerType2 & | weights, |
|
|
double | eps, |
|
|
double | nrmb_correction, |
|
|
std::string | error_norm = "residual", |
|
|
double | res_fac = 1., |
|
|
unsigned | q = 1 ) |
|
inline |
Tridiagonalization of A using Lanczos method.
Tridiagonalize \(A\) using Lanczos algorithm with a residual or universal stopping criterion on the function \( f(x) \)
- Parameters
-
| f | Unary function (ignored if error_norm is not "universal") |
- Parameters
-
| A | A self-adjoint, positive definit matrix |
| b | The initial vector that starts orthogonalization. |
- Attention
- If tridiagonalisation is used to determine Eigenvalues,
b should be initialised with random values to make misconvergence and breakdown due to beta = 0 less likely.
- Parameters
-
| weights | Weights that define the scalar product in which A is self-adjoint and in which the error norm is computed. |
| eps | relative accuracy of residual |
| nrmb_correction | the absolute error C in units of eps to be respected (ignored for "compute_*_EV" error norms) |
| error_norm | Either "residual" or "universal" or "compute_extreme_EV" or "compute_max_EV". In the latter two cases the iteration terminates if the extreme ( = maximum and minimum) /maximum Eigenvalue(s) of the tridiagonal matrix converge i.e. \( | \lambda_{\max,
i} - \lambda_{\max, i-1}|< \epsilon \lambda_{\max,i}\) (analogous for minimum). In the Eigenvalue cases nrmb_correction, res_fac and q are ignored. In the Eigenvalue case the algorithm trows if breakdown occurs. The error conditions is checked every iteration for the first 10 iterations and then only every 10th iteration. |
| res_fac | factor \( \tau\) in the "universal" and "resisdual" error_norm that is multiplied to the norm of the residual. Used to account for specific matrix function and operator in the convergence criterium |
| q | The q-number in the "universal" error_norm (ignored otherwise) |
- See also
- https://www.jstor.org/stable/2007471?seq=1
- Returns
- returns the tridiagonal matrix T. Note that \( T = (MV)^T A V \). The number of iterations is given by
T.num_rows
1.12.0