A Restart Code for the Evaluation
of Matrix Functions
About
This Matlab code computes restarted Krylov approximations to f(A)b, where A is a square, large and sparse matrix, f is a matrix function (such that f(A) is defined) and b is a vector.A simple error bound (indicator in the non-Hermitian case) is provided and it is also possible to perform thick restarts or harmonic restarts. Simply call
[f,out] = funm_kryl (A,b,param)For the two implementations of this restarted Krylov method we refer to [1] and [5].
This Matlab code is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
Download
Download the latest version of funm_kryl here:Download funm_kryl (about 20 kB, version 20090409)To use funm_kryl extract the zip-archive on your computer and change to the directory funm_kryl in Matlab. Run the example files to see how it works.
Input parameter
The structure param stores all relevant parameters for the algorithm to run. A default setup can be generated byparam = param_init ()The fields of param can now be modified as required. Note that not every parameter configuration makes sense, e.g., you cannot use param.H_full==false if param.function is not a structure of type parfrac. Therefore you should rerun
param = param_init (param)which will check whether your settings are valid. If necessary, param_init returns a modified valid param structure.
| Field of param | Type | Values | Examples |
| function | function handle | a handle to a Matlab function which takes a square matrix argument H and returns f(H) | param.function = @sinm; param.function = @myfunm; |
|---|---|---|---|
| struct | a structure representing a partial fraction, see parfrac() for details | param.function = parfrac('expBA',16); | |
| V_full | boolean | false: (default) do not store the full Krylov basis | --- |
| true: keep the full Krylov basis in memory and return it in out.V_full | --- | ||
| H_full | boolean | false: use recursive scheme to compute updates of Arnoldi approximations (only available if param.function is a structure as returned by parfrac) | --- |
| true: (default) evaluate matrix function for full Hessenberg matrix H_full which grows in size with each restart cycle, out.H_full will be returned | --- | ||
| restart_length | integer >= 1 | if the dimension of the Krylov space exceeds param.restart_length the Krylov approximation will be updated and a new Krylov space is generated | param.restart_length = 20; |
| max_restarts | integer >= 1 | number of cycles of the restarting algorithm (if the number is 1 no restart is done) | param.max_restarts = 10; |
| hermitian | boolean | false: (default) A is not Hermitian (the Arnoldi process is used) | --- |
| true: A is Hermitian (the Lanczos process is used) | --- | ||
| reorth_number | 0 (default), 1 or 2 | number reorthonormalizations in the Gram-Schmidt orthonormalization (has no effect in the Hermitian case) | --- |
| thick | [] | (default) No thick-restarts are performed. | --- |
| function handle | Thick-restarts are performed. This requires a function that reorderes the Schur form U*T*U' of a Hessenberg matrix H with eigenvalues D such that wanted eigenvalues occur in the upper left block of the reordered Schur matrix. The wanted eigenvalues are returned in out.thick_replaced and the unwanted eigenvalues are returned in out.thick_interpol. | param.thick = @thick; | |
| harmonic_target | inf (default) or number | Target for harmonic Ritz extraction after each cycle. The value inf corresponds to the standard Ritz extraction. | --- |
| exact | vector | exact solution of f(A)b; if provided, out.err is a vector of the absolute error in the 2-norm of the Krylov approximation for each restart cycle | param.exact = expm(A)*b; param.exact = []; |
| bound | boolean | compute lower and upper error bounds (in non-Hermitan case these may be useful as error indicators) | --- |
| stopping_accuracy | positive number | the restart iteration stops if the upper error bound or the absolute error of the Krylov approximation falls below this treshold (default accuracy is 1e-12) | --- |
| min_decay | positive number | the restart iteration stops if the linear convergence rate of the upper error bound or the absolute error of the Krylov approximation is greater than this number (default linear decay rate is 0.95) | --- |
| waitbar | boolean | show waitbar | --- |
Outputs
[f,out] = funm_kryl (A,b,param) report(param.out) % show simple reportThe restarted Krylov approxiation is returned in f.
out is a structure with the following fields
| Field | Type | Values |
| out.update | vector | norm of the updates of the Krylov approximations for each restart cycle |
| out.err | vector | absolute error of the Krylov approximations for each restart cycle (only available if param.exact is non-empty) |
| out.time | vector | computation time (i.e., Arnoldi/Lanczos decomposition + update of the approximation) for each restart cycle (sum(out.time) is the overall computation time) |
| out.H_full | matrix | block-bidiagonal upper Hessenberg matrix, only returned if param.H_full == 1 |
| out.V_full | matrix | the basis vectors of the full Krylov space, blockwise orthonormal, only returned if param.V_full == 1 |
| out.bound_1, out.bound_2 | vector | lower and upper error estimators, only if param.bound == 1 (requires eigenvalue calculation of a matrix of size (param.restart_length times param.restart_length) in each restart cycle). For A Hermitian these bounds are based on Gauss-Radau and Gauss-Lobatto quadrature, respectively, and would be correct for A's smallest and largest eigenvalue to be known. Since we use the smallest and largest Ritz values instead, these bounds improve for later cycles. These bounds are based on an extension of an idea given in [8] and [7]. |
| out.stop_condition | string | string that tells which stop condition was satisfied at termination of the algorithm |
| out.thick_replaced, out.thick_interpol | cell array | if thick-restarts are used these arrays contain the wanted resp. unwanted eigenvalues which occured during the restart method. |
Please note: You may also call [f,out,param_out] = funm_kryl (A,b,param). Then the output structure param_out is the same as the input structure param if all the parameter settings were valid. Otherwise param_out is the corrected version of the input-param-structure which was then actually used by the algorithm.
Rational approximations
r = parfrac (type,degree,param1) param.parfrac = r x = 3; eval_parfrac(r,x) % test evaluation at point xThe structure r is a representation of a partial fraction having degree poles. The following types are implemented, further details may be found in [9]:
| Type | Available degrees | Meaning | param1 |
| expCF | 2,4,...,16(?) | quasi-best rational approximation to the exponential function on (-inf,0] via Caratheodory-Fejer | n.a. |
| expWP | 1,2,... | rational approximation to the exponential function on (-inf,0] with poles on the Weidemann parabolic contour | n.a. |
| expTW | 2,4,... | rational approximation to the exponential function on (-inf,0] with poles on the Talbot-Weidemann cotangent contour | n.a. |
| expWH | 1,2,... | rational approximation to the exponential function on (-inf,0] with poles on the Weideman hyperbolic contour | n.a. |
| expBA | 2,3,...,18 | best rational approximation to the exponential function on (-inf,0] (Cody-Meinardus-Varga approximation), table lookup with coefficients given in [3], [4] | n.a. |
| cfud | 1,2,...,8 | best rational approximation to function defined by param1 on the unit disk (CF-method) | handle to function analytic on unit disk, e.g., param1=@exp, param1=@cos |
| sqrtZ | --- | sqrtZ has been removed. Use sqrtinvZ to approximate x/sqrt(x). | --- |
| sqrtinvZ | 1,2,...,50 | Zolotarovs relative best rational approximation to 1/sqrt(x) on the interval [1,param1] | right-end of approximation interval, e.g., param1=1e3 |
| signZ | 1,2,...,50 | Zolotarovs best rational approximation to sign(x) on the set [-param1,param1]\[-1,1] | end-point of approximation interval, e.g., param1=1e3 |
| inv | 1 | the function 1/z (solution of linear systems of equations) | n.a. |
References
[1] M. Afanasjew, M. Eiermann, O. G. Ernst, and S. Güttel. Implementation of a Restarted Krylov Subspace Method for the Evaluation of Matrix Functions. To appear, 2007.[2] M. Afanasjew, M. Eiermann, O. G. Ernst, and S. Güttel. On the steepest descent method for matrix functions. In preparation, 2007.
[3] A. J. Carpenter, A. Ruttan, and R. S. Varga. Extended numerical computations on the '1/9' conjecture in rational approximation theory. In P. R. Graves-Morris, E. B. Saff, and R. S. Varga, editors, Rational Approximation and Interpolation, Proceedings, Tampa, Florida 1983, number 1105 in Lecture Notes in Mathematics, pages 383-411, Heidelberg, 1984. Springer-Verlag.
[4] W. J. Cody, G. Meinardus, and R. S. Varga. Chebyshev rational approximation to exp(-x) in [0,+inf) and applications to heat-conduction problems. J. Approx. Theory, 2:50-65, 1969.
[5] M. Eiermann and O. G. Ernst. A restarted Krylov subspace method for the evaluation of matrix functions. SIAM J. Numer. Anal., 44:2481-2504, 2006.
[6] M. Hochbruck and C. Lubich. On Krylov subspace approximations to the matrix exponential operator. SIAM J. Numer. Anal., 34(5):1911-1925, October 1997.
[7] B. Philippe and R. B. Sidje. Transient solutions of Markov processes by Krylov subspaces. Research Report, RR-1989, INRIA, 1995.
[8] Y. Saad. Analysis of some Krylov subspace approximations to the operator. SIAM J. Numer. Anal., 29(1):209-228, February 1992.
[9] L. N. Trefethen, J. A. C. Weideman, and T. Schmelzer. Talbot quadratures and rational approximations. BIT, 46:653-670, 2006.