# Variational procedure and generalized Lanczos recursion for small-amplitude classical oscillations

## Abstract

Variational procedure is developed that yields lowest frequencies of small-amplitude oscillations of classical Hamiltonian systems. Genuine Lanczos recursion is generalized to treat related non-Hermitian eigenvalue problems.

Normal modes $`\xi `$ and frequencies $`\omega `$ of small oscillations of a classical system near the equilibrium are determined by the secular equation 

$$\omega ^2M\xi =K\xi ,$$
(1)   
where $`M`$ and $`K`$ are $`N\times N`$ symmetric positive definite matrices of mass coefficients and spring constants respectively. In many applications the number $`N`$ of degrees of freedom is large, while only a few lowest frequencies are of interest . Equation (1) represents a problem more complex than a regular symmetric eigenvalue problem, unless $`M`$ or $`K`$ is diagonal.

Equation (1) can be transformed into the Hamiltonian form by introducing canonical momentum $`\eta =\omega M\xi `$:

$$K\xi =\omega \eta ,T\eta =\omega \xi ,$$
(2)   
where $`T=M^1`$. Thus, the frequencies of the normal modes are the eigenvalues of a $`2N\times 2N`$ matrix

$$\left(\begin{array}{cc}0& T\\ K& 0\end{array}\right)$$
(3)   
The spectrum of this matrix consists of pairs $`\pm \omega `$, since $`(\xi ,\eta )`$ is also a solution of (2) that corresponds to $`\omega `$. The lowest frequency $`\omega _{\mathrm{min}}`$ is the lowest positive eigenvalue of the matrix (3).

Although the eigenvalues of the matrix (3) are always real, the matrix itself is non-Hermitian, unless $`K=T`$. Therefore, its diagonalization poses a formidable task. The major problem is that no general minimum principle exists that yields eigenvalues of arbitrary diagonalizable non-Hermitian matrices. This does not allow to formulate a variational procedure similar to the Rayleigh-Ritz procedure for Hermitian matrices. If $`K=T`$, the matrix (3) is Hermitian, and its positive eigenvalues coincide with those of $`K`$ and $`T`$.

As known from quantum mechanics, the lowest eigenvalue $`ϵ_{\mathrm{min}}`$ of a Hermitian matrix $`H`$ can be obtained from the minimum principle

$$ϵ_{\mathrm{min}}=\underset{\{\psi \}}{\mathrm{min}}\frac{(\psi H\psi )}{(\psi \psi )}.$$
(4)   
The minimum is to be searched over all vectors $`\psi `$. The Ritz variational procedure is an approximation when the set $`\{\psi \}`$ in (4) is restricted to some subspace $`𝒦`$ of dimension $`n<N`$. The best approximation to $`ϵ_{\mathrm{min}}`$ in the sense of (4) is obtained as the lowest eigenvalue of the $`n\times n`$ Rayleigh matrix $`\stackrel{~}{H}`$, obtained by projection of $`H`$ onto $`𝒦`$.

The special paired structure of the matrix (3) makes it possible to generalize (4) such as to yield $`\omega _{\mathrm{min}}`$. In fact,

$$\omega _{\mathrm{min}}=\underset{\{\xi ,\eta \}}{\mathrm{min}}\frac{(\xi K\xi )+(\eta T\eta )}{2\left|(\xi \eta )\right|}.$$
(5)   
The minimum is to be searched over all possible phase space configurations $`\{\xi ,\eta \}`$. Before providing the proof to this equation let me point out some of its features.

First, it states that $`\omega _{\mathrm{min}}`$ is the minimum harmonic part of the total energy, $`(\xi K\xi )/2+(\eta T\eta )/2`$, over the phase space configurations normalized by $`(\xi \eta )=1`$. Since $`K`$ and $`T`$ are both positive definite, the right-hand side is strictly positive and so is $`\omega _{\mathrm{min}}`$. Second, equation (5) is symmetric in $`K`$ and $`T`$, according to the nature of the problem. When $`K=T`$ the minimum is achieved at $`\xi =\eta `$, and (5) becomes the same as (4).

Note that the functional in (5) has no maximum, since the denominator can be made arbitrarily small. The global minimum, however, always exists. This is not obvious, since a set of all pairs of vectors with $`(\xi \eta )=1`$ is not compact. Indeed, say, any vector orthogonal to $`\eta `$ can be added to $`\xi `$, making $`|\xi |`$ arbitrarily large. However, the functional in (5) grows indefinitely in this case, such that the global minimum is achieved at finite $`|\xi |`$ and $`|\eta |`$.

Variation of (5) with respect to $`\xi `$ and $`\eta `$ yields equations (2). Thus, the solutions of (2) are the stationary points of (5). The global minimum (5), therefore, gives indeed $`\omega _{\mathrm{min}}`$. The singularity in the denominator poses no problem, since it corresponds to infinitely large values of functional, while near the minimum it is analytic.

Minimum principle (5) can, in fact, be obtained from the Thouless minimum principle , derived for non-Hermitian matrices that appear in random phase approximation (RPA). Equation (5) transforms into the Thouless minimum principle by substitution: $`A=(K+T)/2`$, $`B=(KT)/2`$, $`x=(\xi +\eta )/2`$, and $`y=(\xi \eta )/2`$.

Variational procedure similar to the Rayleigh-Ritz procedure can be formulated if coordinates $`\xi `$ and momenta $`\eta `$ in (5) are restricted to some subspaces $`𝒰`$ and $`𝒱`$ of dimension $`n`$, respectively.

Let $`\{\xi _i\}`$ and $`\{\eta _i\}`$ be two sets of vectors that span $`𝒰`$ and $`𝒱`$, such that $`(\xi _i\eta _j)=\delta _{ij}`$. Expanding $`\xi =u_i\xi _i`$, $`\eta =v_i\eta _i`$ and varying (5) with respect to $`u_i`$ and $`v_i`$, we find the latter to obey a $`2n\times 2n`$ eigenvalue equation

$$\left(\begin{array}{cc}0& \stackrel{~}{T}\\ \stackrel{~}{K}& 0\end{array}\right)\left(\begin{array}{c}u\\ v\end{array}\right)=\stackrel{~}{\omega }\left(\begin{array}{c}u\\ v\end{array}\right),$$
(6)   
with $`\stackrel{~}{K}_{ij}=(\xi _iK\xi _j)`$ and $`\stackrel{~}{T}_{ij}=(\eta _iT\eta _j)`$. Equation (6) generalizes Hermitian Rayleigh-Ritz eigenvalue equation for $`\stackrel{~}{H}`$. It has $`2n`$ solutions $`\pm \stackrel{~}{\omega }`$, the lowest positive of which gives the best approximation to $`\omega _{\mathrm{min}}`$ in the sense of equation (5).

Krylov subspace  for the matrix (3) can be constructed by acting with it many times on an arbitrary vector $`(\xi _1,\eta _1)`$:

$$\left(\begin{array}{c}\xi _1\\ \eta _1\end{array}\right),\left(\begin{array}{c}T\eta _1\\ K\xi _1\end{array}\right),\left(\begin{array}{c}TK\xi _1\\ KT\eta _1\end{array}\right),\mathrm{}$$
(7)   
The subspace that spans first $`n`$ vectors of this sequence has the property of approximating an invariant subspace of (3). Thus, it is natural to expand approximation to an eigenvector of (3) as a linear combination of these vectors. In other words, the natural choice for the subspaces $`𝒰`$ and $`𝒱`$ for the variational procedure described above are the subspaces $`𝒰_n`$ and $`𝒱_n`$ that span the upper and lower components of first $`n`$ vectors of (7).

In order to implement the variational procedure, it is necessary to construct a biorthogonal basis $`\{\xi _i,\eta _i\}`$, $`i=1,\mathrm{},n`$ in $`𝒰_n`$ and $`𝒱_n`$ and compute matrix elements of $`\stackrel{~}{K}`$ and $`\stackrel{~}{T}`$. Both tasks can be performed simultaneously using the following recursion:

$`\xi _{i+1}`$ $`=`$ $`\beta _{i+1}^1(T\eta _i\alpha _i\xi _i\beta _i\xi _{i1})`$ (9)   
$`\eta _{i+1}`$ $`=`$ $`\delta _{i+1}^1(K\xi _i\gamma _i\eta _i\delta _i\eta _{i1}).`$ (10)   
The four coefficients $`\alpha _i`$, $`\beta _i`$, $`\gamma _i`$, and $`\delta _i`$ are to be chosen at each step $`i`$ such as to make $`\xi _{i+1}`$ orthogonal to $`\eta _i`$ and $`\eta _{i1}`$, and $`\eta _{i+1}`$ — orthogonal to $`\xi _i`$ and $`\xi _{i1}`$. This appears to be enough to ensure global biorthogonality $`(\xi _i\eta _j)=\delta _{ij}`$.

Indeed, assume biorthogonality to hold up to step $`i`$. Multiplying (9) by $`\eta _j`$, $`j<i1`$, we have $`(\eta _j\xi _{i+1})(\eta _jT\eta _i)=(\eta _iT\eta _j)=0`$ due to Hermiticity of $`T`$ and the fact that $`T\eta _j`$ is a linear combination of all $`\xi _k`$ with $`kj+1<i`$. Thus, the biorthogonality also holds for the step $`i+1`$.

Multiplying (9) by $`\eta _{i1}`$, $`\eta _i`$, and $`\eta _{i+1}`$ and using biorthogonality, we get $`(\xi _i\eta _i)=1`$, $`\stackrel{~}{K}_{ii}=\alpha _i`$, and $`\stackrel{~}{K}_{i,i1}=\stackrel{~}{K}_{i1,i}=\beta _i`$. Similarly, $`\stackrel{~}{T}_{ii}=\gamma _i`$ and $`\stackrel{~}{T}_{i,i1}=\stackrel{~}{T}_{i1,i}=\delta _i`$. All other matrix elements of $`\stackrel{~}{K}`$ and $`\stackrel{~}{T}`$ vanish.

The recursion (Variational procedure and generalized Lanczos recursion for small-amplitude classical oscillations) is a straightforward generalization of the Hermitian Lanczos recursion 

$$\psi _{i+1}=\beta _{i+1}^1(H\psi _i\alpha _i\psi _i\beta _i\psi _{i1})$$
(11)   
applicable to any Hermitian matrix $`H`$. When $`K=T`$ and $`\xi _1=\eta _1`$ both equations (Variational procedure and generalized Lanczos recursion for small-amplitude classical oscillations) coincide with each other and with equation (11), up to the notation.

As in the case of the Hermitian Lanczos algorithm, several lowest frequencies can be found one by one by projecting the $`\xi `$ and $`\eta `$-components of converged eigenvectors out of $`𝒱_n`$ and $`𝒰_n`$ subspaces respectively.

The method was tested on a set of large sparse random matrices of the form (3). Symmetric matrices $`T`$ and $`K`$ were generated to have an average of 40 randomly distributed and randomly positioned matrix elements in each row. Both $`K`$ and $`T`$ were shifted by an appropriate constant to ensure positive definiteness. Figure 1 demonstrates the convergence results for a matrix of the size $`2N=200000`$.

For smaller matrices up to $`2N=2000`$, where it was possible to obtain all eigenvalues with regular methods, the present method has converged to the true lowest frequency in all instances.

In conclusion, the method is proposed that generalizes Rayleigh-Ritz variational procedure and Lanczos recursion to the case of non-Hermitian matrices of the form (3), that determine normal modes and frequencies of small-amplitude oscillations of Hamiltonian systems.

Equations (2) have numerous applications beyond purely mechanical problems. Schroedinger equation in non-orthogonal basis represents a generalized symmetric eigenvalue problem similar to (1). RPA and other time-dependent techniques in nuclear physics and quantum chemistry lead to the equations similar to (2) . At last, eigenvectors of so-called Hamiltonian matrices, to which (3) is a special case, solve the nonlinear algebraic Riccati equation which appears in the theory of stability and optimal control .

I would like to acknowledge numerous enlightening discussions with Vladimir Chernyak during my appointment at the University of Rochester.