---
output:
  pdf_document: default
  html_document: default
---
<!--
  %\VignetteIndexEntry{EfficientMaxEigenpair Vignette} 
  %\VignetteEngine{knitr::knitr}
  %\VignetteEncoding{UTF-8}
-->
---
title: "EfficientMaxEigenpair Vignette"
header-includes:
    - \usepackage{url}
author:
- Mu-Fa Chen
- Assisted by Xiao-Jun Mao
date: '`r Sys.Date()`'
output: html_document
---

`EfficientMaxEigenpair` is a package for computing the maximal eigenpair of the matrices with non-negative off-diagonal elements (or as a dual, the minimal eigenvalue of the matrices with non-positive off-diagonal elements). This vignette is a simple guide to using the package. All the algorithms and examples provided in this vignette are available in both the paper "Efficient initials for computing maximal eigenpair" and "Global algorithms for maximal eigenpair" by Mu-Fa Chen. The papers [Chen (2016)](http://math0.bnu.edu.cn/~chenmf/) and [Chen (2017)](http://math0.bnu.edu.cn/~chenmf/) are now included in the Vol 4, in the middle of the website:
\begin{center}
\url{http://math0.bnu.edu.cn/~chenmf/}
\end{center}
Let us install and require the package `EfficientMaxEigenpair` first.
```{r, eval = TRUE, echo = FALSE}
source("../R/eff.ini.maxeig.general.R")
source("../R/eff.ini.maxeig.tri.R")
source("../R/eff.ini.maxeig.shift.inv.tri.R")
source("../R/eff.ini.seceig.general.R")
source("../R/eff.ini.seceig.tri.R")
source("../R/ray.quot.tri.R")
source("../R/ray.quot.seceig.tri.R")
source("../R/ray.quot.general.R")
source("../R/ray.quot.seceig.general.R")
source("../R/shift.inv.tri.R")
source("../R/tri.sol.R")
source("../R/find_deltak.R")
source("../R/tridiag.R")
```

```{r, eval = FALSE, echo = TRUE}
require(EfficientMaxEigenpair)
```

Before moving to the details, let us illustrate the algorithm by two typical examples.

#### The tridiagonal case (Example 7 in [Chen (2016)](http://math0.bnu.edu.cn/~chenmf/))
The matrices we considered for the tridiagonal case are in the form of
\begin{center}
$Q=\begin{bmatrix} -1 & 1^2 & 0 & 0 & \dots & 0 \\ 1^2 & -5 & 2^2 & 0 & \dots & 0 \\ 0 & 2^2 & -13 & 3^2 & \dots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots & \vdots \\ 0 & \dots & \dots & (n-2)^2 & -(n-2)^2-(n-1)^2 & (n-1)^2 \\ 0 & \dots & \dots & 0 & (n-1)^2 & -(n-1)^2-n^2 \end{bmatrix}$,
\end{center}
where $n$ is the dimension of the matrix. For the infinite case, we have known that the maximal eigenvalue is $-1/4$. We now want to calulate the maximal eigenvalues of matrices in different dimensions by the Rayleigh quotient iteration.

```{r, echo=TRUE, results='hold', tidy=TRUE}
# Define the dimension of matrix to be 100.
nn = 100

# Generate the corresponding tridiagonal matrix.
a = c(1:(nn - 1))^2
b = c(1:(nn - 1))^2
c = rep(0, length(a) + 1)
c[length(a) + 1] = nn^2

# Output the corresponding convergence maximal eigenpairs of the tridiagonal matrix.
eigenpair=eff.ini.maxeig.tri(a, b, c, xi = 7/8)
```

The output here are the approximations $z$ of the maximal eigenvalue of matrix $Q$. However, in what follow, we often write -$z$ instead of $z$ for simplicity.
```{r, echo=TRUE, results='hold', tidy=TRUE}
# The approximating sequence z_0, z_1, z_2 of the maximal eigenvalue.
eigenpair$z
```

The following figure shows the corresponding eigenvetors $v$ up to a sign. To be consistent with the paper's notation, we accept the convention that the first element in vector start with $v[0]$. From the first figure below, one sees that the eigenvectors $-v[1]$ and $v[2]$ are nearly the same. However, if we look at the last parts of the curves, it is clear that $-v[1]$ and $v[2]$ are actually different. Here we plot them in the same figure to compare. The following two figures show the approximating sequence $v[0]$, $v[1]$, $v[2]$(up to a sign) of the maximal eigenvector.

![$v[0]$, $v[1]$, $v[2]$](v012.pdf){width=280px} ![$v[1]$, $v[2]$](v12.pdf){width=280px}

```{r, eval = FALSE, echo=FALSE, results='hold', tidy=TRUE}
plot(eigenpair$v[[1]], type = "l", col = "black", ylab = "", xlab = "")
lines(-eigenpair$v[[2]], col = "blue")
lines(eigenpair$v[[3]], col = "red")
text(50, 0.5, labels = "v[0], -v[1], v[2] on {1,2,...,100}")
text(90, 0.05, labels = "v[0]", col = "black")
text(20, 0.05, labels = "-v[1] and v[2]", col = "black")
```
```{r, eval = FALSE, echo=FALSE, results='hold', tidy=TRUE}
plot(x = c(91:100), -eigenpair$v[[2]][91:100], type = "l", col = "blue", ylab = "", xlab = "")
lines(x = c(91:100), eigenpair$v[[3]][91:100], col = "red")
text(95.5, 0.0018, labels = "-v[1], v[2] on {91,92,...,100}")
text(91.5, 0.0017, labels = "-v[1]", col = "blue")
text(96, 0.0010, labels = "v[2]", col = "red")
```

In the following, because of the space limitation, we do not report the eigenvector results and write -$z$ instead of $z$ for simplicity.
```{r, eval = TRUE, echo=TRUE, results='hold', tidy=TRUE}
# Define the dimension of each matrix.
nn = c(100, 500, 1000, 5000)

zmat = matrix(0, length(nn), 4)
zmat[ ,1] = nn

for (i in 1:length(nn)) {
    # Generate the corresponding tridiagonal matrix for different dimensions.
    a = c(1:(nn[i] - 1))^2
    b = c(1:(nn[i] - 1))^2
    c = rep(0, length(a) + 1)
    c[length(a) + 1] = nn[i]^2

    # Output the corresponding dual case results, i.e,
    # the minimal eigenvalue of the tridiagonal matrix -Q with non-posivie off-diagonal elements.
    zmat[i, -1] = -eff.ini.maxeig.tri(a, b, c, xi = 7/8)$z[1:3]
}

colnames(zmat) = c("Dimension", "-z_0", "-z_1", "-z_2")
zmat
# The approximating sequence -z_0, -z_1, -z_2 of the maximal eigenvalue.
```

#### The general case (Example 20 in [Chen (2016)](http://math0.bnu.edu.cn/~chenmf/))
The matrix we considered for this general case is
```{r, echo=TRUE, results='hold', tidy=TRUE}
# Generate the general matrix A
A = matrix(c(1, 3, 9, 5, 2, 14, 10, 6, 0, 11, 11, 7, 0, 0, 1, 8),
    4, 4)
A
```

This matrix has complex eigenvalues:
```{r, echo=TRUE, results='hold', tidy=TRUE}
# Calculate the corresponding eigenvalues of matrix A
eigen(A)$values
```

We have to be carefull to choice the initial eigenpairs $z_0$ and $v[0]$ for this matrix A. The following is one counterexample with dangerous initials. The approximating sequence converges to the next to maximal eigenvalue rather than the maximal eigenvalue.
```{r, echo=TRUE, results='hold', tidy=TRUE}
# Calulate the approximating sequence of maximal eigenvalue by the Rayleigh quotient iteration.
eff.ini.maxeig.general(A, z0 = "Auto", digit.thresh = 4)$z
```

Here we fixed the problem by using improved algorithm and safer initials which will be illustrated in detail of the algorithm part.
```{r, echo=TRUE, results='hold', tidy=TRUE}
# Calulate the approximating sequence of maximal eigenvalue by the Rayleigh quotient iteration.
eff.ini.maxeig.general(A, xi = 0.65, digit.thresh = 4)$z
```


Having illustrated the examples, we now introduce the package in details. It offers four main algorithms:

* `eff.ini.maxeig.tri()`: calculate the maximal eigenpair for the tridiagonal matrix. The coefficient `xi` is used to form the convex combination of $1/\delta_1$ and $(v_0,Av_0)_\mu$ (in the dual case $(v_0,-Av_0)_\mu$), it should between 0 and 1. The inner product $(,)_\mu$ here is in the space $L^{2}(\mu)$.  The choice `xi=1` is very safe but less effective. When `xi` is a little away from 1, the algorithm becomes more effective. However, when `xi` is closed to 0, the algorithm becomes less effective again and is even dangerous. The "best" choice of `xi` depends on the model. For different models, one may choose different `xi`. In the present tridiagonal case, we suggest a common choice `xi=7/8` (See Example 7 above/below).
    * The precise level used for output results is set to be `digit.thresh = 6` which implies 1e-6 without any special requirement. Same for the following three algorithms and the examples shown in this vignette.
    * This algorithm works not only for tridiagonal matrix, but also for diagonally dominant matrices. In the latter case, one simply picks up the diagonal part from the original one, as illustrated by the last part of Example 17 in the paper [Chen (2016)](http://math0.bnu.edu.cn/~chenmf/). The improved algorithm below Example 10 in the cited paper is recommended.
  
* `eff.ini.maxeig.shift.inv.tri()`: 
    * `eff.ini.maxeig.shift.inv.tri()` is a safe improvement of `eff.ini.maxeig.tri()`. It is simplified based on Section A.4 in [Chen (2017)](http://math0.bnu.edu.cn/~chenmf/). 
    * The construction of the sequences $A_{j}$ and $B_{j}$ in solving the linear equation is not trivial. This is especially meaningful when going to the large scale matrices.
   
* `eff.ini.maxeig.general()`: calculate the maximal eigenpair for the general matrix with non-negative off-diagonal elements. Let the general matrix $A=(a_{ij}:i,j\in E)$.
    * The initial vector $v_0$ is given by the option `v0_tilde`, if `v0_tilde=NULL`, the one computed in Section 4.2 in the paper [Chen (2016)](http://math0.bnu.edu.cn/~chenmf/) is used. Typically, we use two different initial vectors $v_0$, the uniformly distributed one (Section 4.1) and the one computed in Section 4.2, below (10). The first one is safer but less effective and we encourage to use the second one.
    * There are three choices of the initial $z_0$ as the approximation of $\rho(A)$ which is the maximal eigenvalue of the general matrix $A$. The option `z_0="fixed"` corresponding to the choice I, i.e, set $z_0=\sup_{i\in E}A_{i}$, where $A_{i}=\sum_{j\in E}a_{ij}$. The option `z_0="Auto"` corresponding to the choice II, i.e, set $z_0=v_0^{\ast}Av_{0}$. This simpler choice is used in the other ($z_k:k\ge 1$), but it may be dangerous as initial $z_0$ as illustrated by Example 20 above/below. The option `z_0="numeric"` corresponding to the choice III, as mentioned in the first algorithm for diagonally dominant matrices.
    * The option `z0numeric` is only used when the initial $z_0$ is considered based on (11) in the paper [Chen (2016)](http://math0.bnu.edu.cn/~chenmf/).
    * Similarly, there are other options `xi` and `digit.thresh` which are the same as defined in function `eff.ini.maxeig.tri()`. The improved algorithm is recommended. However, here we should choose a smaller `xi`, say 1/3 or 0.65 used in Examples 18-20 below.

* `eff.ini.secondeig.tri()`: calculate the next to maximal eigenpair for the tridiagonal matrix. As mentioned in the cited paper, for simplicity, here we assume that the sums of each row of the input tridiagonal matrix should be 0, i.e, $A_{i}=0$ for all $i\in E$. Similarly, there are other options `xi` and `digit.thresh` which are the same as defined in function `eff.ini.maxeig.tri()`.

* `eff.ini.secondeig.general()`: calculate the next to maximal eigenpair for the general conservative matrix where the conservativity of matrix means that the sums of each row are all 0, i.e, $A_{i}=0$ for all $i\in E$.
    * There are two choices of the initial $z_0$ as the approximation of $\lambda_1(A)$ which is the second largest eigenvalue of the general matrix $A$. The option `z_0="fixed"` corresponding to the approximation to be $z_0\approx\lambda_0(A_1)$, where $A_1$ is an auxiliary matrix $A$-matrix given in below. The option `z_0="Auto"` corresponding to the approximation given in Section 6 in the paper [Chen (2016)](http://math0.bnu.edu.cn/~chenmf/). It may be necessary to use `z_0="fixed"`, especially for large matrices.
    * A large constant $c_1$ is provided to replace the last diagonal element $A[N,N]$ to generate the auxiliary matrix $A$-matrix $A_1$. Similarly, there is an option `digit.thresh` which is the same as defined in function `eff.ini.maxeig.tri()`.

There are two auxiliary functions `tridia()` and `ray.quot()` where:

- `tridiag()`: generate tridiagonal matrix $A$ based on three input vectors.
- `ray.quot.tri()`: rayleigh quotient iteration algorithm for computing the maximal eigenpair of tridigonal matrix $Q$.
- `shift.inv.tri()`: shifted inverse iteration algorithm using $\delta_k$ for computing the maximal eigenpair of tridigonal matrix $Q$.
- `ray.quot.seceig.tri()`: rayleigh quotient iteration algorithm for computing the next to maximal eigenpair of tridigonal matrix $Q$.
- `ray.quot.general()`: rayleigh quotient iteration algorithm for computing the maximal eigenpair of general matrix $A$.
- `ray.quot.seceig.general()`: rayleigh quotient iteration algorithm for computing the next to maximal eigenpair of general matrix $A$.
- `tri.sol()`: solve the linear equation (-$Q$-$z$*$I$)$w$=$v$.
- `delta()`: compute $\delta_k$ for given vector $v$ and matrix $Q$ based on Section A.4 in [Chen (2017)](http://math0.bnu.edu.cn/~chenmf/). 

## Acknowledgement
The research project is supported in part by the National Natural Science Foundation of China (No. 11131003, 11626245, 11771046) and the Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions. The first author acknowledges the assisted one for his hard work in preparing this package in R.

## Examples

We show most of the examples in the paper [Chen (2016)](http://math0.bnu.edu.cn/~chenmf/) in this section. For a convenient comparison, we keep the same subsection names and examples as
the paper [Chen (2016)](http://math0.bnu.edu.cn/~chenmf/).

### Efficient initials. The tridiagonal case.
In this subsection, armed with `eff.ini.maxeig.tri()`, we can calculate the maximal eigenpair for the tridiagonal matrix. Sometimes we will report the minimal eigenvalue of $-Q$ by the convention in the paper [Chen (2016)](http://math0.bnu.edu.cn/~chenmf/).

#### The minimal eigenvalue of -Q Example 7 (Same as in Example 1)
```{r, echo = TRUE, results='hold', tidy=TRUE}
a = c(1:7)^2
b = c(1:7)^2
c = rep(0, length(a) + 1)
c[length(a) + 1] = 8^2
-eff.ini.maxeig.tri(a, b, c, xi = 1)$z
```

#### Example 7 (Ordinary and Improved)
In the tridiagonal case, the improved algorithm is presented below Example 10 in [Chen (2016)](http://math0.bnu.edu.cn/~chenmf/).
```{r, eval = FALSE, echo=TRUE, results='hold', tidy=TRUE}
nn = c(100, 500, 1000, 5000, 7500, 10^4)
for (i in 1:6) {
    a = c(1:(nn[i] - 1))^2
    b = c(1:(nn[i] - 1))^2
    c = rep(0, length(a) + 1)
    c[length(a) + 1] = nn[i]^2

    print(-eff.ini.maxeig.tri(a, b, c, xi = 1)$z)
    print(-eff.ini.maxeig.tri(a, b, c, xi = 7/8)$z)
}
```
```{r, eval = TRUE, echo=FALSE, results='hold', tidy=TRUE}
print(c(0.348549, 0.376437, 0.376383))
print(c(0.387333, 0.376393, 0.376383))
print(c(0.310195, 0.338402, 0.338329))
print(c(0.349147, 0.338342, 0.338329))
print(c(0.299089, 0.327320, 0.327240))
print(c(0.338027, 0.327254, 0.327240))
print(c(0.281156, 0.308623, 0.308529))
print(c(0.319895, 0.308550, 0.308529))
print(c(0.277865, 0.305016, 0.304918))
print(c(0.316529, 0.304942, 0.304918))
print(c(0.275762, 0.302660, 0.302561))
print(c(0.314370, 0.302586, 0.302561))
```
###### With safe improvement
```{r, eval = FALSE, echo=TRUE, results='hold', tidy=TRUE}
nn = c(8, 100, 500, 1000, 5000, 7500, 10^4)
for (i in 1:7) {
    a = c(1:(nn[i] - 1))^2
    b = c(1:(nn[i] - 1))^2
    c = rep(0, length(a) + 1)
    c[length(a) + 1] = nn[i]^2

    print(-eff.ini.maxeig.shift.inv.tri(a, b, c, xi = 1)$z)
}
```
```{r, eval = TRUE, echo=FALSE, results='hold', tidy=TRUE}
print(c(0.485985, 0.524150, 0.525267, 0.525268))
print(c(0.348549, 0.374848, 0.376378, 0.376383))
print(c(0.310195, 0.336860, 0.338320, 0.338329))
print(c(0.299089, 0.325735, 0.327229, 0.327240))
print(c(0.281156, 0.306874, 0.308514, 0.308529))
print(c(0.277865, 0.303213, 0.304903, 0.304918))
print(c(0.275762, 0.300821, 0.302545, 0.302561))
```

#### The maximal eigenvalue of A Example 8 (Due to L.K.Hua)
```{r, echo=TRUE, results='hold', tidy=TRUE}
a = 14/100
b = 40/100
c = c(-25/100 - 40/100, -12/100 - 14/100)
eff.ini.maxeig.tri(a, b, c, xi = 1)$z
eff.ini.maxeig.tri(a, b, c, xi = 7/8)$z
```
###### With safe improvement
```{r, echo=TRUE, results='hold', tidy=TRUE}
eff.ini.maxeig.shift.inv.tri(a, b, c, xi = 1)$z
eff.ini.maxeig.shift.inv.tri(a, b, c, xi = 7/8)$z
```

#### Example 10 (tridiagonal case)
```{r, echo=TRUE, results='hold', tidy=TRUE}
a = c(sqrt(10), sqrt(130)/11)
b = c(11/sqrt(10), 20 * sqrt(130)/143)
c = c(-1 - 11/sqrt(10), -25/11 - sqrt(10) - 20 * sqrt(130)/143, -8/11 -
    sqrt(130)/11)
eff.ini.maxeig.tri(a, b, c, xi = 1, digit.thresh = 5)$z
eff.ini.maxeig.tri(a, b, c, xi = 7/8, digit.thresh = 5)$z
```
###### With safe improvement
```{r, echo=TRUE, results='hold', tidy=TRUE}
eff.ini.maxeig.shift.inv.tri(a, b, c, xi = 1, digit.thresh = 5)$z
eff.ini.maxeig.shift.inv.tri(a, b, c, xi = 7/8, digit.thresh = 5)$z
```

#### Example A3 in [Chen (2017)](http://math0.bnu.edu.cn/~chenmf/)
```{r, echo=TRUE, results='hold', tidy=TRUE}
a = c(0.5142, 0.2115, 0.8442, 0.2347, 0.9837)
b = c(0.9962, 0.1111, 0.1405, 0.7595, 0.0781)
c = c(-2.334-0.9962, -2.6725-0.5142-0.1111, -2.263-0.2115-0.1405, -2.8457-0.8442-0.7595, -2.2257-0.2347-0.0781, -2.1582-0.9837)
eff.ini.maxeig.tri(a, b, c, xi = 1)$z
```
###### With safe improvement
```{r, echo=TRUE, results='hold', tidy=TRUE}
eff.ini.maxeig.shift.inv.tri(a, b, c, xi = 1)$z
```
From the above we can see that using the idea in the paper [Chen (2017)](http://math0.bnu.edu.cn/~chenmf/)...($z_{k}$=$\delta_{k}^{-1}$) may need one or two more steps than the [Chen (2016)](http://math0.bnu.edu.cn/~chenmf/) but [Chen (2017)](http://math0.bnu.edu.cn/~chenmf/)'s improvement is safe and never fall into pitfall, so we can use the new one instead of the old in 2016.


### Efficient initials. The general case.
To get the maximal eigenpair for the general matrix $Q$, there are several choices of the initial vector $v_0$ and the initial $z_0$. In this subsection, we study several combinations of the choices.

#### Fixed Uniformly distributed initial vector $v_0=(1,\dots,1)/\sqrt(N+1)$.
##### Choice I for $z_0=\sup_{i\in E}A_{i}$
###### Example 13
```{r, echo=TRUE, results='hold', tidy=TRUE}
A = matrix(c(1, 1, 3, 2, 2, 2, 3, 1, 1), 3, 3)
eff.ini.maxeig.general(A, v0_tilde = rep(1, dim(A)[1]), z0 = "fixed", digit.thresh = 5)$z
```

###### Example 14
```{r, echo=TRUE, results='hold', tidy=TRUE}
A = t(matrix(seq(1, 16), 4, 4))
eff.ini.maxeig.general(A, v0_tilde = rep(1, dim(A)[1]), z0 = "fixed", digit.thresh = 4)$z
```

###### Example 15
```{r, echo=TRUE, results='hold', tidy=TRUE}
A = matrix(c(1, 3, 9, 5, 2, 14, 10, 6, 0, 11, 11, 7, 0, 0, 1, 8),
    4, 4)
eff.ini.maxeig.general(A, v0_tilde = rep(1, dim(A)[1]), z0 = "fixed", digit.thresh = 4)$z
```

###### Example 16 (Same as Example 1)
```{r, echo=TRUE, results='hold', tidy=TRUE}
a = c(1:7)^2
b = c(1:7)^2
c = rep(0, length(a) + 1)
c[length(a) + 1] = 8^2

N = length(a)
Q = tridiag(b, a, -c(b[1] + c[1], a[1:N - 1] + b[2:N] + c[2:N], a[N] +
    c[N + 1]))

A = 113 * diag(1, (N + 1)) + Q

113 - eff.ini.maxeig.general(A, v0_tilde = rep(1, dim(A)[1]), z0 = "fixed")$z
```

##### Choice II for $z_0=v_0^{\ast}Av_{0}$ for Example 13-15
###### Example 13
```{r, echo=TRUE, results='hold', tidy=TRUE}
A = matrix(c(1, 1, 3, 2, 2, 2, 3, 1, 1), 3, 3)
eff.ini.maxeig.general(A, v0_tilde = rep(1, dim(A)[1]), z0 = "Auto", digit.thresh = 5)$z
```

###### Example 14
```{r, echo=TRUE, results='hold', tidy=TRUE}
A = t(matrix(seq(1, 16), 4, 4))
eff.ini.maxeig.general(A, v0_tilde = rep(1, dim(A)[1]), z0 = "Auto", digit.thresh = 4)$z
```

###### Example 15
```{r, echo=TRUE, results='hold', tidy=TRUE}
A = matrix(c(1, 3, 9, 5, 2, 14, 10, 6, 0, 11, 11, 7, 0, 0, 1, 8),
    4, 4)
eff.ini.maxeig.general(A, v0_tilde = rep(1, dim(A)[1]), z0 = "Auto", digit.thresh = 4)$z
```

##### Choice III for numeric $z_0$
Here we use the combination coefficient `xi=7/8` which is a little different with [Chen (2016)](http://math0.bnu.edu.cn/~chenmf/).

###### Example 17: Symmetrizing matrix $A$
```{r, echo=TRUE, results='hold', tidy=TRUE}
A = matrix(c(1, 3, 9, 5, 2, 14, 10, 6, 0, 11, 11, 7, 0, 0, 1, 8),
    4, 4)
S = (t(A) + A)/2
N = dim(S)[1]
a = diag(S[-1, -N])
b = diag(S[-N, -1])
c = rep(NA, N)
c[1] = -diag(S)[1] - b[1]
c[2:(N - 1)] = -diag(S)[2:(N - 1)] - b[2:(N - 1)] - a[1:(N - 2)]
c[N] = -diag(S)[N] - a[N - 1]

z0ini = eff.ini.maxeig.tri(a, b, c, xi = 7/8)$z[1]
z0ini
eff.ini.maxeig.general(A, v0_tilde = rep(1, dim(A)[1]), z0 = "numeric",
    z0numeric = 28 - z0ini, digit.thresh = 4)$z
```

###### Example 17: Without Symmetrizing matrix $A$
```{r, echo=TRUE, results='hold', tidy=TRUE}
A = matrix(c(1, 3, 9, 5, 2, 14, 10, 6, 0, 11, 11, 7, 0, 0, 1, 8),
    4, 4)
S = A
N = dim(S)[1]
a = diag(S[-1, -N])
b = diag(S[-N, -1])
c[1] = -diag(S)[1] - b[1]
c[2:(N - 1)] = -diag(S)[2:(N - 1)] - b[2:(N - 1)] - a[1:(N - 2)]
c[N] = -diag(S)[N] - a[N - 1]

z0ini = eff.ini.maxeig.tri(a, b, c, xi = 7/8)$z[1]
z0ini
eff.ini.maxeig.general(A, v0_tilde = rep(1, dim(A)[1]), z0 = "numeric",
    z0numeric = 31 - z0ini, digit.thresh = 4)$z
```

####Efficient initial vector $v_0$
The improved algorithm is presented at the end of Section 4 in [Chen (2016)](http://math0.bnu.edu.cn/~chenmf/). The outputs at the last line in Examples 18 and 19 are different from those presented in the cited paper, due to the reason that
a slight different combination was used there.

##### Example 18 (Same as Example 10)
```{r, echo=TRUE, results='hold', tidy=TRUE}
A = matrix(c(1, 1, 3, 2, 2, 2, 3, 1, 1), 3, 3)
eff.ini.maxeig.general(A, z0 = "Auto", digit.thresh = 5)$z
eff.ini.maxeig.general(A, xi = 1, digit.thresh = 5)$z
eff.ini.maxeig.general(A, xi = 1/3, digit.thresh = 5)$z
eff.ini.maxeig.general(A, xi = 0, digit.thresh = 5)$z
```

##### Example 19 (Same as Example 14)
```{r, echo=TRUE, results='hold', tidy=TRUE}
A = t(matrix(seq(1, 16), 4, 4))
eff.ini.maxeig.general(A, z0 = "Auto", digit.thresh = 4)$z
eff.ini.maxeig.general(A, xi = 1, digit.thresh = 4)$z
eff.ini.maxeig.general(A, xi = 1/3, digit.thresh = 4)$z
eff.ini.maxeig.general(A, xi = 0, digit.thresh = 4)$z
```

##### Example 20 (Same as Example 15)
```{r, echo=TRUE, results='hold', tidy=TRUE}
A = matrix(c(1, 3, 9, 5, 2, 14, 10, 6, 0, 11, 11, 7, 0, 0, 1, 8),
    4, 4)
eff.ini.maxeig.general(A, z0 = "Auto", digit.thresh = 4)$z
eff.ini.maxeig.general(A, xi = 1, digit.thresh = 4)$z
eff.ini.maxeig.general(A, xi = 0.65, digit.thresh = 4)$z
eff.ini.maxeig.general(A, xi = 0, digit.thresh = 4)$z
```

##### Example 21
```{r, echo=TRUE, results='hold', tidy=TRUE}
b4 = c(0.01, 1, 100)
digits = c(9, 7, 6)

for (i in 1:3) {
    A = matrix(c(-3, 4, 0, 10, 0, 2, -7, 5, 0, 0, 0, 3, -5, 0, 0,
        1, 0, 0, -16, 11, 0, 0, 0, 6, -11 - b4[i]), 5, 5)
    print(-eff.ini.maxeig.general(A, z0 = "Auto", digit.thresh = digits[i])$z)
}
```

##### Example 22
###### The result given by general algorithm
```{r, echo=TRUE, results='hold', tidy=TRUE}
b4 = c(0.01, 1, 100)
digits = c(9, 7, 6)

for (i in 1:3) {
    A = matrix(c(-5, 3, 0, 0, 0, 5, -7, 2, 0, 0, 0, 4, -3, 10, 0,
        0, 0, 1, -16, 11, 0, 0, 0, 6, -11 - b4[i]), 5, 5)
    print(-eff.ini.maxeig.general(A, z0 = "Auto", digit.thresh = digits[i])$z)
}
```

###### The result given by tri algorithm
```{r, echo=TRUE, results='hold', tidy=TRUE}
b4 = c(0.01, 1, 100)
a = c(3, 2, 10, 11)
b = c(5, 4, 1, 6)
b4 = c(0.01, 1, 100, 10^6)
digits = c(9, 7, 6, 6)

for (i in 1:4) {
    c = c(rep(0, 4), b4[i])
    print(-eff.ini.maxeig.tri(a, b, c, xi = 1, digit.thresh = digits[i])$z)
}
```


#### Example A3 in [Chen (2017)](http://math0.bnu.edu.cn/~chenmf/)
```{r, echo=TRUE, results='hold', tidy=TRUE}
a = c(0.5142, 0.2115, 0.8442, 0.2347, 0.9837)
b = c(0.9962, 0.1111, 0.1405, 0.7595, 0.0781)
c = c(-2.334-0.9962, -2.6725-0.5142-0.1111, -2.263-0.2115-0.1405, -2.8457-0.8442-0.7595, -2.2257-0.2347-0.0781, -2.1582-0.9837)

N = length(a)
A = tridiag(b, a, -c(b[1] + c[1], a[1:N - 1] + b[2:N] + c[2:N], a[N] + c[N + 
        1]))
eff.ini.maxeig.general(A, xi = 1, digit.thresh = 5)$z
eff.ini.maxeig.general(A, xi = 0.18, digit.thresh = 5)$z
eff.ini.maxeig.general(A, xi = 0, digit.thresh = 5)$z
```
From this example, we can see that when `xi=0`, it goes into pitfall. So we
may choose `xi` does not equal to 0.

### The next to the maximal eigenpair. 
### Tridiagonal matrix case.
#### Example 25
```{r, echo=TRUE, results='hold', tidy=TRUE}
a = c(1:7)^2
b = c(1:7)^2

eff.ini.seceig.tri(a, b, xi = 0)$z
eff.ini.seceig.tri(a, b, xi = 1)$z
eff.ini.seceig.tri(a, b, xi = 2/5)$z
```

#### Example 26
```{r, echo=TRUE, results='hold', tidy=TRUE}
a = c(3, 2, 10, 11)
b = c(5, 4, 1, 6)

eff.ini.seceig.tri(a, b, xi = 0, digit.thresh = 5)$z
eff.ini.seceig.tri(a, b, xi = 1, digit.thresh = 5)$z
eff.ini.seceig.tri(a, b, xi = 2/5, digit.thresh = 5)$z
```

### General matrix case.
The results of two choices of the initial $z_0$ are presented in the following.

#### Example 27
```{r, echo=TRUE, results='hold', tidy=TRUE}
Q = matrix(c(-30, 1/5, 11/28, 55/3291, 30, -17, 275/42, 330/1097,
    0, 84/5, -20, 588/1097, 0, 0, 1097/84, -2809/3291), 4, 4)
eff.ini.seceig.general(Q, z0 = "Auto", digit.thresh = 5)$z
eff.ini.seceig.general(Q, z0 = "fixed", digit.thresh = 5)$z
```

#### Example 28
```{r, echo=TRUE, results='hold', tidy=TRUE}
Q = matrix(c(-57, 135/59, 243/91, 351/287, 118/27, -52, 590/91, 118/41,
    91/9, 637/59, -47, 195/41, 1148/27, 2296/59, 492/13, -62/7),
    4, 4)
eff.ini.seceig.general(Q, z0 = "Auto", digit.thresh = 4)$z
eff.ini.seceig.general(Q, z0 = "fixed", digit.thresh = 4)$z
```

# References
[1] M. F. Chen. "Efficient initials for computing maximal eigenpair". In: _Frontiers of Mathematics in China_ 11.6 (2016), pp. 1379-1418.

[2] M. F. Chen. "Global algorithms for maximal eigenpair". In: _Frontiers of Mathematics in China_ 12.5 (2017), pp. 1023-1043.