# Rayleigh quotient
Consider a linear least squares problem. Given a vector $\vec{b} \in \mathrm{R}^m$ and a matrix $A \in \mathrm{R} ^{m\times n}$, the best fit to $\vec{b}$ from the range of $A$ is $A\vec{x}$, where $\vec{x}$ is the linear least square solution and can be found through solving the normal equation

$$\newcommand{\RR}{\mathbb{R}}\newcommand{\ffrac}{\displaystyle \frac}\newcommand{\Tran}[1]{{#1}^{\mathrm{T}}}\Tran{A}A\vec{x} = \Tran{A}\vec{b}$$

Now we let $A$ be a real symmetric $n \times n$ matrix. Given an arbitary vector $\vec{x} \in \RR^n$, which may not be the eigenvector of $A$, we try to best fit that to $A\vec{x}$ among all multiple of $\vec{x}$.

$$\alpha^* = \arg \min_{\alpha\in \RR} \left\| \vec{x}\alpha - A\vec{x} \right\|_2$$

so that the solution $\alpha^*$ satisfies the equation:

$$\Tran{\vec{x}}\vec{x} \alpha^* = \Tran{\vec{x}}A\vec{x} \Longrightarrow \alpha^* = \frac{\Tran{\vec{x}}A\vec{x}} {\Tran{\vec{x}}\vec{x}}$$

And it's obvious that if $\vec{x}$ is indeed an eigenvector of matrix $A$, then

$$\alpha^* = \frac{\Tran{\vec{x}}A\vec{x}} {\Tran{\vec{x}}\vec{x}} = \frac{\Tran{\vec{x}}\lambda\vec{x}} {\Tran{\vec{x}}\vec{x}} = \lambda$$

And in particular, given a real symmetric matrix $A \in \RR^{n \times n}$, we call **Rayleigh quotient** of $A$ at any vector $\vec{x} \in \RR^n$

$$r(\vec{x})= \frac{\Tran{\vec{x}}A\vec{x}} {\Tran{\vec{x}}\vec{x}}$$
***

Then we do Taylor expansion of multiple variables.

$$r\left(\vec{x}\right) = r\left(\vec{q}\right) + \nabla \Tran{r\left(\vec{q}\right)} \left( \vec{x} - \vec{q} \right) + \frac{\Tran{\left( \vec{x} - \vec{q} \right)} \nabla^2 r\left(\tilde{\vec{q}}\right)\left( \vec{x} - \vec{q} \right)} {2}$$

Here $\vec{q}$ is an eigenvector of $A$ corresponding to eigenvalue $\lambda$. $\nabla r\left(\vec{q}\right)$ is the Jacobian of $r\left(\vec{x}\right)$ at $\vec{x} = \vec{q}$; $\nabla^2 r\left(\tilde{\vec{q}}\right)$ is the Hessian matrix of the function $r\left(\vec{x}\right)$ at $\tilde{\vec{q}}$ which is close to $\vec{q}$.

$\nabla r\left(\vec{q}\right)$ is a column vector, and contains the partial drivative of $r\left(\vec{x}\right)$ with respect to the components of $\vec{x}$. From the definition $r(\vec{x})= \ffrac{\Tran{\vec{x}}A\vec{x}} {\Tran{\vec{x}}\vec{x}}$, we have

$$
\frac{\partial r\left(\vec{x}\right)} {\partial x_k} = \frac{\frac{\partial \left(\Tran{\vec{x}}A\vec{x} \right)} {\partial {x_{k}}} \cdot \left(\Tran{\vec{x}}\vec{x}\right) - \left(\Tran{\vec{x}} A \vec{x}\right) \cdot \frac{ \partial \left(\Tran{\vec{x}}\vec{x}\right)} {\partial {x_{k}}}} {\left( \Tran{\vec{x}}\vec{x} \right)^2}
$$

One way to derive the formula for $\ffrac{\partial \left(\Tran{\vec{x}}A\vec{x} \right)} {\partial {x_{k}}}$ is to decomposite it completely, another way is shown below

$$\begin{align}
\frac{\partial \left(\Tran{\vec{x}}A\vec{x} \right)} {\partial {x_{k}}} &= \frac{\partial \left(\Tran{\vec{x}}A \right)} {\partial {x_{k}}} \cdot \vec{x} + \Tran{\vec{x}} \cdot \ffrac{\partial \left(A\vec{x} \right)} {\partial {x_{k}}} \\
&= 2 \cdot \frac{\partial \left(\Tran{\vec{x}}A \right)} {\partial {x_{k}}} \cdot \vec{x} \\
&= 2 \cdot \frac{\partial \left(
\left[
\begin{array}{cccc}
\Tran{\vec{x}} \vec{a}_1^c & \Tran{\vec{x}} \vec{a}_2^c & \cdots & \Tran{\vec{x}} \vec{a}_n^c
\end{array}
\right]
\right)} {\partial {x_{k}}} \cdot \vec{x} \\
&= 2 \cdot \left[
\begin{array}{cccc}
a_{k1} & a_{k2} & \cdots & a_{kn}
\end{array}
\right] \cdot \vec{x} \\
&= 2\cdot \vec{a}_k ^r \cdot \vec{x} \\
&= 2 \cdot k\text{-th row of }A\vec{x}
\end{align}$$

Then we have

$$
\begin{align}
\frac{\partial r\left(\vec{x}\right)} {\partial x_k} &= \frac{\frac{\partial \left(\Tran{\vec{x}}A\vec{x} \right)} {\partial {x_{k}}} \cdot \left(\Tran{\vec{x}}\vec{x}\right) - \left(\Tran{\vec{x}} A \vec{x}\right) \cdot \frac{ \partial \left(\Tran{\vec{x}}\vec{x}\right)} {\partial {x_{k}}}} {\left( \Tran{\vec{x}}\vec{x} \right)^2} \\
&= \frac{2\cdot \left( A \vec{x} \right)_{k\text{-th row}}\left(\Tran{\vec{x}}\vec{x} \right) - \left(\Tran{\vec{x}} A \vec{x}\right)\cdot \left( 2x_k \right)} {\left( \Tran{\vec{x}}\vec{x} \right)^2} \\
&= \frac{2} {\Tran{\vec{x}}\vec{x}} \left( A\vec{x} - r\left(\vec{x}\right) \vec{x} \right)_{k\text{-th row}}
\end{align}$$

So that

$$\nabla r\left( \vec{x} \right) = \Tran{\left[
\begin{array}{cccc}
\frac{\partial r\left(\vec{x}\right)} {\partial x_1}&\frac{\partial r\left(\vec{x}\right)} {\partial x_2} & \cdots & \frac{\partial r\left(\vec{x}\right)} {\partial x_n}
\end{array}
\right]} = \frac{2} {\Tran{\vec{x}}\vec{x}} \left( A\vec{x} - r\left(\vec{x}\right) \vec{x} \right)$$

We can plug in eigenvector $\vec{q}$ for check

$$\nabla r\left( \vec{q} \right) = \frac{2} {\Tran{\vec{q}}\vec{q}} \left( A\vec{q} - \lambda \vec{q} \right) = \vec{0}$$

Then the Taylor expansion of $r\left( \vec{q} \right)$ at $\vec{q}$ becomes

$$\boxed{ r\left( \vec{q} \right) = \lambda + \frac{\Tran{\left( \vec{x} - \vec{q} \right)} \nabla^2 r\left(\tilde{\vec{q}}\right)\left( \vec{x} - \vec{q} \right)} {2} }$$

Assuming that $\left\| \nabla^2 r\left( \tilde{\vec{q}} \right) \right\|$ is bounded by a constant, then we have

$$\left| r\left( \vec{x} \right) - \lambda \right| = O \left( \left\| \vec{x} - \vec{q} \right\|^2 \right)$$

$Conclusion$

Given a symmetric matrix $A$, if $\vec{x}$ is an approximation of an eigenvector of $A$, then its Rayleigh quotient $r\left(\vec{x}\right)$ approximates the corresponding eigenvalue of $A$. Besides, the accuracy is even better than the approximation of the eigenvector.

What's more, if $\vec{x}$ is normalized, $i.e.$, if $\left\|\vec{x}\right\|_2 = 1$, then the Rayleigh quotient of $\vec{x}$ is $r\left( \vec{x} \right) = \Tran{\vec{x}}A \vec{x}$.

# A few properties of eigenvalues
## Inversion
If $A$ is nonsingular and $A\vec{x} = \lambda \vec{x}$, for $\vec{x} \neq \vec{0}$, then $A^{-1} \vec{x} = \ffrac{1} {\lambda} \vec{x}$. In another word, if $\lambda$ is an eigenvalue of $A$, then $\lambda^{-1}$ is the eigenvalue of $A^{-1}$, and they share the same eigenvectors.

## Powers
If $A\vec{x} = \lambda \vec{x}$, then for any integer $k$, we have $A^k \vec{x} = \lambda^k \vec{x}$, $i.e.$,  taking the $k\text{-th}$ power of a matrix also takes the $k\text{-th}$ power of its eigenvalues, while the eigenvectors are the same.

## Shift
If $A\vec{x} = \lambda \vec{x}$, the for any value $\sigma$, $\left( A - \sigma I\right) \vec{x} = \left( \lambda  - \sigma \right) \vec{x}$. Thus, the eigenvalues of the matrix $A - \sigma I$ are shifted from the eigenvalues of
$A$ by $\sigma$, and again their eigenvectors are unchanged.

# Power iteration
Given a real symmetrix $n\times n$ matrix $A$, it has $n$ real eigenvalues $\lambda_1,\lambda_2,\dots, \lambda_n$ and a set of orthonormal eigenvectors $\vec{q}_1, \vec{q}_2, \dots, \vec{q}_n$, which satisfy $A\vec{q}_i = \lambda_i \vec{q}_i$ for $i = 1,2,\dots, n$. And by default, the eigenvalues follow the order $\left|\lambda_1\right| \geq \left|\lambda_2\right| \geq \cdots \geq \left|\lambda_n\right|$.

Since now $\vec{q}_1, \vec{q}_2, \dots, \vec{q}_n$ is also a base of $n\text{-dimensional}$ vector space, so that any vector $\vec{v}^{(0)}$ can be written as 

$$
\begin{array}{rl}
\vec{v}^{(0)}\!\!\! & =\; a_1 \vec{q}_1 + a_2 \vec{q}_2 + \cdots a_n \vec{q}_n \\
A\vec{v}^{(0)}\!\!\! & =\;  A \left( a_1 \vec{q}_1 + a_2 \vec{q}_2 + \cdots a_n \vec{q}_n \right) \\
& =\;  a_1 \lambda_1 \vec{q}_1 + a_2 \lambda_2 \vec{q}_2 + \cdots a_n \lambda_n \vec{q}_n \\
A^2\vec{v}^{(0)}\!\!\! & =\;  A \left( a_1 \lambda_1 \vec{q}_1 + a_2 \lambda_2 \vec{q}_2 + \cdots a_n \lambda_n \vec{q}_n \right)\\
& =\;  a_1 \lambda_1^2 \vec{q}_1 + a_2 \lambda_2^2 \vec{q}_2 + \cdots a_n \lambda_n^2 \vec{q}_n \\
& \;  \vdots \\
A^k\vec{v}^{(0)}\!\!\! &= \;  A \left( a_1 \lambda_1^{k-1} \vec{q}_1 + a_2 \lambda_2^{k-1} \vec{q}_2 + \cdots a_n\lambda_n^{k-1} \vec{q}_n \right) \\
& =\;  a_1 \lambda_1^k \vec{q}_1 + a_2 \lambda_2^k \vec{q}_2 + \cdots a_n \lambda_n^k \vec{q}_n \\
& =\;  a_1 \lambda_1^k \left( \vec{q}_1 + \frac{a_2} {a_1} \left( \frac{\lambda_2} {\lambda_1} \right)^k \vec{q}_2 + \cdots + \frac{a_n} {a_1} \left( \frac{\lambda_n} {\lambda_1}\right)^k \vec{q}_n \right)
\end{array}
$$

Now we need two stronger prerequisites: 
1. $\left|\lambda_1\right| > \left|\lambda_2\right| \geq \cdots \geq \left|\lambda_n\right|$
2. the initial vector $\vec{v}^{(0)}$ has a nonzero component in the direction of eigenvector $\vec{q}_1$.

$$\lim_{k \to \infty} \frac{A^k \vec{v}^{(0)}} {a_1 \lambda_1^k} \to \vec{q}_1$$

And the speed of convergence is determined by 

$$\left\| \frac{A^k \vec{v}^{(0)}} {a_1 \lambda_1^k} - \vec{q}_1 \right\| = O\left( \left| \frac{\lambda_2} {\lambda_1} \right|^k \right)$$

meaning that the convergence will be faster if $\left| \lambda_2 \right|$ is much smaller than $\left| \lambda_1 \right|$, and slower if $\left| \lambda_2 \right|$ is closer to $\left| \lambda_1 \right|$.

In addition, we can also see that $\left\| \ffrac{A^{k-1} \vec{v}^{(0)}} {a_1 \lambda_1^{k-1}} - \vec{q}_1 \right\| = O\left( \left| \ffrac{\lambda_1} {\lambda_2} \right|^{k-1} \right)$, so that

$$\left\| \frac{A^k \vec{v}^{(0)}} {a_1 \lambda_1^k} - \vec{q}_1 \right\| = \left\| \frac{A^{k-1} \vec{v}^{(0)}} {a_1 \lambda_1^{k-1}} - \vec{q}_1 \right\| \cdot O\left( \left| \frac{\lambda_2} {\lambda_1} \right| \right)$$

which implies that *the convergence rate is linear* and is reduced by a factor of $\ffrac{\lambda_2} {\lambda_1}$ after each step. So that we define the sequence to approximate the eigenvector corresponding the largest eigenvalue of $A$.

$$
\begin{align}
\text{Initial } \vec{v}^{(0)}, \left\| \vec{v}^{(0)}\right\| &= 1 \\
\vec{v}^{(1)} &= \frac{A \vec{v}^{(0)}} {\left\| A \vec{v}^{(0)} \right\|} \\
\vec{v}^{(2)} &= \frac{A \vec{v}^{(1)}} {\left\| A \vec{v}^{(1)} \right\|} \\
& \vdots \\
\vec{v}^{(k)} &= \frac{A \vec{v}^{(0)}} {\left\| A \vec{v}^{(0)} \right\|} \\
& \vdots \\
\end{align}
$$

so that $\vec{v}^{(k)}$ converges to the eigenvector $\vec{q}_1$ corresponding the largest eigenvalue of $A$.

$Algorithm$

Given an initial $\vec{v}^{(0)}$ with $\left\| \vec{v}^{(0)}\right\| = 1$

```
for k = 1:n
    v = A v ^{(k − 1)}
    v ^{(k)} = v / \|v\|
    \lambda ^{(k)} = v^{(k)}' A v^{(k)}
end
```

We can check $\left\| A \vec{v}^{(k)} - \lambda^{(k)}\vec{v}^{(k)} \right\|$ so that jump out of the loop when it's sufficiently small.
***
$$\left\| \vec{v}^{(k)} - \vec{q}_1 \right\| = \left\| \vec{v}^{(k-1)} - \vec{q}_1 \right\|\cdot O\left( \left| \frac{\lambda_2} {\lambda_1} \right| \right) \Longrightarrow \left\| \vec{v}^{(k)} - \vec{q}_1 \right\| = O\left( \left| \frac{\lambda_2} {\lambda_1} \right|^{k} \right)$$

$$\left| \lambda^{(k)} - \lambda_1 \right| = \left| \lambda^{(k-1)} - \lambda_1 \right|\cdot O\left( \left| \frac{\lambda_2} {\lambda_1} \right|^{2} \right) \Longrightarrow \left| \lambda^{(k)} - \lambda_1 \right| = O\left( \left| \frac{\lambda_2} {\lambda_1} \right|^{2k} \right)$$

$\dagger$ Here there's a square and it's because the conclusion from the Rayleigh Quotient. Same for the Rayleigh quotient iteration.$\Join$

# Inverse iteration
This is used to approximate an eigenvector corresponding to the eigenvalue closest to an initially specified value.

Given real value $\mu$, the eigenvalues of the matrix $\left( A - \mu I \right)^{-1}$ are simply $\left( \lambda_i - \mu \right)^{-1}$ where $\lambda_i$ are the eigenvalues of $A$. Assuming that $\lambda_j$ is the closest one to $\mu$, and the next one $\lambda_l$, $i.e.$, we have the following relation

$$\left| \frac{1} {\lambda_j - \mu} \right| > \left| \frac{1} {\lambda_l - \mu} \right| \geq \cdots$$

So we apply the power iteration to matrix $\left( A - \mu I \right)^{-1}$ and it will converge to the eigenvector of the matrix $\left( A - \mu I \right)^{-1}$ corresponding to its largest eigenvalue $\frac{1} {\lambda_j- \mu}$, which is the exact eigenvector $\lambda_j$ of $A$.

$Algorithm$

Given an initial $\vec{v}^{(0)}$ with $\left\| \vec{v}^{(0)}\right\| = 1$, and a number $\mu$

```
for k = 1:n
    v = (A - \mu I) ^{-1} v ^{(k − 1)}
    v ^{(k)} = v / \|v\|
    \lambda ^{(k)} = v^{(k)}' A v^{(k)}
end
```
***
And with same logic we have the following

$$\left\| \vec{v}^{(k)} - \vec{q}_j \right\| = \left\| \vec{v}^{(k-1)} - \vec{q}_j \right\|\cdot O\left( \left| \frac{\lambda_j - \mu} {\lambda_l - \mu} \right| \right) \Longrightarrow \left\| \vec{v}^{(k)}-\vec{q}_j \right\| = O\left( \left| \frac{\lambda_j - \mu} {\lambda_l - \mu} \right|^{k} \right)$$

$$\left| \lambda^{(k)} - \lambda_j \right| = \left| \lambda^{(k-1)} - \lambda_j \right|\cdot O\left( \left| \frac{\lambda_j - \mu} {\lambda_l - \mu} \right|^{2} \right) \Longrightarrow \left| \lambda^{(k)} - \lambda_j \right| = O\left( \left| \frac{\lambda_j - \mu} {\lambda_l - \mu} \right|^{2k} \right)$$

This process is more expensive because more calculations are required when solving the linear system, the first step of the iteration, $\vec{v} = \left(A - \mu I \right) ^{-1} \cdot \vec{v} ^{(k − 1)} \Longleftrightarrow \left(A - \mu I \right)\vec{v} = \vec{v} ^{(k − 1)}$. But luckly we only need to implement $LU$or Cholesky  factorization for once the entire time.

# Rayleigh quotient iteration
In the inverse iteration, $\lambda^{(k)}$ converges to $\lambda_j$ and the convergence rate depends on $\mu$. The convergence will become faster if $\mu$ is much closer to $\lambda_j$.

So that we replace $\mu$ with the latest $\lambda^{(k)}$ at each step of the inverse iteration.

$Algorithm$

Given an initial $\vec{v}^{(0)}$ with $\left\| \vec{v}^{(0)}\right\| = 1$, and a number $\mu$

```
\lambda ^{(0)} = v^{(0)}' A v^{(0)}

for k = 1:n
    v = (A - \lambda ^{(k-1)} I) ^{-1} v ^{(k − 1)}
    v ^{(k)} = v / \|v\|
    \lambda ^{(k)} = v^{(k)}' A v^{(k)}
end
```
***
$$\begin{align}
\left\| \vec{v}^{(k)} - \vec{q}_j \right\| &= \left\| \vec{v}^{(k-1)} - \vec{q}_j \right\|\cdot O\left( \left| \frac{\lambda_j - \mu} {\lambda_l - \mu} \right| \right) \\
&= \left\| \vec{v}^{(k-1)} - \vec{q}_j \right\|\cdot O\left( \left| \frac{\lambda_j - \lambda ^{(k-1)}} {\lambda_l - \lambda ^{(k-1)}} \right| \right) \\
&= \left\| \vec{v}^{(k-1)}-\vec{q}_j \right\|\cdot O\left( \left| {\lambda_j- \lambda ^{(k-1)}} \right| \right) \\
&= O\left( \left\| \vec{v}^{(k-1)}-\vec{q}_j \right\|^3 \right)
\end{align}$$

$$\left|\lambda^{k}-\lambda_j\right| = O\left( \left\| \vec{v}^{(k)}-\vec{q}_j \right\|^2 \right) = O\left( \left\| \vec{v}^{(k-1)}-\vec{q}_j \right\|^6\right)$$

And still it's more expensive because at each step you need to apply the $LU$ or Cholesky factorization to solve the linear system.

# Computational cost
If $A$ is an $m \times m$ dense matrix, then the computational cost of the [power iteration](#Power-iteration) for each step is $O(m^2)$ for multiplying a matrix with a vector. For the inverse iteration, one $LU$ factorization of $A-\lambda I$ costs $O(m^3)$ ops, which only needs to be computed once in the inverse iteration, and then in each step of the iteration, the cost is $O(m^2)$ for the forward and backward substitutions. For the Rayleigh quotient iteration, the cost for each step is $O(m^3)$ for solving
the system of linear equations.

However, the computation cost can be greatly reduced, especially for the Rayleigh quotient iteration, if we first reduce the symmetric matrix $A$ to a tridiagonal for by a sequence of symmetric orthogonal transformations which does not change the eigenvalues of $A$.



# Reduction to Tridiagonal form

We know that for a real symmetric matrix $A$, it has an eigenvalue decomposition $A = Q \Lambda Q^{\mathrm{T}}$, where $\Lambda$ is a diagonal matrix which contains the eigenvalues of $A$, and $Q$ is an orthogonal matrix which contains the eigenvectors of $A$.

However in general we can only multiply a sequence of orthogonal matrices to the left and their transpose to the right of $A$ such that the result is a tridiagonal matrix. And such orthogonal matrices can be determined by the
using Householder reflector.

$$\underset{Q^{\mathrm{T}}}{\underbrace{Q_{m-2}\cdots Q_{2}Q_{1}}} \;\; A \;\; \underset{Q}{\underbrace{Q_{1}^{\mathrm{T}}Q_{2}^{\mathrm{T}}\cdots Q_{m-2}^{\mathrm{T}}}} = T = \left[ \begin{array}{ccccc}
\times & \times & & & \\
\times & \times & \times & & \\
& \times & \times & \times & \\
& & \times & \times & \times \\
& & & \times & \times
\end{array}
\right]$$

We can also see that $A$ is similar to $T$ and they have the same eigenvalues.

$Algorithm$

```
for k = 1 to m - 2
    x = A_{k+1:m,k}
    v_{k} = x + sign( x_{1} ) \| x \|_{2} e_{1}
    v_{k} = v_{k} / \| v_{k} \|_{2}
    A_{k+1:m,k:m} = A_{k+1:m,k:m} - 2 v_{k} ( v_{k}' A_{k+1:m,k:m} )
    A_{k:m,k+1:m} = A_{k:m,k+1:m}  2 ( A_{k:m,k+1:m} v_{k} ) v_{k}'
end
```

In general, an eigenvalue problem solver first apply the above algorithm to reduce $A$ into a tridiagonal form $T$. Then an iteration after that. And for the reduction part the cost for each iteration is also $O\left(m^{2}\right)$

eb794246ffe48bf123a8463411ec5219a3dc17f9
ipynb
Jupyter Notebook
Computational/Intro to Numerical Computing/Note_Chap09_Power iteration, inverse iteration, and Rayleigh quotient iteration.ipynb
XavierOwen/Notes
d262a9103b29ee043aa198b475654aabd7a2818d
MIT
2018-11-27T10:31:08.000Z
2019-01-20T03:11:58.000Z
Computational/Intro to Numerical Computing/Note_Chap09_Power iteration, inverse iteration, and Rayleigh quotient iteration.ipynb
XavierOwen/Notes
d262a9103b29ee043aa198b475654aabd7a2818d
MIT
Computational/Intro to Numerical Computing/Note_Chap09_Power iteration, inverse iteration, and Rayleigh quotient iteration.ipynb
XavierOwen/Notes
d262a9103b29ee043aa198b475654aabd7a2818d
MIT
2020-07-14T19:57:23.000Z
2020-07-14T19:57:23.000Z
Qwen/Qwen-72B
1. YES
2. YES