4

What Is the Frank Matrix?

 2 years ago
source link: https://nhigham.com/2022/01/25/what-is-the-frank-matrix/
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.

What Is the Frank Matrix?

The n\times n upper Hessenberg matrix

\notag F_n = \left[\begin{array}{*{6}c} n & n-1 & n-2 & \dots & 2 & 1 \\ n-1 & n-1 & n-2 & \dots & 2 & 1 \\ 0 & n-2 & n-2 & \dots & 2 & 1 \\[-3pt] \vdots & 0 & \ddots & \ddots & \vdots & 1 \\[-3pt] \vdots & \vdots & \dots & 2 & 2 & 1 \\ 0 & 0 & \dots & 0 & 1 & 1 \\ \end{array}\right]

was introduced by Frank in 1958 as a test matrix for eigensolvers.

Determinant

Taking the Laplace expansion about the first column, we obtain \det(F_n) = n\det(F_{n-1}) - (n-1)\det(F_{n-1}) = \det(F_{n-1}), and since \det(F_1) = 1 we have \det(F_n) = 1.

In MATLAB, the computed determinant of the matrix and its transpose can be far from the exact values of 1:

>> n = 20; F = anymatrix('gallery/frank',n); 
>> [det(F), det(F')]
ans =
   1.0000e+00  -1.4320e+01
>> n = 49; F = anymatrix('gallery/frank',n); det(F)
ans =
    -1.406934439401568e+45

This behavior illustrates the sensitivity of the determinant rather than numerical instability in the evaluation and it is very dependent on the arithmetic (different results may be obtained in different releases of MATLAB). The sensitivity can be seen by noting that

\notag \det\bigl(F_n + \epsilon e_1e_n^T\bigr) = 1 + (-1)^{n-1} (n-1)!\epsilon, \qquad (1)

which means that changing the (1,n) element from 1 to 1+\epsilon changes the determinant by (n-1)!\epsilon.

Inverse and Condition Number

It is easily verified that

\notag F_n = \begin{bmatrix} 1 & -1 & & & \\ & 1 & -1 & & \\ & & 1 & \ddots & \\ & & & \ddots & -1\\ & & & & 1 \end{bmatrix}^{-1} \begin{bmatrix} 1 & & & & \\ n-1 & 1 & & & \\ & n-2& 1 & & \\ & & \ddots& \ddots & \\ & & & 1 & 1 \end{bmatrix} \equiv U^{-1}L. \qquad (2)

Hence F_n^{-1} = L^{-1}U is lower Hessenberg with integer entries. This factorization provides another way to see that \det(F_n) = 1.

From (1) we see that F_n + \epsilon e_1e_n^T is singular for \epsilon = (-1)^n/(n-1)!, which implies that

\notag \kappa(F_n) \ge (n-1)! \|F_n\|

for any subordinate matrix norm, showing that the condition number grows very rapidly with n. In fact, this lower bound is within a factor 3.5 of the condition number for the 1-, 2-, and \infty-norms for n\le 20.

Eigenvalues

Denote by p_n(t) = \det(F_n - tI) the characteristic polynomial of F_n. By expanding about the first column one can show that

\label{Fnpn-rec} p_n(t) = (1-t)p_{n-1}(t) - (n-1)t p_{n-2}(t). \qquad (3)

Using (3), one can show by induction that

\notag p_n(t^{-1}) = (-1)^n t^{-n} p_n(t).

This means that the eigenvalues of F_n occur in reciprocal pairs, and hence that \lambda = 1 is an eigenvalue when n is odd. It also follows that p_n is palindromic when n is even and anti-palindromic when n is odd. Examples:

>> charpoly( anymatrix('gallery/frank',6) )
ans =
     1   -21   120  -215   120   -21     1
>> charpoly( anymatrix('gallery/frank',7) )
ans =
     1   -28   231  -665   665  -231    28    -1

Since F_n has nonzero subdiagonal entries, \mathrm{rank}(F_n - t I) \ge n-1 for any t, and hence F_n is nonderogatory, that is, no eigenvalue appears in more than one Jordan block in the Jordan canonical form. It can be shown that the eigenvalues are real and positive and that they can be expressed in terms of the zeros of Hermite polynomials. Furthermore, the eigenvalues are distinct.

If each eigenvalue of an n\times n matrix undergoes a small perturbation then the determinant, being the product of the eigenvalues, undergoes a perturbation up to about n times as large. From (1) we see that a change to F_n of order \epsilon can perturb \det(F_n) by (n-1)!\epsilon, and it follows that some eigenvalues must undergo a large relative perturbation. The condition number of a simple eigenvalue is defined as the reciprocal of the cosine of the angle between its left and right eigenvectors. For the Frank matrix it is the small eigenvalues that are ill conditioned, as shown here for n = 9.

>> n = 9; F = anymatrix('gallery/frank',n);
>> [V,D,cond_evals] = condeig(F); evals = diag(D);
>> [~,k] = sort(evals,'descend');
>> [evals(k) cond_evals(k)]
ans =
   2.2320e+01   1.9916e+00
   1.2193e+01   2.3903e+00
   6.1507e+00   1.5268e+00
   2.6729e+00   3.6210e+00
   1.0000e+00   6.8615e+01
   3.7412e-01   1.5996e+03
   1.6258e-01   1.1907e+04
   8.2016e-02   2.5134e+04
   4.4803e-02   1.4762e+04

Notes

Frank found that F_n “gives our selected procedures difficulties” and that “accuracy was lost in the smaller roots”. The difficulties encountered by Frank’s codes were shown by Wilkinson to be caused by the inherent sensitivity of the eigenvalues to perturbations in the matrix.

References

This is a minimal set of references, which contain further useful references within.

Posted on January 25, 2022January 25, 2022 by Nick HighamPosted in what-is

Post navigation

One thought on “What Is the Frank Matrix?”

  1. Some time ago, I had devised a modified version of the Frank matrix, using the ideas in Varah’s paper (https://doi.org/10.1137/0907056). In particular, using his recipe with the Clement-Kac matrix (MATLAB’s gallery(‘clement’)), the $n\times n$ upper Hessenberg matrix with entries $h_{i,j} = j (n-j)+[i \leq j]$ (where $[p]$ is an Iverson bracket) has determinant $1$ and eigenvalues $\left(\frac{k}{2}+\sqrt{1+\frac{k^2}{4}}\right)^2$, $k=-(n-1), -(n-3), \dots, n-1$. It seems to be even more badly behaved than the original.

Leave a Reply Cancel reply


report

Recommend

About Joyk


Aggregate valuable and interesting links.
Joyk means Joy of geeK