Nonlinear Least Square Regression

The nonlinear least square fit

The nonlinear least square fit (regression) minimizes the target action: \[ S(\theta) = \sum_i w_i (y_i - f(x_i; \theta))^2
\] with positive weights $w_i \ge 0$, for a given set of data $\{ (x_1, y_1), (x_2,y_2),\cdots \}$ and a form function $f(x; \theta) \equiv f(x;\theta_1,\theta_2,\cdots \theta_p)$.

The minima lies in the extrema: \[
0 = \frac{\partial S}{\partial \theta_r} = -2 \sum_i w_i (y_i - f(x_i;\theta))\frac{\partial f(x_i;\theta)}{\partial \theta_r}. \quad r = 1,2,\cdots p \]

Gauss-Newton algorithm

The Newton's method for solving for a local root of a function $\varphi(x)$ is, \[
x_{n+1} = x_n - \frac{\varphi'(x_n)}{\varphi(x_n)}
 \] In this problem, $\varphi(\theta) = S'(\theta)$. Therefore, the Newton's method can be realized as \[
\theta^{(k+1)} = \theta^{(k)} - (H S(\theta^{(k)}))^{-1} \nabla S(\theta^{(k)})
\] where $[H S(\theta)]_{rs} = {\partial^2 S}{\partial \theta_r \partial \theta_s}$ is the Hesse matrix (Jacob matrix) of $S$. It is easy to show \[
\left[ H S \right]_{rs} = 2 \sum_i w_i \frac{\partial f_i}{\partial\theta_r}\frac{\partial f_i}{\partial \theta_s} - 2 \sum_i w_i (y_i - f(x_i)) \frac{\partial^2 f_i}{\partial\theta_r\partial\theta_s}
\] Ignore the second term, as the residue and the second derivative should be small. Then the Gauss-Newton algorithm gives the iteration procedure: \[
\theta^{(k+1)} = \theta^{(k)} + \frac{\mathbf{J}^T\cdot \mathbf{W}\cdot (\mathbf{y}-f(\mathbf{x};\theta^{(k)}) )}{\mathbf{J}^T\cdot\mathbf{W}\cdot\mathbf{J}}
\]

Levenberg-Marquardt algorithm

Levenberg-Marquardt algorithm is one of the most important algorithm based on the Gauss-Newton algorithm. Let $\theta_r, r = 1,2,\cdots p$ be the fitting parameters such as $a,b,c$. Define the Jacobian \[
J_{i,r} = \frac{\partial f(x_i;\theta)}{\partial \theta_r}
\] Define $ W_{ij} = w_i \delta_{ij}$. The iterative procedure is: \[
\mathbf{\theta}^{(k+1)} = \mathbf{\theta}^{(k)} + \alpha \frac{\mathbf{J}^T \cdot \mathbf{W} \cdot (\mathbf{y} - f(\mathbf{x};\mathbf{\theta}^{(k)}))}{\lambda \mathbf{\Lambda} + \mathbf{J}^T \cdot \mathbf{W} \cdot \mathbf{J}}
\] where $\mathbf{\Lambda}$ is a diagonal matrix, either chosen as the identity matrix or the diagonal term of $\mathbf{J}^T \cdot \mathbf{W} \cdot \mathbf{J}$, $0<\alpha<1$ is a small positive number.

Therefore, the solution $\bar\theta$ satisfies equation $\mathbf{J}^T_{\bar\theta} \cdot \mathbf{W} \cdot (\mathbf{y} - f(\mathbf{x};\bar{\mathbf{\theta}})) = 0$ or \[
\sum_i w_i  \frac{\partial f(x_i; \bar\theta)}{\partial \theta_r} (y_i - f(x_i; \bar\theta)) = 0
\] The choice of $\lambda$ is very important. Nonzero $\lambda$ can avoid local minimum but large $\lambda$ is likely to slow down the convergence. A useful strategy is to to choose a nonzero $\lambda$ to start with and decrease $\lambda$ as the iteration converges.

Parameter Error Estimation

It is important to attach estimated error bars/bands to the fitted function (or the fitted parameters). Let's write the data as $ y_i = f(x_i) + \xi_i$, where $\xi_i$ are the residuals. If we assume $\xi_i, i=1,2,3,\cdots n$ are $n$ independent random variables in Normal distribution. The likelihood of the fitted parameters can be written as \[
L(\theta_r, \sigma)
= N_\sigma \exp\left[ -\frac{1}{2\sigma^2}\sum_i w_i (y_i-f(x_i;\theta_r))^2\right] \\
\simeq N'_\sigma \exp\left[-\frac{1}{2\sigma^2} \sum_{rs} (\theta_r-\bar\theta_r) A_{rs} (\theta_s-\bar\theta_s) \right]
\] where $A_{rs} = \sum_i w_i  \frac{\partial f(x_i;\bar\theta)}{\partial\theta_r} \frac{\partial f(x_i;\bar\theta)}{\partial\theta_s} = \left[ \mathbf{J}^T \cdot \mathbf{W} \cdot \mathbf{J} \right]_{rs}
$  and $N_\sigma, N'_\sigma$ are normalization factors.

The variance $\sigma^2$ can be estimated by the estimator \[
\hat{\sigma}^2 = \frac{1}{n-p}\sum_i w_i (y_i - f(x_i; \bar\theta))^2 = \frac{S({\bar\theta})}{n-p}.
 \] The covariance matrix of the parameters is, \[ \text{Cov}(\theta) = \hat\sigma^2 A^{-1}
\] It's also useful to find the prediction band. The width of the band is the variance of function $f(x; \theta)$. The "Delta Method" gives an approximative way to calculate $\text{Var}\left[ f(x; \theta) \right]$ from $\text{Cov}[ \theta ]$. \[
\sigma_f(x) = \text{Var}\left[ f(x; \bar\theta) \right] = \frac{\partial f^T(x;\bar\theta)}{\partial\mathbf{\theta}}\cdot\text{Cov}( \mathbf\theta) \cdot \frac{\partial f(x;\bar\theta)}{\partial\mathbf{\theta}} \\
 = \hat\sigma^2 \sum_{r,s} \frac{\partial f(x;\bar\theta)}{\partial\theta_r}\left[ A^{-1} \right]_{rs}\frac{\partial f(x;\bar\theta)}{\partial\theta_s}
\]
Finally,the distribution of the residuals can be test by the standard $\chi^2$-test.

Example: a power-law fitting

In our example, we will explore fitting some data to the power-law function $f(x) = a + b/x^c $.

Fig. 1: (Color Lines) Black, the theoretical curve; Red, fit with uniform weights; Green fit with weights proportional to $x$

The corresponding covariance matrices for the parameters $a,b,c$.


The diagonal terms gives the variance of the parameters.

On the Parallelism of Symmetric Sparse Matrix Vector Multiplication

Symmetric Sparse Matrix 

Sparse matrices are the matrices populated primarily with zeros.  Iterative linear solvers, such as the Lanczos-Arnoldi procedure, are among the most important computational tools for solving large symmetric (Hermitian) sparse linear problems. The fundamental operation in the iterative methods is the sparse matrix vector multiplication (SpMV).

Sparse matrix can be stored in several formats. The most flexible one is the coordinate storage scheme (COO). In this scheme, one use three arrays (two integral and one scalar):

row(nnz), col(nnz), mat(nnz)

to store the row indices, column indices and the values of the nonzero matrix elements. It's easy to convert a coordinate storage scheme into other schemes such as the compressed column/row scheme (CCS/CRS). For SpMV, coordinate scheme is as efficient as CCS/CRS. It demands $2 n_{nz} + n_{nz}$ memory comparing with CCS/CRS $(n_{nz} + n) + n_{nz}$, where $n$ is the dimensionality of the matrix and $n_{nz}$ is the number of nonzeros, and $n^2 \gg n_{nz} \gg n$ for most useful cases. Only half of the matrix is stored, because the matrix is symmetric (Hermitian).

Serial Algorithm

subroutine SpMV(nnz, row, col, mat, n, x, y, ierr)
implicit none
! A*x + y -> y
! n is the dimension of the vector x, y
! nnz is the number of nonzeros
integer, intent(in) :: nnz, n, row(nnz), col(nnz)
double precision, intent(in) :: mat(nnz), x(n)
double precision, intent(inout) :: y(n)
integer, intent(out) :: ierr
integer :: i

  do i = 1, nnz

    y(row(i)) = y(row(i)) + mat(i)*x(col(i))
    if(row(i) .ne. col(i)) then
      y(col(i)) = y(col(i)) + x(row(i))*conjg(mat(i))
    end if

  end do
  ierr = 0
end subroutine

Parallel Algorithm

Here we only talk about MPI parallelism. In reality, it is most efficient to use both MPI and OpenMP. The OpenMP parallelism can be used in the series part.

2D cyclic distribution

A parallel SpMV depends on the the distribution of the processors across the processors. 1D distribution has a serious load balancing problem. An simple and popular scheme is the 2D cyclic distribution.

Left: A naive 2D distribution. Center: A cyclic 2D distribution. Right: A cyclic 2D distribution with non-square block matrices.

In this scheme, the $n_p$ processors are aligned as a $n_d \times n_d$ 2D grid. For symmetric matrices, only half of the grid is needed. Then, for $n_p$ processors, $ n_d (n_d + 1)/2 \le n_p \implies n_d \le \text{floor}\left( \frac{\sqrt{1+n_p}-1}{2}\right)$. It is useful to assign each processor a 2D coordinate $(I, J)$, where $I,J = 1,2,3,\cdots, n_d$. We will discuss the processor assignment later.

The vector $x$ is stored on the diagonal and the matrix $A$ is distributed across the processor grid. More specifically, vector element $x_i$ is stored in processor $(I, I)$, where $I = \mod(i, n_d) + 1$. Matrix element $A_{ij}$ is distributed to processor $(I, J)$, $I = \mod(i, n_d)+1, J = \mod(j,n_d)+1$, forming a local block matrix $A_{IJ}$.

The local blocks need not be square matrix, except the diagonal blocks. Let $n_I$ be the dimensionality of the $I$-th diagonal blocks. It is easy to see, the row and column dimensionality of the block $A_{IJ}$ is $n_I$ and $n_J$, respectively. \[ n_I = \text{floor}\left(\frac{n-1}{n_d}\right) + \vartheta\left(n-\text{floor}\left(\frac{n-1}{n_d}\right)\cdot n_d - I \ge 0\right) \]

Processor Assignment: Diagonal Major (DM) Scheme

In the DM scheme, the processors are assigned to the matrix blocks according to their diagonal coordinate:

Matrix block $A_{IJ}$ is on processor $ p = F(I-J-1)+J-1$, where $F(x) = (x+1) (2 n_d - x)/2$. In practice, one knows the local process id $p$ and wants to determine which matrix block is stored on it: \[
\begin{split}
I &= 1+p+1+G(p+1)-F(G(p+1)-1); \\
J &= p+1-F(G(p+1)-1).
\end{split}
\] where $G(x) = \text{floor}\left(n_d - \sqrt{\left(n_d+\frac{1}{2}\right)^2-2x} - \frac{1}{2}\right)$.

To store the node information, we introduce a FORTRAN module (a struct is a similar structure in C):

module nodeinfo
implicit none

integer :: nproc, myrank, myrow, mycol, ndiag, idiag;
integer :: proc_comm, row_comm, col_comm;

contains

subroutine initnodeinfo
use mpi
implicit none
integer :: ierr

! determine ndiag, nproc and rank
  call mpi_comm_size(mpi_comm_world, nproc, ierr)
  call mpi_comm_rank(mpi_comm_world, myrank, ierr)
  ndiag = floor(0.5*sqrt(1.+8.*nproc)-0.5)
  
  color = 0
  if(ndiag*(ndiag+1)/2>rank) then
    color = 1
  endif

  call mpi_comm_split(mpi_comm_world, color, myrank, proc_comm, ierr)
  
  call mpi_comm_size(proc_comm, nproc, ierr)
  call mpi_comm_rank(proc_comm, myrank, ierr)

! determine myrow and mycol i.e. (I,J) the block ID
  mycol = 1+rank-F(G(myrank,ndiag),ndiag)
  myrow = G(myrank,ndiag)+mycol+1
  idiag = myrow - mycol

! create row_comm and col_comm, diag_comm
  call mpi_comm_split(proc_comm, myrow, abs(idiag), row_comm, ierr)
  call mpi_comm_split(proc_comm, mycol, abs(idiag), col_comm, ierr)
  call mpi_comm_split(proc_comm, idiag, myrank, diag_comm, ierr)
  
end subroutine
  
integer function F(x, nd)
implicit none
integer, intent(in) :: x, nd

  F = (x+1)*(2*nd-x)/2

end function

integer function G(x, nd)
implicit none
integer, intent(in) :: x, nd

  G = floor(nd-0.5-sqrt((nd+0.5)**2-2*x)-0.5)

end function

integer function nloc(n, nd, i)
implicit none
integer, intent(in) :: n, nd, i

  nloc = (n-1)/nd
  if(n.ge.((n-1)/nd)*nd+i) then
    nloc = nloc + 1
  end if

end function

end module

The problem of the DM scheme is that the number of processors in each row & column communicator is not the same. Indeed, it is not well balanced: some rows & columns have $n_d$ (ndiag) processors, some of them have only 1. This problem can be solved by using the so-called Balanced Column Major scheme.

Processor Assignment: Balanced Column Major (BCM) Scheme

In the BCM scheme, each row or column has $\frac{n_d+1}{2}$ blocks.

The corresponding matrix blocks distribution on processor $p$ is,
if $ p\le \text{ceiling}\left(\frac{n_d+1}{2}\right)\left[ n_d+1-\text{ceiling}\left(\frac{n_d+1}{2}\right)\right]-1$ then \[ \begin{split}
J &= \text{floor}\left(\frac{p}{\text{ceiling}\left(\frac{n_d+1}{2}\right)}\right)+1; \\ I &= \mod\left( p, \text{ceiling}\left(\frac{n_d+1}{2}\right)\right)+J;
\end{split} \]
else,
\[ \begin{split} J &= n_d+2-\text{ceiling}\left(\frac{n_d+1}{2}\right)+\text{floor}\left( \frac{q}{\text{floor}\left(\frac{n_d+1}{2}\right)}\right); \\ I &= \mod\left( \mod\left(q, \text{floor}\left(\frac{n_d+1}{2}\right)\right) +J-1, n_d\right)+1; \end{split} \] where $q = p-\text{ceiling}\left(\frac{n_d+1}{2}\right)\left[n_d+1-\text{ceiling}\left(\frac{n_d+1}{2}\right)\right]$.

The nodeinfo.mod can be modified accordingly.

The Parallel SpMV Algorithm Based on the 2D cyclic Distribution

$A = L + D + L^t$, where $L$ is a strict lower triangular matrix and $D$ is a diagonal matrix. $T = L+D$.
\[ A x + y \to y \Rightarrow T x + (x^t L)^t + y \to y \]

subroutine SpMV(nnz, row, col, mat, ncol, x, yt, nrow, y, xt, ierr)
use nodeinfo
use mpi
implicit none
! A*x + y -> y
! A = L + D + L^t, where L is a strict lower triangular matrix, 
! D is a diagonal matrix, L^t is the transpose of L.
! Let T = L+D; A = T + L^t.  
! A*x = T*x + L^t*x = T*x + (x^t*L)^t
! 
! x/y is the local vector. Only x/y on the diagonal are initialized/output
! xt/yt is the local vector of the transpose vector of x/y
! note their local dimensionalities nrow, ncol are not the same as the original vector.
integer :: nnz, ncol, nrow, ierr, row(nnz), col(nnz);
double precision :: x(ncol), yt(ncol), y(nrow), xt(nrow), mat(nnz);
intent(in) :: nnz, ncol, nrow, row, col, mat;

! (L+D)*x + y -> y
  call mpi_bcast(x, ncol, mpi_real8, 0, col_comm, ierr)
  
! local SpMV
  if(idiag.eq.0) then
    xt = y
  end if
  do i = 1, nnz
    xt( 1+(row(i)-1)/ndiag ) = xt( 1+(row(i)-1)/ndiag ) + mat(i)*x( 1+(col(i)-1)/ndiag )
  end do

  call mpi_reduce(xt, y, nrow, mpi_real8, mpi_sum, 0, row_comm, ierr)

! x^t * L -> y^t
  if(idiag.eq.0) then
    xt = x
    x = y
  end if
  call mpi_bcast(xt, nrow, mpi_real8, 0, row_comm, ierr)
  
! local SpMV
  do i = 1, nnz
    if( row(i) .eq. col(i) ) then
      x( 1+(col(i)-1)/ndiag ) = x( 1+(col(i)-1)/ndiag ) + xt( 1+(row(i)-1)/ndiag )*conjg(mat(i))
    end if
  end do

  call mpi_reduce(x, yt, ncol, mpi_real8, mpi_sum, 0, col_comm, ierr)

  if(idiag.eq.0) then
    y = y + yt
    x = xt
  endif

end subroutine

Sleep Sort

The ingenious sleep sort algorithm makes use of the system (hardware) scheduler to sort a list of data. The original version was implemented with shell scipt:

#!/bin/bash
function f() {
    sleep "$1"
    echo "$1"
}
while [ -n "$1" ]
do
    f "$1" &
    shift
done
wait

example:

Of course, it can be implemented with any language that has sleep/schedule functions. A couple of examples have been shown in this wiki.

This algorithm is not as hilarious as it seems. It only takes O(n) in the idea case (ignore the time cost in scheduling and output). This is faster than fastest comparing sort algorithm such as quick sort, which takes O(n log n). The sleep time is proportional to the hardware time. It can be scaled down  dramatically. For example we can set the sleeping time to be the cpu clock period.

People have discussed the implementation of scheduling and argued that if the scheduling is achieved via priority queue it will eventually take O(n log n). But if we have sufficient "timers/schedulers", this algorithm is equivalent to bead sort.

on the algorithms of binomial coefficients

factorial based algorithms: 
The most naive definition of binomial coefficients is: \[ { n \choose m } = \frac{n!}{m! (n-m)!}, \] where $n!$ is factorial of $n$. However this formula is not a practical algorithm. The reason is the value of factorial increases extremely fast. For example if we want to calculate ${13 \choose 6 }$. Using above algorithm, $13! = 6227020800 > 2^{32} = 4294967296$ overflows the typical unsigned 4-byte integer. On the other hand, ${13 \choose 6 }= 1716$. In addition, the time complexity is $T(n, m) = O( n )$.

A more practical variation of this formula is: \[ { n \choose m } = \prod_{k=1}^m \frac{n-m+k}{k} = \frac{n}{m} \cdot \frac{n-1}{m-1} \cdots \frac{n-m+1}{1}. \] Also notice that $ { n \choose m } = { n \choose n - m }$. One can choose the smaller one to do computing: \[ { n \choose m } = \prod_{k=1}^{\bar{m}} \frac{n-\bar{m}+k}{k} \] where $\bar{m} = \min\{ m, n-m \}$. It can be shown the time complexity is only $T(n, m) = O( \bar{ m } )$.

recursive algorithms:
A recursive algorithm based on above algorithm is, \[ { n \choose m } = \frac{n}{m} \cdot { n - 1 \choose m - 1} \] with initial values $ {k \choose 0} = { 1 \choose l } = 1 \quad k, l \in \mathbb{N}, l \le 1$. $T(n,m) = O(\bar{m})$. Another recursive algorithm, known as Pascal's triangle (Meru-prastaara, Khayyam triangle, Yang Hui's triangle, Tartaglia's triangle) is, \[ {n \choose m} = { n - 1 \choose m - 1} + { n - 1 \choose m }\] with the same initial values above. $T(n,m) = O( n^2 )$.

falling factorial and gamma function based algorithm for generalized binomial coefficients:
Another thing one often encounters is the generalized binomial coefficients. For instance, the generalized binomial theorem states: \[ (x + y)^\alpha = \sum_{k=0}^\infty {\alpha \choose k} x^k y^{\alpha - k} \qquad ( x, y \in \mathbb{C}, | x | < | y |, \alpha \in \mathbb{R} ) \] where $ { \alpha \choose k } = \frac{\alpha \cdot (\alpha - 1) \cdots ( \alpha - k + 1)}{ k !} \equiv \frac{(\alpha)_k}{k!} $. Symbol $ (\cdot )_k$ is termed falling factorial. The algorithm for computing falling factorial fashioned binomial coefficients is the same as the one mentioned above.

The most general definition of binomial coefficients is via Gamma function or equivalent limit of some Gamma function: \[ { \alpha \choose \beta } \equiv \frac{\Gamma(\alpha + 1)}{\Gamma(\beta + 1) \Gamma(\alpha - \beta + 1)} \] with $\alpha, \beta \in \mathbb{C} $. There are several issues for using gamma function in practical computing. First of all, they are computationally complex. Therefore, only if $m$ in ${n \choose m}$ is not an integer (which is rare in practice), Gamma function definition could be used. In that case, one actually only need the initial values for $\Gamma(s) \quad (0 \le s \le 1)$. Calculation of gamma function can be obtained by Lanczos approximation:
\[
\Gamma(z) = \sqrt{2 \pi} (z + g + \frac{1}{2})^{z + \frac{1}{2}}e^{- (z + g + \frac{1}{2})} A_g(z)
\] where, $g$ is an arbitrary constant subject to restriction $\text{Re}( z + g + \frac{1}{2}) > 0$. The reflection formula \[ \Gamma(z) \Gamma(1-z) = \frac{\pi}{\sin(\pi z)} \] can be used to achieve this. \[ A_g(z) = c_0 + \sum_{k=1}^\infty \frac{c_k}{z + k} \] where $c_k$ are pre-calculable constants. GNU Scientific Library choose $g = 7$, and the first 9 coefficients are given such that 15 correct digits are guaranteed:

$c_0 = 0.99999999999980993$,
$c_1 = 676.5203681218851$,
$c_2 = -1259.1392167224028$,
$c_3 = 771.32342877765313$,
$c_4 = -176.61502916214059$,
$c_5 = 12.507343278686905$,
$c_6 = -0.13857109526572012$,
$c_7 = 9.9843695780195716\text{E}-6$,
$c_8 = 1.5056327351493116\text{E}-7 $.

Stirling's approximation and Spouge's approximation may also be used to calculate gamma function numerically, under different circumstance.

The second issue of gamma function definition, also reviewed in Lanczos approximation, is singularities at negative integers. However, binomial coefficients do not suffer from singularities. M. J. Kronenburg  (arXiv:1105.3689v1 [math.CO], 2011) extended the definition of binomial coefficients to arbitrary integer coefficients,  adopted by Mathematica 8. Let $n, m$ be integers, \[
{ n \choose m } = \left\{ \begin{array}{l l l}
& \frac{n!}{m! (n-m)!}                      &\qquad  (n \ge 0, 0\le m \le n );                    \\
& (-1)^m { -n + m - 1 \choose m }    &\qquad  (n < 0, m \ge 0 );                             \\
& (-1)^{n-m} { -m-1 \choose n - m} &\qquad  (n < 0, m \le n  );                             \\
& 0                                                     &\qquad   \text{ others }.                                \\
\end{array}
\right.
\]

logarithm method:
${ n \choose m }$ is bounded in, \[ \left(\frac{n}{k} \right) ^k \le { n \choose k } \le \left(\frac{n \text{e}}{k} \right)^k. \] This still may be a large number. In actuality, binomial coefficients emerging in intermediate processes may be very large, yet several binomial coefficients cancels out and result a relative small number at the end of the day. Thus it is inappropriate to evaluate binomial coefficients directly. The solution is to calculate \[ \log {n \choose m } = \sum_{k=1}^\bar{m} \log (n-\bar{m}+k) - \log(k) \] for each binomial coefficients, and add/subtract them depending on multiplication/division demand. At the end of the calculation, take exponential. Generalized binomial coefficients may not be positive. When using logarithm method, the absolute value and the sign should be considered separately.  

QR algorithm as a sophistical power iteration

The QR algorithm is one of the most popular diagonalization algorithms. In this post, we compare it with several variations of the power iteration. Throughout this post, we only consider positive-defined hermitian matrices. If a matrix is not positively defined, one can always apply a shift: $A \to A + s I, \quad s = |\min\{\lambda_i, i = 1, 2, \cdots \}| $.

QR algorithm

Recall the QR decomposition (QR factorization) matrix $A$ factorizes $A$ into the matrix product of a orthogonal matrix $Q$ and a upper  triangle matrix $R$, i.e. $A = QR$.  QR decomposition can be seen as the reorthogonalization of the column vectors of $A$. Namely, column vectors $a_i \quad (i = 1, \cdots, n)$ of $A = [a_1, a_2, \cdots, a_n]$ are orthogonalized to $Q = [q_1, q_2, \cdots, q_n]$ such that $a_i \in \text{span}(q_1, q_2, \cdots, q_i)$.  Based on QR decomposition, QR algorithm to diagonalize a matrix $A$ reads:

1. $ A_0 = A $;
2. $ A_{k-1} \to Q_k R_k $ (QR decomposition);
3. $ A_k = R_k Q_k $.

Note that after procedures (2) and (3) $ A_k = Q^\dagger_k A_{k-1} Q_k $. Hence $A_k \simeq A_{k-1} \simeq A_{k-2} \simeq \cdots A_0 = A$. There is a theorem to guarantee that under certain conditions, the sequence converges to diagonal matrix. We won't discuss the theorem here.

Power iteration

Power iteration is fairly simple. Choose a random vector $v = \sum_i \alpha_i x_i$ where $x_i \quad (i = 1, 2, \cdots, n)$ is eigenvector of $A$. Then
\[
A^k v = \sum_i \alpha_i \lambda_i^k x_i.
\] Therefore, if $k$ is sufficiently large, $\lambda_m \equiv \max \{ \lambda_i, i = 1, 2, \cdots, n \}$ dominates, so $A^k v \to \lambda_m^k x_m$. In practice, the resulted vector each step is normalized.

1. $ v_0 = \text{ random unit vector }$;
2. $ v_k = A \cdot v_{k-1} $;
3. $ v_k \leftarrow v_k / \| v_k \|$.

There are two issues in the power iteration. First of all, if the largest eigen-pair (eigenvalue and eigenvector) is degenerate, the iterated vector can only converge to a vector of the eigen-space. Secondly,  it can only generate the largest eigen-pair .

simultaneous power iteration

In order to overcome the first issue, more than one vector can be used simultaneously in iteration. The method is known as simultaneous power iteration. The algorithm can be described as following:

1. $ V_0 = [ v^{(0)}_1, v^{(0)}_2, \cdots, v^{(0)}_m ] \quad v^{(0)}_i \text{ is random vector }, i = 1, 2, \cdots, m; m \le n$;
2. $ V_k = A \cdot V_{k-1} $;
3. $ v^{(k)}_i \leftarrow v^{(k)}_i / \| v^{(k)}_i \|  \quad i = 1, 2, \cdots m$.

The resulted vectors can be used to construct the degenerate eigen-space.

orthogonal power iteration

Simultaneous power iteration does not resolve the second issue. In principle, after obtaining the first eigenvector, it can be subtracted from the entire space. Then perform the power iteration for the rest subspace, the second largest eigen-pair can be obtained. However, this method is never practical. Because any numerical error may add the subtracted eigenvector back. In order to keep the calculated eigen-vector off, we can re-orthogonalize the vectors after each iteration. Notice that QR decomposition as described above can be used for the orthogonalization. This iteration scheme is known as orthogonal power iteration. The algorithm is sketched below:

1. $ V_0 $ is a $n \times n$ random matrix ;
2. $ V_{k-1} \to Q_k R_k $ ( QR decomposition) ;
3. $ V_k = A Q_k $.

Power iteration vs. QR algorithm

Now let's compare orthogonal power iteration and QR algorithm. In orthogonal power iteration, define $ A_k \equiv Q_k^\dagger V_k =Q_k^\dagger A Q_k$, $U_k \equiv Q_{k-1}^\dagger Q_k$. It can be seen

• $U_k$ is orthogonal;
• $A_{k-1} = Q_{k-1}^{\dagger} Q_k R_k = U_k R_k$;
• $A_k = Q_k^\dagger Q_{k-1} A_{k-1} Q_{k-1}^\dagger Q_k = U_k^\dagger A_{k-1} U_k$. 

In fact, $A_{k-1} = U_k R_k$ is a QR decomposition of $A_{k-1}$. That is to say, in orthogonal power iteration, $V_{k-1}$ is QR-factorized, whereas in QR algorithm, $A_{k-1}$ is QR-factorized. These two methods are equivalent. The operations on $A_k$ yields a sequence identical to those in QR algorithm, if we further set $Q_0 = I \implies V_0 = A$. The benefit of QR algorithm is that matrices $U_k$ and $A_k$ are exactly what we need.

This comparison shows that QR algorithm is actually equivalent to some variation of the power iteration. Therefore we could expect it has the similar performance (convergence rate, memory consumption etc.) with the power iteration.