22 Functions of a vector

In quantum chemistry, one often views multivariable functions as functions of one or more vectors. The archetypal example is the Rayleigh quotient in a non-orthonormal basis. Let \(\mathbf H\) and \(\mathbf S\) be Hermitian \(n\times n\) matrices, with \(\mathbf S\) positive definite. Think of \(\mathbf H\) as the matrix of a one-electron Hamiltonian in a non-orthonormal AO basis \(\{\chi_\mu\}_{\mu=1}^n\), and \(\mathbf S\) as the overlap matrix. Their elements are

\[ H_{\mu\nu}=\langle \chi_\mu|\hat H|\chi_\nu\rangle, \qquad S_{\mu\nu}=\langle \chi_\mu|\chi_\nu\rangle. \]

For

\[ \psi=\sum_{\mu=1}^n c_\mu\chi_\mu, \qquad c_\mu\in\mathbb C, \]

we collect the coefficients in the vector \(\mathbf c\in\mathbb C^n\). The expectation value of \(\hat H\) is

\[ \langle \hat H\rangle=\frac{\langle\psi|\hat H|\psi\rangle}{\langle\psi|\psi\rangle}=\frac{\mathbf c^\dagger \mathbf H\mathbf c}{\mathbf c^\dagger \mathbf S\mathbf c}. \]

Recall that \(\mathbf c^\dagger\) is the Hermitian transpose (conjugate transpose) of \(\mathbf c\). The last expression is the Rayleigh quotient

\[ R(\mathbf c)=\frac{\mathbf c^\dagger \mathbf H\mathbf c}{\mathbf c^\dagger \mathbf S\mathbf c}. \]

Since \(\mathbf H\) is Hermitian and \(\mathbf S\) is positive definite, \(R(\mathbf c)\) is real for every nonzero \(\mathbf c\). Notice also that multiplication of \(\mathbf c\) by a nonzero scalar does not change the quotient:

\[ R(a\mathbf c)=R(\mathbf c), \qquad a\ne0. \]

Thus the Rayleigh quotient depends on the direction represented by \(\mathbf c\), but not on its overall normalization or phase.

The gradient of a vector-dependent function

A function of a real vector \(\mathbf x=(x_1,\ldots,x_n)^T\) is simply a multivariable function written using vector notation,

\[ f:\mathbb R^n\to\mathbb R, \qquad f(\mathbf x)=f(x_1,\ldots,x_n). \]

Its gradient is the vector of partial derivatives

\[ \nabla_{\mathbf x}f=\begin{pmatrix}\partial f/\partial x_1\\ \vdots\\ \partial f/\partial x_n\end{pmatrix}. \]

For example, if

\[ f(\mathbf x)=\mathbf a^T\mathbf x=\sum_i a_i x_i, \]

then

\[ \nabla_{\mathbf x}f=\mathbf a. \]

A particularly important case is a quadratic form. For a real symmetric matrix \(\mathbf A\),

\[ f(\mathbf x)=\mathbf x^T \mathbf A\mathbf x \]

has gradient

\[ \nabla_{\mathbf x}f=2\mathbf A\mathbf x. \]

One way of seeing this is to write the quadratic form with indices,

\[ f(\mathbf x)=\sum_{ij}x_iA_{ij}x_j, \]

and differentiate with respect to one component \(x_k\). The variable occurs twice, giving two contributions. Symmetry of \(\mathbf A\) makes the two contributions equal.

In quantum chemistry the vectors are often complex. A real-valued function such as \(R(\mathbf c)\) then depends on both the coefficients \(c_i\) and their complex conjugates \(c_i^*\). A convenient notation is to temporarily treat \(\mathbf c\) and \(\mathbf c^*\) as independent variables. For example,

\[ f(\mathbf c,\mathbf c^*)=\mathbf c^\dagger \mathbf A\mathbf c=\sum_{ij}c_i^*A_{ij}c_j. \]

Differentiating with respect to \(c_k^*\) gives

\[ \frac{\partial f}{\partial c_k^*}=\sum_j A_{kj}c_j, \]

or, in vector notation,

\[ \frac{\partial}{\partial\mathbf c^*}(\mathbf c^\dagger \mathbf A\mathbf c)=\mathbf A\mathbf c. \]

Similarly,

\[ \frac{\partial}{\partial\mathbf c}(\mathbf c^\dagger \mathbf A\mathbf c)=\mathbf A^T\mathbf c^*. \]

This way of differentiating with respect to a complex variable and its conjugate is closely related to Wirtinger differentiation. For the elementary applications here, the main practical rule is simple: when differentiating with respect to \(\mathbf c^*\), regard \(\mathbf c\) as constant, and vice versa.

We can now differentiate the Rayleigh quotient. Using the quotient rule,

\[ \frac{\partial R}{\partial\mathbf c^*}=\frac{(\mathbf H\mathbf c)(\mathbf c^\dagger \mathbf S\mathbf c)-(\mathbf c^\dagger \mathbf H\mathbf c)(\mathbf S\mathbf c)}{(\mathbf c^\dagger \mathbf S\mathbf c)^2}. \]

Since \(R(\mathbf c)=(\mathbf c^\dagger \mathbf H\mathbf c)/(\mathbf c^\dagger \mathbf S\mathbf c)\), this can be written more transparently as

\[ \boxed{\frac{\partial R}{\partial\mathbf c^*}=\frac{(\mathbf H-R\mathbf S)\mathbf c}{\mathbf c^\dagger \mathbf S\mathbf c}.} \]

At a stationary point the derivative vanishes. Therefore

\[ (\mathbf H-R\mathbf S)\mathbf c=0, \]

\[ \boxed{\mathbf H\mathbf c=R\,\mathbf S\mathbf c.} \]

Thus stationary points of the Rayleigh quotient are generalized eigenvectors of \(\mathbf H\) with respect to \(\mathbf S\), and the stationary value of the quotient is the corresponding generalized eigenvalue. This is one of the fundamental variational ideas underlying quantum chemistry.

The Hessian of a vector-dependent function

The gradient describes the first-order change of a scalar function. The Hessian describes its second-order change. For a real-valued function \(f(\mathbf x)\), the Hessian is the matrix

\[ \mathbf H_f(\mathbf x)=\nabla_{\mathbf x}^2f, \qquad (\mathbf H_f)_{ij}=\frac{\partial^2f}{\partial x_i\partial x_j}. \]

For example, for the quadratic function

\[ f(\mathbf x)=\frac12\mathbf x^T \mathbf A\mathbf x+\mathbf b^T\mathbf x+c, \]

where \(\mathbf A\) is real and symmetric,

\[ \nabla f=\mathbf A\mathbf x+\mathbf b, \qquad \mathbf H_f=\mathbf A. \]

The Hessian enters the multivariable Taylor expansion. For a small displacement \(\Delta\mathbf x\),

\[ f(\mathbf x+\Delta\mathbf x)=f(\mathbf x)+\nabla f(\mathbf x)^T\Delta\mathbf x+\frac12\Delta\mathbf x^T \mathbf H_f(\mathbf x)\Delta\mathbf x+O(\|\Delta\mathbf x\|^3). \]

At a stationary point the gradient is zero, so the leading change is usually quadratic:

\[ f(\mathbf x+\Delta\mathbf x)-f(\mathbf x)\approx\frac12\Delta\mathbf x^T \mathbf H_f\Delta\mathbf x. \]

The signs of the Hessian eigenvalues therefore tell us about the local shape of the function. Positive eigenvalues correspond to directions in which the function curves upward, negative eigenvalues to directions in which it curves downward, and zero eigenvalues indicate flat directions that require more care.

For functions of complex vectors, the complete second-derivative structure can be written using derivatives with respect to both \(\mathbf c\) and \(\mathbf c^*\). We will not develop the full formalism here. The important idea is the same as in the real case: near a stationary vector, the second-order change tells us how sensitive the function is to small changes in the vector. For the Rayleigh quotient, this connects energy curvature to gaps between eigenvalues.

Exercises

Exercise 22.1 (From components to a gradient) Let \(\mathbf x=(x_1,x_2)^T\) and

\[ f(\mathbf x)=x_1^2+3x_2^2. \]

Write \(f\) without vector notation, as a function \(f(x_1,x_2)\).
Compute \(\partial f/\partial x_1\).
Compute \(\partial f/\partial x_2\).
Collect the two partial derivatives into the gradient \(\nabla_{\mathbf x}f\).
Evaluate the gradient at \(\mathbf x=(1,2)^T\).

Exercise 22.2 (Differentiating a quadratic form) Consider

\[ \mathbf A=\begin{pmatrix}2&0\\0&5\end{pmatrix}, \qquad \mathbf x=\begin{pmatrix}x_1\\x_2\end{pmatrix}, \qquad f(\mathbf x)=\mathbf x^T \mathbf A\mathbf x. \]

Multiply out \(\mathbf A\mathbf x\).
Multiply out \(\mathbf x^T\mathbf A\mathbf x\) and write \(f\) in terms of \(x_1\) and \(x_2\).
Differentiate the expression from part b with respect to \(x_1\) and \(x_2\).
Verify explicitly that

\[ \nabla_{\mathbf x}f=2\mathbf A\mathbf x. \]
Use a short Python script to evaluate \(f\) and its gradient at \(\mathbf x=(1,-1)^T\).

Exercise 22.3 (A first Rayleigh quotient) Let

\[ \mathbf H=\begin{pmatrix}1&0\\0&3\end{pmatrix}, \qquad \mathbf S=\mathbf I, \qquad \mathbf c=\begin{pmatrix}c_1\\c_2\end{pmatrix}, \]

where \(c_1\) and \(c_2\) are real.

Show that

\[ R(\mathbf c)=\frac{c_1^2+3c_2^2}{c_1^2+c_2^2}. \]
Evaluate \(R\) for \(\mathbf c=(1,0)^T\), \((0,1)^T\), and \((1,1)^T\).
Verify directly that multiplying \(\mathbf c\) by \(2\) does not change \(R\).
The stationary vectors satisfy \(\mathbf H\mathbf c=R\mathbf S\mathbf c\). Insert \(\mathbf c=(1,0)^T\) and determine \(R\).
Repeat part d for \(\mathbf c=(0,1)^T\). What familiar matrix problem have you recovered?

Exercise 22.4 (Treating \(c\) and \(c^*\) as independent variables) When differentiating expressions involving a complex number \(c\) and its complex conjugate \(c^*\), it is often useful to formally treat \(c\) and \(c^*\) as independent variables.

This means that when differentiating with respect to \(c\), we regard \(c^*\) as a constant, and when differentiating with respect to \(c^*\), we regard \(c\) as a constant.

Using this rule, compute \[ \frac{\partial c}{\partial c}, \qquad \frac{\partial c^*}{\partial c}. \]
Similarly, compute \[ \frac{\partial c}{\partial c^*}, \qquad \frac{\partial c^*}{\partial c^*}. \]
Consider the function \[ f(c,c^*)=c^*c. \] Compute \[ \frac{\partial f}{\partial c} \] by treating \(c^*\) as a constant.
Compute \[ \frac{\partial f}{\partial c^*} \] by treating \(c\) as a constant.
Now consider \[ g(c,c^*)=c^*Ac, \] where \(A\) is a constant number. Compute \(\partial g/\partial c\) and \(\partial g/\partial c^*\).
Let \(c=x+iy\). Check that \(c^*c=x^2+y^2\). Thus, although we treated \(c\) and \(c^*\) as independent when differentiating, are they actually independent complex numbers?

Exercise 22.5 (From the Rayleigh quotient to an eigenvalue problem) Let

\[ R(\mathbf c)=\frac{\mathbf c^\dagger \mathbf H\mathbf c}{\mathbf c^\dagger \mathbf S\mathbf c}, \]

where \(\mathbf H\) and \(\mathbf S\) are Hermitian and \(\mathbf S\) is positive definite.

Define

\[ N(\mathbf c)=\mathbf c^\dagger \mathbf H\mathbf c, \qquad D(\mathbf c)=\mathbf c^\dagger \mathbf S\mathbf c. \]

Show that

\[ \frac{\partial N}{\partial\mathbf c^*}=\mathbf H\mathbf c, \qquad \frac{\partial D}{\partial\mathbf c^*}=\mathbf S\mathbf c. \]
Use the quotient rule to show that

\[ \frac{\partial R}{\partial\mathbf c^*}=\frac{(\mathbf H-R\mathbf S)\mathbf c}{\mathbf c^\dagger \mathbf S\mathbf c}. \]
Show that a stationary point satisfies

\[ \mathbf H\mathbf c=R\mathbf S\mathbf c. \]
Suppose \(\mathbf c\) satisfies \(\mathbf H\mathbf c=\varepsilon\mathbf S\mathbf c\). Show directly that \(R(\mathbf c)=\varepsilon\).
Explain in words why minimizing the Rayleigh quotient is closely connected to finding the lowest-energy eigenvector.

Exercise 22.6 (Bonding and antibonding orbitals of \(\mathrm{H}_2^+\)) Consider a minimal basis for \(\mathrm{H}_2^+\) consisting of two normalized hydrogen \(1s\) orbitals, \(\chi_A\) and \(\chi_B\). In this basis, suppose

\[ \mathbf H=\begin{pmatrix}\alpha&\beta\\\beta&\alpha\end{pmatrix}, \qquad \mathbf S=\begin{pmatrix}1&s\\s&1\end{pmatrix}, \]

where all quantities are real and \(0<s<1\).

Consider the two coefficient vectors

\[ \mathbf c_+=\begin{pmatrix}1\\1\end{pmatrix}, \qquad \mathbf c_-=\begin{pmatrix}1\\-1\end{pmatrix}. \]

Write the corresponding molecular orbitals in terms of \(\chi_A\) and \(\chi_B\). Which one has the same sign on the two atoms, and which one changes sign?
Show that \(\mathbf c_+\) and \(\mathbf c_-\) satisfy the generalized eigenvalue equation \(\mathbf H\mathbf c=\varepsilon\mathbf S\mathbf c\).
Show that their energies are

\[ \varepsilon_+=\frac{\alpha+\beta}{1+s}, \qquad \varepsilon_-=\frac{\alpha-\beta}{1-s}. \]
Take \(\alpha=-0.50\), \(\beta=-0.20\), and \(s=0.30\) in atomic units. Compute the two energies. Which orbital is lower in energy?
Normalize \(\mathbf c_+\) and \(\mathbf c_-\) with respect to the overlap matrix, so that \(\mathbf c^T \mathbf S\mathbf c=1\).
Use Python and scipy.linalg.eigh(H, S) to calculate the generalized eigenvalues and eigenvectors numerically, and compare with your results.

Exercise 22.7 (Exploring an energy surface) Let

\[ \mathbf H=\begin{pmatrix}1&0\\0&4\end{pmatrix}, \qquad \mathbf S=\mathbf I. \]

Because the Rayleigh quotient is independent of the length of \(\mathbf c\), restrict the coefficient vector to the unit circle,

\[ \mathbf c(\theta)=\begin{pmatrix}\cos\theta\\\sin\theta\end{pmatrix}. \]

Insert \(\mathbf c(\theta)\) into the Rayleigh quotient and show that

\[ R(\theta)=\cos^2\theta+4\sin^2\theta. \]
Differentiate \(R(\theta)\) and find its stationary points for \(0\leq\theta<2\pi\).
Identify which stationary points are minima and which are maxima. Relate the corresponding vectors to the eigenvectors of \(\mathbf H\).
Write a short Python script that plots \(R(\theta)\) for \(0\leq\theta\leq2\pi\). Mark or print the minimum and maximum values.
What does this example suggest about the relationship between the smallest/largest eigenvalues and the minimum/maximum of the Rayleigh quotient?

Exercise 22.8 (Second-order variation of the Rayleigh quotient) Let \(\mathbf H\) be a real symmetric matrix with orthonormal eigenvectors \(\mathbf u_1,\ldots,\mathbf u_n\) and eigenvalues

\[ \lambda_1\leq\lambda_2\leq\cdots\leq\lambda_n. \]

For simplicity take \(\mathbf S=\mathbf I\). Consider a normalized eigenvector \(\mathbf u_k\) and a small displacement in a perpendicular direction,

\[ \mathbf c=\mathbf u_k+\delta\mathbf v, \qquad \mathbf u_k^T\mathbf v=0, \]

where \(|\delta|\ll1\) and \(\|\mathbf v\|=1\).

Show that

\[ \mathbf c^T\mathbf c=1+\delta^2. \]
Use \(\mathbf H\mathbf u_k=\lambda_k\mathbf u_k\) and \(\mathbf u_k^T\mathbf v=0\) to show that

\[ \mathbf c^T \mathbf H\mathbf c=\lambda_k+\delta^2\mathbf v^T \mathbf H\mathbf v. \]
Hence show that

\[ R(\mathbf c)=\frac{\lambda_k+\delta^2\mathbf v^T\mathbf H\mathbf v}{1+\delta^2}. \]
Use \((1+\delta^2)^{-1}=1-\delta^2+O(\delta^4)\) to obtain the expansion

\[ R(\mathbf c)=\lambda_k+\delta^2(\mathbf v^T\mathbf H\mathbf v-\lambda_k)+O(\delta^4). \]

Why is there no term proportional to \(\delta\)?
Now choose \(\mathbf v=\mathbf u_j\) with \(j\ne k\). Show that

\[ R(\mathbf u_k+\delta\mathbf u_j)=\lambda_k+\delta^2(\lambda_j-\lambda_k)+O(\delta^4). \]
Use this result to explain why the lowest eigenvector is a minimum of the Rayleigh quotient, the highest eigenvector is a maximum, and an eigenvector in the middle generally has both uphill and downhill directions.
Choose a diagonal \(3\times3\) matrix with three different eigenvalues and verify the second-order formula numerically with a short Python script. Try several decreasing values of \(\delta\) and examine how rapidly the approximation improves.