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) and Chen (2017) 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.
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")
require(EfficientMaxEigenpair)
Before moving to the details, let us illustrate the algorithm by two typical examples.
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.
# 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.
# 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}
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")
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.
# 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 matrix we considered for this general case is
# 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:
# 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.
# 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.
# 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).
digit.thresh = 6
which implies 1e-6 without any special requirement. Same for the following three algorithms and the examples shown in this vignette.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). 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)$.
v0_tilde
, if v0_tilde=NULL
, the one computed in Section 4.2 in the paper Chen (2016) 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.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.z0numeric
is only used when the initial $z_0$ is considered based on (11) in the paper Chen (2016).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$.
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). It may be necessary to use z_0="fixed"
, especially for large matrices.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). 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.
We show most of the examples in the paper Chen (2016) in this section. For a convenient comparison, we keep the same subsection names and examples as the paper Chen (2016).
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).
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
In the tridiagonal case, the improved algorithm is presented below Example 10 in Chen (2016).
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) }
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))
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) }
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))
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
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
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
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
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
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)...($z_{k}$=$\delta_{k}^{-1}$) may need one or two more steps than the Chen (2016) but Chen (2017)'s improvement is safe and never fall into pitfall, so we can use the new one instead of the old in 2016.
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.
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
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
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
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
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
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
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
Here we use the combination coefficient xi=7/8
which is a little different with Chen (2016).
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
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
The improved algorithm is presented at the end of Section 4 in Chen (2016). 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.
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
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
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
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) }
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) }
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) }
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.
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
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
The results of two choices of the initial $z_0$ are presented in the following.
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
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
[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.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.