home » rbk » Xprojects » linalg

Chebyshev.cpp File Reference

Source code for Chebyshev polynomial approximations. More...

#include "linalg.h"
#include "quantum.h"
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_errno.h>

Include dependency graph for Chebyshev.cpp:


Functions
void	GenerateChebyshevMatrix (itype N, DenseMatrix< complex > &M)
	Computes an NxN matrix, where $M_{m,n}$ is the coefficient of in the 'th Chebyshev polynomial.
void	ResizeChebyshevMatrix (itype N, DenseMatrix< complex > &M)
	Efficiently resizes an existing Chebyshev polynomial matrix to NxN, using the coefficients already computed.
void	GenerateChebyshevEvolutionVector (double time, itype N, DenseVector< complex > &V)
	Computes an N-element vector, where is the coefficient of the 'th Chebyshev polynomial in the expansion of $e^{-ix}$ .
int	ChebyshevPrep (double timestep, double prec, DenseVector< complex > &PowerCoeffs, DenseVector< double > &CoeffsBoundSeq)
	Precomputes the coefficients needed to evolve a state using ChebyshevStep(), for time `timestep`, with accuracy `prec`.

Detailed Description

Source code for Chebyshev polynomial approximations.

A standard and nontrivial task in quantum simulations is to compute the unitary evolution of a state vector:

$\left|\psi\right\rangle \longrightarrow \left|\psi(t)\right\rangle = e^{-i\mathbf{H}t}\left|\psi\right\rangle$

Matrix exponentiation is nontrivial, especially if the matrix in question is very large. However, computing $\left|\psi(t)\right\rangle$ does not actually require computing $e^{-i\mathbf{H}t}$ , or even $\mathbf{H}$ . It can be written as a power series in $\mathbf{H}$ :

$e^{-i\mathbf{H}t}\left|\psi\right\rangle \approx \left(\mathbf{1} -i\mathbf{H}t - \frac{t^2}{2}\mathbf{H}^2\ldots\right) \left|\psi\right\rangle$

This requires only that we compute $\mathbf{H}^n\left|\psi\right\rangle$ , which in turn can be done recursively if we can compute $\left|\psi\right\rangle \longrightarrow \mathbf{H}\left|\psi\right\rangle$ .

The Taylor expansion given above turns out to converge rather slowly; the Taylor series is optimized for local convergence at $ t=0 $ , at the cost of all larger $ t $ . Other power series expansions exist, and have better global convergence properties. The power series with the best global convergence properties over a fixed intervals is the Chebyshev series:

The th Chebyshev polynomial is an th order polynomial in , defined on $x \in [-1\ldots 1]$ .
Chebyshev polynomials are orthonormal with respect to the integration kernel $\kappa(x)\mathrm{d}x = \frac{\mathrm{d}x}{\sqrt{1-x^2}}$ , so letting $f(x) \cdot g(x) \equiv \int_{-1}^1{f(x)g(x)\kappa(x)\mathrm{d}x}$ , we have $T_n(x) \cdot T_m(x) = \delta_{nm}$ .
Any function can be expanded as
$f(x) \approx \tilde{f}_N(x) = \sum_{n=0}^N{\left[f(x) \cdot T_n(x)\right] T_n(x)}$
$\tilde{f}_N(x)$ is almost as close to , by $L_\infty$ distance (on the interval $x \in [-1\ldots 1]$ !) as the optimal 'th order polynomial approximation.

This makes Chebyshev polynomials ideally suited for quantum state evolution. When they are generalized to matrices, the argument $ x $ becomes a matrix, which typically has many eigenvalues. We want the approximation of $ f(x) $ to be accurate for all the eigenvalues. This motivates the uniform convergence of the Chebyshev polynomial expansion. It is necessary to bound the highest and lowest eigenvalues of $\mathbf{H}$ , and then rescale energy and time so that $||\mathbf{\tilde{H}}||_\infty < 1$ , but once this is done, computing the correct power series is straightforward. From the rescaled time $\tilde{t}$ , it is possible to identify how many coefficients will be needed to approximate the matrix exponential well; beyond a critical value, the error declines exponentially.

Representing $e^{-i\mathbf{H}t}$ as a power series in $\mathbf{H}$ requires first computing the Chebyshev polynomials themselves. The method GenerateChebyshevMatrix() does this via a simple recursion relation, storing the coefficient of $ x^m $ in $ T_n(x) $ as an element $M_{mn}$ in a (triangular) matrix. If more coefficients are needed later, ResizeChebyshevMatrix() can compute just the additional ones needed.

The second step involves computing the expansion of $e^{-i\tilde{t}x}$ in Chebyshev polynomials. The method GenerateChebyshevEvolutionVector() does this using Bessel functions (borrowed from the Gnu Scientific Library), and stores the expansion in a complex vector. The coefficients of $\mathbf{\tilde{H}}^n$ in the approximated matrix exponential can then be computed by matrix-vector multiplication.

For a particular Hamiltonian and time, ChebyshevPrep() first determines what order approximation will be required to bound the error by some very small number, then computes the coefficients. A single step of duration $ t=1 $ can then be accomplished by ChebyshevStep(). The entire process of evolving for time $ t $ (including both ChebyshevPrep() and ChebyshevStep()), can be done in one step using ChebyshevEvolve().

However, taking a single step of large duration can lead to substantially increased error -- not because of error in the Chebyshev polynomials, but because the coefficients become so large for high-order approximations that they amplify rounding errors. This could very likely be taken care of by a more sophisticated method of evaluating polynomials, but this would require keeping additional temporary vectors (e.g., $\mathbf{\tilde{H}}^n\left|\psi\right\rangle$ for several values of $ n $ ), which could be memory-intensive for extremely large vectors. As an alternative, ChebyshevStepEvolve() will automatically break up large timesteps into several smaller jumps of bounded size.

While the Chebyshev polynomial method is (unlike split-operator or Cayley methods) not explicitly norm-preserving, its accuracy is so absurdly good (errors are typically $O(10^{-13})$ in double precision) that unitarity is not a concern.

Generated on Wed Jun 14 22:25:26 2006 for linalg by

1.4.4

Chebyshev.cpp File Reference

Functions

Detailed Description